#include <math.h>
#include <stdlib.h>
#include <stdlib.h>
#include <fstream.h>
#include <iostream.h>
#include <unistd.h>


void panic (char *str); 

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

    if (argc < 3) {
        cout << "Info file format" << endl ;
        cout << "# of points in first gauss. # in second" << endl;
        cout << "dimensionality" << endl ;    
        cout << "sign: -1, 0 or 1" << endl ; // sign of ALL variables.
        cout << "sigma_11 sigma_21 center_11 center_21" << endl; 
        cout << "......." << endl ; 
        cout << "sigma_1n sigma_2n center_1n center_2n" << endl;
        cout << "(first index refers to number of distribution, " << endl ; 
        cout << "second - to that of coordinate)" << endl;
        panic ("gauss <infofile> <output file>");
    }

// reading data, opening files

    ifstream    infofile (argv[1]); 
    ofstream    outfile (argv[2]); 

    int     n1, n2, d, sign;
    double  c1[20], c2[20], sigma1[20], sigma2[20];
    
    infofile >> n1 >> n2 >> d >> sign; 
    for (int cnt=0 ; cnt < d ; cnt++) 
       infofile >> sigma1[cnt] >> sigma2[cnt] >> c1[cnt] >> c2[cnt];

// -------------loop over points  

    double   coord[20];  
    double   euclisum; 

// seeding rand
    int SEED; 
    system ("date +%M%S > dummy") ; 
    ifstream dum ("dummy" ) ; 
    dum >> SEED;     
    system ("rm dummy") ; 
    srand (SEED); 

// ------------- NOISE

    for (int npts1 = 0 ; npts1 < n1 ; ) {
       euclisum=0; 
       for (cnt = 0 ; cnt < d ; cnt++ ) { 
          switch (sign) {  
             case 0:
                 coord[cnt] = c1[cnt] + (rand()-16384.0)/1024.0;  
                 break; 
             case -1:
                 coord[cnt] = c1[cnt] - (rand())/2048.0;
                 break ; 
             case 1: 
                 coord[cnt] = c1[cnt] + (rand())/2048.0;
   	  }     

          euclisum -= (coord[cnt]-c1[cnt])*(coord[cnt]-c1[cnt]) /
                                             (sigma1[cnt] * sigma1[cnt]) ;
       }
       if ((rand()+1)/32768.0 < exp (euclisum) * 0.5) {
           npts1++; 
           for (cnt = 0 ; cnt < d ; cnt++ )  outfile << coord[cnt] << " " ;
           outfile << 0 << endl; 
       } 
    } 

// ++++++++++"SIGNAL" 

    for (int cnt2 = 0 ; cnt2 < n2 ; ) { 
       euclisum=0; 
       for (cnt = 0 ; cnt < d ; cnt++ ) { 
          switch (sign) {  
             case 0:
                 coord[cnt] = c2[cnt] + (rand()-16384.0)/1024.0;  
                 break; 
             case -1:
                 coord[cnt] = c2[cnt] - (rand())/2048.0;
                 break ; 
             case 1: 
                 coord[cnt] = c2[cnt] + (rand())/2048.0;
   	  }     
          euclisum -= (coord[cnt]-c2[cnt])*(coord[cnt]-c2[cnt]) / 
                                         (sigma2[cnt] * sigma2[cnt]) ; 
       }
       if ((rand()+1)/32768.0 < exp (euclisum) * 0.5) {
           cnt2++; 
           for (cnt = 0 ; cnt < d ; cnt++ )  outfile << coord[cnt] << " " ;
           outfile << 1 << endl; 
       } 
    } 
}

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