//////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// Questa macro per ROOT (http://root.cern.ch) decide l'accoppiamento
// ottimale tra moduli usando il metodo di "simulated annealing".
//
// Andrea Giammanco, 11/6/2005
//
// (Anche il file di input e' online: http://cern.ch/andrea.giammanco/ds/modulevdepl.txt)
// (Un esempio di file di output generato da questo programma: http://cern.ch/andrea.giammanco/ds/out.txt)
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

//gROOT->Reset();

const Int_t nSteps=3000; // numero iterazioni per l'annealing
const Double_t TI=1000.0; // "temperatura" iniziale
const Double_t TF=1.0; // "temperatura" finale
const Double_t K=1.0; // analogo della K di Boltzmann (non serve...)

//void ds(Int_t strategy=1, Char_t input1[40]="moduli_stereolr_25_6_2005.txt", Char_t input2[40]="moduli_rphi_25_6_2005.txt", Char_t out[40]="out.txt", Double_t cut=0.3, Int_t verbosity=0)
void ds(Int_t strategy=1,  Char_t out[40]="out.txt", Double_t cut=0.2, Int_t verbosity=0)
{
  /*
  sceglibassi(180,80.,"moduli_rphi_25_6_2005.txt","moduli_rphi_180.txt",0);
  sceglibassi(90,80.,"moduli_stereol_25_6_2005.txt","moduli_stereol_90.txt",0);
  sceglibassi(90,80.,"moduli_stereor_25_6_2005.txt","moduli_stereor_90.txt",0);
  */
  //  unisci("moduli_stereol_v230_4_7_2005.txt","moduli_stereor_v230_4_7_2005.txt","moduli_stereolr_v230_4_7_2005.txt",0);
  //  run(strategy,"moduli_rphi_4_7_2005.txt","moduli_stereolr_v230_4_7_2005.txt",out,cut,verbosity);
  Int_t n1=conta("moduli_stereolr_12_7_2005.txt", verbosity);
  Int_t n2=conta("moduli_rphi_12_7_2005.txt", verbosity);
  Double_t media1=media(n1,"moduli_stereolr_12_7_2005.txt", verbosity);
  Double_t media2=media(n2,"moduli_rphi_12_7_2005.txt", verbosity);
  Int_t nmax=min(n1,n2);
  /*
  if (n1>n2) { // il n.1 dev'essere il piu' piccolo
    swap( input1, input2 );
    swap( media1, media2 );
    swap( n1, n2 );
  }
  Char_t middle1[40];
  Char_t middle2[40];
  Char_t middle2bis[40];
  sprintf(middle1,"%s_reordered",input1);
  riordina(nmax,input1,middle1,verbosity);
  sprintf(middle2,"%s_closer",input2);
  */
  sceglivicini(nmax,media1,"moduli_rphi_12_7_2005.txt","moduli_rphi_temp.txt");
  //  sprintf(middle2bis,"%s_reordered",input2);
  riordina(nmax,"moduli_rphi_temp.txt","moduli_rphi_reord.txt",verbosity);
  riordina(nmax,"moduli_stereolr_12_7_2005.txt","moduli_stereolr_reord.txt",verbosity);
  run(strategy,"moduli_rphi_reord.txt","moduli_stereolr_reord.txt",out,cut,verbosity);
}


Int_t conta(Char_t input[40]="moduli_stereolr_25_6_2005.txt", Int_t verbosity=0, Double_t cut=0.2)
{
  Int_t n=0;
  TString name[2000];
  Char_t temporary[14];
  Double_t y[2000];
  ifstream in;
  in.open(input);
  do {
    in >> temporary >> y[n];
    name[n]=temporary;
    if (!in.good()) break;
    if (verbosity==1) cout << "n = " << n << ", x = " << name[n] << ", y = " << y[n] << endl;
    n++;
  } while (1);
  if (verbosity==1) cout << "Numero moduli: " << n << endl;
  
  Double_t mean=0;
  for (Int_t i=0;i<n;i++) {
    mean+=y[i];
  }
  mean=mean/n;
  if (verbosity==1) cout << "Media: " << mean << endl;

  Double_t sigma=0;
  for (Int_t i=0;i<n;i++) {
    sigma+=(y[i]-mean)**2;
  }
  sigma=sqrt(sigma/n);
  Double_t r=100.*sigma/mean;
  if (verbosity==1) cout << "Sigma: " << sigma << " (" << r << " %)" << endl;

  if (verbosity==1) cout << endl << "I seguenti moduli si discostano del " << cut*100 << " % dalla media:" << endl;
  for (Int_t i=0;i<n;i++) {
    if (fabs(y[i]-mean)/mean>cut) {
      if (verbosity==1) cout << name[i] << " V = " << y[i] << endl;
    }
  }

  return n;
}

