#include "libscl.h"
#include "kronprd.h"
#include "spline_interpolator.h"

#include "var1EXX.h"
#include "lprior.h"
#include "crra_mf.h"

#include <map>

using namespace std;
using namespace scl;

const bool debug = true;

// Raw sim is DCFdatayp.dat, M=5, n=86

// Col 1: Year: 1930 -- 2015
// Col 2: log of extracted mrs
// Col 3: growth of Real per capita GDP
// Col 4: GDP/Corporate profits
// Col 5: growth of Corporate profits

// Data for use is sim_orig, n=86, M=3

// Row 1: log of extracted mrs
// Row 2: log GDP/Corporate profits
// Row 3: log growth of Corporate profits

REAL psi2phi(const realmat& z, const spline_interpolator& Fz1, const
             spline_interpolator& qnorm)
{
  INTEGER d= z.size();
  if (z.ncol() != 1) error("Error,  psi2phi, z not a vector");

  REAL z1 = z[1];
  REAL u1 = Fz1(z1);
  if (u1 < 0.0001) {
    u1 = 0.0001;
  }
  else if (u1 > 0.9999) {
    u1 = 0.9999;
  }
  REAL x1 = qnorm(u1);
  REAL dz1 = qnorm.derivative(u1)*Fz1.derivative(z1);
  dz1 = (dz1 < 0.0 ? -dz1 : dz1);
  realmat z_trf = z;
  z_trf[1] = x1;  

  REAL sse = 0.0;
  for (INTEGER k=1; k<=d; ++k) sse += pow(z_trf[k],2);
  const REAL logtwopi = log(8.0*atan(1.0));
  REAL log_likelihood_trf = -0.5*sse - 0.5*REAL(d)*logtwopi + log(dz1);
  return log_likelihood_trf;
}

