// usage bestcuts <Data file name> <Info file name> <Output File Name> 

//    Info File format : 
//  eff.boundary 1 ... eff.boundary 4
//  eff step 1 ... eff step 4         { these define points e/p curve density }
//  # of coords 
//  # of lines
//  # of signal events
//  suppression (double) 
//  var. 1 lower end of upper cut {  total of  } 
//  var. 1 upper end of upper cut {  4 lines   } 
//  var. 1 upper cut initial step {  for each  } 
//  var. 1 upper cut final step   {  cut       } 
//  the same for var.1 lower cut  {  total of 2 cuts on each variable  }
//  the same for all vars.
//  for the program not to make steps on certain cut - set final step
//  greater than initial and initial step greater than boundaries of the cut.


#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>
#include "linklist.h" 
#include <stdio.h>
#include <unistd.h>

// global variables

struct  best
{
    double   efficiency; 
    double   purity; 
    double   *vector; 
    best () {
       vector = new double[100]; 
    }    
}   bests[5000];           // maximal number of points on e/p graph 5000


int     ncoord, nlines, ntotal, SoughtIdt;   // basic run parameters

double  ulower[100] ;      // maximal # of variables = 100 
double  llower[100] ;      // min lower boundary limit
double  uupper[100] ;      // max upper boundary limit 
double  lupper[100] ;      // min upper boundary limit 

double  lstep[100] ;       // lower boundary step 
double  ustep[100] ;       // upper boundary step
double  lminstep[100];     // lower boundary min. step
double  uminstep[100];     // upper ----- // -----  
double  cup[100], clow[100]; 

double  **Data ;           // the array for holding data

double  mark[4]; 
double  purstep[3];       

int     d = 0 ;            // number of points on e/p graph
double  suppression;       // noise suppression

ofstream   OutFile; 


// functions 

void          ReadData ( char* DataName ) ; 
inline void   TheLoop ( int c ) ;  // c - recursion counter 
void          panic (char*) ; 


inline INT (double val)
{
    return    val < 0 ? (int) (val - 0.5) : int (val + 0.5) ;
}





int main (int argc, char* argv[])
{
    char* DataName = new char[100]; 
    char* InfoName = new char[100];
    char* OutName = new char[100];
    char *com     = new char[100];

    ifstream  InfoFile; 

    unsigned  int c = 0 ;     // counter of recursion

    if (argc != 4) panic ("bestcuts <Data file> <Info file> <Out file> "); 
    DataName = argv[1];      // name of data file
    InfoName = argv[2];      // name of info file 
    OutName = argv[3];       // name of output file    

    InfoFile.open (InfoName, ios::nocreate ) ; 

// reading run parameters from info file.
   
    InfoFile >> mark[0] >> mark[1] >> mark[2] >> mark[3]; 
    InfoFile >> purstep[0] >> purstep[1] >> purstep[2]; 
    
    int     cnt1; 
    double  dnt2;
  
    for (cnt1 = 0; cnt1 <= 3; cnt1++)
    { 
      for (dnt2 = mark[cnt1] ; dnt2 <= mark[cnt1+1] ; dnt2 += purstep[cnt1])
      {
            bests[d].purity = dnt2 ; 
            bests[d++].efficiency = 0 ; 
      }
    }  
    d-- ;  

    InfoFile >> ncoord ;       // number of coords in data file, not incl. idt
    if (ncoord > 100)        panic ("Too many coordinates. Maximal = 100");  
    InfoFile >> nlines ;       // number of lines in DataFile
    if (ncoord*nlines > 5000000)  panic ("Too big sample"); // max. RAM
    InfoFile >> ntotal;        // total number of sought particles
    InfoFile >> SoughtIdt;     // sought idt
    InfoFile >> suppression;   // suppression 


    for (cnt1 = 0; cnt1 < ncoord; cnt1++ )  // setting borders 
    {
          InfoFile >> llower[cnt1];
          InfoFile >> ulower[cnt1]; 
          InfoFile >> lstep[cnt1];
          InfoFile >> lminstep[cnt1]; 
          InfoFile >> lupper[cnt1];   
          InfoFile >> uupper[cnt1]; 
          InfoFile >> ustep[cnt1];  
          InfoFile >> uminstep[cnt1]; 
    } 

    ReadData (DataName);         // setting the array Data
    OutFile.open (OutName) ; 
    TheLoop (0);       
    for (cnt1 = 0 ; cnt1 < d; cnt1++ ) 
    {   OutFile << bests[cnt1].efficiency << " " << bests[cnt1].purity 
               << endl; 
       cout << bests[cnt1].efficiency << " " << bests[cnt1].purity 
               << endl;}
    OutFile.close (); 

    while (1) {
          for (cnt1 = 0 ; cnt1 < ncoord ; cnt1++ ) // setting new step
          { 
             if (ustep[cnt1] > uminstep[cnt1]*1.01)  ustep[cnt1] /= 2; 
             if (lstep[cnt1] > lstep[cnt1]*1.01)     lstep[cnt1] /= 2; 
          }   
          for (cnt1 = 0 ; cnt1 < d ; cnt1++ )  // improving all the points
       	  {     			       // with new step
             if (bests[cnt1].efficiency < 0.001)       continue; 
             if (ustep[cnt1] >= 0.99*uminstep[cnt1]) { 
                 ulower[cnt1] = bests[cnt1].vector[cnt1*2] + 2*lstep[cnt1];
                 llower[cnt1] = bests[cnt1].vector[cnt1*2] - 2*lstep[cnt1];
             }
             if (lstep[cnt1] >= 0.99*lminstep[cnt1]) {
                 uupper[cnt1] = bests[cnt1].vector[cnt1*2+1] + 2*ustep[cnt1];
                 lupper[cnt1] = bests[cnt1].vector[cnt1*2+1] - 2*ustep[cnt1];
             } 
             TheLoop(0); 
   	  }
          OutFile.open (OutName) ;
          for ( cnt1 = 0 ; cnt1 < d ; cnt1++ )  
              OutFile << bests[cnt1].efficiency << " "  
		      << bests[cnt1].purity << endl;
          OutFile.close (); 
          cout << "file created" << endl; 
          sprintf (com, "cp file1 dummy"); 
          system (com); 
    } 
    return 1; 
}