Double_t media(Int_t n, Char_t input[40]="moduli_stereolr_25_6_2005.txt", Int_t verbosity=0)
{
  TString name[2000];
  Char_t temporary[14];
  Double_t y[2000];
  ifstream in;
  in.open(input);
  for (Int_t i=0;i<n;i++) {
    in >> temporary >> y[i];
    name[i]=temporary;
    if (verbosity==1) cout << "i = " << i << ", x = " << name[i] << ", y = " << y[i] << endl;
  }
  
  Double_t mean=0;
  for (Int_t i=0;i<n;i++) {
    mean+=y[i];
  }
  mean=mean/n;
  if (verbosity==1) cout << "Media: " << mean << endl;

  return mean;
}

void riordina(Int_t n, Char_t input[40]="moduli_stereolr_25_6_2005.txt",Char_t out[40]="moduli_stereolr_reordered.txt", Int_t verbosity=0)
{
  TString name[2000];
  Char_t temporary[14];
  Double_t y[2000];
  ifstream in;
  in.open(input);
  for (Int_t i=0;i<n;i++) {
    in >> temporary >> y[i];
    name[i]=temporary;
    if (verbosity==1) cout << "i = " << i << ", x = " << name[i] << ", y = " << y[i] << endl;
  }
 
    // uso l'algoritmo "bubblesort" (http://www.aihorizon.com/essays/basiccs/lists/sorting/bubble.htm) per ordinare in Vdepl decrescente
  for ( int i = 0; i < n-1; ++i ) {
      for ( int j = 1; j < n-i; ++j ) {
	if ( y[j-1] < y[j] ) {
	  swap( y[j-1], y[j] );
	  swap( name[j-1], name[j] );
	}
      }
    }

  // output su file:
  ofstream output(out);
  if (! output) {
    cout << "Error opening output file" << endl;
  }
  for (Int_t i=0;i<n;i++) {
    output << name[i] << " " << y[i] << endl;
  }

  output.close();
  cout << "File riordinato: " << out << endl;

}

void unisci(Char_t input1[40]="moduli_stereol_25_6_2005.txt", Char_t input2[40]="moduli_stereor_25_6_2005.txt", Char_t out[40]="moduli_stereolr.txt", Int_t verbosity=1)
{
  TString name1[2000];
  TString name2[2000];
  Char_t temporary[14];
  Double_t y1[2000];
  Double_t y2[2000];
  ifstream in1;
  in1.open(input1);
  Int_t n1=0;
  do {
    in1 >> temporary >> y1[n1];
    name1[n1]=temporary;
    if (!in1.good()) break;
    if (verbosity==1) cout << "n1 = " << n1 << ", x1 = " << name1[n1] << ", y1 = " << y1[n1] << endl;
    n1++;
  } while (1);
  if (verbosity==1) cout << "Numero moduli: " << n1 << endl;
  ifstream in2;
  in2.open(input2);
  Int_t n2=0;
  do {
    in2 >> temporary >> y2[n2];
    name2[n2]=temporary;
    if (!in2.good()) break;
    if (verbosity==1) cout << "n2 = " << n2 << ", x2 = " << name2[n2] << ", y2 = " << y2[n2] << endl;
    n2++;
  } while (1);
  if (verbosity==1) cout << "Numero moduli: " << n2 << endl;


  // output su file:
  ofstream output(out);
  if (! output) {
    cout << "Error opening output file" << endl;
  }
  for (Int_t i=0;i<n1;i++) {
    output << name1[i] << " " << y1[i] << endl;
  }
  for (Int_t i=0;i<n2;i++) {
    output << name2[i] << " " << y2[i] << endl;
  }

  output.close();
  cout << "File somma: " << out << endl;

}

