#undef   QUALITY_CHECKS

#include "scltypes.h"
#include "sclfuncs.h"

#if defined QUALITY_CHECKS
#include "libscl.h"
#endif

using namespace std;
using namespace scl;

int main()
{

#if defined QUALITY_CHECKS

  static INT_32BIT ix;

  ix = 421005;

  realmat probit;
  if (!vecread("probit.dat",probit)) error("Error, vecread failed");

  INTEGER ntrials = 1000000000;
  INTEGER nbins = probit.nrow();
  scl::realmat bins(nbins,1,0.0);

  REAL sum = 0.0;
  REAL sumsqr = 0.0;
  REAL true_mean = 0.0;
  REAL true_var = 1.0;

  INTEGER nex = 7;
  realmat exprobit(nex,1);
  realmat exexpect(nex,1);
  realmat excum(nex,1);
  realmat exbins(nex,1,0.0);

  exprobit[1] = -4.0;        excum[1] = 1.0 - 0.9999683288;
  exprobit[2] = -3.5;        excum[2] = 1.0 - 0.9997673709; 
  exprobit[3] = -3.0;        excum[3] = 1.0 - 0.998650101968370;
  exprobit[4] =  3.0;        excum[4] = 0.998650101968370;
  exprobit[5] =  3.5;        excum[5] = 0.9997673709;
  exprobit[6] =  4.0;        excum[6] = 0.9999683288;
  exprobit[7] =  REAL_MAX;   excum[7] = 1.0;

  exexpect[1] = excum[1];
  exexpect[2] = excum[2] - excum[1];
  exexpect[3] = excum[3] - excum[2];
  exexpect[4] = excum[4] - excum[3];
  exexpect[5] = excum[5] - excum[4];
  exexpect[6] = excum[6] - excum[5];
  exexpect[7] = excum[7] - excum[6];

  //cout << exprobit << excum << exexpect;

  for (INTEGER t=1; t<=ntrials; ++t) {
    REAL x = unsk(ix);
    sum += x;
    sumsqr += pow(x-true_mean,2);
    for (INTEGER i=1; i<=nbins; ++i) {
      if (x <= probit[i]) { bins[i] += 1.0; break; }
    }
    for (INTEGER i=1; i<=nex; ++i) {
      if (x <= exprobit[i]) { exbins[i] += 1.0; break; }
    }
  }
  REAL mean = sum/REAL(ntrials);
  REAL var = sumsqr/REAL(ntrials);

  std::cout << '\n';
  std::cout << "Mean is " << mean << ". "
       << "Variance is " << var << ". "
       << "Sample size is " << ntrials << ". ";
  std::cout << '\n';
  std::cout << "The mean should be " << true_mean << " and the variance "
       << true_var << ". ";
  std::cout << '\n';

  REAL expected = REAL(ntrials)/REAL(nbins);
  REAL chisq = 0.0;
  for (INTEGER i=1; i<=nbins; ++i) {
    chisq += pow(bins[i] - expected,2);
  }
  chisq /= expected;
  INTEGER df = nbins - 1;
  std::cout << '\n';
  std::cout << "Chi square is " << chisq
       << " on " << df << " degrees freedom."
       << " The z statistic is "
       << (chisq - REAL(df))/sqrt(2.0*REAL(df)) << '\n';
  std::cout << "Remember that this test is one sided. "
       << "Reject for large positive values only. " << '\n';

  chisq = 0.0;
  for (INTEGER i=1; i<=nex; ++i) {
    exexpect[i] *= REAL(ntrials);
    chisq += pow(exbins[i]-exexpect[i],2)/exexpect[i];
  }
  df = nex - 1;
  std::cout << '\n';
  std::cout << "Chi square is " << chisq
       << " on " << df << " degrees freedom."
       << " The z statistic is "
       << (chisq - REAL(df))/sqrt(2.0*REAL(df)) << '\n';
  std::cout << "Remember that this test is one sided. "
       << "Reject for large positive values only. " << '\n';

  std::cout << exbins << exexpect;
  for (INTEGER i=1; i<=nex; ++i) {
    exbins[i] -= exexpect[i];
    exbins[i] /= sqrt(exexpect[i]);
  }
  std::cout << exbins;
  std::cout << '\n';

#else

  static INT_32BIT ix;
  INTEGER i;
  REAL    x;
  
  ix = 15;
  for (i=0; i < 50; i++) {
    x=unsk(&ix);
    cout << x << "\n";
  }

  cout << endl;

  ix = 15;
  for (i=0; i < 50; i++) {
    x=unsk(ix);
    cout << x << "\n";
  }

#endif

  exit(0);
}
