// ------------------------------------------------------------------
//  linkmeth.c - methods operating on linked lists
// ------------------------------------------------------------------

#include "skiplist.h"
#include "neighbours.h"
#include <unistd.h>
#include <sys/wait.h>
#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>

// ------------------------------------------------------------------
//  1.  For ordered lists
// ------------------------------------------------------------------


// FINDING BY KEY 

template <class T> Node <T> *LinkedList<T>::Find (T val, double TheKey) 
{
    Node <T> *current_node = first[0]->next;

    while ( TheKey > current_node->key && current_node != last[0] ) {
            if ( ! current_node ) Panic () ;
       	    current_node = current_node->next;
    }    

    if (current_node->key != TheKey) return NULL;
    else                             return current_node;
}



// INSERTION BY VALUE

template <class T>
void LinkedList<T>::Insert (T val)
{
   Node <T> *insert_node;

   if (chlen && length == MaxLength && val > MaxValue) return;

   insert_node = new Node<T> (val);
   current = first[0]->next; 
  
   while ( val > current->value && current != last[0] ) {
           if ( ! current )    
                Panic ("Null pointer in LinkedList::Insert(<T>)"); 
           current = current->next;
   }

   insert_node->prev = current->prev;
   insert_node->next = current;
   current->prev->next = insert_node;
   current->prev = insert_node;
   
   if (chlen) {  
        if (length == MaxLength)  delete last[0]->prev;  
        else                      length++;  
        MaxValue = last[0]->prev->value; 
   } 
   return ; 
}


// INSERTION BY KEY

template <class T>
void  LinkedList<T>::Insert (T val, double TheKey)
{
   Node <T>  *insert_node;

   if (chlen && length == MaxLength && TheKey > MaxKey)  return;

// going down, while searchig the place to insert

   current = first[nlevels-1]->next;
   for (int levelcount = nlevels-1 ; levelcount >= 0 ; levelcount-- ) {
      while (TheKey > current->key && current != last[levelcount]) {
            if ( ! current )    Panic ("No node in list"); 
	    current = current->next;
      }
      if (current->lower)  current = current->prev->lower->next; 
   }

// going up, if there is a need to insert a node there.
 
   Node <T> *prev_node = current->prev;
   Node <T> *next_node = current->prev->next; 
   Node <T> *lower_node = NULL;
 
   insert_node  = new Node <T> (val, TheKey);    

   for (int level = 0 ; ; level++ ) {     

// inserting new node ; 
       insert_node->lower = lower_node; 
       if (lower_node) lower_node->upper = insert_node; 
       insert_node->next = next_node;
       insert_node->prev = prev_node; 
       prev_node->next = insert_node;
       next_node->prev = insert_node;   
       insert_node->level = level;  

       if (    !prev_node->upper && !next_node->upper  ) { 
           lower_node = insert_node ; 
           next_node = next_node->next->upper ;  
           prev_node = prev_node->prev->upper ; 
           insert_node = new Node <T> (insert_node->value, insert_node->key); 
           continue; 
       }

       else if ( !next_node->upper && next_node->next && 
                 !next_node->next->upper                   ) { 
           insert_node = new Node <T> (next_node->value, next_node->key); 
           lower_node = next_node ; 
           next_node = next_node->next->next->upper;
           prev_node = prev_node->upper; 
           continue; 
       } 
  
       else if ( !prev_node->upper && prev_node->prev && 
                 !prev_node->prev->upper                   ) {
           insert_node = new Node <T> (prev_node->value, prev_node->key);
           lower_node = prev_node; 
           next_node = next_node->upper; 
           prev_node = prev_node->prev->prev->upper; 
           continue; 
       } 

       else { 
           if ( nlevels < level )  nlevels++ ;       
           break;       
       } 
   } 
    
// deleting the last node if MaxLenngth number of entries reached
 
   if (chlen) {
               if (length == MaxLength)   delete last[0]->prev; 
               else                       length++; 
               MaxKey = last[0]->prev->key; 
   } 

   return ; 
}

// DELETION BY KEY 

