// Implementation of Latin Square Algorithm 5/31/2003 by MJ
//      1. do "gcc -o latin latin.c" and "latin 3 2 1"
//      2. changed check condition
//      3. added clock
// Code description

// The algorithms for building valid L-G and L-L schedules
// both use a subroutine implementing the Birkhoff-von
// Neumann decomposition theorem using edge coloring in a regular
// bipartite graph. For a regular bipartite graph in which all nodes
// have a degree $k$, this subroutine distinguishes two cases. If $k$
// is even, it will use Euler splits in order to split the
// graph into two sub-graphs of degree $k/2$, and will recursively
// continue $n$ times until the degree $k/2^n$ of the $2^n$ resulting
// sub-graphs is odd. If the degree $k/2^n$ is equal to $1$, each of
// the sub-graphs is a perfect match, and thus the subroutine stops
// because th edge coloring is done. Otherwise, if the degree $l$ of
// each sub-graph is odd and different from $1$, the subroutine
// computes a perfect match in each sub-graph using Ford-Fulkerson.
// It then removes the match from the sub-graph, obtaining a regular
// bipartite sub-graph of even degree $l-1$. It can then apply the
// same subroutine to this sub-graph using Euler splits, since $l-1$
// is even, and so on iteratively. The algorithm clearly converges
// since the degree always decreases \cite{cole,schrijver}.



#include <stdio.h>
#include <math.h>
#include <stdlib.h> //isaac
#include <time.h> //isaac

//Basic Definitions for Max flow algorithm
#define WHITE 0
#define GRAY 1
#define BLACK 2
#define MAX_NODES 2000 
#define oo 1000000000

//basic definitions for english square
#define MAX_TIME_UNITS 3000   //maximum number of time units
#define MAX_LINE_CARDS 3000   //maximum number of line cards
#define MAX_GROUPS 52       //maximum number of line card groups

//Declarations for max-flow algorithm
int n, group;  // number of nodes
int e;  // number of edges
int capacity[MAX_NODES][MAX_NODES]; // capacity matrix
int flow[MAX_NODES][MAX_NODES];     // flow matrix
int color[MAX_NODES]; // needed for breadth-first search               
int pred[MAX_NODES];  // array to store augmenting path

//declaration of the latin square algorithm
int NN;                           //number of time units
int G;                            //number of groups 
int L[MAX_GROUPS];                //number of line cards in each group 
int begin[MAX_GROUPS];            //beginning of index of a group
int end[MAX_GROUPS];              //end of index of a group
int M[MAX_GROUPS][MAX_GROUPS];    //working array of M
int P[MAX_GROUPS][MAX_GROUPS];    //working array of P,Q,R,S,a,b,linecard
int Q[MAX_GROUPS][MAX_GROUPS]; 
int R[MAX_GROUPS][MAX_GROUPS]; 
int S[MAX_GROUPS][MAX_GROUPS][MAX_TIME_UNITS]; 
int a[MAX_GROUPS], b[MAX_GROUPS]; 
int linecard[MAX_LINE_CARDS];
int letter[MAX_LINE_CARDS];       //quick lookup of a line card group id
int row_max[MAX_LINE_CARDS];  
int col_max[MAX_LINE_CARDS];
int es[MAX_LINE_CARDS][MAX_LINE_CARDS];                    //english square
int neighbour[MAX_LINE_CARDS][MAX_LINE_CARDS];                    //english square
int final[MAX_LINE_CARDS][MAX_LINE_CARDS];                 //final output array of the latin square
//int check[MAX_LINE_CARDS][MAX_LINE_CARDS][MAX_LINE_CARDS]; //to check if the resulted array satisfy the MEMS requirement
int XX[MAX_LINE_CARDS][MAX_LINE_CARDS];                            //working array
int ZZ[MAX_LINE_CARDS][MAX_LINE_CARDS];                            //working array
clock_t mytime; //isaac
double cpu_time_used; //isaac

int W[640][640], W_8[640][640][2], W_4[640][640][4], W_2[640][640][8], W_1[640][640][16]; 
int D[640];

