/*-----------------------------------------------------------------------------

Copyright (C) 2005, 2006, 2008.

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA   

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-----------------------------------------------------------------------------*/
#include "libsnp.h"
#include "libgsm.h"
#include "gsm.h"

using namespace scl;
using namespace libgsm;
using namespace gsm;
using namespace std;

bool gsm::datblock::read_data(realmat& data)
{
   string pathname = "";
   return this->read_data(pathname,data);
}

bool gsm::datblock::read_data(const string& pathname, realmat& data)
{
  data.resize(M, n);

  ifstream* dat_ifs_ptr;

  if (dsn[0] == '/') {
    dat_ifs_ptr = new ifstream(dsn.c_str());
    if (dat_ifs_ptr == 0 || !*dat_ifs_ptr) return false;
  }
  else {
    string filename = pathname + dsn;
    dat_ifs_ptr = new ifstream(filename.c_str());
    if (dat_ifs_ptr == 0 || !*dat_ifs_ptr) return false;
  }
  
  ifstream& dat_ifs = *dat_ifs_ptr;

  INTEGER max=0; 
  for (INTEGER i=1; i<=fields.size(); ++i) {
    max = fields[i] > max ? fields[i] : max ;
  }
  intvec idx(max,0);
  for (INTEGER i=1; i<=fields.size(); ++i) idx[fields[i]] = i;
  
  string discard;

  for (INTEGER t=1; t<=n; ++t) {
    for (INTEGER j=1; j<=max; ++j) {
      if (idx[j]==0) {
        dat_ifs >> discard;
      }
      else {
        dat_ifs >> data(idx[j],t);
      }
    }
    getline(dat_ifs,discard);
  }

  return dat_ifs.good();
}


void gsm::sci_summary::update(const realmat& sci_parm_chain, 
           const realmat& sci_func_chain, 
           const realmat& stat_parm_chain, 
           const realmat& stat_func_chain,
           const sci_val& sci_mode)
{
  mode = sci_mode;  // Assumed called with overall mode from implied_map

  INTEGER len_sci_parm = sci_parm_chain.nrow();
  INTEGER len_sci_func = sci_func_chain.nrow();
  INTEGER len_stat_parm = stat_parm_chain.nrow();
  INTEGER len_stat_func = stat_func_chain.nrow();
  INTEGER n = sci_parm_chain.ncol();

  if (N == 0) {
    sci_parm_sum.resize(len_sci_parm,1,0.0);
    sci_func_sum.resize(len_sci_func,1,0.0);
    stat_parm_sum.resize(len_stat_parm,1,0.0);
    stat_func_sum.resize(len_stat_func,1,0.0);

    sci_parm_ctr.resize(len_sci_parm,1,0.0);
    sci_func_ctr.resize(len_sci_func,1,0.0);
    stat_parm_ctr.resize(len_stat_parm,1,0.0);
    stat_func_ctr.resize(len_stat_func,1,0.0);

    sci_parm_sse.resize(len_sci_parm,len_sci_parm,0.0);
    sci_func_sse.resize(len_sci_func,len_sci_func,0.0);
    stat_parm_sse.resize(len_stat_parm,len_stat_parm,0.0);
    stat_func_sse.resize(len_stat_func,len_stat_func,0.0);
  }

  if (N == 0) { 
    for (INTEGER t=1; t<=n; ++t) {
      for (INTEGER i=1; i<=len_sci_parm; ++i) {
        sci_parm_ctr[i] += sci_parm_chain(i,t);
      }
      for (INTEGER i=1; i<=len_sci_func; ++i) {
        sci_func_ctr[i] += sci_func_chain(i,t);
      }
      for (INTEGER i=1; i<=len_stat_parm; ++i) {
        stat_parm_ctr[i] += stat_parm_chain(i,t);
      }
      for (INTEGER i=1; i<=len_stat_func; ++i) {
        stat_func_ctr[i] += stat_func_chain(i,t);
      }
    }
    sci_parm_ctr = sci_parm_ctr/n;
    sci_func_ctr = sci_func_ctr/n;
    stat_parm_ctr = stat_parm_ctr/n;
    stat_func_ctr = stat_func_ctr/n;
  }

  for (INTEGER t=1; t<=n; ++t) {
    for (INTEGER i=1; i<=len_sci_parm; ++i) {
      sci_parm_sum[i] += (sci_parm_chain(i,t) - sci_parm_ctr[i]);
    }
    for (INTEGER i=1; i<=len_sci_func; ++i) {
      sci_func_sum[i] += (sci_func_chain(i,t) - sci_func_ctr[i]);
    }
    for (INTEGER i=1; i<=len_stat_parm; ++i) {
      stat_parm_sum[i] += (stat_parm_chain(i,t) - stat_parm_ctr[i]);
    }
    for (INTEGER i=1; i<=len_stat_func; ++i) {
      stat_func_sum[i] += (stat_func_chain(i,t) - stat_func_ctr[i]);
    }
  }

  for (INTEGER t=1; t<=n; ++t) {
    for (INTEGER j=1; j<=len_sci_parm; ++j) {
      for (INTEGER i=1; i<=len_sci_parm; ++i) {
        sci_parm_sse(i,j) += (sci_parm_chain(i,t) - sci_parm_ctr[i])
                               *(sci_parm_chain(j,t) - sci_parm_ctr[j]);
      }
    }
    for (INTEGER j=1; j<=len_sci_func; ++j) {
      for (INTEGER i=1; i<=len_sci_func; ++i) {
        sci_func_sse(i,j) += (sci_func_chain(i,t) - sci_func_ctr[i])
                               *(sci_func_chain(j,t) - sci_func_ctr[j]);
      }  
    }
    for (INTEGER j=1; j<=len_stat_parm; ++j) {
      for (INTEGER i=1; i<=len_stat_parm; ++i) {
        stat_parm_sse(i,j) += (stat_parm_chain(i,t) - stat_parm_ctr[i])
                                 *(stat_parm_chain(j,t) - stat_parm_ctr[j]);
      }
    }
    for (INTEGER j=1; j<=len_stat_func; ++j) {
      for (INTEGER i=1; i<=len_stat_func; ++i) {
        stat_func_sse(i,j) += (stat_func_chain(i,t) - stat_func_ctr[i])
                                 *(stat_func_chain(j,t) - stat_func_ctr[j]);
      }
    }
  }
    
  N += n;
}    

