/* *************************************************
   optimal design for drug screening, no covariates
   Model: 	see forwsim.h
   design:	see forwsim.h

   main:	main driver
   init:	read in initialization file "init.dsg"

   dsg_forward:	forward simulation
   dsg_opt:	optimization
   dsg_u:	observed utility

   mcmc:	posterior estimates
   mcmc_ms:	posterior moments
   
   Version 7:   - allows hyperparameters for theta generated from a random distrib.
                - allows pooling when estimating posterior moments via Empirical Bayes
                - sd can be set to log or original scale
                - decision rule reparameterized as 3 points
                - saves estimated type I&II error probabilities, posterior mean & sd, terminal decisions
                - allows to use real theta instead of posterior expectation to evaluate phase III utility
                - allows to drop phase III term from utility
                - terminal decision is to stop if m<=th0. Phase III sample size calculated under H1:theta=m
		- accepts general discrete distribution for number of treatments arising over time
		- allows to include competition for patients by seting max nb patients available at any time point

   Version for R: - init.dsg cannot have any comments, basically it has to be formatted as R function forwsim outputs init.dsg
                  - option verbose is disabled so no progress info is displayed on the screen.

 ************************************************* */


#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include "cforwsim.h"
#include "cstat.h"

struct parameter th;
struct data dta;
struct design dsg;
struct mcpars mc;

int K;

void cforwsim()
{
  init();
  dsg_forw(1);
  Uhat_scale(1);
  Uhat_write(1);
  /* finish();  */
}



/* *************************************************
   init
 ************************************************* */

void init(){
  init_read();
  init_memalloc();
  init_pars();
}

void init_read(){
  FILE *ifile;
  int N;
  double lo,hi,s0;


  /* control file */
  ifile = openIn("init.dsg");	/* open for input */

  /* dimensions */
  scanInt(" dta-N ", &dta.N);
  N = dta.N;

  /* hyperpriors */
  scanDouble("th-u ",&th.u);   /* prior pars for th_i  */
  scanDouble("th-v ",&th.v);   
  scanDouble("th-au ",&th.au); /* hyperpriors */
  scanDouble("th-bu ",&th.bu);
  scanDouble("th-av ",&th.av);
  scanDouble("th-bv ",&th.bv);

  /* Utility function to use, use of hyperpriors, empirical Bayes pooling, SD scale */
  scanInt("dsg-uth ", &dsg.uth);
  if ((dsg.uth<0) | (dsg.uth>2)) err_msg("init","invalid uth (should be between 0 and 2)\n",0,0,0);
  scanInt("dsg-hyper ", &dsg.hyper);
  if ((dsg.hyper<0) | (dsg.hyper>1)) err_msg("init","invalid hyper (should be 0 or 1)\n",0,0,0);
  scanInt("dsg-pool ", &dsg.pool);
  if ((dsg.pool<0) | (dsg.pool>1)) err_msg("init","invalid pool (should be 0 or 1)\n",0,0,0);
  scanInt("dsg-logsd", &dsg.logsd);
  if ((dsg.logsd<0) | (dsg.logsd>1)) err_msg("init","invalid logsd (should be 0 or 1)\n",0,0,0);

  /* limit of available patients at any given time point */
  scanInt("dsg-patlimit ", &dsg.patlimit);
  
  /* initial number of treatments */
  scanInt(" dsg-kini ", &dsg.kini);
  if (dsg.kini<1) err_msg("init","invalid k-ini (should be >=1)\n",0,0,0);
  
  /* design paramters */
  scanDouble("dsg-c0 ", &dsg.c0);
  scanDouble("dsg-c1 ", &dsg.c1);
  scanDouble("dsg-c2 ", &dsg.c2);
  
  scanDouble("dsg-p0 ", &dsg.p0);
  scanDouble("dsg-a ", &dsg.a);   /* alpha */
  scanDouble("dsg-b ", &dsg.b);   /* beta */
  scanDouble("dsg-s0 ", &dsg.s0); /* lower bound on log(s.d) or s.d. */
  s0 = dsg.s0;
  if (dsg.logsd==1 && ((s0> -2) | (s0 < -4)) )
    err_msg("init","invalid s0 (in log scale it should be -4<s0< -2).\n",0,0,0);
  if (dsg.logsd==0 && ((s0> .15) | (s0 < .01)) )
    err_msg("init","invalid s0 (in original scale it should be .01<s0< .15).\n",0,0,0);
  scanDouble("dsg-s1 ", &dsg.s1);
  if (s0>=dsg.s1)
      err_msg("init","invalid s1 (should be >s0)\n",0,0,0);

  /* mcmc tuning pars */
  scanInt("mc-nsim ", &mc.nsim);  /* MC sample size for forward sim */
  scanInt("mc-nt ",   &mc.nt);    /* n periods in forward sim */
  scanInt("mc-Npk ", &mc.Npk);    /* max nb of new treats per period */
  K = mc.Npk*mc.nt;               // can't have more than mc.Npk new trtmt's per period
  mc.pk= dvector(0,mc.Npk+1);
  scanDoubleArray("mc-pk", mc.pk, mc.Npk+1); /* pk */

  fscanf(ifile," uhat-b0grid %lf %lf %d ", &lo, &hi, &dsg.nb0);
  dsg.b0grid = dvector(0,dsg.nb0);
  if ( (lo<0.1) | (hi<lo) | (hi>0.9) )
    err_msg("init","invalid b0 grid (b0>0 required).\n",0,0,0);
  grid (lo,hi,dsg.nb0,dsg.b0grid);
  fscanf(ifile," uhat-b1grid %lf %lf %d ", &lo, &hi, &dsg.nb1);
  dsg.b1grid = dvector(0,dsg.nb1);
  if ( (lo<0) | (hi<lo) )
    err_msg("init","invalid b1 grid (b1>0 required).\n",0,0,0);
  grid (lo,hi,dsg.nb1,dsg.b1grid);
  fscanf(ifile," uhat-b2grid %lf %lf %d ", &lo, &hi, &dsg.nb2);
  dsg.b2grid = dvector(0,dsg.nb2);
  if ( (lo<0) | (hi<lo) )
    err_msg("init","invalid b2 grid (b2>0 required).\n",0,0,0);
  grid (lo,hi,dsg.nb2,dsg.b2grid);
  fscanf(ifile," uhat-bsave %lf %lf %lf ",&dsg.b0,&dsg.b1,&dsg.b2);

  fclose(ifile);
}