template <class T>
Bool LinkedList<T>::Delete (T dummyvalue, double TheKey)
{
    Node <T> *current_node = Find (dummyvalue, TheKey);

    if (current_node == NULL) return FALSE;

    ~current_node();
    if (chlen)  {
         MaxKey = last[0]->prev->key; 
         length--; 
    } 
    return TRUE; 
}




// ------------------------------------------------------------------
//  2.  For UNordered lists
// ------------------------------------------------------------------

template <class T>
Node <T> *LinkedList<T>::Find (T value)     // goes through all the nodes
{
    Node <T> *current_node = first[0]->next;

    while (value != current_node->value && current_node != last[0]) {
            if ( ! current_node ) Panic ("No node in LinkedList::Find(..)") ;
       	    current_node = current_node->next;
    }    

    if (current_node->value != value) return NULL;
    else                             return current_node;
}


template <class T>
void LinkedList<T>::UnOrdInsert (T val)
{
   Node <T> *insert_node;
   
   if (chlen && length = MaxLength) return;  // if the length is already max
   if (!last[0] || !last[0]->prev ) 
              Panic ("No last node in LinkedList::Find(..)");

   insert_node = new Node <T> (val); 
   insert_node->prev = last[0]->prev;
   insert_node->next = last[0];
   last[0]->prev->next = insert_node;
   last[0]->prev = insert_node;
   current = last[0]->prev; 

   return ; 
}

template <class T>
Bool   LinkedList<T>::Delete (T value)
{
    Node <T> *current_node = Find (value);

    if (current_node == NULL) return FALSE;

    delete current_node;
    length--; 
    if (chlen) {
          MaxKey = last[0]->prev->key; 
          length--; 
    } 
    return TRUE; 
}






// ------------------------------------------------------------------
//  2.  METHODS for both types of lists)
//      1. Bool Isempty
//      2. void Panic
//      3. void Start // 3-6 motin methods
//      4. void End
//      5. T * NextValue
//      6. T * PrevValue
// ------------------------------------------------------------------


template <class T>
Bool LinkedList<T>::Isempty (void)
{
     if  ( first[0]->next == last[0] )  return  TRUE;
     else                         return  FALSE;
}

template <class T>
void  LinkedList<T>::Panic ( const char* str)
{ 
     cout << str << endl ; 
     exit (0) ;
}

template <class T>
Bool  LinkedList<T>::Start (void)
{
     if ( Isempty () ) return FALSE; 
     current = first[0]->next;
     return TRUE ; 
}


template <class T>
Bool  LinkedList<T>::End (void)
{
     if ( Isempty () ) return FALSE; 
     current = last[0]->prev;
     return TRUE;     
}

template <class T>
Bool  LinkedList<T>::Forward (void)  
{
     if ( Isempty () || current->next == last[0] ) return FALSE ; 
     current = current->next; 
     return TRUE; 
}

template <class T>
Bool  LinkedList<T>::Backward (void)
{
    if ( Isempty || current->prev == first[0] ) return FALSE;
    current = current->prev ; 
    return TRUE;
}

template <class T>
T  LinkedList<T>::GetValue (void)
{   
     return current->value; 
}


template <class T>
double  LinkedList<T>::GetKey (void)
{
     return current->key; 
}

template <class T>
void LinkedList<T>::Emptify (void)
{
   for ( current = first[0]->next;
         current != last[0];
         current = current->next)  
      if (current->prev != first[0] )              delete current->prev;

   if (current->prev != first[0] && current->prev) delete current->prev;

   current = first[0];
   length = 0; nlevels = 1; 
}
//---------------------------------------------------------------------
// SRL algorithm, working version 1.02. 29.12.1997
// feature - forks out m processes, each SRL-ing 10000 pts of data. Waits 
// till they finish. Then forks another m processes, if sample is not 
// finished yet
// (c) Eugene Prosso. WIS, 1997-98
// Scientific supervisor and permanent source of good advices
// Ehud Duchovni, WIS
//---------------------------------------------------------------------

//---------------------------------------------------------------------
// input - input file with format
// "data" file name
// "MC" file name
// output file name 
// id # of sought particle
// flag (1 or 0) 
// number of different numbers of NN
// # of nn_1
// ....
// # of nn_last 
//_____________________________________________________________________