//=============================================================================
//=============================================================================
int main(int argc, char **argv){ 

   int i,j,k, sum, temp, x, y, ii,jj,lll, card, c;
   int neig_c[MAX_LINE_CARDS],temp1;

   //the following is the input for the algorithm 
   //    NN  -  total number if line cards
   //    G   -  number of line card groups
   //    L[i]-  number of line cards in each group

   G=argc-1;
	printf("G = %d\n", G);
   for (i= 0; i < G; i++) {
    	L[i]=atoi(argv[i+1]);	
    	printf("L[%i] = %d\n",i,L[i]);
	}

   NN=0;
   for(i=0; i<G; i++){
	NN=NN+L[i];
   }
   
   //initilization of M[][]
   for(i=0; i<G; i++){
      for(j=0; j<G; j++){
         M[i][j] = L[i]*L[j];
      }
   }

   for(i=0; i<NN; i++){
         linecard[i] = i+1;
   }
   
   //inplementation of the algorithm --- part 1: colum to column
   for(k=NN; k>0; k--){

   	for(i=0; i<G; i++){
      	   for(j=0; j<G; j++){
              P[i][j] = M[i][j]/k;
           }
        }

   	for(i=0; i<G; i++){
      	   for(j=0; j<G; j++){
              Q[i][j] = M[i][j]-P[i][j]*k;
           }
        }

   	for(i=0; i<G; i++){
           sum = 0;
      	   for(j=0; j<G; j++){
              sum = sum + Q[i][j];
           }
           a[i]=sum/k;
        }
        
   	for(j=0; j<G; j++){
           sum = 0;
      	   for(i=0; i<G; i++){
              sum = sum + Q[i][j];
           }
           b[j]=sum/k;
        }

       n = 2*G+2;
       e = G + G*G + G;
       
       // initialize empty capacity matrix 
       for (i=0; i<n; i++) {
           for (j=0; j<n; j++) {
	       capacity[i][j] = 0;
	   }
       }
       
       // read edge capacities
       for (i=0; i<G; i++) {
           capacity[0][i+1] = a[i];
       }

       for (i=1; i<=G; i++) {
           for (j=1+G; j<=G+G; j++) {
               if(Q[i-1][j-1-G]>0){
	           capacity[i][j] = 1;
               }
           }
       }

       for (i=G+1; i<=G+G; i++) {
           capacity[i][G+G+1] = b[i-G-1];
       }

       temp = max_flow(0,n-1);


       for (i=0; i<G; i++) {
           for (j=0; j<G; j++) {
	       R[i][j] = flow[i+1][j+1+G];
           }
       }

       for (i=0; i<G; i++) {
           for (j=0; j<G; j++) {
	       S[i][j][k] = P[i][j] + R[i][j];
	       M[i][j] = M[i][j] - S[i][j][k];
           }
       }

   }

   //implementation of the algorithm --- part 2: from the english square to a Deiscrete Programming Formulation
   sum = 0;
   for(i=0;i<G;i++){
      begin[i]=sum;
      end[i]=sum+L[i];
      sum = end[i];

      for(j=begin[i]; j<end[i]; j++){
         letter[j]=i;
      }
   }

   lll=0;

   for(x=0; x<G; x++){ 

       printf("group %d\n", x);

       for (ii=0; ii<NN; ii++) {
           for (jj=0; jj<NN; jj++) {
                XX[ii][jj]=0;
           }
       }

       for(ii=0; ii<NN+NN+1; ii++){
           neig_c[ii]=0;
       }


       for(y=0; y<G; y++){ 
           for(j=0; j<NN; j++){
           col_max[j]= S[x][y][j+1];
           }
           temp=0;
           for (jj=0; jj<NN; jj++) {
                for (ii=begin[y]; ii<end[y]; ii++) {
                     if(col_max[jj]>0){
                         XX[temp%L[y]+begin[y]][jj]=1;
                         temp1=temp%L[y]+begin[y];
                         neighbour[temp1+1][neig_c[temp1]]=jj + NN +1;
                         neig_c[temp1]++;
                         neighbour[jj+NN+1][neig_c[jj+NN]]=temp1 + 1;
                         neig_c[jj+NN]++;
                         temp++;
                         col_max[jj]--;
                     }
                }
           }
       }

if(L[x]==16){
      temp=euler_partition_16();
      for(k=0; k<L[x]; k++){
           for (i=0; i<NN; i++) {
              for (j=0; j<NN; j++) {
                 if(W_1[i][j][k]==1){
                    es[lll][j] = letter[i];
                 }
              }
           }
           lll++;  //next line of english square
      }
}
else if(L[x]==8){
      temp=euler_partition_8();
      for(k=0; k<L[x]; k++){
           for (i=0; i<NN; i++) {
              for (j=0; j<NN; j++) {
                 if(W_1[i][j][k]==1){
                    es[lll][j] = letter[i];
                 }
              }
           }
           lll++;  //next line of english square
      }
}
else if(L[x]==4){
      temp=euler_partition_4();
      for(k=0; k<L[x]; k++){
           for (i=0; i<NN; i++) {
              for (j=0; j<NN; j++) {
                 if(W_1[i][j][k]==1){
                    es[lll][j] = letter[i];
                 }
              }
           }
           lll++;  //next line of english square
      }
}
else if(L[x]==2){
      temp=euler_partition_2();
      for(k=0; k<L[x]; k++){
           for (i=0; i<NN; i++) {
              for (j=0; j<NN; j++) {
                 if(W_1[i][j][k]==1){
                    es[lll][j] = letter[i];
                 }
              }
           }
           lll++;  //next line of english square
      }
}
else{
      n = 2*NN+2;
      e = NN + NN*NN + NN;

      for(k=0; k<L[x]; k++){

           for (i=0; i<2*NN; i++) {
               for (j=0; j<2*NN; j++) {
                    capacity[i][j] = 0;
               }
           }

           for (i=0; i<NN; i++) {
               capacity[0][i+1] = 1;
           }

           for (i=1; i<=NN; i++) {
               c=0;
               for (j=1+NN; j<=NN+NN; j++) {
                   if(XX[i-1][j-1-NN]>0){
                       capacity[i][j] = 1;
                       neighbour[i][c]=j;
                       c++;
                   }
                   else{
                       capacity[i][j] = 0;
                   }
               }
           }


           for (j=1+NN; j<=NN+NN; j++) {
              c=0;
              for (i=1; i<=NN; i++) {
                  if(XX[i-1][j-1-NN]>0){
                  neighbour[j][c] = i;
                  c++;
                  }
              }
           }




           for (i=NN+1; i<=NN+NN; i++) {
               capacity[i][NN+NN+1] = 1;
           }

           group = L[x]-k;
           temp = new_max_flow(0,n-1);

           for(i=0; i<NN; i++){
              for(j=0; j<NN; j++) {
                 if(flow[i+1][j+1+NN]==1){
                    ZZ[i][k]=j;
                    break;
                 }
              }
           }


           for (i=0; i<NN; i++) {
              for (j=0; j<NN; j++) {
                 if(ZZ[i][k]==j){
                    XX[i][j]--;
                 }
              }
           }

          //compress
           for (i=0; i<NN; i++) {
              for (j=0; j<NN; j++) {
                 if(ZZ[i][k]==j){
                    es[lll][j] = letter[i];
                 }
              }
           }

           for (i=0; i<NN; i++) {
              for (j=0; j<NN; j++) {
                 ZZ[i][j] = 0;
              }
           }

           lll++;  //next line of english square
      }
}
}

   //inplementation of the algorithm --- part 3: from the english square to a latin square

   n = 2*NN+2;
   e = NN + NN*NN + NN;

   card=0;
   for(x=0; x<G; x++){
          printf("group -> %d\n", x);

if(L[x]==16){
          for (i=0; i<NN; i++) {
             for (j=0; j<NN; j++) {
                if(es[i][j]==x){
                   XX[i][j] = 1;
                }
                else{
                   XX[i][j] = 0;
                }
             }
          }

          temp=euler_partition_16();

          for(k=0; k<L[x]; k++){

             for (i=0; i<NN; i++) {
                for (j=0; j<NN; j++) {
                    if(W_1[i][j][k]==1){
                       final[i][j]=card;
                       break;
                    }
                }
             }
        
             card++;
          }
}//if=else
else if(L[x]==8){
          for (i=0; i<NN; i++) {
             for (j=0; j<NN; j++) {
                if(es[i][j]==x){
                   XX[i][j] = 1;
                }
                else{
                   XX[i][j] = 0;
                }
             }
          }

          temp=euler_partition_8();

          for(k=0; k<L[x]; k++){

             for (i=0; i<NN; i++) {
                for (j=0; j<NN; j++) {
                    if(W_1[i][j][k]==1){
                       final[i][j]=card;
                       break;
                    }
                }
             }
        
             card++;
          }
}//if=else
else if(L[x]==4){
          for (i=0; i<NN; i++) {
             for (j=0; j<NN; j++) {
                if(es[i][j]==x){
                   XX[i][j] = 1;
                }
                else{
                   XX[i][j] = 0;
                }
             }
          }

          temp=euler_partition_4();

          for(k=0; k<L[x]; k++){

             for (i=0; i<NN; i++) {
                for (j=0; j<NN; j++) {
                    if(W_1[i][j][k]==1){
                       final[i][j]=card;
                       break;
                    }
                }
             }
        
             card++;
          }
}//if=else88
else if(L[x]==2){
          for (i=0; i<NN; i++) {
             for (j=0; j<NN; j++) {
                if(es[i][j]==x){
                   XX[i][j] = 1;
                }
                else{
                   XX[i][j] = 0;
                }
             }
          }

          temp=euler_partition_2();

          for(k=0; k<L[x]; k++){

             for (i=0; i<NN; i++) {
                for (j=0; j<NN; j++) {
                    if(W_1[i][j][k]==1){
                       final[i][j]=card;
                       break;
                    }
                }
             }
        
             card++;
          }
}//if=else88
else{

      for(k=0; k<L[x]; k++){
          for (i=0; i<2*NN; i++) {
             for (j=0; j<2*NN; j++) {
                capacity[i][j] = 0;
             }
          }

          // read edge capacities
          for (i=0; i<NN; i++) {
             capacity[0][i+1] = 1;
          }

          for (i=1; i<=NN; i++) {
             c=0;
             for (j=1+NN; j<=NN+NN; j++) {
                if(es[i-1][j-1-NN]==x){
                   capacity[i][j] = 1;
                   neighbour[i][c] = j;
                   c++;
                }
                else{
                   capacity[i][j] = 0;
                }
             }
          }


           for (j=1+NN; j<=NN+NN; j++) {
              c=0;
              for (i=1; i<=NN; i++) {
                  if(es[i-1][j-1-NN]==x){
                  neighbour[j][c] = i;
                  c++;
                  }
              }
           }

          for (i=NN+1; i<=NN+NN; i++) {
             capacity[i][NN+NN+1] = 1;
          }

          group = L[x] - k;
          temp = new_max_flow(0,n-1);

          for (i=0; i<NN; i++) {
             for (j=0; j<NN; j++) {
                 if(flow[i+1][j+1+NN]==1){
                    ZZ[i][k]=j;
                    final[i][j]=card;
                    break;
                 }

             }
          }

          for (i=0; i<NN; i++) {
             for (j=0; j<NN; j++) {
                 if(ZZ[i][k]==j){
                    es[i][j]--;
                 }
             }
          }

          card++;


      }


      for (i=0; i<NN; i++) {
         for (j=0; j<NN; j++) {
            ZZ[i][j] = 0;
         }
      }


}


   }

      printf("\n\nOutput of the algorithm -----    \n");
      for (i=0; i<NN; i++) {
          for (j=0; j<NN; j++) {
             printf("%d  ", final[i][j]);
          }
          printf("\n");
      }

} 