void init_memalloc(){
  
  /* data */
  dta.t0 = ivector(0,K);
  dta.t1 = ivector(0,K);

  
  th.th  = dvector(0,K);

  /* grids */
  dsg.Ugrid = darray3(dsg.nb0,dsg.nb1,dsg.nb2);
  dsg.pBgrid = darray3(dsg.nb0,dsg.nb1,dsg.nb2);
  dsg.n3grid = darray3(dsg.nb0,dsg.nb1,dsg.nb2);
  dsg.nrecogrid = darray3(dsg.nb0,dsg.nb1,dsg.nb2);
  dsg.ntorecogrid = darray3(dsg.nb0,dsg.nb1,dsg.nb2);
  dsg.nrejgrid = darray3(dsg.nb0,dsg.nb1,dsg.nb2);
  dsg.dtgrid = darray3(dsg.nb0,dsg.nb1,dsg.nb2);
  dsg.errIgrid = darray3(dsg.nb0,dsg.nb1,dsg.nb2);
  dsg.errIIgrid = darray3(dsg.nb0,dsg.nb1,dsg.nb2);
  dsg.msegrid = darray3(dsg.nb0,dsg.nb1,dsg.nb2);

  dsg.ngoodgrid = iarray3(dsg.nb0,dsg.nb1,dsg.nb2);
  dsg.nbadgrid = iarray3(dsg.nb0,dsg.nb1,dsg.nb2);
  
}

void init_pars(){
  /* init par VECTORS and latent vars */
  long is1=71374549;
  long is2=90820;

  /* r.v. */
  setseed(is1,is2);

  /* phase III */
  dsg.za = invcdfn(1-dsg.a/2,0,1.0);
  dsg.zb = invcdfn(1-dsg.b,0,1.0);


}


/* *************************************************
   dsg_forw
 ************************************************* */