bool gsm::sci_summary::write(string prefix)
{
  if (N == 0) error("Error, sci_summary, write, N == 0");

  INTEGER len_sci_parm = sci_parm_sum.nrow();
  INTEGER len_sci_func = sci_func_sum.nrow();
  INTEGER len_stat_parm = stat_parm_sum.nrow();
  INTEGER len_stat_func = stat_func_sum.nrow();

  realmat sci_parm_mean = sci_parm_sum/N;
  realmat sci_func_mean = sci_func_sum/N;
  realmat stat_parm_mean = stat_parm_sum/N;
  realmat stat_func_mean = stat_func_sum/N;

  realmat sci_parm_var = sci_parm_sse/N;
  realmat sci_func_var = sci_func_sse/N;
  realmat stat_parm_var = stat_parm_sse/N;
  realmat stat_func_var = stat_func_sse/N;

  for (INTEGER j=1; j<=len_sci_parm; ++j) {
    for (INTEGER i=1; i<=len_sci_parm; ++i) {
      sci_parm_var(i,j) -= sci_parm_mean[i]*sci_parm_mean[j];
    }
  }
  for (INTEGER j=1; j<=len_sci_func; ++j) {
    for (INTEGER i=1; i<=len_sci_func; ++i) {
      sci_func_var(i,j) -= sci_func_mean[i]*sci_func_mean[j];
    }  
  }
  for (INTEGER j=1; j<=len_stat_parm; ++j) {
    for (INTEGER i=1; i<=len_stat_parm; ++i) {
      stat_parm_var(i,j) -= stat_parm_mean[i]*stat_parm_mean[j];
    }
  }
  for (INTEGER j=1; j<=len_stat_func; ++j) {
    for (INTEGER i=1; i<=len_stat_func; ++i) {
      stat_func_var(i,j) -= stat_func_mean[i]*stat_func_mean[j];
    }
  }
    
  REAL scale = REAL(N)/REAL(N-1);
  sci_parm_var = scale*sci_parm_var;
  sci_func_var = scale*sci_func_var;
  stat_parm_var = scale*stat_parm_var;
  stat_func_var = scale*stat_func_var;

  sci_parm_mean += sci_parm_ctr;
  sci_func_mean += sci_func_ctr;
  stat_parm_mean += stat_parm_ctr;
  stat_func_mean += stat_func_ctr;

  realmat sci_parm_sdev(len_sci_parm,1);
  realmat sci_func_sdev(len_sci_func,1);
  realmat stat_parm_sdev(len_stat_parm,1);
  realmat stat_func_sdev(len_stat_func,1);
  
  for (INTEGER i=1; i<=len_sci_parm; ++i) {
    sci_parm_sdev[i] = sqrt(sci_parm_var(i,i));
  }
  for (INTEGER i=1; i<=len_sci_func; ++i) {
    sci_func_sdev[i] = sqrt(sci_func_var(i,i));
  }
  for (INTEGER i=1; i<=len_stat_parm; ++i) {
    stat_parm_sdev[i] = sqrt(stat_parm_var(i,i));
  }
  for (INTEGER i=1; i<=len_stat_func; ++i) {
    stat_func_sdev[i] = sqrt(stat_func_var(i,i));
  }

  prefix += ".summary.";
  string filename;

  filename = prefix + "sci_parm_mode.dat";
  vecwrite(filename.c_str(),mode.sci_mod_parm);
  filename = prefix + "sci_func_mode.dat";
  vecwrite(filename.c_str(),mode.sci_mod_func);
  filename = prefix + "stat_parm_mode.dat";
  vecwrite(filename.c_str(),mode.stat_mod_parm);
  filename = prefix + "stat_func_mode.dat";
  vecwrite(filename.c_str(),mode.stat_mod_func);

  filename = prefix + "sci_parm_mean.dat";
  vecwrite(filename.c_str(),sci_parm_mean);
  filename = prefix + "sci_func_mean.dat";
  vecwrite(filename.c_str(),sci_func_mean);
  filename = prefix + "stat_parm_mean.dat";
  vecwrite(filename.c_str(),stat_parm_mean);
  filename = prefix + "stat_func_mean.dat";
  vecwrite(filename.c_str(),stat_func_mean);

  filename = prefix + "sci_parm_var.dat";
  vecwrite(filename.c_str(),sci_parm_var);
  filename = prefix + "sci_func_var.dat";
  vecwrite(filename.c_str(),sci_func_var);
  filename = prefix + "stat_parm_var.dat";
  vecwrite(filename.c_str(),stat_parm_var);
  filename = prefix + "stat_func_var.dat";
  vecwrite(filename.c_str(),stat_func_var);

  filename = prefix + "sci.dat";
  ofstream fout(filename.c_str());
  if (!fout) return false;

  fout << "     "
       << "                mode"
       << "                mean"
       << "                sdev"  << '\n'; 

  fout << "sci_parm" << '\n';
  for (INTEGER i=1; i<=len_sci_parm; ++i) {
    fout << fmt('d',5,i);
    fout << fmt('f',20,8, mode.sci_mod_parm[i]);
    fout << fmt('f',20,8, sci_parm_mean[i]);
    fout << fmt('f',20,8, sci_parm_sdev[i]) << '\n';
  }

  fout << "sci_func" << '\n';
  for (INTEGER i=1; i<=len_sci_func; ++i) {
    fout << fmt('d',5,i);
    fout << fmt('f',20,8, mode.sci_mod_func[i]);
    fout << fmt('f',20,8, sci_func_mean[i]);
    fout << fmt('f',20,8, sci_func_sdev[i]) << '\n';
  }

  fout << "stat_parm" << '\n';
  for (INTEGER i=1; i<=len_stat_parm; ++i) {
    fout << fmt('d',5,i);
    fout << fmt('f',20,8, mode.stat_mod_parm[i]);
    fout << fmt('f',20,8, stat_parm_mean[i]);
    fout << fmt('f',20,8, stat_parm_sdev[i]) << '\n';
  }

  fout << "stat_func" << '\n';
  for (INTEGER i=1; i<=len_stat_func; ++i) {
    fout << fmt('d',5,i);
    fout << fmt('f',20,8, mode.stat_mod_func[i]);
    fout << fmt('f',20,8, stat_func_mean[i]);
    fout << fmt('f',20,8, stat_func_sdev[i]) << '\n';
  }

  fout << '\n';
  fout << "log sci_mod_prior       = " 
       << fmt('f',20,8,mode.sci_mod_prior.log_den) << '\n';
  fout << "log stat_mod_prior      = " 
       << fmt('f',20,8,mode.stat_mod_prior.log_den) << '\n';
  fout << "log stat_mod_likelihood = " 
       << fmt('f',20,8,mode.stat_mod_logl.log_den) << '\n';
  if (est_blk.sci_draw_from_prior) {
    fout << "Draws are from prior, no posterior printed" << '\n';
  }
  else {
    REAL log_post = mode.sci_mod_prior.log_den;
    log_post += mode.stat_mod_prior.log_den;
    log_post += mode.stat_mod_logl.log_den;
    fout << "log sci_mod_posterior   = " << fmt('f',20,8,log_post) << '\n';
  }

  fout << '\n';
  fout << '\n';
  fout << "     "
       << "                mode"
       << "                mean"
       << "                sdev"  << '\n'; 

  fout << "sci_parm" << '\n';
  for (INTEGER i=1; i<=len_sci_parm; ++i) {
    fout << fmt('d',5,i);
    fout << fmt('e',20,10, mode.sci_mod_parm[i]);
    fout << fmt('e',20,10, sci_parm_mean[i]);
    fout << fmt('e',20,10, sci_parm_sdev[i]) << '\n';
  }

  fout << "sci_func" << '\n';
  for (INTEGER i=1; i<=len_sci_func; ++i) {
    fout << fmt('d',5,i);
    fout << fmt('e',20,10, mode.sci_mod_func[i]);
    fout << fmt('e',20,10, sci_func_mean[i]);
    fout << fmt('e',20,10, sci_func_sdev[i]) << '\n';
  }

  fout << "stat_parm" << '\n';
  for (INTEGER i=1; i<=len_stat_parm; ++i) {
    fout << fmt('d',5,i);
    fout << fmt('e',20,10, mode.stat_mod_parm[i]);
    fout << fmt('e',20,10, stat_parm_mean[i]);
    fout << fmt('e',20,10, stat_parm_sdev[i]) << '\n';
  }

  fout << "stat_func" << '\n';
  for (INTEGER i=1; i<=len_stat_func; ++i) {
    fout << fmt('d',5,i);
    fout << fmt('e',20,10, mode.stat_mod_func[i]);
    fout << fmt('e',20,10, stat_func_mean[i]);
    fout << fmt('e',20,10, stat_func_sdev[i]) << '\n';
  }

  fout << '\n';
  fout << "log sci_mod_prior       = " 
       << fmt('e',20,10,mode.sci_mod_prior.log_den) << '\n';
  fout << "log stat_mod_prior      = " 
       << fmt('e',20,10,mode.stat_mod_prior.log_den) << '\n';
  fout << "log stat_mod_likelihood = " 
       << fmt('e',20,10,mode.stat_mod_logl.log_den) << '\n';
  if (est_blk.sci_draw_from_prior) {
    fout << "Draws are from prior, no posterior printed" << '\n';
  }
  else {
    REAL log_post = mode.sci_mod_prior.log_den;
    log_post += mode.stat_mod_prior.log_den;
    log_post += mode.stat_mod_logl.log_den;
    fout << "log sci_mod_posterior   = " << fmt('e',20,10,log_post) << '\n';
  }

  fout.clear();
  fout.close();
  return true;
}