//==========================================


int min (int x, int y) {
    return x<y ? x : y;  // returns minimum of x and y
}


//A Queue for Breadth-First Search
int head,tail;
int q[MAX_NODES+2];

void enqueue (int x) {
    q[tail] = x;
    tail++;
    color[x] = GRAY;
}

int dequeue () {
    int x = q[head];
    head++;
    color[x] = BLACK;
    return x;
}

//Breadth-First Search for an augmenting path
int bfs (int start, int target) {
    int aa,bb,num_nodes, extra;
    int u,v;
    for (u=0; u<n; u++) {
	color[u] = WHITE;
    }   
    head = tail = 0;
    enqueue(start);
    pred[start] = -1;
    while (head!=tail) {
	u = dequeue();
        // Search all adjacent white nodes v. If the capacity
        // from u to v in the residual network is positive,
        // enqueue v.
        num_nodes=(n-2)/2;
        if(u==0){
           aa=1;
           bb=num_nodes+1;
        }
        else if(u>=1 && u<=num_nodes){
           aa=num_nodes+1;
           bb=num_nodes+num_nodes+1;
        }
        else if(u>=num_nodes+1 && u<=num_nodes+num_nodes+1){
           aa=1;
           bb=num_nodes+1;
           extra =1;
        }
        else{
           aa=0;
           bb=n;
        }
           
	for (v=aa; v<bb; v++) {
	    if (color[v]==WHITE && capacity[u][v]-flow[u][v]>0) {
		enqueue(v);
		pred[v] = u;
	    }
	}
       
        if(extra==1){
            v=n-1;
	    if (color[v]==WHITE && capacity[u][v]-flow[u][v]>0) {
		enqueue(v);
		pred[v] = u;
	    }

           extra = 0;
        }
 
    }
    // If the color of the target node is black now,
    // it means that we reached it.
    return color[target]==BLACK;
}