void ReadInput ( char* InputFileName,
                 char *DataName, char* MCName, 
                 char *OutName,  int &nnns, int nns[], int &flag, 
                 int &idt  ) ;

void CheckInput ( char* DataName, char* MCName,
                  long &ndata, long &nMC, int &d ) ; 

void ReadArrays ( char* DataName, char* MCName,
                  double** &Data, double** &MC, 
                  int ndata, int nMC, int d );

void Forker ( double **Data, double **MC, 
                      char *OutName,
                      int ndata, int nMC,
                      int d, int nnns, int nns[], int flag,
                      int idt                              ); 

void FindNeighbours ( double** Data,
                      double** MC,
                      char * OutName,
                      int ndata1, int ndata2, int nMC,
                      int d, int nnns, int nns[], int flag,
                      int SoughtIdt              ) ; 

inline double EucliDistance ( double *vect1, 
                              double *vect2,
                              int d               ) ; 

inline void CalcProbability ( LinkedList <sosed> &list, 
                              int SoughtIdt, int nnns, int nns[], int flag,
                              double probability[] ) ; 

void panic (char* str) ; 

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


// ----------------------
//   MAIN
// ----------------------

int main (int  argc, char* argv[])
{
   int     idt;               /* idt being sought */
   int     d;                 /* number of variables */
   long     ndata, nMC ;       /* number of lines in data and MC */

   char*   DataName = new char[30];       // data file name */
   char*   MCName = new char[30];         // MC file name */
   char*   OutName = new char[30];        // output file name */
   double  **Data, **MC; 
   int     nns[20];          // number of neighbours to consider
   int     nnns;             // total number of different nns
   int     flag;             // prob = sqrt(noise)/signal  OR signal/total 

   if (argc < 2) panic ("neighbours <Input File>") ;

   ReadInput ( argv[1],                   // reads command-line arguments
               DataName, MCName, 
               OutName,  nnns, nns, flag,
               idt                    );   
   CheckInput ( DataName, MCName,         // half-checks compatibility 
                ndata, nMC, d   );        // of data & MC files
   ReadArrays ( DataName, MCName,
                Data, MC, ndata, nMC, d );
   
   Forker ( Data, MC, OutName, ndata, 
            nMC, d, nnns, nns, flag,
            idt) ; 
 
   return 1;
}


// forks the process into m, each 
void Forker (double **Data, double **MC, 
                      char *OutName,
                      int ndata, int nMC,
                      int d, int nnns, int nns[], int flag,
                      int idt                              )

{
   char *command = new char[200]; 
   char *OutName1 = new char[200];


   if (flag == 1 || flag == 0) {
       FindNeighbours ( Data, MC, OutName,     // finds nns[i] nearest
                        0, ndata, nMC, d,      // MCnbrs for each data pt.,
                        nnns, nns, flag, idt ); // prints results in OF[i]
       return; 
   } 

   pid_t  pid; 
   int ndata1, ndata2; 
   int m = 0; 
   int status;
   int waited = 0;  
   long  counter, cnt1, cnt2; 

   for (counter  = 1 ; counter <= flag ; counter++) { 
      if ((pid = fork()) < 0) panic ("Unable to fork. Exit");  
      if (pid == 0)  { 
          ndata1 = (ndata * (counter-1)) / flag ; 
          ndata2 = (ndata * counter) / flag ;
          sprintf (OutName1, "%s%d", OutName, m+1); 
          FindNeighbours ( Data, MC, OutName1,    
                           ndata1, ndata2, nMC, d,
                           nnns, nns, flag, idt   );
          exit (0); 
      }

      m++ ; 

      if (m == 1) {
	 m = 0 ; 
         for ( cnt1 = 0; cnt1 < 1 ; cnt1++)    wait (&status);
         waited++; 
         for ( int cnt1 = 0 ; cnt1 < nnns ; cnt1++ )  
         { 
             sprintf (command, "cat ") ; 
             for ( int cnt2 = 1; cnt2 <= 1 ; cnt2++ )  
                 sprintf (command, "%s%s%d.%d ", command, OutName, 
                                                 cnt2, nns[cnt1] ); 
             sprintf (command, "%s >> %s.%d", command, OutName, nns[cnt1]);
             system (command); 
         } 
      }  
   }
   for (counter = 0 ; counter < flag-waited ; counter++ )   
            if (pid > 0) wait (&status); 

   sleep (240); 

   for ( int cnt1 = 0 ; cnt1 < nnns ; cnt1++ )  
   { 
     sprintf (command, "cat ") ; 
     for ( int cnt2 = 1; cnt2 <= flag-waited ; cnt2++ )  
        sprintf (command, "%s%s%d.%d ", command, OutName, cnt2, nns[cnt1]); 
     sprintf (command, "%s >> %s.%d", command, OutName, nns[cnt1]);
     system (command); 
   } 
}