void dsg_forw(int on_grid)
{ /* forward simulation */
    int i,i0,i1,i2,j,jact,t,
        trtlimit,  //maximum number of treatments at any given time point
        *dk,   /* k[t]-k[t-1] */
        *kt,   // k[t]
        *ktavail, // k[t] discounting treatments for which no patients have been available yet
        **Y,
        *ycum,    //cumulative number of successes
        *trtact,  //list of active treatments
        nact,     //number of active treatments
        navail, //number of active treatments for which data can be collected
        ntot,     //cumulative overall number of patients
        *term;     //terminal decisions
    double a=0,b=0,*m,*s, b0=0, b1=0, b2=0, u, v, *weight;
    FILE *mcfile, *thfile, *trajfile, *termfile;
    double muemp, varemp, *themp;

  dk = ivector(0,mc.nt);
  kt = ivector(0,mc.nt);
  ktavail= ivector(0,mc.nt);

  mcfile = openOut("mc.dsg");
  thfile = openOut("th.dsg");
  trajfile = openOut("traj.dsg");
  termfile = openOut("term.sim");
  
  Uhat_init(on_grid);
  for(i=0; i<mc.nsim; i++){    // loop over forw simulations

    multinomial(mc.nt, mc.Npk+1, mc.pk, dk);   // simulate dk, kt, t0 (start with k.ini treatments)
    dk[0] = dsg.kini;          
    for (j=0;j<dsg.kini;j++) { dta.t0[j] = 0; }
    for (t=1,kt[0]=dsg.kini;t<mc.nt;t++) {
        kt[t] = kt[t-1] + dk[t];
        for (j=kt[t-1];j<kt[t];j++) { dta.t0[j] = t; }
    }

    if (dsg.hyper==1) {       // simulate u~Gamma(th.au,th.bu), v~Gamma(th.av,th.bv), th ~ Be(u,v)
        u = gengam(th.bu,th.au);
        v = gengam(th.bv,th.av);
    } else {                  // used fixed u,v
        u = th.u;
        v = th.v;
    }
    for(j=0;j<kt[mc.nt-1];j++) { th.th[j] = betadev(u,v); }  //simulate parameter values
    Y = imatrix(0,mc.nt,0,kt[mc.nt-1]);        // simulate data
    for (t=0;t<mc.nt;t++) { for (j=0;j<kt[t];j++) { Y[t][j] = binomial(dta.N,th.th[j]); } }

    trtact = ivector(0,kt[mc.nt-1]);  // loop over grid on (b0,b1,b2)
    m = dvector(0,kt[mc.nt-1]);
    s = dvector(0,kt[mc.nt-1]);
    ycum = ivector(0,kt[mc.nt-1]);
    themp = dvector(0,kt[mc.nt-1]);
    weight = dvector(0,kt[mc.nt-1]);
    term = ivector(0,kt[mc.nt-1]);
    for (i0=0;i0<dsg.nb0;i0++) {
        for (i1=(dsg.nb1-1);(i1>=0) && (dsg.b0grid[i0]<=dsg.b1grid[i1]);i1--) {
            for (i2=0;(i2<dsg.nb2)&&(dsg.b0grid[i0]>=dsg.b2grid[i2]);i2++) {
                nact = dsg.kini;
                for (j=0;j<dsg.kini;j++) { trtact[j] = j; }
                for (j=dsg.kini;j<kt[mc.nt-1];j++) { trtact[j] = -1; }
                for (j=0;j<kt[mc.nt-1];j++) { ycum[j] = weight[j]= m[j] = s[j] = 0; themp[j]= 0; }
		for (j=0;j<mc.nt;j++) { ktavail[j]= 0; }
                ntot = dsg.dt = dsg.U = dsg.pB = dsg.n3 = dsg.nreco = dsg.mse = 0;
                dsg.errI = dsg.errII = dsg.ngood= dsg.nbad= 0;
                
                for (t=0;t<mc.nt;t++) {      // loop over time
		  if (t>0) ktavail[t]= ktavail[t-1];
		  if (dsg.patlimit>0) {      // if there's a limit on the nb of patients available
		    trtlimit= dsg.patlimit/dta.N;
		    if (trtlimit<nact) { navail= trtlimit; } else { navail= nact; }
		  } else { navail= nact; }   // if there's no limit on the nb of patients available
                  for(j=0;j<navail;j++) {      //update nb successes & theta estimate for active treats
                   jact = trtact[j];
                   ycum[jact] += Y[t][jact];
		   if (weight[jact]==0) {
		     ktavail[t]+= 1;
		     if (th.th[jact]<dsg.p0) { dsg.nbad++; } else { dsg.ngood++; }
		   }
		   weight[jact]+= dta.N;
                   themp[jact] = (0.0 + ycum[jact])/(weight[jact]);
                   ntot += dta.N;
                  }  // end j
                  if (navail>0) {               
                    if (dsg.pool==1) {       //empirical Bayes estimates
                      if (ktavail[t]>4) {             //only compute them if we have >4 treatments with data
                        muemp = wmeanx(themp,kt[t]-1,weight);
                        varemp = wvarx(themp,kt[t]-1,weight);
                        a = max_xy(.5,muemp*ntot*(muemp*(1-muemp)-varemp)/(varemp*ntot-muemp*(1-muemp)*ktavail[t]));
                        b = max_xy(.5,(1-muemp)*ntot*(muemp*(1-muemp)-varemp)/(varemp*ntot-muemp*(1-muemp)*ktavail[t]));
                        if (a+b>10) { a= 10*a/(a+b); b= 10 - a; } //don't allow prior sample size >10
                      } else {
                        a = th.au/th.bu;
                        b = th.av/th.bv;
                      }
                    } else {                 //use fixed prior
                        a = th.u;
                        b = th.v;
                    }
                  }
                  for (j=0;j<nact;j++) {     //post. moments, stopping decision & utility evaluation
                    jact = trtact[j];
                    m[jact] = (a+ycum[jact])/(a+b+weight[jact]);
                    s[jact] = sqrt((a+ycum[jact])*(b+weight[jact]-ycum[jact])/(a+b+weight[jact]+1))/(a+b+weight[jact]);
                    if (dsg.logsd==1) { if (s[jact]<0.018) { s[jact] = log(0.018); } else { s[jact] = log(s[jact]); } }
                    b0 = dsg.b0grid[i0]; b1 = dsg.b1grid[i1]; b2 = dsg.b2grid[i2];
                    if ((m[jact]<b0+(s[jact]-dsg.s0)*(b2-b0)/(dsg.s1-dsg.s0)) || (m[jact]>b0+(s[jact]-dsg.s0)*(b1-b0)/(dsg.s1-dsg.s0))) {
                      trtact[j] = -1;
                      dta.t1[jact] = t;
                      dsg.dt += weight[jact];
                      dsg.mse += pow(m[jact]-th.th[jact],2);
                      Uhat_b(m[jact],s[jact],weight[jact],th.th[jact],&dsg.U,&dsg.pB,&dsg.n3,&dsg.nreco,&dsg.errI,&dsg.errII,&dsg.term,b0,b1,b2);
                      term[jact] = dsg.term;
                    }
                  }  // end j
                  for(j=0,jact=0;j<nact;j++) {                              // remove -1's from trtact
                    if (trtact[j]!=-1) { trtact[jact] = trtact[j]; jact++; }
                  }
                  for (j=jact;j<nact;j++) { trtact[j] = -1; }
                  nact = jact;
                  for (j=kt[t],jact=nact;j<kt[t+1];j++,jact++) {            // add new active treatments
                      trtact[jact] = j;
                      nact++;
                  }
                  if (i==0 && i0==dsg.b0 && i1==dsg.b1 && i2==dsg.b2) {    // save trajectories for 1 grid point
                      for (j=0;j<kt[t];j++) { fprintf(trajfile,"%d %d %lf %lf %lf \n",t,j,m[j],s[j],th.th[j]); }
                  }
                } // end t
                                

                for (j=0;j<nact;j++) {               // utility of treats active at t=time horizon
                    jact = trtact[j];
                    dta.t1[jact] = mc.nt-1;
		    if (weight[jact]>=1) {           // only take decisions for treatments for which we collected some data
                      dsg.dt +=  weight[jact];
                      dsg.mse += pow(m[jact]-th.th[jact],2);
                      Uhat_b(m[jact],s[jact],weight[jact],th.th[jact],&dsg.U,&dsg.pB,&dsg.n3,&dsg.nreco,&dsg.errI,&dsg.errII,&dsg.term,b0,b1,b2);
                      term[jact] = dsg.term;
		    } else { term[jact]= -1; }
                }

                dsg.Ugrid[i0][i1][i2] += dsg.U;
                dsg.nrecogrid[i0][i1][i2] += dsg.nreco;
                dsg.ntorecogrid[i0][i1][i2] += dsg.dt/dsg.nreco;
                dsg.pBgrid[i0][i1][i2] += dsg.pB;
                dsg.n3grid[i0][i1][i2] += dsg.n3;
                dsg.dtgrid[i0][i1][i2] += dsg.dt;
                dsg.errIgrid[i0][i1][i2] += dsg.errI;
                dsg.errIIgrid[i0][i1][i2] += dsg.errII;
                dsg.nrejgrid[i0][i1][i2] += ktavail[mc.nt-1]-dsg.nreco;
                dsg.msegrid[i0][i1][i2] += dsg.mse/ktavail[mc.nt-1];
		dsg.ngoodgrid[i0][i1][i2] += dsg.ngood;
		dsg.nbadgrid[i0][i1][i2] += dsg.nbad;
                if (i0==dsg.b0 && i1==dsg.b1 && i2==dsg.b2) {
                  for (j=0;j<kt[mc.nt-1];j++) { fprintf(termfile,"%lf %lf %lf %d %d %lf %lf %lf %d %d %d\n",dsg.b0grid[i0],dsg.b1grid[i1],dsg.b2grid[i2],i,j,m[j],s[j],th.th[j],term[j],dta.t0[j],dta.t1[j]); }
                }
            }  // end i2
        }  // end i1
    }  // end i0
    
    fprintf(mcfile,"%d %d \n",i,kt[mc.nt-1]);
    fwriteDoubleArray(thfile,th.th,1,kt[mc.nt-1]);

    free_imatrix(Y,0,mc.nt,0,kt[mc.nt-1]);
    free_ivector(trtact,0,kt[mc.nt-1]);
    free_dvector(m,0,kt[mc.nt-1]);
    free_dvector(s,0,kt[mc.nt-1]);
    free_ivector(ycum,0,kt[mc.nt-1]);
    free_dvector(themp,0,kt[mc.nt-1]);
    free_dvector(weight,0,kt[mc.nt-1]);
    free_ivector(term,0,kt[mc.nt-1]);
  }  // end i
  
  free_ivector(dk,0,mc.nt);
  free_ivector(kt,0,mc.nt);
  fclose(mcfile);
  fclose(thfile);
  fclose(trajfile);
  fclose(termfile);
  
}