void sceglibassi(Int_t nmax, Double_t safetycut=80., Char_t input[40]="moduli_stereolr_25_6_2005.txt",Char_t out[40]="moduli_stereolr_reordered.txt", Int_t verbosity=1)
{
  Int_t n=0;
  TString name[2000];
  Char_t temporary[14];
  Double_t y[2000];
  ifstream in;
  in.open(input);
  do {
    in >> temporary >> y[n];
    name[n]=temporary;
    if (!in.good()) break;
    if (verbosity==1) cout << "n = " << n << ", x = " << name[n] << ", y = " << y[n] << endl;
    n++;
  } while (1);
  if (verbosity==1) cout << "Numero moduli: " << n << endl;

  if (n<nmax) {nmax=n;}

    // uso l'algoritmo "bubblesort" (http://www.aihorizon.com/essays/basiccs/lists/sorting/bubble.htm) per ordinare in Vdepl crescente
  for ( int i = 0; i < n-1; ++i ) {
      for ( int j = 1; j < n-i; ++j ) {
	if ( y[j-1] > y[j] ) {
	  swap( y[j-1], y[j] );
	  swap( name[j-1], name[j] );
	}
      }
    }


  // output su file:
  ofstream output(out);
  if (! output) {
    cout << "Error opening output file" << endl;
  }
  for (Int_t i=0;i<nmax;i++) {
    if (y[i]>safetycut) {
      output << name[i] << " " << y[i] << endl;
    } else {
      cout << name[i] << " " << y[i] << ": sotto il taglio Vdepl > " << safetycut << endl;
      nmax++;
    }
  }

  output.close();
  cout << "File riordinato: " << out << endl;

}


void sceglivicini(Int_t n1, Double_t valore, Char_t input[40]="moduli_stereolr_25_6_2005.txt",Char_t out[40]="moduli_stereolr_closer.txt", Int_t verbosity=0)
{
  Int_t n2=0;
  TString name[2000];
  Char_t temporary[14];
  Double_t y[2000];
  ifstream in;
  in.open(input);
  do {
    in >> temporary >> y[n2];
    name[n2]=temporary;
    if (!in.good()) break;
    if (verbosity==1) cout << "n = " << n2 << ", x = " << name[n2] << ", y = " << y[n2] << endl;
    n2++;
  } while (1);
  if (verbosity==1) cout << "Numero moduli: " << n2 << endl;

  // selezione dei piu' vicini al valore scelto:
  // uso l'algoritmo "bubblesort" (http://www.aihorizon.com/essays/basiccs/lists/sorting/bubble.htm) per ordinare in |Vdepl-valore|
  for ( int i = 0; i < n2-1; ++i ) {
      for ( int j = 1; j < n2-i; ++j ) {
	if ( fabs(y[j-1]-valore) > fabs(y[j]-valore) ) {
	  swap( y[j-1], y[j] );
	  swap( name[j-1], name[j] );
	}
      }
    }

  // output su file:
  ofstream output(out);
  if (! output) {
    cout << "Error opening output file" << endl;
  }
  for (Int_t i=0;i<n1;i++) {
    output << name[i] << " " << y[i] << endl;
  }

  output.close();
  cout << "File dopo selezione: " << out << endl;

}


