#enter the following
nburn<-100;  #burn in iterations
nsave<-1;  #1 if a machine readable output is to be written, 0 otherwise
apmout<-"momout1"; #name of the machine readable output if nsave=1
hfile<-"sim";  #input data file
outfile<-"momout"; #name of the output file
jf<-c(1,2); #functions of interest


#do not write anything below
NG<-100;
LG<-1000;
ntaper<-c(4,8,15);

simp <-
function(x,y,lag)
{
   if (lag==0)
      sum(x*y)
   else
      sum(x[-(0:lag)]*y[-((length(y)-lag+1): (length(y)+1))])
}

gen<-
function(x,y)
{
   z_1:length(x);
   for (i in 1:length(x))
      z[i]_simp(x,y,i-1);
   z
}        


#script starts here
if (nburn<0) nburn_0; 
cat(file=outfile,'\n');

nrun_0;
in1_scan(hfile);
niter<-in1[1];
nf<-in1[2];
if (length(in1) != (niter*(4+nf) + 2))
{
  cat(file=outfile,append=T, 
      "unexpected end of simulation input file reached at line ",
      (niter*nf +2),"\n");
  q();
};
in1_in1[-(1:2)];
dim(in1)_c(nf+4,niter);     #reshape in1

errcheck<-
function(x) {
  err<-0;
  if (x <= 0) 
  { 
    cat(file=outfile,append=T, 'function of interest negative; skipped\n'); 
    err<-1;
  }
  else
  { 
     if (x > nf)
     {
        cat(file=outfile,append=T,'function of interest greater than ',jf);
        cat(file=outfile,append=T,'limit of , skipped',nf);
        err<-1;
     }
     else if (nburn>niter)
     {
        cat(file=outfile,append=T,'burn in requested but only ',nburn);
        cat(file=outfile,append=T,'avaliable', niter);
        q();
     }
     else if (nf>LG)
     {
        cat(file=outfile,append=T,'file contains functions per record, ',nf);
        cat(file=outfile,append=T,'exceeds storage capacity of ', LG);
        q();
      } 
  }
    err
} #errcheck 

  #if (!err)
mom<-
function(x,nrun)
  {
     nuse<-niter-nburn;
     if (nuse<NG)
     {
        cat(file=outfile, append=T,'minimum analysis iterations are required; ',NG);
        cat(file=outfile, append=T,' avaliable \n', nuse);
        q();
     }
     ns<-floor(nuse/NG);
     nuse<-ns*NG;
     if (nrun==1) 
     {
        cat(file=outfile, append=T,'Posterior moment analysis \n');
        cat(file=outfile, append=T,'Total iterations available: ',niter, '\n');
        cat(file=outfile, append=T, '        Burn-in iterations: ',nburn,'\n');
        cat(file=outfile, append=T, '       Analysis iterations: ',nuse,'\n');                   
     } 
   # form sufficiency statistics for all the subsequent work
wtlog_in1[2,];              #extract wtlogs 
f_in1[4+x,];               #extract function of interest
#skip the burn in iterations:
#f_f[-(1:nburn)]; wtlog_wtlog[-(1:nburn)];
f_f[((nburn+1):(nburn+nuse))]; wtlog_wtlog[((nburn+1):(nburn+nuse))];
wtlog[wtlog>-700]_exp(wtlog[wtlog>-700]);
wtlog[wtlog<=-700]_0;
dgd_wtlog;  #can be removed, need changes below
dgn_wtlog*f;
dgdn_f*wtlog^2;
dgdd_wtlog^2;
dgnn_dgn^2;
dgvar_wtlog*f^2;

#form posterior moment and standard deviation, an print
tn_sum(dgn) 
td_sum(wtlog); 
eg<-tn/td

tvar_sum(dgvar); 
varg<-tvar/td -eg^2
sdg<--1; if (varg>0) sdg<-sqrt(varg); 

cat(file=outfile,append=T,'\n\nFunction ',x,'Posterior mean:', eg,"\n");
cat(file=outfile,append=T,'             Posterior st dev:',sdg, "\n");
cat(file=outfile,append=T, '      Accuracy of approximation:');
cat(file=outfile,append=T, '        Numerical    Relative\n');
cat(file=outfile,append=T, '                                ');
cat(file=outfile,append=T, '        Standard     Numerical\n');
cat(file=outfile,append=T, 'Method--                        ');
cat(file=outfile,append=T, '        Error        Efficiency\n');
                                                                              

if (nsave)
  cat(file=apmout,append=T,x,eg,nuse,sdg,"\n");
#numerical standard error assuming no serial correlation

tnn_sum(dgnn); tdn_sum(dgdn); tdd_sum(dgdd);
varnum<-(tnn-2*eg*tdn+tdd*eg^2)/(td^2);
sdnum<--1; if (varnum>0) sdnum<-sqrt(varnum);  
#header2(file=outfile);
cat(file=outfile,append=T,'Assuming no serial correlation:    '); 
cat(file=outfile,append=T, sdnum, varg/(nuse*varnum), "\n");
if (nsave)
  cat(file=apmout,append=T,sdnum,"  ");

#get autocovariance of grouped means
barn<-tn/nuse;
bard<-td/nuse;

dim(dgd)_c(ns,NG); dim(dgn)_c(ns,NG);
dim(dgdn)_c(ns,NG); dim(dgdd)_c(ns,NG);
dim(dgnn)_c(ns,NG); dim(dgvar)_c(ns,NG);
a_rep(1,ns);
cd_(a %*% dgd)/ns-bard; cn_(a %*% dgn)/ns-barn;
#cdn_(a %*% dgdn)/ns; cdd_(a %*% dgdd)/ns;
#cnn_(a %*% dgnn)/ns; cvar_(a %*% dgvar)/ns;

ann_gen(cn,cn); add_gen(cd,cd); and_gen(cn,cd); adn_gen(cd,cn);
#to save memory, rnn is ann and so on
rnn_ann/NG; rdd_add/NG; rnd_and/NG; rdn_adn/NG;

# numerical standard error with tapered autocovariance functions
for (mm in 1:3)
{
   m<-ntaper[mm];
   am<-m;
   snn<-rnn[1]; sdd<-rdd[1]; snd<-rnd[1];
   for (lag in 1:(m-1))
   {
       att<-1-lag/am;
       snn<-snn+2*att*rnn[lag+1];
       sdd<-sdd+2*att*rdd[lag+1];
       snd<-snd+att*(rdn[lag+1] + rnd[lag+1]); 
    } #lag
        varnum<-ns*nuse*(snn-2*eg*snd+sdd*eg^2)/(td^2);
        sdnum<--1; if (varnum>0) sdnum<-sqrt(varnum); 
        c<-'%';
        cat(file=outfile,append=T,'Autocovariance function tapered to ', m,c);
        cat(file=outfile,append=T,sdnum, varg/(nuse*varnum),'\n'); 
        if (nsave)
           cat(file=apmout,append=T,sdnum,"  ");
      } #mm
     if (nsave)
        cat(file=apmout,append=T,"\n");
} #mom 


for (i in 1:length(jf))
 if (!errcheck(jf[i]))
{ 
  nrun<-nrun+1;
  mom(jf[i],nrun);
}



if (nrun==0)
{
  cat(file=outfile,append=T, 
      'error on input for function number, program stops'); 
  q();
}
#clean up
remove(objects());
q();