/*  Transforms program arguments into corresponding variables
*/

void ReadInput ( char *FileName,
                 char *DataName, char* MCName, 
                 char *OutName,  int &nnns, int nns[],
                 int &flag, int &idt                   )
{
    ifstream   InFile (FileName) ;
    if (!InFile) panic ("No input file" );
  
    InFile >> DataName ;
    InFile >> MCName;
    InFile >> OutName;
    InFile >> idt ;
    InFile >> flag ;      // alfa = sqrt(noise)/signal  OR  signal/total  
    InFile >> nnns; 
    for (int cnt1 = 0; cnt1 < nnns; cnt1++ )  {
        if (InFile.eof()) panic ("Unexpected end of input file"); 
        InFile >> nns[cnt1]; 
    } 
    InFile.close(); 
}


/*  1. compares number of columns in MC and data files. 
       (#column(MC) - #coulmn(data) should be 1 )
    2. counts number of variables and number of points in data and MC
*/

void CheckInput ( char* DataName, char* MCName, 
                  long &ndata, long &nMC, int &d )
{
    int      nwords, nlines;
    double   DataNcolumns, MCNcolumns;
    char      *com = (char*) malloc (512);  // command
    int       cnt ;  /* counter */

    for (cnt = 0 ; cnt++ <= 1; ) 
    {
            if  (cnt == 1)   sprintf (com, "wc %s > dummyfile", DataName); 
            else             sprintf (com, "wc %s > dummyfile", MCName); 
   	    system (com);    

            ifstream  dummy ("dummyfile");
	    dummy >> nlines >> nwords; 
            if     (cnt == 1)  {
                       DataNcolumns = (double) nwords / nlines ;   
                       ndata = nlines ;
            } else { 
                      MCNcolumns = (double) nwords / nlines ; 
                      nMC = nlines ; 
            }
	    dummy.close();
    }

    if ( fabs ( DataNcolumns - INT (DataNcolumns) > 0.0001 )     ||
         fabs ( MCNcolumns - INT (MCNcolumns + 0.5) > 0.0001 )  ||
         INT (MCNcolumns) - INT (DataNcolumns) != 1  )
             panic ("Bad input files or unmatched number of variables\n") ; 

    d = INT (DataNcolumns) ; 
}


// ------------------------------------------------------------------------
// in the simplest and most naive version of this program 
// this function just reads the data and MC files 
// INTO RAM!!!
// ------------------------------------------------------------------------


void   ReadArrays ( char* DataName, char* MCName,
             double** &Data, double** &MC, 
             int  ndata, int  nMC, int d )
{
      ifstream  DataFile (DataName), MCFile (MCName); 
      long      cnt1; 

      if (! DataFile || ! MCFile )    panic (" No input Data or MC file");  

      Data = new double* [ndata]; 

      MC = new double* [nMC];


// reading from data file to data array

      for ( cnt1 = 0 ; cnt1 < ndata; cnt1++ ) 
      { 
          Data[cnt1] = new double [d]; 
          for (int cnt2 = 0; cnt2 < d; cnt2++ ) { 
               DataFile >> Data[cnt1][cnt2]; 
       	       if ( DataFile.eof() )
                  panic ( "Unexpected data EOF in ReadArrays" );
	  } 
      } 

// reading from MC file to MC array

      for ( cnt1 = 0 ; cnt1 < nMC; cnt1++ ) 
      { 
          MC[cnt1] = new double [d+1]; 
          for (int cnt2 = 0; cnt2 <= d; cnt2++ ) {
              MCFile >> MC[cnt1][cnt2]; 
	      if (MCFile.eof()) 
                    panic( "Unexpected MC EOF in ReadArrays\n");
          } 
      } 

      DataFile.close (); 
      MCFile.close (); 
}