// the main loop 


inline void TheLoop ( int c )  // c - recursion counter 
{
   unsigned   int   ntrue = 0, nfalse = 0 ; 
   double     efficiency, purity ;    
// loops over upper and lower cut on each variable

   int cnt1, cnt2; 
   Bool match; 

   for ( cup[c] = lupper[c] ; cup[c] <= uupper[c] ; cup[c] += ustep[c] ) 
   { 
       for ( clow[c] = llower[c]; clow[c] <= ulower[c]; clow[c] += lstep[c] )  
       {
          if ( c < ncoord-1 )  TheLoop (c+1) ;  // next recursion step 
          // else 
          for ( cnt1 = 0 ; cnt1 < nlines ; cnt1++ ) // loop over points
          {
              match = TRUE;                // does a point match cuts ?  
              for ( cnt2 = 0; cnt2 < ncoord; cnt2++) //loop over coords
              {
   	           if ( Data[cnt1][cnt2] > cup[cnt2] || 
    		        Data[cnt1][cnt2] < clow [cnt2]   )    match = FALSE;  
              }
              if (  match && abs(INT(Data[cnt1][ncoord]))
                    == SoughtIdt)                          ntrue++; 
              else if (match)                              nfalse++;  
           }
           // ------------------------ writing output
           efficiency = ntrue / (double) ntotal;
           purity = ntrue / ((double) ntrue + suppression*nfalse) ;
           for (cnt1 = 0 ; purity > bests[cnt1].purity && cnt1 < d ; cnt1++);  
           if (efficiency > bests[cnt1].efficiency) {    
                  bests[cnt1].efficiency = efficiency;  
                  for ( cnt2 = 0 ; cnt2 < ncoord ; cnt2++ ) {
                      bests[cnt1].vector[cnt2] = clow[cnt2];
                      bests[cnt1].vector[cnt2+1] = cup[cnt2]; 
		  }
	   }
           ntrue = nfalse = 0; // wrongly and truly found # of pts set 0
       }
    }
} 



// reading data file into memory.

void     ReadData (char* DataName)
{
      ifstream DataFile(DataName); 
   
      if (! DataFile) panic ("Unable to open data file"); 

      Data = new double* [nlines];
      if (!Data) panic ("memory allocation error"); 

      for ( unsigned long cnt1 = 0 ; cnt1 < nlines; cnt1++ ) 
      {
          Data[cnt1] = new double [ncoord+1];

          for (int cnt2 = 0; cnt2 <= ncoord; cnt2++ ) {
               DataFile >> Data[cnt1][cnt2];
               if ( DataFile.eof() )
                  panic ( "Unexpected data EOF in FindNeighbours" );
          }
      }

      DataFile.close (); 
} 

void  panic (char *phrase)
{
    cout << phrase << endl ;
    exit (0); 
}