void gsm::assess_summary::update(const realmat& assess_stat_parm_chain, 
           const realmat& assess_stat_func_chain,
           const stat_val& stat_mode)
{
  bool draw_from_posterior = !est_blk.stat_draw_from_prior;

  if (N == 0) {
    mode = stat_mode;
  }
  else {
    REAL log_pi_old = mode.stat_mod_prior.log_den;
    log_pi_old += mode.assess_prior.log_den;
    if (draw_from_posterior) {
      log_pi_old += mode.stat_mod_logl.log_den;
    }
    REAL log_pi_new = stat_mode.stat_mod_prior.log_den;
    log_pi_new += stat_mode.assess_prior.log_den;
    if (draw_from_posterior) {
      log_pi_new += stat_mode.stat_mod_logl.log_den;
    }
    if (log_pi_new > log_pi_old) {
      mode = stat_mode;
    }
  }

  INTEGER len_assess_stat_parm = assess_stat_parm_chain.nrow();
  INTEGER len_assess_stat_func = assess_stat_func_chain.nrow();
  INTEGER n = assess_stat_parm_chain.ncol();

  if (N == 0) {
    assess_stat_parm_sum.resize(len_assess_stat_parm,1,0.0);
    assess_stat_func_sum.resize(len_assess_stat_func,1,0.0);

    assess_stat_parm_ctr.resize(len_assess_stat_parm,1,0.0);
    assess_stat_func_ctr.resize(len_assess_stat_func,1,0.0);

    assess_stat_parm_sse.resize(len_assess_stat_parm,len_assess_stat_parm,0.0);
    assess_stat_func_sse.resize(len_assess_stat_func,len_assess_stat_func,0.0);
  }

  if (N == 0) { 
    for (INTEGER t=1; t<=n; ++t) {
      for (INTEGER i=1; i<=len_assess_stat_parm; ++i) {
        assess_stat_parm_ctr[i] += assess_stat_parm_chain(i,t);
      }
      for (INTEGER i=1; i<=len_assess_stat_func; ++i) {
        assess_stat_func_ctr[i] += assess_stat_func_chain(i,t);
      }
    }
    assess_stat_parm_ctr = assess_stat_parm_ctr/n;
    assess_stat_func_ctr = assess_stat_func_ctr/n;
  }

  for (INTEGER t=1; t<=n; ++t) {
    for (INTEGER i=1; i<=len_assess_stat_parm; ++i) {
      assess_stat_parm_sum[i] 
        += (assess_stat_parm_chain(i,t) - assess_stat_parm_ctr[i]);
    }
    for (INTEGER i=1; i<=len_assess_stat_func; ++i) {
      assess_stat_func_sum[i] 
        += (assess_stat_func_chain(i,t) - assess_stat_func_ctr[i]);
    }
  }

  for (INTEGER t=1; t<=n; ++t) {
    for (INTEGER j=1; j<=len_assess_stat_parm; ++j) {
      for (INTEGER i=1; i<=len_assess_stat_parm; ++i) {
        assess_stat_parm_sse(i,j) 
          += (assess_stat_parm_chain(i,t) - assess_stat_parm_ctr[i])
               *(assess_stat_parm_chain(j,t) - assess_stat_parm_ctr[j]);
      }
    }
    for (INTEGER j=1; j<=len_assess_stat_func; ++j) {
      for (INTEGER i=1; i<=len_assess_stat_func; ++i) {
        assess_stat_func_sse(i,j) 
          += (assess_stat_func_chain(i,t) - assess_stat_func_ctr[i])
               *(assess_stat_func_chain(j,t) - assess_stat_func_ctr[j]);
      }
    }
  }
    
  N += n;
}    