// ------------------------------------------------------------------------
// The major procedure - finding nn-s for each point &
// printing the results.
// Complexity = O (D * M * d) 
// D - # of data points, M - # of MC points, d - dimensionality  
// ------------------------------------------------------------------------

void FindNeighbours ( double **Data, 
                      double **MC,
                      char *OutName, 
                      int ndata1, int ndata2, int nMC, 
                      int d, int nnns, int nns[], int flag, 
                      int SoughtIdt       )
{

    LinkedList <sosed> nearest(nns[nnns-1]);     // nn l. list for

    ifstream   DataFile, MCFile;
    ofstream   OutputFile[200];

    int         cnt1, cnt2, cnt;                 // counters
    double      prob[200];                       // returned probability 
    sosed       TempNeighbour;     
    char*       OName = new char[100]; 
//cout << "poehali-2" <<endl; 
    if (nnns > 200)  panic (" Too much output nn-s" ) ; 

    for (cnt = 0; cnt < nnns ; cnt++ ) {
            sprintf (OName, "%s.%d", OutName, nns[cnt]) ;       
            OutputFile[cnt].open (OName);
    } 
// cout << "poehali-3" <<endl; 
    for (cnt1 = ndata1 ; cnt1 < ndata2; cnt1++ ) // loop on all data 
    {
        for (cnt2 = 0; cnt2 < nMC; cnt2++)       // loop on all MC points
	{
	   TempNeighbour.distance = EucliDistance (Data[cnt1], MC[cnt2], d);
                                         // trying to put TempNbr. on list
           if ( cnt2 < nns[nnns-1] || 
                TempNeighbour.distance < nearest.MaxKey )
           {
                   TempNeighbour.id = INT (MC[cnt2][d]) ;  
                   nearest.Insert (TempNeighbour, TempNeighbour.distance);
// cout << "poehali-4" <<endl;
           }
	}

        CalcProbability (nearest, SoughtIdt, nnns, nns, flag, prob) ;
        for (cnt = 0; cnt < nnns ; cnt++ ) 
              OutputFile[cnt] << prob[cnt] << endl;   // writing output

        nearest.Emptify ();
    }
    for (cnt = 0; cnt < nnns ; cnt++ )  OutputFile[cnt].close ();
}


// gets distance between 2 linked lists of doubles; # of coords is chosen to 
// be the minimal of 2  

inline double EucliDistance ( double* vect1, double* vect2, int d ) 
{
    double   SquareSum = 0 ; 

    for ( int c = 0 ; c < d ; c++ )   
          SquareSum += (vect1[c]-vect2[c]) * (vect1[c]-vect2[c]) ;  

    return SquareSum;
}



// calculating probability for a data point to belong to chosen group
// prob = (number of nn with chosen ID) / (total number of nn) 

inline void CalcProbability ( LinkedList <sosed> &list, 
                              int SoughtIdt, int nnns, int nns[],
                              int flag, double probability[]) 
{
    sosed      ANeighbour;
    int        cnt1, cnt2 = 0 ;
    char       *com = new char[50];      
    int        ngood = 0;  

    if ( ! list.Start () )  panic ("Empty list in CalcProbability\n");

    for (cnt1 = 0; cnt1 < nnns; cnt1++ ) 
    { 
       while (cnt2++ < nns[cnt1]) {
          ANeighbour = list.GetValue() ; 
          if ( ! list.Forward () ) panic ("Bad list in CalcProbability"); 
          if ( abs (ANeighbour.id) == SoughtIdt )  ngood++;  
       } 
       if (flag>0)   probability[cnt1] = ngood / (double) nns[cnt1]; 
       else          probability[cnt1] = sqrt(nns[cnt1]-ngood) / ngood;  
    } 
}


//  panic

void panic (char* str)
{
    cerr  << str << endl ; 
    exit (0);
}