//*************************************************************
//	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[]) {

  int h, i, j, k, l, m, n, n2, n3, n4, matrix_size, le, *d, iter,  **xb, **x;
  double **a, *u, *uu, *uut, *v, *vv, *vvt, *t, *w, lb, lbb, ub, ubb, cf, cv, gap, **c, *y ;
  double splx, splt, stp, lambda, ynorma, ynorma2, wh;
  float lef;

  ifstream data(argv[1]);

  data>>n;;

  n2 = n*n;
  n3 = n2*n;
  n4 = n3*n;
  
//  alocacao dinamica : 
  a = new double* [n];
  for(i=0; i<n; i++) a[i] = new double [n]; 
  v = new double [n3];
  u = new double [n3]; 
  uu = new double [n3];
  vv = new double [n3];
  d = new int [n4];
  x = new int* [n];
  for(i=0; i<n; i++) x[i] = new int [n]; 
  xb = new int* [n];
  for(i=0; i<n; i++) xb[i] = new int [n]; 
  c = new double* [n];
  for(i=0; i<n; i++) c[i] = new double [n];

//  leitura de dados :
  for(i=0; i<n; i++){
   for(j=0; j<n; j++){
     data>>a[i][j];
   }
  }

for(k=0; k<n; k++){
 for(i=0; i<n; i++){
  for(l=0; l<n; l++){
   for(j=0; j<n; j++){
    data>>d[k*n3+i*n2+l*n+j];
   }
  }
 }
}

//  inicializando upperbound :
  ub = ubb = 1.e6;

  // solucao inicial:
  lb = splx_lap(n, a, x);
   printf("\n");
   for(k=0; k<n; k++){
        for(i=0; i<n; i++){
	     printf(" %d ", x[k][i]);
        }
	printf("\n");
   }
   printf("\n");

  lbb = lb;
  printf(" Limite Inferior Inicial: %10.5f \n\n", lbb);
  cf = lb;
  cv = solve_slave(n, d, x, v, u);
   printf(" u: --- v: \n");
   for(k=0;k<n;k++){
       for(l=0;l<n;l++){
           for(i=0;i<n;i++){
		printf(" %10.5f     %10.5f \n",u[k*n2+l*n+i], v[k*n2+l*n+i] );
	   }
       }
   }

  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; 
   for(k=0;k<n;k++){
       for(l=0;l<n;l++){
           for(i=0;i<n;i++){
            uu[n3*(h-1)+k*n2+l*n+i] = u[k*n2+l*n+i];
	    vv[n3*(h-1)+k*n2+l*n+i] = v[k*n2+l*n+i];
	   }
       }
   }