void run(Int_t strategy=1, Char_t input1[40]="moduli_stereolr_25_6_2005.txt", Char_t input2[40]="moduli_rphi_25_6_2005.txt", Char_t out[40]="out.txt", Double_t cut=0.2, Int_t verbosity=0)
{
  const Int_t nModules=2000;
  Double_t x1[nModules];
  Double_t x2[nModules];
  TString name1[nModules];
  TString name2[nModules];
  Char_t temporary[14];

  ifstream in1;
  in1.open(input1);
  Int_t n1=0;
  do {
    in1 >> temporary >> x1[n1];
    name1[n1]=temporary;
    if (!in1.good()) break;
    if (verbosity==1) cout << "n1 = " << n1 << ", x1 = " << name1[n1] << ", V = " << x1[n1] << endl;
    n1++;
  } while (1);
  cout << "Numero moduli file 1: " << n1 << endl;
  ifstream in2;
  in2.open(input2);
  Int_t n2=0;
  do {
    in2 >> temporary >> x2[n2];
    name2[n2]=temporary;
    if (!in2.good()) break;
    if (verbosity==1) cout << "n2 = " << n2 << ", x2 = " << name2[n2] << ", V = " << x2[n2] << endl;
    n2++;
  } while (1);
  cout << "Numero moduli file 2: " << n2 << endl;

  Int_t nmax=min(n1,n2);

  // array delle distanze tra coppie:
  Double_t dpair[nModules][nModules];
  for (Int_t i=0;i<nmax;i++) {
    for (Int_t j=0;j<nmax;j++) {
      //      dpair[i][j]=fabs(x1[i]-x2[j]);
      //      dpair[i][j]=fabs(x1[i]-x2[j])/(x1[i]+x2[j]);
      dpair[i][j]=fabs(x1[i]-x2[j])/min(x1[i],x2[j]);
      if (i<4 && j<4) cout << "i = " << i << ", j = " << j << ", dpair = " << dpair[i][j] << endl;
    }
  }


  // ricerca di un buon punto di partenza per la minimizzazione:
  Int_t jmin[nModules];
  if (strategy==1) {
    // parte dal piu' vicino (in Vdepl) al primo modulo, e cosi' via
    // (NOTA: lento!!!)
    Double_t min[nModules];
    for (Int_t i=0;i<nmax;i++) {
      min[i]=9999;
    }
    for (Int_t i=0;i<nmax;i++) {
      for (Int_t j=1;j<nmax;j++) {
	bool utilizzabile=true;
	if (j!=i) {
	  for (Int_t m=0;m<i;m++) { // evitare di usare indici gia' assegnati
	    if (j==jmin[m]) { 
	      utilizzabile=false; 
	      m=i;
	    }
	  }
	  if (utilizzabile && dpair[i][j]<min[i]) {
	    min[i]=dpair[i][j];
	    jmin[i]=j;
	  }
	}
      }
    }
  } else if (strategy==2) {
    // uso l'algoritmo "bubblesort" (http://www.aihorizon.com/essays/basiccs/lists/sorting/bubble.htm) per ordinare in Vdepl entrambi gli insiemi di moduli, e poi faccio l'associazione 1 a 1
    for ( int i = 0; i < nmax-1; ++i ) {
      for ( int j = 1; j < nmax-i; ++j ) {
	if ( x1[j-1] < x1[j] ) {
	  swap( x1[j-1], x1[j] );
	  swap( name1[j-1], name1[j] );
	}
	if ( x2[j-1] < x2[j] ) {
	  swap( x2[j-1], x2[j] );
	  swap( name2[j-1], name2[j] );
	}
      }
    }
    for (Int_t i=0;i<nmax;i++) {
      jmin[i]=i;
      cout << name1[i] << " V = " << x1[i] << " - " << name2[i] << " V = " << x2[i] << endl;
    }
  } else if (strategy==3) {
    // come strategy==1, ma prima ordino gli array (usando bubblesort) non per Vdepl, ma per distanza dalla media di Vdepl sull'unione dei due insiemi

  } else if (strategy==0) {
    // un array di partenza qualsiasi, che potrebbe in generale essere troppo lontano da qualunque minimo per dare una convergenza sensata
    for (Int_t i=0;i<nmax;i++) {
      jmin[i]=nmax-i-1;
    }
  }

  // usare questo minimo come punto di partenza per la ricerca di altri minimi
  cout << endl << "jmin = ";
  for (Int_t i=0;i<nmax;i++) { cout << jmin[i] << " "; }
  cout << endl;

  // "simulated annealing":
  Int_t jcurr[nModules]; // soluzione corrente
  Int_t jnew[nModules]; // nuova soluzione
  Double_t dmin[nModules]; // array di distanze 1-jmin[1],2-jmin[2],...,N-jmin[N]
  Double_t dcurr[nModules];
  Double_t dnew[nModules];
  if (verbosity==1) cout << "dmin = ";
  for (Int_t i=0;i<nmax;i++) {
    dmin[i]=dpair[i][jmin[i]];
    if (verbosity==1) cout << dmin[i] << " ";
  }
  Double_t vmin=sommadistanze(dmin,nmax);
  if (verbosity==1) cout << ", somma: " << vmin << endl << endl;
  Double_t vcurr=vmin;
  Double_t vnew=vmin;
  copiaarray(jmin,jcurr,nmax);
  copiaarray(jmin,jnew,nmax);
  Double_t Tt=TI; // temperatura all'istante t
  for (Int_t t=0; t<nSteps; t++) { // mega-ciclo di annealing
    //copio la sol. corrente nella sol. nuova
    copiaarray(jcurr,jnew,nmax);
    //perturbo la nuova (evitando di cadere in 0,1,2,...,N)
    perturba(jnew,nmax);
    if (verbosity==1) cout << "jnew = ";
    for (Int_t i=0;i<nmax;i++) { if (verbosity==1) cout << jnew[i] << " "; }
    if (verbosity==1) cout << endl;
    //calcolo la somma della sol. nuova
    if (verbosity==1) cout << "dnew = ";
    vnew=sommadistanze(dnew,nmax);
    if (verbosity==1) cout << ", somma: " << vnew << endl;
    //se la nuova e' migliore dell'ottimo (e tutti i moduli sono piu' vicini del vincolo), aggiorno l'ottimo
    if (vnew<vmin) {
      vmin=vnew;
      copiaarray(jnew,jmin,nmax);
      if (verbosity==1) cout << " *** nuovo minimo ***" << endl;
    }
    //se la nuova e' migliore della corrente, oppure e' peggiore ma la routine casuale tengo() decide di tenerlo, aggiorno la corrente
    if (vnew<vcurr || tengo(vnew,vcurr,Tt)) {
      vcurr=vnew;
      copiaarray(jnew,jcurr,nmax);
    }
    //raffreddo il sistema
    Tt=raffredda(t);
    if (verbosity==1) cout << ", temperatura: " << Tt << endl << endl;;
  }

  // conclusione:
  cout << "---------------" << endl;
  cout << "Miglior soluzione trovata:" << endl;
  cout << endl << "jmin = ";
  for (Int_t i=0;i<nmax;i++) { cout << jmin[i] << " "; }
  cout << endl;
  cout << "dmin = ";
  for (Int_t i=0;i<nmax;i++) {
    dmin[i]=dpair[i][jmin[i]];
    cout << dmin[i] << " ";
  }
  Double_t vmin=sommadistanze(dmin,nmax);
  cout << ", somma: " << vmin << endl << endl;

  // output su file:
  ofstream output(out);
  if (! output) {
    cout << "Error opening output file" << endl;
  }
  Int_t counter=0;
  for (Int_t i=0;i<nmax;i++) {
    output << "Il modulo n." << name1[i] << " va accoppiato con il n." << name2[jmin[i]];
    output << "  (" << x1[i] << " - " << x2[jmin[i]] << ")";
    //    Double_t discrep=2*fabs(x1[i]-x2[jmin[i]])/(x1[i]+x2[jmin[i]]);
    Double_t discrep=fabs(x1[i]-x2[jmin[i]])/min(x1[i],x2[jmin[i]]);
    if (discrep>cut) {
      output << "  *";
      counter++;
    }
    output << endl;
  }
  output << "L'asterisco significa accoppiamento troppo distante in V_depl (piu` del " << 100*cut << "%)" << endl;
  output << "Questo capita il " << 100*ratio(counter,nmax) << " % delle volte." << endl;

  output.close();
  cout << "File di output con gli accoppiamenti: " << out << endl;

}


