//*************************************************************
//	resolvendo assignments com o método húngaro
//	gilberto miranda jr. - 2006
//*************************************************************
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include "lap.h"
#include "lap.cpp"
#include "splx_lap.cpp"
#include "subproNEW01.h"

using namespace std;

typedef int boolean;

int main(int argc, char *argv[]) {

  FILE *fp1;
  int h, i, j, k, l, m, n, n2, n3, n4, le, *d, iter, *f, **x;
  double **a, *u, *v, *w, *r, *z, lb, lbb, ub, ubb, cf, cv, gap, **c, *y;
  double splx, splt, stp, lambda, ynorma, ynorma2, wh;
  float lef;
  fp1 = fopen(argv[1], "r");

  fscanf(fp1, "%d", &n);

  n2 = n*n;
  n3 = n2*n;
  n4 = n3*n;
  
//  alocacao dinamica : 
  a = (double **)malloc(n*sizeof(double*));
  for(i=0; i<n; i++) a[i] = (double *)malloc(n*sizeof(double)); 
  v = (double *)malloc(n3*sizeof(double));
  u = (double *)malloc(n3*sizeof(double));
  r = (double *)malloc(n3*sizeof(double));
  d = (int *)malloc(n4*sizeof(int));
  f = (int *)malloc(n4*sizeof(int));
  x = (int **)malloc(n*sizeof(int*));
  for(i=0; i<n; i++) x[i] = (int *)malloc(n*sizeof(int)); 
  c = (double **)malloc(n*sizeof(double*));
  for(i=0; i<n; i++) c[i] = (double *)malloc(n*sizeof(double)); 

//  leitura de dados :
  for(i=0; i<n; i++){
   for(j=0; j<n; j++){
     fscanf(fp1, "%f", &lef);
     a[i][j] = lef;
   }
  }

for(k=0; k<n; k++){
 for(i=0; i<n; i++){
  for(l=0; l<n; l++){
   for(j=0; j<n; j++){
    fscanf(fp1, "%d", &le);
    d[k*n3+i*n2+l*n+j] = le;
   }
  }
 }
}

  // inicializando upperbound :
  ub = ubb = 10000.0;

  // solucao inicial:
  lb = splx_lap(n, a, x);
  lbb = lb;
  int iternmelhora =0;

ubb = 0;
for(k=0; k<n; k++){
 for(i=0; i<n; i++){
  ubb+= a[k][i]*x[k][i];
  for(l=(k+1); l<n; l++){
   for(j=0; j<n; j++){
    ubb += (d[k*n3+i*n2+l*n+j]+d[l*n3+j*n2+k*n+i])*x[l][j]*x[k][i];
   }
  }

 }
}


  printf(" Limite Inferior Inicial: %10.5f \n\n", lbb);
  printf(" Limite Superior Inicial: %10.5f \n\n", ubb);
  gap = 100*(ubb-lbb)/ubb;
  printf(" Gap-otimalidade Inicial: %10.5f \n\n", gap);
 
   w = (double*)malloc(n3*sizeof(double));
   y = (double*)malloc(n3*sizeof(double));
   z = (double*)malloc(n3*sizeof(double));

   // inicializando vetor de precos:

    for (i=0; i<n3; i++) w[i] = 0.0;
    for (i=0; i<n3; i++) y[i] = 0.0;
    for (i=0; i<n3; i++) z[i] = 0.0;
    for (i=0; i<n3; i++) v[i] = 0.0;
    for (i=0; i<n3; i++) u[i] = 0.0;
    for (i=0; i<n3; i++) r[i] = 0.0;
    for (i=0; i<n4; i++) f[i] = 0;
    lambda = 0.0001;

    // disparando método do subgradiente :
    iter = 0;
    ynorma = 1;
    //while(ynorma > 1.e-12){
      while(1){
     //inicializando matriz c :
      for(k=0; k<n; k++){
       for(i=0; i<n; i++){
        c[k][i] = a[k][i];
               for(j=0;j<n;j++){
	       c[k][i] -= r[k*n2+i*n+j];
	       }

	       for(j=0;j<n;j++){
	       c[k][i] -= r[k*n2+j*n+i];
	       }
       }
      }
     // subproblema lagrangeano em x :
     // montando matriz c:
     	    for(k=0;k<n;k++){
	       for(l=(k+1);l<n;l++){

               for(j=0;j<n;j++){
	       c[l][j] -= v[k*n2+l*n+j];
	       }

	       for(i=0;i<n;i++){
	       c[k][i] -= u[k*n2+l*n+i];
	       }
 
               } 
	   }   
 
       splx = splx_lap(n, c, x);
       // subproblema lagrangeano em f :
         splt = 0.0;
	for(k=0; k<n; k++){
 		for(i=0; i<n; i++){
		  for(l=(k+1); l<n; l++){
   			for(j=0; j<n; j++){
			if(i!=j){
			 double del = d[k*n3+i*n2+l*n+j]+d[l*n3+j*n2+k*n+i]+v[k*n2+l*n+j]+u[k*n2+l*n+i]+r[k*n2+i*n+j]+r[l*n2+j*n+i];
			if( del>0.0) {splt+= 0.0; f[k*n3+i*n2+l*n+j]=0;}
			else {splt-= del; f[k*n3+i*n2+l*n+j]=1;}
			}
			}
  		  }
 		}
	}

       // calculando lowerbound lagrangeano
       lb = splx + splt;  

       // determinando subgradientes :

         for(k=0; k<n; k++){
          for(l=(k+1); l<n; l++){
           for(j=0; j<n; j++){
            y[k*n2+l*n+j] = 0;
	    for(i=0; i<n; i++){
             if(i!=j)y[k*n2+l*n+j] += f[k*n3+i*n2+l*n+j];
            }
            y[k*n2+l*n+j] -= x[l][j];
           }
          }
         }

         for(k=0; k<n; k++){
          for(l=(k+1); l<n; l++){
           for(i=0; i<n; i++){
            w[k*n2+l*n+i] = 0;
	    for(j=0; j<n; j++){
             if(i!=j)w[k*n2+l*n+i] += f[k*n3+i*n2+l*n+j];
            }
            w[k*n2+l*n+i] -= x[k][i];
           }
          }
         }

         for(k=0; k<n; k++){
          for(i=0; i<n; i++){
           for(j=0; j<n; j++){
           if(i!=j){
            z[k*n2+i*n+j] = 0;
	    for(l=(k+1); l<n; l++) z[k*n2+i*n+j] += f[k*n3+i*n2+l*n+j];
 	    for(l=0; l<k; l++) z[l*n2+j*n+i] += f[l*n3+j*n2+k*n+i];
            z[k*n2+i*n+j] -= x[k][i];
           }
           }
          }
         }
   

       // determinando norma-2 do subgradiente :                
        ynorma2 = 0.0;
        for(j=0; j<n; j++){
         for(k=0; k<n; k++){
          for(l=(k+1); l<n; l++){
            ynorma2 += y[k*n2+l*n+j] * y[k*n2+l*n+j];
          }
         }
        }
        double wnorma2 = 0.0;
        for(j=0; j<n; j++){
         for(k=0; k<n; k++){
          for(l=(k+1); l<n; l++){
            wnorma2 += w[k*n2+l*n+j] * w[k*n2+l*n+j];
          }
         }
        }
        double znorma2 = 0.0;
        for(k=0; k<n; k++){
         for(i=0; i<n; i++){
          for(j=0; j<n; j++){
            znorma2 += z[k*n2+i*n+j] * z[k*n2+i*n+j];
          }
         }
        }

        ynorma = sqrt(ynorma2);                
        
// determinando passo do subgradiente :
        double stpy = lambda*(ubb - lb)/ynorma2;
        double stpw = lambda*(ubb - lb)/wnorma2;
        double stpz = lambda*(ubb - lb)/znorma2;
        //atualizando variáveis duais :
        
           for(k=0; k<n; k++){
            for(l=(k+1); l<n; l++){
		for(j=0; j<n; j++){            
                  v[k*n2+l*n+j] = v[k*n2+l*n+j] + stpy*y[k*n2+l*n+j];
	          //if(v[k*n2+l*n+j] <0.0)v[k*n2+l*n+j] = 0.0;
                }
            }
           }

        
           for(k=0; k<n; k++){
            for(l=(k+1); l<n; l++){
	     for(j=0; j<n; j++){
             u[k*n2+l*n+j] = u[k*n2+l*n+j] + stpw*w[k*n2+l*n+j];
	     //if(u[k*n2+l*n+j] <0.0)u[k*n2+l*n+j] = 0.0;	
             }
            }
           }

    
           for(k=0; k<n; k++){
            for(i=0; i<n; i++){
	     for(j=0; j<n; j++){           
             if(i!=j)r[k*n2+i*n+j] = r[k*n2+i*n+j] + stpz*z[k*n2+i*n+j];
	     //if(r[k*n2+i*n+j] <0.0)r[k*n2+i*n+j] = 0.0;
             }
            }
           }
        

        iter ++;
	printf(" Iteracao: %d --- Lowerbound: %10.5f --- \n\n\n", iter, lb);
       // impressao intermediária :
       if(lb >lbb){
        lbb = lb;
	printf(" Melhorou!!! Iteracao: %d --- Lowerbound: %10.5f --- Upperbound: %10.5f --- Lambda: %10.8f \n\n\n", iter, lbb, ubb,lambda);
        iternmelhora = 0;
        }
        else iternmelhora++;
        if(iternmelhora == 150000){
	 iternmelhora  = 0;
         lambda = 0.5 * lambda;
	printf(" Reduziu!!! Iteracao: %d --- Lowerbound: %10.5f --- Upperbound: %10.5f --- Lambda: %10.8f \n\n\n", iter, lbb, ubb,lambda);
        }
        if(lambda<1.e-8) break; 
ub = 0;
for(k=0; k<n; k++){
 for(i=0; i<n; i++){
  ub+= a[k][i]*x[k][i];
  for(l=(k+1); l<n; l++){
   for(j=0; j<n; j++){
    ub += (d[k*n3+i*n2+l*n+j]+d[l*n3+j*n2+k*n+i])*x[l][j]*x[k][i];
   }
  }
 }
}
if(ub<ubb) ubb = ub;

    } // end while : fim do método do subgradiente


  return 0;
}