// loop de "benders" :  
  //while((ubb - lbb/ubb)<0.01){
 while(h<17){

   w = new double [n2*h];
   y = new double [n2*h];
   t = new double [n2];

   // 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;
    int iternmelhora = 0;
    ynorma = 1;
    while(ynorma > 1.e-12){
    //while(iter<=8850){
     //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*m+k*n+l]*vv[n3*m+k*n2+l*n+j];
	        }
	        for(i=0;i<n;i++){
	         c[k][i] -= w[n2*m+k*n+l]*uu[n3*m+k*n2+l*n+i];
	        }
 
               }
	    }
     }
     splx = splx_lap(n, c, x);
     cv = solve_slave(n, d, x, v, u);
     cf = 0.0;
     for(k=0; k<n; k++){
        for(i=0; i<n; i++){
             cf += a[k][i]*x[k][i]; 
         }
     }
     printf("\n");
     ub = cf + cv;
     if (ub<ubb) ubb = ub;
     //printf(" Limite Inferior : %10.5f \n\n", lbb);
     //printf(" Limite Superior : %10.5f \n\n", ub);

       // 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*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; 
       int imprimir = 0; 
       if (lb>lbb) 
	{ 
	   lbb = lb; 
	   imprimir = 1;
	   iternmelhora = 0;
	   for(k=0; k<n; k++){
              for(i=0; i<n; i++){
               xb[k][i] = x[k][i];
              }
           }
	}
	else iternmelhora++;
       // determinando subgradiente :
       for(j=0; j<h; j++){
        for(k=0; k<n; k++){
          for(l=(k+1); l<n; l++){
            y[n2*j+k*n+l] = 0.0;
            for(i=0; i<n; i++){
            y[n2*j+k*n+l]+= vv[n3*j+k*n2+l*n+i]*x[l][i];
            }
            for(i=0; i<n; i++){
            y[n2*j+k*n+l]-= uu[n3*j+k*n2+l*n+i]*x[k][i];
            }
            y[n2*j+k*n+l] -= t[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*j+k*n+l] * y[n2*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++){
           for(k=0; k<n; k++){
            for(l=(k+1); l<n; l++){
             w[n2*j+k*n+l] = w[n2*j+k*n+l] + stp*y[n2*j+k*n+l];
             if(w[n2*j+k*n+l]<0.0) w[n2*j+k*n+l]=0.0;
            }
           }
        }
        iter ++;
       // impressao intermediária :
       //if(imprimir && (iter%100)==0) printf("  Iteracao: %d --- Lowerbound: %10.5f --- NormaSG: %10.5f \n", iter, lbb, ynorma);
       printf("  Iteracao: %d --- Lowerbound: %10.5f --- Lambda: %15.10f \n", iter, lb, lambda);
      if (iternmelhora==300) {
           lambda = 0.1*lambda;
           //cout<<"diminuindo passo : "<<iter<<" "<<lambda<<"\n";
           iternmelhora = 0;
           //printf("  Iteracao: %d --- Lowerbound: %10.5f --- NormaSG: %10.5f \n", iter, lbb, ynorma);
       }
      if(lambda < 1.e-5) break;
    } // end while : fim do método do subgradiente

      /*cout<<"Imprimindo variáveis duais:\n";
        for(j=0; j<h; j++){
           for(k=0; k<n; k++){
            for(l=(k+1); l<n; l++){
             cout<<" "<<w[n2*j+k*n+l];
            }
	    cout<<"\n";
           }

        }*/

    h++;
    uut = new double [n3*h];
    vvt = new double [n3*h];

	// preenchendo uut e vvt:
	  for(j=0; j<(h-1); j++){
   	    for(k=0;k<n;k++) {
	       for(l=0;l<n;l++){
                for(i=0;i<n;i++){
                uut[n3*j+k*n2+l*n+i] = uu[n3*j+k*n2+l*n+i];
		vvt[n3*j+k*n2+l*n+i] = vv[n3*j+k*n2+l*n+i];
		}
               }
            }
	  }
	delete[] uu;
	delete[] vv;

	//aumentando uu e vv :
	uu = new double [n3*h];
	vv = new double [n3*h];

	// preenchendo uu e vv aumentados :
	  for(j=0; j<(h-1); j++){
   	    for(k=0;k<n;k++) {
	       for(l=0;l<n;l++){
                for(i=0;i<n;i++){
                uu[n3*j+k*n2+l*n+i] = uut[n3*j+k*n2+l*n+i];
		vv[n3*j+k*n2+l*n+i] = vvt[n3*j+k*n2+l*n+i];
		}
               }
            }
	  }
        //
	delete[] uut;
	delete[] vvt;

   cv = solve_slave(n, d, xb, v, u);
   printf(" u: --- v: \n");
   for(k=0;k<n;k++){
       for(l=0;l<n;l++){
           for(i=0;i<n;i++){
		printf(" %10.5f     %10.5f \n",u[k*n2+l*n+i], v[k*n2+l*n+i] );
	   }
       }
   }

   cf = 0.0;
   printf("\n");
   for(k=0; k<n; k++){
        for(i=0; i<n; i++){
             cf += a[k][i]*xb[k][i]; 
	     printf(" %d ", xb[k][i]);
        }
	printf("\n");
   }
   printf("\n");
   ub = cf + cv;
   if (ub<ubb) ubb = ub;
   printf("\n\n Iteracao Benders: %d \n\n", h);
   printf(" Limite Inferior : %10.5f \n\n", lbb);
   printf(" Limite Superior : %10.5f \n\n", ubb);
   printf(" Limite Superior : %10.5f \n\n", ub);
   gap = 100*(ubb-lbb)/ubb;
   printf(" Gap-otimalidade  : %10.5f \n\n", gap);
   // atualizando uu e vv :

   for(k=0;k<n;k++){
       for(l=0;l<n;l++){
           for(i=0;i<n;i++){
            uu[n3*(h-1)+k*n2+l*n+i] = u[k*n2+l*n+i];
	    vv[n3*(h-1)+k*n2+l*n+i] = v[k*n2+l*n+i];
	   }
       }
   }
   /*printf(" uu: --- vv: \n");
	  for(j=0; j<h; j++){
   	    for(k=0;k<n;k++) {
	       for(l=0;l<n;l++){
                for(i=0;i<n;i++){
		printf(" %10.5f     %10.5f \n",uu[n3*j+k*n2+l*n+i], vv[n3*j+k*n2+l*n+i] );
		}
               }
            }
	    printf("\n Novo h! : \n");
	  }
*/
 getchar();
 } //end while: fim do método de Benders
  return 0;
}