int euler_partition_16(void){

   int ii,jj, kk,d,kkk,x;
   int neigh[2*640][32], neigh_c[2*640];

   for (ii=0; ii<NN; ii++) {
      D[ii]=16;
      for (jj=0; jj<NN; jj++) {
          W[ii][jj]=XX[ii][jj];
          for(kk=0; kk<2;  kk++) W_8[ii][jj][kk]=0;
          for(kk=0; kk<4;  kk++) W_4[ii][jj][kk]=0;
          for(kk=0; kk<8;  kk++) W_2[ii][jj][kk]=0;
          for(kk=0; kk<16; kk++) W_1[ii][jj][kk]=0;
      }
   }

//=========    16 to 8  ===============

   
   for (ii=0; ii<2*NN; ii++) neigh_c[ii]=0;

   for (ii=0; ii<NN; ii++) {
      for (jj=0; jj<NN; jj++) {
          if(XX[ii][jj]==1){
              neigh[ii][neigh_c[ii]]=jj+NN;
              neigh_c[ii]++;
              neigh[jj+NN][neigh_c[jj+NN]]=ii;
              neigh_c[jj+NN]++;
          }
      }
   }


   d=NN*16;

   while(d>0){

         ii=0;
         while( neigh_c[ii] == 0) ii++;

         while(neigh_c[ii]!=0){    // use the last edge in the edge list

           
            x=0;
            while(neigh[ii][x]==-1 && x<16) x++;
            jj=neigh[ii][x];
            W_8[ii][jj-NN][0]=1;
            neigh[ii][x]=-1;    //remove the edge
            neigh_c[ii]--;
            x=0;
            while(neigh[jj][x]!=ii && x<16) x++;
            neigh[jj][x]=-1;
            neigh_c[jj]--;
            d--;

            x=0;
            while(neigh[jj][x]==-1 && x<16) x++;
            ii=neigh[jj][x];
            W_8[ii][jj-NN][1]=1;
            neigh[jj][x]=-1;    //remove the edge
            neigh_c[jj]--;
            x=0;
            while(neigh[ii][x]!=jj && x<16) x++;
            neigh[ii][x]=-1;
            neigh_c[ii]--;
            d--;
         }
    }



//=========    8 to 4  ===============
for(kkk=0; kkk<2; kkk++){
   d=NN*8;

   for (ii=0; ii<2*NN; ii++) neigh_c[ii]=0;

   for (ii=0; ii<NN; ii++) {
      for (jj=0; jj<NN; jj++) {
          if(W_8[ii][jj][kkk]==1){
              neigh[ii][neigh_c[ii]]=jj+NN;
              neigh_c[ii]++;
              neigh[jj+NN][neigh_c[jj+NN]]=ii;
              neigh_c[jj+NN]++;
          }
      }
   }


   while(d>0){

         ii=0;
         while( neigh_c[ii] == 0) ii++;

         while(neigh_c[ii]!=0){    // use the last edge in the edge list

           
            x=0;
            while(neigh[ii][x]==-1 && x<8) x++;
            jj=neigh[ii][x];
            W_4[ii][jj-NN][2*kkk]=1;
            neigh[ii][x]=-1;    //remove the edge
            neigh_c[ii]--;
            x=0;
            while(neigh[jj][x]!=ii && x<8) x++;
            neigh[jj][x]=-1;
            neigh_c[jj]--;
            d--;

            x=0;
            while(neigh[jj][x]==-1 && x<8) x++;
            ii=neigh[jj][x];
            W_4[ii][jj-NN][2*kkk+1]=1;
            neigh[jj][x]=-1;    //remove the edge
            neigh_c[jj]--;
            x=0;
            while(neigh[ii][x]!=jj && x<8) x++;
            neigh[ii][x]=-1;
            neigh_c[ii]--;
            d--;
         }
    }


}

//=========    4 to 2  ===============

for(kkk=0; kkk<4; kkk++){
   d=NN*4;


   for (ii=0; ii<2*NN; ii++) neigh_c[ii]=0;

   for (ii=0; ii<NN; ii++) {
      for (jj=0; jj<NN; jj++) {
          if(W_4[ii][jj][kkk]==1){
              neigh[ii][neigh_c[ii]]=jj+NN;
              neigh_c[ii]++;
              neigh[jj+NN][neigh_c[jj+NN]]=ii;
              neigh_c[jj+NN]++;
          }
      }
   }


   while(d>0){

         ii=0;
         while( neigh_c[ii] == 0) ii++;

         while(neigh_c[ii]!=0){    // use the last edge in the edge list

           
            x=0;
            while(neigh[ii][x]==-1 && x<4) x++;
            jj=neigh[ii][x];
            W_2[ii][jj-NN][2*kkk]=1;
            neigh[ii][x]=-1;    //remove the edge
            neigh_c[ii]--;
            x=0;
            while(neigh[jj][x]!=ii && x<4) x++;
            neigh[jj][x]=-1;
            neigh_c[jj]--;
            d--;

            x=0;
            while(neigh[jj][x]==-1 && x<4) x++;
            ii=neigh[jj][x];
            W_2[ii][jj-NN][2*kkk+1]=1;
            neigh[jj][x]=-1;    //remove the edge
            neigh_c[jj]--;
            x=0;
            while(neigh[ii][x]!=jj && x<4) x++;
            neigh[ii][x]=-1;
            neigh_c[ii]--;
            d--;
         }
    }



}

//=========    2 to 1  ===============

for(kkk=0; kkk<8; kkk++){
   d=NN*2;


   for (ii=0; ii<2*NN; ii++) neigh_c[ii]=0;

   for (ii=0; ii<NN; ii++) {
      for (jj=0; jj<NN; jj++) {
          if(W_2[ii][jj][kkk]==1){
              neigh[ii][neigh_c[ii]]=jj+NN;
              neigh_c[ii]++;
              neigh[jj+NN][neigh_c[jj+NN]]=ii;
              neigh_c[jj+NN]++;
          }
      }
   }


   while(d>0){

         ii=0;
         while( neigh_c[ii] == 0) ii++;

         while(neigh_c[ii]!=0){    // use the last edge in the edge list

           
            x=0;
            while(neigh[ii][x]==-1 && x<2) x++;
            jj=neigh[ii][x];
            W_1[ii][jj-NN][2*kkk]=1;
            neigh[ii][x]=-1;    //remove the edge
            neigh_c[ii]--;
            x=0;
            while(neigh[jj][x]!=ii && x<2) x++;
            neigh[jj][x]=-1;
            neigh_c[jj]--;
            d--;

            x=0;
            while(neigh[jj][x]==-1 && x<2) x++;
            ii=neigh[jj][x];
            W_1[ii][jj-NN][2*kkk+1]=1;
            neigh[jj][x]=-1;    //remove the edge
            neigh_c[jj]--;
            x=0;
            while(neigh[ii][x]!=jj && x<2) x++;
            neigh[ii][x]=-1;
            neigh_c[ii]--;
            d--;
         }
    }
}

   return 1;

}