/* *************************************************
   U
 ************************************************* */


void Uhat_init(int on_grid)
{
  int i0;

  /* initialize at 0 */
  if (on_grid){
    for(i0=0;i0<dsg.nb0;i0++){
      R_zero(dsg.Ugrid[i0],dsg.nb1,dsg.nb2);
      R_zero(dsg.pBgrid[i0],dsg.nb1,dsg.nb2);
      R_zero(dsg.n3grid[i0],dsg.nb1,dsg.nb2);
      R_zero(dsg.nrecogrid[i0],dsg.nb1,dsg.nb2);
      R_zero(dsg.ntorecogrid[i0],dsg.nb1,dsg.nb2);
      R_zero(dsg.nrejgrid[i0],dsg.nb1,dsg.nb2);
      R_zero(dsg.errIgrid[i0],dsg.nb1,dsg.nb2);
      R_zero(dsg.errIIgrid[i0],dsg.nb1,dsg.nb2);
      R_zero(dsg.msegrid[i0],dsg.nb1,dsg.nb2);
      R_izero(dsg.ngoodgrid[i0],dsg.nb1,dsg.nb2);
      R_izero(dsg.nbadgrid[i0],dsg.nb1,dsg.nb2);
    }
  } 
  dsg.U = dsg.pB = dsg.n3 = dsg.nreco = dsg.ntoreco = dsg.nrej = dsg.errI = dsg.errII = dsg.mse = dsg.ngood = dsg.nbad = 0;
}