Double_t raffredda(Int_t t)
{
  return TI*pow(TF/TI,(double)t / (double)nSteps);
}

bool tengo(Double_t nuovo, Double_t vecchio, Double_t temp)
{
  Double_t e = exp((nuovo-vecchio)/(K*temp));
  Double_t r= gRandom->Rndm(1);
  return r<e;
}

void perturba(Int_t *iarray, Int_t n)
{
  //estrae due posizioni a caso:
  Int_t a,b;
  do {
    a = random(n);
    b = random(n);
  }  while(a==b);
  //swap delle due posizioni:
  swap(iarray[a],iarray[b]);
}

inline int swap( int &x, int &y )
{
  int z = x;
  x = y;
  y = z;
}
inline int swap( double &x, double &y )
{
  double z = x;
  x = y;
  y = z;
}
inline int swap( TString &x, TString &y )
{
  TString z = x;
  x = y;
  y = z;
}

Int_t random(Int_t nmax) // estrae un numero intero tra 0 e nmax-1
{
  Float_t r;
  r= gRandom->Rndm(1)*nmax;
  Int_t i=r;
  //  cout << i << endl;
  return i;
}

Double_t sommadistanze(Double_t *darray, Int_t n) // la quantita' da minimizzare nell'annealing
{
  Double_t sum=0.;
  for (Int_t i=0; i<n; i++) {
    sum+=darray[i];
  }
  return sum/n;
}

void copiaarray(Int_t *iarray1, Int_t *iarray2, Int_t n)
{
  for (Int_t i=0;i<n;i++) { iarray2[i]=iarray1[i]; }
}


Double_t min(Int_t a, Int_t b)
{
  Double_t res=a;
  if (b<a) res=b;

  return res;
}

void bubblesort( int array[], double x[], int n )
{
  for ( int i = 0; i < n-1; ++i )
    for ( int j = 1; j < n-i; ++j )
      if ( x[j-1] < x[j] ) swap( array[j-1], array[j] );
}

Double_t ratio(Int_t n1,Int_t n2)
{
  Double_t x1=1.*n1;
  Double_t x2=1.*n2;
  Double_t r;
  if (x2!=0){
    r=x1/x2;
  } else {
    cout << "**** dividing by zero!!!!!!!" << endl;
    r=-99999.;
  }
  return r;
}