int euler_partition_8(void){

   int ii,jj, kk,d,kkk,x;
   int neigh[2*640][32], neigh_c[2*640];

   for (ii=0; ii<NN; ii++) {
      D[ii]=8;
      for (jj=0; jj<NN; jj++) {
          W[ii][jj]=XX[ii][jj];
          for(kk=0; kk<4;  kk++) W_4[ii][jj][kk]=0;
          for(kk=0; kk<8;  kk++) W_2[ii][jj][kk]=0;
          for(kk=0; kk<16; kk++) W_1[ii][jj][kk]=0;
      }
   }

//=========    8 to 4  ===============
   d=NN*8;

   for (ii=0; ii<2*NN; ii++) neigh_c[ii]=0;

   for (ii=0; ii<NN; ii++) {
      for (jj=0; jj<NN; jj++) {
          if(W[ii][jj]==1){
              neigh[ii][neigh_c[ii]]=jj+NN;
              neigh_c[ii]++;
              neigh[jj+NN][neigh_c[jj+NN]]=ii;
              neigh_c[jj+NN]++;
          }
      }
   }


   while(d>0){

         ii=0;
         while( neigh_c[ii] == 0) ii++;

         while(neigh_c[ii]!=0){    // use the last edge in the edge list

           
            x=0;

            while(neigh[ii][x]==-1 && x<8) x++;
            jj=neigh[ii][x];
            W_4[ii][jj-NN][0]=1;
            neigh[ii][x]=-1;    //remove the edge
            neigh_c[ii]--;
            x=0;
            while(neigh[jj][x]!=ii && x<8) x++;
            neigh[jj][x]=-1;
            neigh_c[jj]--;
            d--;

            x=0;
            while(neigh[jj][x]==-1 && x<8) x++;
            ii=neigh[jj][x];
            W_4[ii][jj-NN][1]=1;
            neigh[jj][x]=-1;    //remove the edge
            neigh_c[jj]--;
            x=0;
            while(neigh[ii][x]!=jj && x<8) x++;
            neigh[ii][x]=-1;
            neigh_c[ii]--;
            d--;
         }
    }


//=========    4 to 2  ===============

for(kkk=0; kkk<2; kkk++){
   d=NN*4;

   for (ii=0; ii<2*NN; ii++) neigh_c[ii]=0;

   for (ii=0; ii<NN; ii++) {
      for (jj=0; jj<NN; jj++) {
          if(W_4[ii][jj][kkk]==1){
              neigh[ii][neigh_c[ii]]=jj+NN;
              neigh_c[ii]++;
              neigh[jj+NN][neigh_c[jj+NN]]=ii;
              neigh_c[jj+NN]++;
          }
      }
   }


   while(d>0){

         ii=0;
         while( neigh_c[ii] == 0) ii++;

         while(neigh_c[ii]!=0){    // use the last edge in the edge list

           
            x=0;
            while(neigh[ii][x]==-1 && x<4) x++;
            jj=neigh[ii][x];
            W_2[ii][jj-NN][2*kkk]=1;
            neigh[ii][x]=-1;    //remove the edge
            neigh_c[ii]--;
            x=0;
            while(neigh[jj][x]!=ii && x<4) x++;
            neigh[jj][x]=-1;
            neigh_c[jj]--;
            d--;

            x=0;
            while(neigh[jj][x]==-1 && x<4) x++;
            ii=neigh[jj][x];
            W_2[ii][jj-NN][2*kkk+1]=1;
            neigh[jj][x]=-1;    //remove the edge
            neigh_c[jj]--;
            x=0;
            while(neigh[ii][x]!=jj && x<4) x++;
            neigh[ii][x]=-1;
            neigh_c[ii]--;
            d--;
         }
    }



}

//=========    2 to 1  ===============

for(kkk=0; kkk<4; kkk++){
   d=NN*2;


   for (ii=0; ii<2*NN; ii++) neigh_c[ii]=0;

   for (ii=0; ii<NN; ii++) {
      for (jj=0; jj<NN; jj++) {
          if(W_2[ii][jj][kkk]==1){
              neigh[ii][neigh_c[ii]]=jj+NN;
              neigh_c[ii]++;
              neigh[jj+NN][neigh_c[jj+NN]]=ii;
              neigh_c[jj+NN]++;
          }
      }
   }


   while(d>0){

         ii=0;
         while( neigh_c[ii] == 0) ii++;

         while(neigh_c[ii]!=0){    // use the last edge in the edge list

           
            x=0;
            while(neigh[ii][x]==-1 && x<2) x++;
            jj=neigh[ii][x];
            W_1[ii][jj-NN][2*kkk]=1;
            neigh[ii][x]=-1;    //remove the edge
            neigh_c[ii]--;
            x=0;
            while(neigh[jj][x]!=ii && x<2) x++;
            neigh[jj][x]=-1;
            neigh_c[jj]--;
            d--;


            x=0;
            while(neigh[jj][x]==-1 && x<2) x++;
            ii=neigh[jj][x];
            W_1[ii][jj-NN][2*kkk+1]=1;
            neigh[jj][x]=-1;    //remove the edge
            neigh_c[jj]--;
            x=0;
            while(neigh[ii][x]!=jj && x<2) x++;
            neigh[ii][x]=-1;
            neigh_c[ii]--;
            d--;
         }
    }
}

   return 1;

}