void Uhat_scale(int on_grid)
{ int i0,i1,i2;
  double r; 

  /* get average observered utility */
  if (on_grid){
    for(i0=0;i0<dsg.nb0;i0++){
      for(i1=0;i1<dsg.nb1;i1++) // scale by total # nreco
	for(i2=0;i2<dsg.nb2;i2++){
          dsg.errIgrid[i0][i1][i2] *= 1.0/dsg.nbadgrid[i0][i1][i2];
          dsg.errIIgrid[i0][i1][i2] *= 1.0/dsg.ngoodgrid[i0][i1][i2];
          if (dsg.nrecogrid[i0][i1][i2]>0){
	    r = 1.0/dsg.nrecogrid[i0][i1][i2];   // scale by nsim
	    dsg.pBgrid[i0][i1][i2] *= r;
	    dsg.n3grid[i0][i1][i2] *= r;
          } else {
	    dsg.pBgrid[i0][i1][i2] = 0;
	    dsg.n3grid[i0][i1][i2] = 0;
          } 
          dsg.dtgrid[i0][i1][i2] *= 1.0/(dsg.nrecogrid[i0][i1][i2]+dsg.nrejgrid[i0][i1][i2]);
        }
      r = 1.0/mc.nsim;          // scale by nsim
      rA(r,dsg.Ugrid[i0],dsg.Ugrid[i0],        dsg.nb1,dsg.nb2);
      rA(r,dsg.nrecogrid[i0],dsg.nrecogrid[i0],dsg.nb1,dsg.nb2);
      rA(r,dsg.ntorecogrid[i0],dsg.ntorecogrid[i0],dsg.nb1,dsg.nb2);
      rA(r,dsg.nrejgrid[i0],dsg.nrejgrid[i0],dsg.nb1,dsg.nb2);
      rA(r,dsg.msegrid[i0],dsg.msegrid[i0],dsg.nb1,dsg.nb2);
    }// i0
  } // on_grid
}


