#include <iostream>
#include <fstream>
#include <math.h>
#include "libscl.h"

       
using namespace scl;
using namespace std; 

int main()
{
   INTEGER ier, nvar, lsim; 
   ier  = 0;
   INT_32BIT jseed;
   jseed = 19291011;

   nvar = 3;
   lsim = 5000;
   REAL b[1+nvar];    
   REAL r12, r13, r23;
   REAL sx[nvar]; 
   REAL z1, z2, z3, logit;
   REAL zbar = 0.0, ybar = 0.0, u;
   
   REAL z[lsim], x1[lsim], x2[lsim], x3[lsim];
   INTEGER y[lsim];

   b[0] = 1.00; b[1] = 1.5; b[2] = 2.1; b[3] = 0.75;
   
   r12 = 0.00; r13=0.00;
   r23 = 0.95;
   
   sx[0] = 1.00;
   sx[1] = 3.45;
   sx[2] = 2.66;
  
  for (int i=0; i<lsim; i++) {
   z1 = unsk(jseed);
   z2 = unsk(jseed);
   z3 = unsk(jseed);
   x3[i] = z3;
   x2[i] = sqrt(1.0 - r23*r23)*z2 + r23*z3;
   x1[i] = sqrt(1.0 - r12*r12 - r13*r13)*z1 + r12*z2 + r13*z3;

   x1[i] *= sx[0];
   x2[i] *= sx[1];
   x3[i] *= sx[2];
      
   z[i] = b[0] + b[1]*x1[i] + b[2]*x2[i] + b[3]*x3[i];
   zbar = zbar + z[i];
   
   y[i] = 0;
   
   logit = exp(z[i])/(1.0 + exp(z[i]));
   
   u = ran(jseed);
   
  if ( u < logit ) {y[i] = 1;}

   ybar += y[i];
  }
  
  zbar /= lsim;
  ybar /= lsim;
  
   ofstream logitsim("logitsim.out");
   if (! logitsim)  {  
       cout << "Error opening logitsim.out file" << endl;  return -1;
   }
   
   logitsim << "beta: " << 
    "  " << fmt('f',8,4,b[0]) << 
    "  " << fmt('f',8,4,b[1]) <<
    "  " << fmt('f',8,4,b[2]) <<
    "  " << fmt('f',8,4,b[3]) << endl;
    
    logitsim << "r12 = " << fmt('f',8,4,r12) << "  " <<
                "r13 = " << fmt('f',8,4,r13) << "  " <<
                "r23 = " << fmt('f',8,4,r23) << endl;
                
    logitsim << "sx1 = " << fmt('f',8,4,sx[0]) << "  " <<
                "sx2 = " << fmt('f',8,4,sx[1]) << "  " <<
                "sx3 = " << fmt('f',8,4,sx[2]) << "  " << endl;
                
    logitsim << endl;
    
    logitsim << "zbar = " << fmt('f',8,4,zbar) << "  ybar = " << fmt('f',8,4,ybar) << endl << endl;
    
    
    
    for (int i=0; i<lsim; i++) {
       logitsim << fmt('i',2,y[i]   ) << " " <<
                   fmt('f',8,4,x1[i]) << " " <<
                   fmt('f',8,4,x2[i]) << " " <<
                   fmt('f',8,4,x3[i]) << " " << endl;
    }      
    

   logitsim.close();
   return 0;
/*end*/
}

 /*
      sampxout << fmt('i',12,k) << " " 
         << fmt('e',26,17,p[k])   << " " 
         << fmt('e',26,17,v[k])   << " "
         << fmt('e',26,17,levy[k])<< " " 
         << fmt('i',10,njmp) << endl;        
 */
   
//   cout << i << " " << x1[i] << " " << x2[i] << " " << x3[i] <<endl;