int euler_partition_4(void){

   int ii,jj, kk,d,kkk,x;
   int neigh[2*640][32], neigh_c[2*640];

   for (ii=0; ii<NN; ii++) {
      D[ii]=4;
      for (jj=0; jj<NN; jj++) {
          W[ii][jj]=XX[ii][jj];
          for(kk=0; kk<8;  kk++) W_2[ii][jj][kk]=0;
          for(kk=0; kk<16; kk++) W_1[ii][jj][kk]=0;
      }
   }

//=========    4 to 2  ===============

   d=NN*4;

   for (ii=0; ii<2*NN; ii++) neigh_c[ii]=0;

   for (ii=0; ii<NN; ii++) {
      for (jj=0; jj<NN; jj++) {
          if(W[ii][jj]==1){
              neigh[ii][neigh_c[ii]]=jj+NN;
              neigh_c[ii]++;
              neigh[jj+NN][neigh_c[jj+NN]]=ii;
              neigh_c[jj+NN]++;
          }
      }
   }


   while(d>0){

         ii=0;
         while( neigh_c[ii] == 0) ii++;

         while(neigh_c[ii]!=0){    // use the last edge in the edge list

           
            x=0;
            while(neigh[ii][x]==-1 && x<4) x++;
            jj=neigh[ii][x];
            W_2[ii][jj-NN][0]=1;
            neigh[ii][x]=-1;    //remove the edge
            neigh_c[ii]--;
            x=0;
            while(neigh[jj][x]!=ii && x<4) x++;
            neigh[jj][x]=-1;
            neigh_c[jj]--;
            d--;

            x=0;
            while(neigh[jj][x]==-1 && x<4) x++;
            ii=neigh[jj][x];
            W_2[ii][jj-NN][1]=1;
            neigh[jj][x]=-1;    //remove the edge
            neigh_c[jj]--;
            x=0;
            while(neigh[ii][x]!=jj && x<4) x++;
            neigh[ii][x]=-1;
            neigh_c[ii]--;
            d--;
         }
    }

//=========    2 to 1  ===============

for(kkk=0; kkk<2; kkk++){
   d=NN*2;


   for (ii=0; ii<2*NN; ii++) neigh_c[ii]=0;

   for (ii=0; ii<NN; ii++) {
      for (jj=0; jj<NN; jj++) {
          if(W_2[ii][jj][kkk]==1){
              neigh[ii][neigh_c[ii]]=jj+NN;
              neigh_c[ii]++;
              neigh[jj+NN][neigh_c[jj+NN]]=ii;
              neigh_c[jj+NN]++;
          }
      }
   }


   while(d>0){

         ii=0;
         while( neigh_c[ii] == 0) ii++;

         while(neigh_c[ii]!=0){    // use the last edge in the edge list

           
            x=0;
            while(neigh[ii][x]==-1 && x<2) x++;
            jj=neigh[ii][x];
            W_1[ii][jj-NN][2*kkk]=1;
            neigh[ii][x]=-1;    //remove the edge
            neigh_c[ii]--;
            x=0;
            while(neigh[jj][x]!=ii && x<2) x++;
            neigh[jj][x]=-1;
            neigh_c[jj]--;
            d--;


            x=0;
            while(neigh[jj][x]==-1 && x<2) x++;
            ii=neigh[jj][x];
            W_1[ii][jj-NN][2*kkk+1]=1;
            neigh[jj][x]=-1;    //remove the edge
            neigh_c[jj]--;
            x=0;
            while(neigh[ii][x]!=jj && x<2) x++;
            neigh[ii][x]=-1;
            neigh_c[ii]--;
            d--;
         }
    }
}

   return 1;

}