int main(int argc, char** argp, char** envp)
{
  const INTEGER n_obs = 50;

  const INTEGER mf_lags = 1;
  const INTEGER HAC_lags = 0;
  const INTEGER n_prior_reps = 1000;
  const INTEGER ngrid = 11;

  const INTEGER L = 1;
  const INTEGER M = 3;
  const INTEGER years = 30;
  const INTEGER p_rho = M + M*M + (M*M+M)/2 + 2;

  const INTEGER n_raw_items = 5;  //year, mrs, gdp, gdp_cp, cp
  const INTEGER n_raw_obs   = 86; //Annual 1930 to 2015

  //const INTEGER loc_dat_yr     = 1;  
  const INTEGER loc_dat_mrs    = 2;
  //const INTEGER loc_dat_gdp    = 3;
  const INTEGER loc_dat_gdp_cp = 4;
  const INTEGER loc_dat_cp     = 5;

  const INTEGER loc_var_mrs    = 1;
  const INTEGER loc_var_gdp_cp = 2;
  const INTEGER loc_var_cp     = 3;

  const INTEGER mrs_pos = loc_var_mrs;
  const INTEGER cf_pos = loc_var_cp;

  const REAL beta_lo = 0.8;
  const REAL beta_hi = 0.99;

  const REAL gamma_lo = 0.5;
  const REAL gamma_hi = 100.0;

  const REAL beta_prior_mean = 0.95319;
  const REAL beta_prior_sdev = 0.054691;
  const REAL gamma_prior_mean = 24.503;
  const REAL gamma_prior_sdev = 28.243;

  const REAL beta_inc = 0.01;
  const REAL gamma_inc = 3.0;

  const REAL minus_log_root_two_pi = -0.5*(M_LN2 + log(M_PI));

  realmat sim_raw;
  INTEGER ct = n_raw_items*n_raw_obs; 

  if (ct!=readtable("DCFdatayp.dat",sim_raw)) error("Error, readtable failed");

  realmat sim_orig(M,n_raw_obs);
  for (INTEGER j=1; j<=n_raw_obs; ++j) {
    sim_orig(loc_var_mrs,j)    = sim_raw(j,loc_dat_mrs);
    sim_orig(loc_var_gdp_cp,j) = sim_raw(j,loc_dat_gdp_cp);
    sim_orig(loc_var_cp,j)     = sim_raw(j,loc_dat_cp);
  }

  realmat y0(M,1);
  for (INTEGER i=1; i<=M; ++i) y0[i] = sim_orig(i,n_raw_obs);

  if (debug) {
    cout << '\n';
    cout << "First 6 observations: mrs,gdp_cp,cp, n_obs = " << n_obs;
    cout << sim_orig("",seq(1,6));
    cout << '\n';
    cout << "Last 6 observations " << sim_orig("",seq(n_raw_obs-5,n_raw_obs));
    cout << '\n';
    cout << "y0 " << y0 << '\n';
    cout << '\n';
    cout.flush();
  }
  
  realmat rho_orig, hess_orig;

  if(!vecread("dcf.rho_mode.dat",rho_orig)) error("Error, vecread failed");
  if(!vecread("dcf.V_hat_hess.dat",hess_orig)) error("Error, vecread failed");

  realmat R_orig = hess_orig;
  if(factor(R_orig)!=0) error("Error, factor failed");

  if (debug) {
    cout << '\n';
    cout << "rho_orig = " << rho_orig << '\n';
    cout << "hess_orig = " << hess_orig << '\n';
    //cout << "R_orig = " << R_orig << '\n';
    cout << '\n';
    cout.flush();
  } 

  realmat rho = rho_orig;

  realmat b0, B, S, R, theta;
  frac_rho(M, L, rho, b0, B, S, R, theta);

  REAL beta = theta[1];
  REAL gamma = theta[2];

  if (debug) {
    cout << '\n';
    cout << "b0 = " << b0 << '\n';
    cout << "B = " << B << '\n';
    cout << "R = " << R << '\n';
    cout << "S = " << S << '\n';
    cout << "theta = " << theta << '\n';
    cout.flush();
  }

  string filename = "Z_" + fmt('d',4,n_obs)('0') + "_knots.txt";
  realmat Fz_knots;
  if (!readtable(filename.c_str(),Fz_knots)) error("Error, " + filename);
  realmat z1_knots = Fz_knots("",1);
  realmat pz1_knots = Fz_knots("",2);
  spline_interpolator Fz1(z1_knots,pz1_knots);

  if (debug) {
    cout << '\n';
    cout << "\t Fz1 constructed from " << filename << '\n';
    cout.flush();
  }

  filename = "normal_knots.txt";
  realmat normal_knots;
  if (!readtable(filename.c_str(),normal_knots)) error("Error, " + filename);
  realmat xn_knots = normal_knots("",1);
  realmat pn_knots = normal_knots("",2);
  spline_interpolator qnorm(pn_knots,xn_knots);

  if (debug) {
    cout << '\n';
    cout << "\t qnorm constructed from " << filename << '\n';
    cout.flush();
  }

  realmat ecf, pvcf, pv1, cumecf, dcf, yld;
  realmat ylag = y0;

  dcfyld(b0,B,S,ylag,mrs_pos,cf_pos,years,ecf,pvcf,pv1,cumecf,dcf,yld);
  REAL price_lag = dcf[years];

  INT_32BIT seed = 893475;

  realmat sim(M,n_obs);
  realmat data(2,n_obs);

  realmat y(M,1), u(M,1);

  filename = "lsrlcg.txt";
  ofstream fout;
  fout.open(filename.c_str());
  if (!fout) error("Error, cannot open " + filename);

  for (INTEGER t=1; t<=5000; ++t) {
    for (INTEGER i=1; i<=M; ++i) u[i] = unsk(seed);
    y = b0 + B*ylag + R*u;
    dcfyld(b0,B,S,y,mrs_pos,cf_pos,years,ecf,pvcf,pv1,cumecf,dcf,yld);
    REAL price = dcf[years];
    REAL lsr = log(price+y[cf_pos]) - log(price_lag);
    REAL lcg = (log(beta) - y[mrs_pos])/gamma;
    fout << fmt('f',22,16,lsr) << fmt('f',22,16,lcg) << '\n';
    if (t<=n_obs) {
      for (INTEGER i=1; i<=M; ++i) sim(i,t) = y[i];
      data(1,t) = lsr;
      data(2,t) = lcg;
    }
    ylag = y;
    price_lag = price;
  }

  fout.clear(); fout.close();

  if (debug) {
    cout << '\n';
    cout << "First 6 sim: mrs,gdp_cp,cp " << sim("",seq(1,6));
    cout << '\n';
    cout << "Last 6 sim: mrs,gdp_cp,cp " << sim("",seq(n_obs-5,n_obs));
    cout << '\n';
    cout << "First 6 data: lsr, lcg " << data("",seq(1,6));
    cout << '\n';
    cout << "Last 6 data: lsr, lcg " << data("",seq(n_obs-5,n_obs));
    cout << '\n';
    cout.flush();
  }

  bool success;
  bool rv = true;

  mle::crra_moment_function crra_mf;

  success = crra_mf.set_data(&data);
  if (!success) rv = false;

  success = crra_mf.set_sample_size(n_obs);
  if (!success) rv = false;

  success = crra_mf.set_L(mf_lags);
  if (!success) rv = false;

  success = crra_mf.set_theta(theta);
  if (!success) rv = false;

  // INTEGER T0 = crra_mf.get_T0();

  INTEGER d = crra_mf.get_d();

  gmm crra_gmm(&crra_mf, mf_lags, &data, n_obs, HAC_lags);

  crra_gmm.set_correct_W_for_mean(true);
  crra_gmm.set_warning_messages(false);
  crra_gmm.set_regularize_W(true);

  success = crra_gmm.set_data(&data);
  if (!success) rv = false;

  realmat z;
  denval likelihood = crra_gmm.likelihood(theta,z);

  if (debug) {
    cout << '\n';
    cout << "d " << d << '\n';
    cout << '\n';
    cout << "likelihood " 
         << "("<< likelihood.positive <<", "<< likelihood.log_den << ")"<<'\n';
    cout << '\n';
    cout << "z " << z << '\n';
    cout << '\n';
  }

  if (p_rho != rho.size()) error("Error, p_rho !=rho.size()");

  realmat e(p_rho,1);

  realmat Z_save(n_prior_reps,d);
  realmat rho_save(n_prior_reps,p_rho);

  realmat rho_left(p_rho,1,REAL_MAX);
  realmat rho_rite(p_rho,1,-REAL_MAX);

  realmat xgrid(ngrid,1);       // beta
  realmat ygrid(ngrid,1);       // gamma
  realmat lgrid(ngrid,ngrid);   // likelihood
  realmat lgrid_alt(ngrid,ngrid); 
  realmat tgrid(ngrid,ngrid);   // transformed likelihood
  realmat pgrid(ngrid,ngrid);   // prior
  realmat agrid(ngrid,ngrid);   // adj
  realmat dgrid(ngrid,ngrid);   // fabs(lgrid) - fabs(agrid)
  realmat lgrid_inc(ngrid,ngrid); // likelihood increment
  realmat tgrid_inc(ngrid,ngrid); // transformed likelihood increment
  realmat lgrid_save(ngrid,ngrid);
  realmat lgrid_alt_save(ngrid,ngrid);
  realmat tgrid_save(ngrid,ngrid);
  realmat pgrid_save(ngrid,ngrid);
  realmat agrid_save(ngrid,ngrid);
  realmat lgrid_inc_save(ngrid,ngrid);
  realmat tgrid_inc_save(ngrid,ngrid);

  for (INTEGER i=1; i<=ngrid; ++i) {
    REAL delta = REAL(i-1)/REAL(ngrid-1);
    xgrid[i] = beta_lo + (beta_hi-beta_lo)*delta;
    ygrid[i] = gamma_lo + (gamma_hi-gamma_lo)*delta;
  }

  REAL min_adj = REAL_MAX;
  REAL max_adj = -REAL_MAX;
  REAL max_dif = -REAL_MAX;

  REAL tmax_ij = -REAL_MAX;
  REAL tmax = -REAL_MAX;

  INTEGER W_error_count = 0;
  INTEGER W_success_count = 0;

  //INTEGER T_error_count = 0;
  //INTEGER T_success_count = 0;

  INTEGER Z_error_count = 0;

  for (INTEGER rep=1; rep<=n_prior_reps; ++rep) {

    bool bad_draw = false;
    do {
      for (INTEGER i=1; i<=p_rho; ++i) e[i] = unsk(seed);
      rho = rho_orig + R_orig*e;

      frac_rho(M, L, rho, b0, B, S, R, theta);
      beta = theta[1];
      gamma = theta[2];

      bad_draw = false;
      if ( (beta  <= 0.8) || (0.99 <= beta  ) ) bad_draw = true;
      if ( (gamma <= 0.0) || (100.0 <= gamma) ) bad_draw = true;

    } while (bad_draw);

    if (debug) if (rep%100 == 0) {
      cout << rep << '\n';
      cout.flush();
    }

    for (INTEGER i=1; i<=p_rho; ++i) {
      rho_save(rep,i) = rho[i];
      if (rho[i] < rho_left[i]) rho_left[i] = rho[i];
      if (rho[i] > rho_rite[i]) rho_rite[i] = rho[i];
    }

    ylag = y0;
    dcfyld(b0,B,S,ylag,mrs_pos,cf_pos,years,ecf,pvcf,pv1,cumecf,dcf,yld);
    REAL price_lag = dcf[years];

    for (INTEGER t=1; t<=n_obs; ++t) {
      for (INTEGER i=1; i<=M; ++i) u[i] = unsk(seed);
      y = b0 + B*ylag + R*u;
      for (INTEGER i=1; i<=M; ++i) sim(i,t) = y[i];
      dcfyld(b0,B,S,y,mrs_pos,cf_pos,years,ecf,pvcf,pv1,cumecf,dcf,yld);
      REAL price = dcf[years];
      data(1,t) = log(price+y[cf_pos]) - log(price_lag);
      data(2,t) = (log(beta) - y[mrs_pos])/gamma;
      ylag = y;
      price_lag = price;
    }

    success = crra_gmm.set_data(&data);
    if (!success) rv = false;

    INTEGER ier = 0;

    for (INTEGER i=1; i<=ngrid; ++i) {
      for (INTEGER j=1; j<=ngrid; ++j) {

	beta  = theta[1] = xgrid[i] + beta_inc;
        gamma = theta[2] = ygrid[j] + gamma_inc;

        lgrid_inc(i,j)  = crra_gmm.likelihood(theta,z).log_den;
	tgrid_inc(i,j)  = psi2phi(z,Fz1,qnorm);

	if (crra_gmm.get_W_numerr()) ier = crra_gmm.get_W_numerr();

        beta  = theta[1] = xgrid[i];
        gamma = theta[2] = ygrid[j];

        REAL log_likelihood = crra_gmm.likelihood(theta,z).log_den;

	if (crra_gmm.get_W_numerr()) ier = crra_gmm.get_W_numerr();

	REAL log_likelihood_trf = psi2phi(z,Fz1,qnorm);

        /*
        if (debug) {
          cerr << fmt('f',10,5,beta) << fmt('f',10,5,gamma) << fmt('f',10,5,u1) 
               << fmt('f',10,5,x1) << fmt('f',10,5,dz1) 
               << fmt('f',10,5,log_likelihood_trf) << '\n';
        }
        */

        REAL z_beta = (beta-beta_prior_mean)/beta_prior_sdev;
        REAL z_gamma = (gamma-gamma_prior_mean)/gamma_prior_sdev;
        REAL e_beta 
             = minus_log_root_two_pi-log(beta_prior_sdev)-0.5*pow(z_beta,2);
        REAL e_gamma 
             = minus_log_root_two_pi-log(gamma_prior_sdev)-0.5*pow(z_gamma,2);
        REAL log_pe_prior = e_beta + e_gamma;

        const REAL sf = 1.0 - exp(1.0);

        REAL tz1 = tanh(z[1]/4.0);
        REAL dZ1;
        if (-1.0 < tz1 && tz1 < 1.0) {
          REAL top = -4.0*(1.0 - tz1);
          REAL bot = 1 - pow(tz1,2);
          dZ1 = top/bot;
        }
        else if (tz1 >= 1.0){
          dZ1 = -2.0;
        }
        else {
          dZ1 = -REAL_MAX;
          ++Z_error_count;
        }
        REAL det = dZ1*sf*sf;
        REAL adj = det < 0.0 ? -det : det;
        REAL log_adj = log(adj);
        
        if (adj != adj) {
          //cerr << beta <<' '<< gamma <<' '<< z[1]
          //     <<' '<< s1 <<' '<< esb <<' '<< adj << '\n';
          cerr << z[1] <<' '<< tanh(z[1]/4.0) << '\n';
        }

        lgrid(i,j) = log_likelihood;       // likelihood
        tgrid(i,j) = log_likelihood_trf;   // transformed likelihood
        pgrid(i,j) = log_pe_prior;         // prior
        agrid(i,j) = log_adj;              // adj

        min_adj = log_adj < min_adj ? log_adj : min_adj;
        max_adj = log_adj > max_adj ? log_adj : max_adj;
        REAL tval = fabs(lgrid(i,j) - tgrid(i,j));
        tmax_ij = (tval > tmax_ij ? tval : tmax_ij);

        if (ier == 0) {
          ++W_success_count;
        }
        else {
          ++W_error_count;
          /*
          if (debug) {
            cerr << "W_numerr = " << fmt('d',3,ier)
                 << fmt('f',8,4,beta) << fmt('f',9,4,gamma) 
                 << fmt('f',7,2,lgrid(i,j)) << fmt('f',7,2,agrid(i,j)) <<'\n';
          }
          */
        }

      }
    }
    if (tmax_ij > tmax) {
      tmax = tmax_ij;
      lgrid_alt_save = lgrid;
      tgrid_save = tgrid;
      lgrid_inc_save = lgrid_inc;
      tgrid_inc_save = tgrid_inc;
      filename = "transform_rho_" + fmt('d',4,n_obs)('0') + ".txt";
      INTEGER cnt = writetable(filename.c_str(),rho,20,16);
      if (cnt == 0) warn("Warning, writetable failed, " + filename);
    }
    REAL dif = max_adj-min_adj;
    if (dif > max_dif && ier == 0) {
      max_dif = dif;
      lgrid_save = lgrid;
      pgrid_save = pgrid;
      agrid_save = agrid;
    }
    ier = 0;
  }

  lgrid = lgrid_save;
  lgrid_alt = lgrid_alt_save;
  tgrid = tgrid_save;
  pgrid = pgrid_save;
  agrid = agrid_save;
  lgrid_inc = lgrid_inc_save;
  tgrid_inc = tgrid_inc_save;

  for (INTEGER i=1; i<=ngrid*ngrid; ++i) {
     dgrid[i] = fabs(lgrid[i]) - fabs(agrid[i]);
     dgrid[i] *= M_LOG10E;
  } 

  cout << "prior_range"<< cbind(rho_left("19:20",1),rho_rite("19:20",1)) <<'\n';
  cout << "xgrid " << xgrid << '\n';
  cout << "ygrid " << ygrid << '\n';
  cout << "lgrid " << lgrid << '\n';
  cout << "pgrid " << pgrid << '\n';
  cout << "agrid " << agrid << '\n';
  cout << "dgrid " << dgrid << '\n';
  cout << "lgrid_alt " << lgrid_alt << '\n';
  cout << "lgrid_inc " << lgrid_inc << '\n';
  cout << "tgrid " << tgrid << '\n';
  cout << "tgrid_inc " << tgrid_inc << '\n';
  //cout << "lgrid - tgrid " << lgrid - tgrid << '\n';

  REAL dmin = REAL_MAX;
  for (INTEGER i=1; i<=ngrid*ngrid; ++i) {
    if (IsFinite(dgrid[i])) {
      dmin = (dgrid[i] < dmin ? dgrid[i] : dmin);
    }
  }

  cout << '\n';
  cout << "dmin = " << dmin << '\n';
  cout << "tmax = " << tmax << '\n';

  filename = "contour_" + fmt('d',4,n_obs)('0') + ".tex";
  fout.open(filename.c_str());
  if (!fout) error("Error, cannot open " + filename);

  fout << "\\multicolumn{" << ngrid+1 << "}{c}{Log Likelihood}\\\\" << '\n';
  fout << "$\\gamma/\\beta$&";
  for (INTEGER j=1; j<=ngrid; ++j) fout << "& " <<fmt('f',7,2,xgrid[j]) <<' ';
  fout << " \\\\" << '\n';
  for (INTEGER i=1; i<=ngrid; ++i) {
    fout << fmt('f',6,2,ygrid[i]) << " &";
    for (INTEGER j=1; j<=ngrid; ++j) fout <<"& "<<fmt('f',7,2,lgrid(i,j)) <<' ';
    fout << " \\\\" << '\n';
  }

  fout << " \\\\" << '\n';

  fout << "\\multicolumn{" << ngrid+1 << "}{c}{Log Adjustment}\\\\" << '\n';
  fout << "$\\gamma/\\beta$&";
  for (INTEGER j=1; j<=ngrid; ++j) fout << "& " <<fmt('f',7,2,xgrid[j]) <<' ';
  fout << " \\\\" << '\n';
  for (INTEGER i=1; i<=ngrid; ++i) {
    fout << fmt('f',6,2,ygrid[i]) << " &";
    for (INTEGER j=1; j<=ngrid; ++j) fout <<"& "<<fmt('f',7,2,agrid(i,j)) <<' ';
    fout << " \\\\" << '\n';
  }

  fout.clear(); fout.close();

  filename = "transform_" + fmt('d',4,n_obs)('0') + ".tex";
  fout.open(filename.c_str());
  if (!fout) error("Error, cannot open " + filename);

  fout << "\\multicolumn{" << ngrid+1 << "}{c}{Log Likelihood}\\\\" << '\n';
  fout << "$\\gamma/\\beta$&";
  for (INTEGER j=1; j<=ngrid; ++j) fout << "& " <<fmt('f',7,2,xgrid[j]) <<' ';
  fout << " \\\\" << '\n';
  for (INTEGER i=1; i<=ngrid; ++i) {
    fout << fmt('f',6,2,ygrid[i]) << " &";
    for (INTEGER j=1; j<=ngrid; ++j) {
      fout <<"& "<<fmt('f',7,2,lgrid_alt(i,j))<<' ';
    }
    fout << " \\\\" << '\n';
  }

  fout << " \\\\" << '\n';

  fout << 
    "\\multicolumn{" <<ngrid+1<< "}{c}{Log Transformed Likelihood}\\\\"<<'\n';
  fout << "$\\gamma/\\beta$&";
  for (INTEGER j=1; j<=ngrid; ++j) fout << "& " <<fmt('f',7,2,xgrid[j]) <<' ';
  fout << " \\\\" << '\n';
  for (INTEGER i=1; i<=ngrid; ++i) {
    fout << fmt('f',6,2,ygrid[i]) << " &";
    for (INTEGER j=1; j<=ngrid; ++j) fout <<"& "<<fmt('f',7,2,tgrid(i,j)) <<' ';
    fout << " \\\\" << '\n';
  }

  fout.clear(); fout.close();

  filename = "reject_" + fmt('d',4,n_obs)('0') + ".tex";
  fout.open(filename.c_str());
  if (!fout) error("Error, cannot open " + filename);

  realmat rej(ngrid,ngrid);

  for (INTEGER j=1; j<=ngrid; ++j) {
    for (INTEGER i=1; i<=ngrid; ++i) {
      REAL diff = lgrid(i,j) - lgrid_inc(i,j);
      if (diff > 0.0) {
        rej(i,j) = 1.0;
      }
      else {
        rej(i,j) = exp(diff);
      }
    }
  }

  fout << "\\multicolumn{" << ngrid+1 << "}{c}{Log Likelihood}\\\\" << '\n';
  fout << "$\\gamma/\\beta$&";
  for (INTEGER j=1; j<=ngrid; ++j) fout << "& " <<fmt('f',7,2,xgrid[j]) <<' ';
  fout << " \\\\" << '\n';
  for (INTEGER i=1; i<=ngrid; ++i) {
    fout << fmt('f',6,2,ygrid[i]) << " &";
    for (INTEGER j=1; j<=ngrid; ++j) fout <<"& "<<fmt('f',4,2,rej(i,j)) <<' ';
    fout << " \\\\" << '\n';
  }
  fout << " \\\\" << '\n';

  for (INTEGER j=1; j<=ngrid; ++j) {
    for (INTEGER i=1; i<=ngrid; ++i) {
      REAL diff = tgrid(i,j) - tgrid_inc(i,j);
      if (diff > 0.0) {
        rej(i,j) = 1.0;
      }
      else {
        rej(i,j) = exp(diff);
      }
//cerr << diff << '\n';
    }
  }

  fout << 
    "\\multicolumn{" <<ngrid+1<< "}{c}{Log Transformed Likelihood}\\\\"<<'\n';
  fout << "$\\gamma/\\beta$&";
  for (INTEGER j=1; j<=ngrid; ++j) fout << "& " <<fmt('f',7,2,xgrid[j]) <<' ';
  fout << " \\\\" << '\n';
  for (INTEGER i=1; i<=ngrid; ++i) {
    fout << fmt('f',6,2,ygrid[i]) << " &";
    for (INTEGER j=1; j<=ngrid; ++j) fout <<"& "<<fmt('f',4,2,rej(i,j)) <<' ';
    fout << " \\\\" << '\n';
  }

  fout.clear(); fout.close();

  cout << "W_success_count = " << W_success_count << '\n';
  cout << "W_error_count = " << W_error_count << '\n';
  cout << "Z_error_count = " << Z_error_count << '\n';

  cout << '\n';
  cout << (rv ? "success" : "failure") << '\n';

  return 0;
}