void Uhat_traj(double **M, double **S, int kmax, int i, FILE *trajfile)
{ /* saves trajectories to file for later plotting */
  int j, t;

  /* evaluate utility on grid over dsgs */
  for(j=0;j<kmax;j++)
    for(t=dta.t0[j];t<mc.nt;t++)
      fprintf(trajfile,"%d %d %d %lf %lf %lf\n",
	      i,j,t,M[t][j],S[t][j],th.th[j]);
}

void Uhat_b(double m,double logs,double weight,double th,double *U,double *pB,double *n,double *nreco,int *errI,int *errII,int *term,double b0,double b1,double b2)
{ 
  int n3j;
  double pBj,uterm,s;
  
    if (dsg.uth==2) {    //phaseIII term dropped from utility
        uterm = 0;
        if (m<b0+(logs-dsg.s0)*(b2-b0)/(dsg.s1-dsg.s0)) { *term = 0; if (th>=dsg.p0) { *errII += 1; } }
        else if (m>b0+(logs-dsg.s0)*(b1-b0)/(dsg.s1-dsg.s0)) { *term=1; *nreco += 1; if (th<dsg.p0) { *errI += 1; } }
        else {
            if (m<dsg.p0) { *term=0; if (th>=dsg.p0) { *errII += 1; } }
            else { *term=1; *nreco+=1; if (th<dsg.p0) { *errI += 1; } }
        }
    } else {                //phaseIII term not dropped from utility
        if (dsg.logsd==1) { s = exp(logs); } else { s = logs; }
        uterm = uobs(m,s,&pBj,&n3j,th,term);
        if (*term == 1){ // go to phase III
              *nreco += 1;
              *pB += pBj;
              *n += 1.0*n3j;
              if (th<dsg.p0) { *errI += 1; }
        } else {
            uterm = 0;
            if (th>=dsg.p0) { *errII += 1; }
        }
    }
    *U += uterm -weight*dsg.c1;
}