int euler_partition_2(void){

   int ii,jj, kk,d,kkk,x;
   int neigh[2*640][32], neigh_c[2*640];

   for (ii=0; ii<NN; ii++) {
      D[ii]=2;
      for (jj=0; jj<NN; jj++) {
          W[ii][jj]=XX[ii][jj];
          for(kk=0; kk<8;  kk++) W_2[ii][jj][kk]=0;
          for(kk=0; kk<16; kk++) W_1[ii][jj][kk]=0;
      }
   }

//=========    2 to 1  ===============

   d=NN*2;

   for (ii=0; ii<2*NN; ii++) neigh_c[ii]=0;

   for (ii=0; ii<NN; ii++) {
      for (jj=0; jj<NN; jj++) {
          if(W[ii][jj]==1){
              neigh[ii][neigh_c[ii]]=jj+NN;
              neigh_c[ii]++;
              neigh[jj+NN][neigh_c[jj+NN]]=ii;
              neigh_c[jj+NN]++;
          }
      }
   }


   while(d>0){

         ii=0;
         while( neigh_c[ii] == 0) ii++;

         while(neigh_c[ii]!=0){    // use the last edge in the edge list

           
            x=0;
            while(neigh[ii][x]==-1 && x<2) x++;
            jj=neigh[ii][x];
            W_1[ii][jj-NN][0]=1;
            neigh[ii][x]=-1;    //remove the edge
            neigh_c[ii]--;
            x=0;
            while(neigh[jj][x]!=ii && x<2) x++;
            neigh[jj][x]=-1;
            neigh_c[jj]--;
            d--;


            x=0;
            while(neigh[jj][x]==-1 && x<2) x++;
            ii=neigh[jj][x];
            W_1[ii][jj-NN][1]=1;
            neigh[jj][x]=-1;    //remove the edge
            neigh_c[jj]--;
            x=0;
            while(neigh[ii][x]!=jj && x<2) x++;
            neigh[ii][x]=-1;
            neigh_c[ii]--;
            d--;
         }
    }

   return 1;

}