bool gsm::assess_summary::write(string prefix)
{
  if (N == 0) error("Error, assess_summary, write, N == 0");

  INTEGER len_assess_stat_parm = assess_stat_parm_sum.nrow();
  INTEGER len_assess_stat_func = assess_stat_func_sum.nrow();

  realmat assess_stat_parm_mean = assess_stat_parm_sum/N;
  realmat assess_stat_func_mean = assess_stat_func_sum/N;

  realmat assess_stat_parm_var = assess_stat_parm_sse/N;
  realmat assess_stat_func_var = assess_stat_func_sse/N;

  for (INTEGER j=1; j<=len_assess_stat_parm; ++j) {
    for (INTEGER i=1; i<=len_assess_stat_parm; ++i) {
      assess_stat_parm_var(i,j) 
        -= assess_stat_parm_mean[i]*assess_stat_parm_mean[j];
    }
  }
  for (INTEGER j=1; j<=len_assess_stat_func; ++j) {
    for (INTEGER i=1; i<=len_assess_stat_func; ++i) {
      assess_stat_func_var(i,j) 
        -= assess_stat_func_mean[i]*assess_stat_func_mean[j];
    }
  }
    
  REAL scale = REAL(N)/REAL(N-1);
  assess_stat_parm_var = scale*assess_stat_parm_var;
  assess_stat_func_var = scale*assess_stat_func_var;

  assess_stat_parm_mean += assess_stat_parm_ctr;
  assess_stat_func_mean += assess_stat_func_ctr;

  realmat assess_stat_parm_sdev(len_assess_stat_parm,1);
  realmat assess_stat_func_sdev(len_assess_stat_func,1);
  
  for (INTEGER i=1; i<=len_assess_stat_parm; ++i) {
    assess_stat_parm_sdev[i] = sqrt(assess_stat_parm_var(i,i));
  }
  for (INTEGER i=1; i<=len_assess_stat_func; ++i) {
    assess_stat_func_sdev[i] = sqrt(assess_stat_func_var(i,i));
  }

  prefix += ".summary.";
  string filename;

  filename = prefix + "as_st_parm_mode.dat";
  vecwrite(filename.c_str(),mode.stat_mod_parm);
  filename = prefix + "as_st_func_mode.dat";
  vecwrite(filename.c_str(),mode.stat_mod_func);

  filename = prefix + "as_st_parm_mean.dat";
  vecwrite(filename.c_str(),assess_stat_parm_mean);
  filename = prefix + "as_st_func_mean.dat";
  vecwrite(filename.c_str(),assess_stat_func_mean);

  filename = prefix + "as_st_parm_var.dat";
  vecwrite(filename.c_str(),assess_stat_parm_var);
  filename = prefix + "as_st_func_var.dat";
  vecwrite(filename.c_str(),assess_stat_func_var);

  filename = prefix + "as_st.dat";
  ofstream fout(filename.c_str());
  if (!fout) return false;

  fout << "     "
       << "                mode"
       << "                mean"
       << "                sdev"  << '\n'; 

  fout << "assess_stat_parm" << '\n';
  for (INTEGER i=1; i<=len_assess_stat_parm; ++i) {
    fout << fmt('d',5,i);
    fout << fmt('f',20,8, mode.stat_mod_parm[i]);
    fout << fmt('f',20,8, assess_stat_parm_mean[i]);
    fout << fmt('f',20,8, assess_stat_parm_sdev[i]) << '\n';
  }

  fout << "assess_stat_func" << '\n';
  for (INTEGER i=1; i<=len_assess_stat_func; ++i) {
    fout << fmt('d',5,i);
    fout << fmt('f',20,8, mode.stat_mod_func[i]);
    fout << fmt('f',20,8, assess_stat_func_mean[i]);
    fout << fmt('f',20,8, assess_stat_func_sdev[i]) << '\n';
  }

  fout << '\n';
  fout << "log stat_mod_prior      = " 
       << fmt('f',20,8,mode.stat_mod_prior.log_den) << '\n';
  fout << "log assess_prior        = " 
       << fmt('f',20,8,mode.assess_prior.log_den) << '\n';
  fout << "log stat_mod_likelihood = " 
       << fmt('f',20,8,mode.stat_mod_logl.log_den) << '\n';
  if (est_blk.stat_draw_from_prior) {
    fout << "Draws are from prior, no posterior printed" << '\n';
  }
  else {
    REAL log_post = mode.stat_mod_prior.log_den;
    log_post += mode.assess_prior.log_den;
    log_post += mode.stat_mod_logl.log_den;
    fout << "log stat_mod_posterior  = " << fmt('f',20,8,log_post) << '\n';
  }

  fout << '\n';
  fout << '\n';
  fout << "     "
       << "                mode"
       << "                mean"
       << "                sdev"  << '\n'; 

  fout << "assess_stat_parm" << '\n';
  for (INTEGER i=1; i<=len_assess_stat_parm; ++i) {
    fout << fmt('d',5,i);
    fout << fmt('e',20,10, mode.stat_mod_parm[i]);
    fout << fmt('e',20,10, assess_stat_parm_mean[i]);
    fout << fmt('e',20,10, assess_stat_parm_sdev[i]) << '\n';
  }

  fout << "assess_stat_func" << '\n';
  for (INTEGER i=1; i<=len_assess_stat_func; ++i) {
    fout << fmt('d',5,i);
    fout << fmt('e',20,10, mode.stat_mod_func[i]);
    fout << fmt('e',20,10, assess_stat_func_mean[i]);
    fout << fmt('e',20,10, assess_stat_func_sdev[i]) << '\n';
  }

  fout << '\n';
  fout << "log stat_mod_prior      = " 
       << fmt('e',20,10,mode.stat_mod_prior.log_den) << '\n';
  fout << "log assess_prior        = " 
       << fmt('e',20,10,mode.assess_prior.log_den) << '\n';
  fout << "log stat_mod_likelihood = " 
       << fmt('e',20,10,mode.stat_mod_logl.log_den) << '\n';
  if (est_blk.stat_draw_from_prior) {
    fout << "Draws are from prior, no posterior printed" << '\n';
  }
  else {
    REAL log_post = mode.stat_mod_prior.log_den;
    log_post += mode.assess_prior.log_den;
    log_post += mode.stat_mod_logl.log_den;
    fout << "log stat_mod_posterior  = " << fmt('e',20,10,log_post) << '\n';
  }

  fout.clear();
  fout.close();
  return true;
}