void Uhat_write(int on_grid)
{
  int i0,n0,n1,n2;
  FILE *Ufile,*b0file, *b1file, *b2file, *nrfile, *ntorfile, *pBfile, *n3file, *dtfile, *errIfile, *errIIfile, *msefile; 

  if (!on_grid)
    return;

  n0 = dsg.nb0;
  n1 = dsg.nb1;
  n2 = dsg.nb2;
  Ufile = openOut("U.sim");
  b0file = openOut("b0.sim");
  b1file = openOut("b1.sim");
  b2file = openOut("b2.sim");
  nrfile = openOut("nreco.sim");
  ntorfile = openOut("ntoreco.sim");
  pBfile = openOut("pB.sim");
  n3file = openOut("n3.sim");
  dtfile = openOut("dt.sim");
  errIfile = openOut("errI.sim");
  errIIfile = openOut("errII.sim");
  msefile = openOut("mse.sim");
  
  fwriteDoubleArray(b0file,dsg.b0grid,1,n0);
  fwriteDoubleArray(b1file,dsg.b1grid,1,  n1);
  fwriteDoubleArray(b2file,dsg.b2grid,1,n2);
  
  for(i0=0;i0<n0;i0++){
    fwriteDoubleMatrix2(Ufile,dsg.Ugrid[i0],n1,n2);
    fwriteDoubleMatrix2(nrfile,dsg.nrecogrid[i0],n1,n2);
    fwriteDoubleMatrix2(ntorfile,dsg.ntorecogrid[i0],n1,n2);
    fwriteDoubleMatrix2(pBfile,dsg.pBgrid[i0],n1,n2);
    fwriteDoubleMatrix2(n3file,dsg.n3grid[i0],n1,n2);
    fwriteDoubleMatrix2(dtfile,dsg.dtgrid[i0],n1,n2);
    fwriteDoubleMatrix2(errIfile,dsg.errIgrid[i0],n1,n2);
    fwriteDoubleMatrix2(errIIfile,dsg.errIIgrid[i0],n1,n2);
    fwriteDoubleMatrix2(msefile,dsg.msegrid[i0],n1,n2);
  }

  fclose(Ufile);
  fclose(b0file);
  fclose(b1file);
  fclose(b2file);
  fclose(nrfile);
  fclose(ntorfile);
  fclose(pBfile);
  fclose(n3file);
  fclose(dtfile);
  fclose(errIfile);
  fclose(errIIfile);
  fclose(msefile);
}
    

/* *************************************************
   uobs
 ************************************************* */
double uobs(double m, double s, double *pB, int *N3, double th, int *term)
{ // returns expected utility for treatment with final moments (mi,si), true success prob th
  int n;
  double Edp,u2,mbar,Fcut,mtstat,sdtstat;

  if (m <= dsg.p0)
    u2 = *pB = *N3 = *term = 0;
  else {
    n = n2(m,s);
    if (n>10000) {       //if phase III sample size is too large
        u2 = *pB = *N3 = *term = 0;
    } else {            //if phase III sample size is not too large
        mbar = 0.5*(dsg.p0+m);
        Fcut = dsg.za*sqrt(2*mbar*(1-mbar)/n);
        mtstat = m - dsg.p0;
        sdtstat = sqrt(m*(1-m)/n+s*s*(1.0-1.0/n)+dsg.p0*(1-dsg.p0)/n);
        *pB = 1.0 - cdfnorm(Fcut,mtstat,sdtstat);

        if (dsg.uth==0) {
            Edp = mtstat + sdtstat*exp(-pow((Fcut-mtstat)/(sdtstat*sqrt(2)),2))/(sqrt(2*M_PI)*(*pB));
            if (Edp>1) Edp = 1;   //difference in proportions cannot be greater than 1
            u2  = (-n)*dsg.c1 + (*pB)*(Edp*dsg.c2) - dsg.c0;
            if (u2>0) { *term=1; } else { *term=0; }
        }
        else {
            Edp = th - dsg.p0;
            u2  = (-n)*dsg.c1 + (*pB)*(Edp*dsg.c2) - dsg.c0;
            if (n*dsg.c1<(*pB)*((m-dsg.p0)*dsg.c2) - dsg.c0) { *term=1; } else { *term=0; }
        }

        *N3  = n;
    }
  }
  return(u2);
}

int n2(double m, double s)
{ // sample size for phase III
  // (zb*sqrt(m*(1-m)+th0*(1-th0))-za*sqrt(2*mbar*(1-mbar)))^2 / (th0-m)^2
  // ,where mbar=0.5*(dsg.p0+m)
  int n;
  double mbar;
  
  mbar = 0.5*(dsg.p0+m);
  n = pow((dsg.zb*sqrt(m*(1-m)+dsg.p0*(1-dsg.p0))+dsg.za*sqrt(2*mbar*(1-mbar)))/(dsg.p0-m), 2);
  return(n);
}
  
