//*************************************************************
//	resolvendo assignments com o método húngaro
//	gilberto miranda jr. - 2006
//*************************************************************
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "hungarian.h"
#include "hungarian.c"
#include "subproNEW01.h"
#include "assfunc.c"

double trunc5(double n){
  unsigned long int c = floor(n*100);
  double r = c/100.0;
  return r;
}

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

  FILE *fp1;
  int h, i, j, k, l, m, n, n2, n3, n4, le, *d, iter;
  double **a, *u, *v, *t, *w, lb, lbb, ub, ubb, **x, 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));
  d = (int *)malloc(n4*sizeof(int));
  x = (double **)malloc(n*sizeof(double*));
  for(i=0; i<n; i++) x[i] = (double *)malloc(n*sizeof(double)); 
  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 = 1.e6;

  // solucao inicial:
  lb = assign_a(n, a, x);
  lbb = lb;
  printf(" Limite Inferior Inicial: %10.5f \n\n", lbb);
  cf = lb;
  cv = solve_slave(n, d, x, v, u);
  ub = cf + cv;
  if (ub<ubb) ubb = ub;
  printf(" Limite Superior Inicial: %10.5f \n\n", ubb);
  gap = 100*(ubb-lbb)/ubb;
  printf(" Gap-otimalidade Inicial: %10.5f \n\n", gap);
  h = 1; 
// loop de "benders" :  
//  while((ubb - lbb/ubb)<0.01){

   w = (double*)malloc(n2*h*sizeof(double));
   y = (double*)malloc(n2*h*sizeof(double));
   t = (double*)malloc(n2*sizeof(double));
   // inicializando vetor de precos:
    for (i=0; i<(n2*h); i++) w[i] = 0.0;                 
    for (i=0; i<(n2*h); i++) y[i] = 0.0;                   
    lambda = 0.1;

    // disparando método do subgradiente :
    iter = 0;
    ynorma = 1;              
    while(ynorma > 1.e-12){                    
     //inicializando matriz c :
      for(k=0; k<n; k++){
       for(i=0; i<n; i++){
        c[k][i] = a[k][i];       
       }
      }
     // subproblema lagrangeano em x :
     // montando matriz c:
     for(m=0;m<h;m++){            
	    for(k=0;k<n;k++){
	       for(l=(k+1);l<n;l++){
                    
           for(j=0;j<n;j++){
	       c[l][j] += w[n2*h*m+k*n+l]*v[n2*h*m+k*n2+l*n+j];
	       }
	       for(i=0;i<n;i++){
	       c[k][i] -= w[n2*h*m+k*n+l]*u[n2*h*m+k*n2+l*n+i];
	       }
 
           } 
	    }   
     }   
     for(k=0; k<n; k++){
        for(i=0; i<n; i++){      
          c[k][i] = trunc5(c[k][i]);       
          printf(" %15.8f ",c[k][i]);           
        }
        printf("\n");
     }
     printf("\n");
     splx = assign_a(n, c, x);
       // subproblema lagrangeano em t :
         splt = 0.0;
           for (i=0; i<n2; i++) t[i] = 0.0;
             for(k=0; k<n; k++){
                  for(l=(k+1); l<n; l++){
                    wh = 0.0;
                    for(j=0; j<h; j++) wh += w[n2*h*j+k*n+l];
                      if((1.0-wh) > -1.e-12){ splt += 0.0; t[k*n+l] = 0.0; }
                      else { splt += -1.e6; t[k*n+l] = (1.e6/(1.0-wh)); }
             }
           }
       // calculando lowerbound lagrangeano
       lb = splx + splt;  
       // determinando subgradiente :
       for(j=0; j<h; j++){
        //printf(" Imprimido subgradiente: \n");
        for(k=0; k<n; k++){
          for(l=(k+1); l<n; l++){
            y[n2*h*j+k*n+l] = 0.0;
            for(i=0; i<n; i++){
            y[n2*h*j+k*n+l]+= v[n2*h*j+k*n2+l*n+i]*x[l][i];
            //printf("v: %10.5f --- x: %10.5f \n", v[n2*h*j+l*n2+k*n+i], x[l][i]);
            }
            for(i=0; i<n; i++){
            y[n2*h*j+k*n+l]-= u[n2*h*j+k*n2+l*n+i]*x[k][i];
            //printf("u: %10.5f --- x: %10.5f \n", u[n2*h*j+k*n2+l*n+i], x[k][i]);
            }
            y[n2*h*j+k*n+l] -= t[k*n+l];        
            //printf("y: %10.5f \n", y[n2*h*j+k*n+l]);                 
            }
          }
        }
        // determinando norma-2 do subgradiente :                
        ynorma2 = 0.0;
        for(j=0; j<h; j++){
         for(k=0; k<n; k++){
          for(l=(k+1); l<n; l++){
            ynorma2 += y[n2*h*j+k*n+l] * y[n2*h*j+k*n+l];
          }
         }
        }
        ynorma = sqrt(ynorma2);                
        // determinando passo do subgradiente :
        stp = lambda*(ubb - lb)/ynorma2;
        //atualizando variáveis duais :
        for(j=0; j<h; j++){
          printf(" Imprimindo variaveis duais: \n");       
           for(k=0; k<n; k++){
            for(l=(k+1); l<n; l++){            
             w[n2*h*j+k*n+l] = w[n2*h*j+k*n+l] + stp*y[n2*h*j+k*n+l];
             if(w[n2*h*j+k*n+l]<0.0) w[n2*h*j+k*n+l]=0.0;
             printf(" %5.3f ", w[n2*h*j+k*n+l]);
            }
            printf("\n");
           }
        }
        iter ++;
       // impressao intermediária :
       printf("  Iteracao: %d --- Lowerbound: %10.5f --- NormaSG: %10.5f \n\n\n", iter, lb, ynorma);             
  
    } // end while : fim do método do subgradiente

                       
//  } //end while  
  return 0;
}