//Ford-Fulkerson Algorithm
int max_flow (int source, int sink) {
    int i,j,u;
    // Initialize empty flow.
    int max_flow = 0;
    for (i=0; i<n; i++) {
	for (j=0; j<n; j++) {
	    flow[i][j] = 0;
	}
    }
    // While there exists an augmenting path,
    // increment the flow along this path.
    while (bfs(source,sink)) {
        // Determine the amount by which we can increment the flow.
	int increment = oo;
	for (u=n-1; pred[u]>=0; u=pred[u]) {
	    increment = min(increment,capacity[pred[u]][u]-flow[pred[u]][u]);
	}
        // Now increment the flow.
	for (u=n-1; pred[u]>=0; u=pred[u]) {
	    flow[pred[u]][u] += increment;
	    flow[u][pred[u]] -= increment;
	}
	max_flow += increment;
    }
    // No augmenting path anymore. We are done.
    return max_flow;
}

int new_bfs (int start, int target) {
    int i,j,aa,bb,num_nodes, special;
    int u,v;

    for (u=0; u<n; u++) {
	color[u] = WHITE;
    }   
    head = tail = 0;
    enqueue(start);
    pred[start] = -1;
    while (head!=tail) {
	u = dequeue();
        // Search all adjacent white nodes v. If the capacity
        // from u to v in the residual network is positive,
        // enqueue v.
        num_nodes=(n-2)/2;
        if(u==0){
           aa=1;
           bb=num_nodes+1;
           special=0;
        }
        else if(u>=1 && u<=num_nodes){
           special=1;
        }
        else if(u>=num_nodes+1 && u<num_nodes+num_nodes+1){
           special=2;
        }
        else{
           aa=0;
           bb=n;
           special=0;
        }
           
        if(special == 0){
	   for (v=aa; v<bb; v++) {
	       if (color[v]==WHITE && capacity[u][v]-flow[u][v]>0) {
	        	enqueue(v);
		        pred[v] = u;
	       }
	   }
           
        }
        else if(special == 1){
	   for (i=0; i<group; i++) {
               v=neighbour[u][i];
	       if (color[v]==WHITE && capacity[u][v]-flow[u][v]>0) {
	        	enqueue(v);
		        pred[v] = u;
	       }
	   }
           special=0;
        }
        else if(special == 2){
	   for (i=0; i<group; i++) {
               v=neighbour[u][i];
	       if (color[v]==WHITE && capacity[u][v]-flow[u][v]>0) {
	        	enqueue(v);
		        pred[v] = u;
	       }
	   }


           v=n-1;
	   if (color[v]==WHITE && capacity[u][v]-flow[u][v]>0) {
	       	enqueue(v);
	        pred[v] = u;
	   }


           special=0;
        }
        else{
           printf("error\n");
        }
 
    }
    // If the color of the target node is black now,
    // it means that we reached it.
    return color[target]==BLACK;
}

int new_max_flow (int source, int sink) {
    int i,j,u;
    // Initialize empty flow.
    int max_flow = 0;
    for (i=0; i<n; i++) {
	for (j=0; j<n; j++) {
	    flow[i][j] = 0;
	}
    }
    // While there exists an augmenting path,
    // increment the flow along this path.
    while (new_bfs(source,sink)) {
        // Determine the amount by which we can increment the flow.
	int increment = oo;
	for (u=n-1; pred[u]>=0; u=pred[u]) {
	    increment = min(increment,capacity[pred[u]][u]-flow[pred[u]][u]);
	}
        // Now increment the flow.
	for (u=n-1; pred[u]>=0; u=pred[u]) {
	    flow[pred[u]][u] += increment;
	    flow[u][pred[u]] -= increment;
	}
	max_flow += increment;
    }
    // No augmenting path anymore. We are done.
    return max_flow;
}