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

Copyright (C) 2005, 2006, 2008, 2009.

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 "pathname.h"

#include <typeinfo>
#include "mpi.h"

#include "libscl.h"
#include "libgsm.h"
#include "gsmusr.h"
#include "gsm.h"
#include "snp_base.h"

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

namespace {

  int my_rank;
  int no_procs;

  void mpi_error (string msg) {
    cerr << msg << endl; MPI_Abort(MPI_COMM_WORLD, my_rank);
  }
  void mpi_warn (string msg) {
    cerr << msg << endl;
  }
}

int main(int argc, char** argp, char** envp)
{
  MPI_Init(&argc, &argp);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &no_procs);

  LIB_ERROR_HANDLER_PTR previous_error = set_lib_error_handler(&mpi_error);
  LIB_WARN_HANDLER_PTR previous_warn = set_lib_warn_handler(&mpi_warn);

  string pathname = "/tmp/";

  if (my_rank == 0) {
    pathname = string(PATHNAME) + string("/");
  }

  if (my_rank == 0) { // Check some type assumptions

    REAL x1; double x2;
    if (typeid(x1) != typeid(x2) ) {
      error("Error, a REAL is not a double");
      return 1;
    }
    
    INTEGER i1; int i2;   
    if (typeid(i1) != typeid(i2)) {
      error("Error, an INTEGER is not an int");
      return 1;
    }

    /*
    Cannot check these assumptions:
    MPI_DOUBLE is double
    MPI_INT    is int
    MPI__CHAR  is char
    */

  }

  ostream* detail_ptr = &cout;

  string parmfile, prefix;
 
  vector<string> pfvec; // Parmfile as a vector of strings

  vecstrbuf pfbuf;         // Parmfile as char array of null
                        // terminated strings stored end to end

  int  dim[4];          // dim[0] size, dim[1] rows, dim[2] cols, dim[3] extra

  dim[0] = dim[1] = dim[2] = dim[3] = 0;

  if (my_rank == 0) {

    string filename = pathname + "control.dat";

    ifstream ctrl_ifs(filename.c_str());
    if (!ctrl_ifs) error("Error, gsm, control.dat open failed");

    ctrl_ifs >> parmfile >> prefix;

    filename = pathname + prefix + ".detail.dat";

    detail_ptr = new(nothrow) ofstream(filename.c_str());
    if ( (detail_ptr==0) || (!*detail_ptr) ) 
      error("Error, gsm, detail.dat open failed");

    filename = pathname + parmfile;

    ifstream pf_ifs(filename.c_str());
    if (!pf_ifs) error("Error, gsm, cannot open " + parmfile);

    string line;

    while (getline(pf_ifs, line)) pfvec.push_back(line);

    vecstrbuf send_buf(pfvec);
    
    dim[0] = send_buf.size();
    dim[1] = send_buf.get_rows();
    dim[2] = send_buf.get_cols();

    pfbuf = send_buf;
  }

  MPI_Bcast(&dim,4,MPI_INT,0,MPI_COMM_WORLD);
  if (my_rank != 0) pfbuf.resize(dim[1],dim[2]);      
  MPI_Bcast(pfbuf.get_ptr(),dim[0],MPI_CHAR,0,MPI_COMM_WORLD);  
  if (my_rank != 0) pfvec = pfbuf.get_vec();       

  if (my_rank != 0) { // Must suppress print if my_rank != 0
    snp::keyword kw;
    vector<string>::iterator kw_ptr;
    string header=kw.set_keyword("ESTIMATION DESCRIPTION");
    kw_ptr = find_if(pfvec.begin(), pfvec.end(), kw);       
    if (kw_ptr == pfvec.end() || kw_ptr + 3 >= pfvec.end())
      error("Error, gsm, " + header);
    kw_ptr += 3;
    *kw_ptr = "0 ";
  }

  ostream& detail = *detail_ptr;

  gsmblocks pf;
  pf.set_blks(pfvec, detail);

  estblock est_blk = pf.get_estblock();
  datblock dat_blk = pf.get_datblock();
  sciblock sci_blk = pf.get_sciblock();
  statblock stat_blk = pf.get_statblock();

  realmat data;

  if (my_rank == 0) {
  
    if (!dat_blk.read_data(pathname,data)) {
      error("Error, gsm, cannot read data, dsn = " + dat_blk.dsn);
    }

    if (est_blk.print) {
      detail << starbox("/First 12 observations//");
      detail << data("",seq(1,12));
      detail << starbox("/Last 12 observations//");
      detail << data("",seq(dat_blk.n-11, dat_blk.n));
      detail.flush();
    }
        
    dim[0] = data.size();
    dim[1] = data.get_rows();
    dim[2] = data.get_cols();
  }

  MPI_Bcast(&dim,4,MPI_INT,0,MPI_COMM_WORLD);
  if (my_rank != 0) data.resize(dim[1],dim[2]);      
  MPI_Bcast(data.get_x(),dim[0],MPI_DOUBLE,0,MPI_COMM_WORLD);  

  vector<string> stat_pfvec;
  vecstrbuf stat_pfbuf;
  
  if (my_rank == 0) {
    if (est_blk.print) {
      detail << starbox("/Begin instantiation of stat_mod//");
      detail.flush();
    }
  }
  
  if (stat_blk.is_stat_parmfile) {

    if (my_rank == 0) {
      string filename;
  
      if (stat_blk.stat_parmfile[0] == '/') {
        filename = stat_blk.stat_parmfile;
      } 
      else {
        filename = pathname + stat_blk.stat_parmfile;
      }

      ifstream stat_pf_ifs(filename.c_str());
      if (!stat_pf_ifs) error("Error, gsm, cannot open " + filename);
  
      string line;
  
      while (getline(stat_pf_ifs, line)) stat_pfvec.push_back(line);
  
      if (est_blk.print) {
        detail << starbox("/stat_pfvec//") << '\n';
        vector<string>::const_iterator itr = stat_pfvec.begin();
        while (itr != stat_pfvec.end()) detail << *itr++ << '\n';
        detail.flush(); 
      }

      vecstrbuf send_buf(stat_pfvec);
      
      dim[0] = send_buf.size();
      dim[1] = send_buf.get_rows();
      dim[2] = send_buf.get_cols();

      stat_pfbuf = send_buf;
    }

    MPI_Bcast(&dim,4,MPI_INT,0,MPI_COMM_WORLD);
    if (my_rank != 0) stat_pfbuf.resize(dim[1],dim[2]);
    MPI_Bcast(stat_pfbuf.get_ptr(),dim[0],MPI_CHAR,0,MPI_COMM_WORLD);
    if (my_rank != 0) stat_pfvec = stat_pfbuf.get_vec();
  } 
  
  vector<string> stat_alvec = pf.get_stat_alvec();

  if (my_rank == 0) {
    if (est_blk.print) {
      detail << starbox("/stat_alvec//") << '\n';
      vector<string>::const_iterator itr = stat_alvec.begin();
      while (itr != stat_alvec.end()) detail << *itr++ << '\n';
      detail.flush(); 
    }
  
    if (est_blk.print) {
      detail << starbox("/stat_mod constructor now controls print//");
      detail.flush();
    }
  }

  stat_mod_type stat_mod(&data, stat_pfvec, stat_alvec, detail);

  if (stat_blk.len_stat_parm != stat_mod.len_parm()) {
    fmt lpf('d',5,stat_blk.len_stat_parm);
    fmt lsm('d',5,stat_mod.len_parm());
    string msg = "Error, gsm, parmfile len_parm is";
    msg += lpf.get_ostr();
    msg += ", but stat_mod returns";
    msg += lsm.get_ostr();
    error(msg);
  }

  if (stat_blk.len_stat_func != stat_mod.len_func()) {
    fmt lpf('d',5,stat_blk.len_stat_func);
    fmt lsm('d',5,stat_mod.len_func());
    string msg = "Error, gsm, parmfile len_func is";
    msg += lpf.get_ostr();
    msg += ", but stat_mod returns";
    msg += lsm.get_ostr();
    error(msg);
  }

  stat_mod.set_parm(pf.get_stat_parm());
  stat_mod.set_data_ptr(&data);

  if (my_rank == 0) {

    if (est_blk.print) {
      detail << starbox("/gsm now controls print//");
      detail.flush();
    }

    if (est_blk.print) {
      realmat parm;
      stat_mod.get_parm(parm);
      detail << parm << '\n';
      detail << starbox("/Log likelihood at the data//") << '\n';
      detail << stat_mod.loglikelihood() << '\n';
      realmat sim, func;
      stat_mod.gen_sim(sim, func);
      detail << starbox("/First 12 stat_mod simulations//");
      detail << sim("",seq(1,12));
      detail << starbox("/Last 12 stat_mod simulations//");
      intvec idx = seq(sim.get_cols()-11, sim.get_cols());
      detail << sim("",idx);
      detail << starbox("/stat_mod func//");
      detail << func << '\n';
      detail.flush();
    }
  
    if (est_blk.print) {
      detail << starbox("/End instantiation of stat_mod//");
      detail.flush();
    }
        
    if (est_blk.print) {
      detail << starbox("/Begin instantiation of sci_mod//");
      detail.flush();
    }
  }

  vector<string> sci_pfvec;
  vecstrbuf sci_pfbuf;

  if (sci_blk.is_sci_parmfile) {

    if (my_rank == 0) {
      string filename;

      if (sci_blk.sci_parmfile[0] == '/') {
        filename = sci_blk.sci_parmfile;
      } 
      else {
        filename = pathname + sci_blk.sci_parmfile;
      }

      ifstream sci_pf_ifs(filename.c_str());
      if (!sci_pf_ifs) error("Error, gsm, cannot open " + filename);

      string line;

      while (getline(sci_pf_ifs, line)) sci_pfvec.push_back(line);

      if (est_blk.print) {
        detail << starbox("/sci_pfvec//") << '\n';
        vector<string>::const_iterator itr = sci_pfvec.begin();
        while (itr != sci_pfvec.end()) detail << *itr++ << '\n';
        detail.flush(); 
      }

      vecstrbuf send_buf(sci_pfvec);
      
      dim[0] = send_buf.size();
      dim[1] = send_buf.get_rows();
      dim[2] = send_buf.get_cols();

      sci_pfbuf = send_buf;
    }

    MPI_Bcast(&dim,4,MPI_INT,0,MPI_COMM_WORLD);
    if (my_rank != 0) sci_pfbuf.resize(dim[1],dim[2]);
    MPI_Bcast(sci_pfbuf.get_ptr(),dim[0],MPI_CHAR,0,MPI_COMM_WORLD);
    if (my_rank != 0) sci_pfvec = sci_pfbuf.get_vec();
  }

  vector<string> sci_alvec = pf.get_sci_alvec();

  if (my_rank == 0) {
    if (est_blk.print) {
      detail << starbox("/sci_alvec//") << '\n';
      vector<string>::const_iterator itr = sci_alvec.begin();
      while (itr != sci_alvec.end()) detail << *itr++ << '\n';
      detail.flush(); 
    }
  
    if (est_blk.print) {
      detail << starbox("/sci_mod constructor now controls print//");
      detail.flush();
    }
  }

  sci_mod_type sci_mod(&data, sci_pfvec, sci_alvec, detail);

  if (sci_blk.len_sci_parm != sci_mod.len_parm()) {
    fmt lpf('d',5,sci_blk.len_sci_parm);
    fmt lsm('d',5,sci_mod.len_parm());
    string msg = "Error, gsm, parmfile len_parm is";
    msg += lpf.get_ostr();
    msg += ", but sci_mod returns";
    msg += lsm.get_ostr();
    error(msg);
  }

  if (sci_blk.len_sci_func != sci_mod.len_func()) {
    fmt lpf('d',5,sci_blk.len_sci_func);
    fmt lsm('d',5,sci_mod.len_func());
    string msg = "Error, gsm, parmfile len_func is";
    msg += lpf.get_ostr();
    msg += ", but sci_mod returns";
    msg += lsm.get_ostr();
    error(msg);
  }

  sci_summary sci_sum(est_blk);
  assess_summary assess_sum(est_blk);

  sci_mod.set_parm(pf.get_sci_parm());
  
  if (my_rank == 0) {

    if (est_blk.print) {
      detail << starbox("/gsm now controls print//");
      detail.flush();
    }
  
    if (est_blk.print) {
      detail << starbox("/Parameter values//");
      realmat parm;
      sci_mod.get_parm(parm);
      detail << parm << '\n';
      realmat sim, func;
      sci_mod.gen_sim(sim, func);
      detail << starbox("/First 12 sci_mod simulations//");
      detail << sim("",seq(1,12));
      detail << starbox("/Last 12 sci_mod simulations//");
      intvec idx = seq(sim.get_cols()-11, sim.get_cols());
      detail << sim("",idx);
      detail << starbox("/sci_mod func//");
      detail << func << '\n';
      detail.flush();
    }
  
    if (est_blk.print) {
      detail << starbox("/End instantiation of sci_mod//");
      detail.flush();
    }
  
    if (est_blk.print) {
      detail << starbox("/Instantiate the proposal densities//");
      detail.flush();
    }
  }

  group_move  stat_mod_prop(pf.get_stat_prop_def());
  grid_group_move  sci_mod_prop(pf.get_sci_prop_def());

  if (my_rank == 0) {
    if (est_blk.print) {
      detail << starbox("/sci_mod proposal density//");
      sci_mod_prop.write_proposal(detail);
      detail << starbox("/Instantiate the implied map//");
      detail.flush();
    }
  }

  implied_map im;
  im.set_draw_from_posterior(!est_blk.sci_draw_from_prior);

  implied_map_buffer imbuf;

  dim[0] = dim[1] = dim[2] = dim[3] = 0;

  if (my_rank == 0) {
    string filename = pathname + prefix + ".implied_map.dat";

    if (im.read_map(filename.c_str())) {
      if (im.empty()) error("Error, gsm, implied_map.dat empty");

      implied_map_buffer send_buf(im);

      imbuf = send_buf;

      dim[0] = imbuf.size();
      dim[1] = imbuf.get_rows();
      dim[2] = imbuf.get_cols();
    }
    else {
      if (est_blk.run_sci_chain) {
        warn("Warning, gsm, could not read implied map, will compute afresh");
      }
      else if (est_blk.run_stat_chain) {
        error("Error, gsm, could not read implied map, you must run sci_chain");
      }
    }
  }

  MPI_Bcast(&dim,4,MPI_INT,0,MPI_COMM_WORLD);
  if (dim[0] > 0) {
    if (my_rank != 0) imbuf.resize(dim[1],dim[2]);
    MPI_Bcast(imbuf.get_ptr(),dim[0],MPI_DOUBLE,0,MPI_COMM_WORLD);
  }

  if (my_rank != 0) {
    im.clear(); 
    for (INTEGER j=1; j<=dim[2]; ++j) {
      sci_val sv = unpack_sci_val(imbuf.get_col(j));
      im.set_val(sv);
    } 
  }

  if (my_rank == 0) {
    if (est_blk.print) {
      detail << starbox("/Instantiate assess_prior//"); 
      detail.flush();
    }
  }

  assess_prior ap(stat_mod.len_parm());
  ap.set_kappa(est_blk.kappa); 

  realmat sigma;
  INTEGER nobs = 0;

  dim[0] = dim[1] = dim[2] = 0;

  if (my_rank == 0) {
    string filename = pathname + prefix + ".assess_sigma.dat";

    if (ap.read_sigma(filename.c_str())) {

      if (stat_blk.len_stat_parm != ap.get_sigma().get_rows()) {
        error("Error, gsm, sigma read in has wrong dimensions");
      }

      sigma = ap.get_sigma();
      nobs = ap.get_n();

      dim[0] = sigma.size();
      dim[1] = sigma.get_rows(); 
      dim[2] = sigma.get_cols(); 
      dim[3] = nobs; 
    }
    else {
      if (est_blk.run_sci_chain) {
        warn("Warning, gsm, could not read assess sigma, will compute afresh");
      }
      else if (est_blk.run_stat_chain) {
        error("Error, gsm, could not read sigma, you must run sci_chain");
      }
    }
  }

  MPI_Bcast(&dim,4,MPI_INT,0,MPI_COMM_WORLD);
  if (dim[0] > 0) {
    if (my_rank != 0) {
      sigma.resize(dim[1],dim[2]);
      nobs = dim[3];
    }
    MPI_Bcast(sigma.get_x(),dim[0],MPI_DOUBLE,0,MPI_COMM_WORLD);
    if (my_rank != 0) ap.update_sigma(sigma, nobs);
  }

  realmat sci_parm_chain; 
  realmat sci_func_chain;
  realmat sci_prior_chain;
    
  realmat sci_rej;

  realmat sci_sub_parm_chain;
  realmat sci_sub_logl_chain;
  realmat sci_sub_rej_chain;

  realmat stat_parm_chain; 
  realmat stat_func_chain;
  realmat stat_prior_chain; 
  realmat stat_logl_chain;
  realmat stat_sub_logl_chain;

  realmat assess_stat_parm_chain;
  realmat assess_stat_func_chain;
  realmat assess_stat_prior_chain;
  realmat assess_stat_logl_chain;

  realmat assess_prior_chain;

  realmat reject;

  if (my_rank == 0) {

    INTEGER count = 1;

    while (count < no_procs) {

      int  src;
      int  tag;
      int  i;
      char number[5];
      char sender[5];

      MPI_Status status;

      MPI_Recv(&i,1,MPI_INT,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&status);

      src = status.MPI_SOURCE;
      tag = status.MPI_TAG;

      sprintf(number,"%03d",i);
      sprintf(sender,"%03d",src);

      string id = string(sender) + string(".") + string(number);

      REAL* buf;
      string filename;

      if (tag == 50) {

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        sci_parm_chain.resize(dim[1],dim[2]);
        buf = sci_parm_chain.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".sci_parm." + id;
        vecwrite(filename.c_str(),sci_parm_chain);
          
        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        sci_func_chain.resize(dim[1],dim[2]);
        buf = sci_func_chain.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".sci_func." + id;
        vecwrite(filename.c_str(),sci_func_chain);
          
        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        sci_prior_chain.resize(dim[1],dim[2]);
        buf = sci_prior_chain.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".sci_prior." + id;
        vecwrite(filename.c_str(),sci_prior_chain);

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        stat_parm_chain.resize(dim[1],dim[2]);
        buf = stat_parm_chain.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".stat_parm." + id;
        vecwrite(filename.c_str(),stat_parm_chain);

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        stat_func_chain.resize(dim[1],dim[2]);
        buf = stat_func_chain.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".stat_func." + id;
        vecwrite(filename.c_str(),stat_func_chain);

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        stat_prior_chain.resize(dim[1],dim[2]);
        buf = stat_prior_chain.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".stat_prior." + id;
        vecwrite(filename.c_str(),stat_prior_chain);

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        stat_logl_chain.resize(dim[1],dim[2]);
        buf = stat_logl_chain.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".stat_logl." + id;
        vecwrite(filename.c_str(),stat_logl_chain);

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        reject.resize(dim[1],dim[2]);
        buf = reject.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".sci_rej." + id;
        vecwrite(filename.c_str(),reject);

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        imbuf.resize(dim[1],dim[2]);
        buf = imbuf.get_ptr();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        implied_map recv_map;
        for (INTEGER j=1; j<=imbuf.get_cols(); ++j) {
          sci_val sv = unpack_sci_val(imbuf.get_col(j));
          recv_map.set_val(sv);
        }

        im.update_map(recv_map);
        filename = pathname + prefix + ".implied_map.new";
        im.write_map(filename.c_str());

        implied_map_buffer send_buf(im);
        imbuf = send_buf;
        dim[0] = imbuf.size();
        dim[1] = imbuf.get_rows();
        dim[2] = imbuf.get_cols();
        MPI_Send(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD);
        buf = imbuf.get_ptr();
        MPI_Send(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD);
          
        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        realmat sigma(dim[1],dim[2]);
        buf = sigma.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        
        ap.update_sigma(sigma,dim[3]);
        filename = pathname + prefix + ".assess_sigma.new";
        ap.write_sigma(filename.c_str());

        INTEGER len_chain = sci_parm_chain.get_cols();
        pf.set_sci_parm(sci_parm_chain("",len_chain));
        pf.set_stat_parm(stat_parm_chain("",len_chain));
        string pre = pathname + prefix;
        string post = id + ".end";
        pf.write_blks(parmfile,pre,post);

        sci_val mode = im.get_mode();
        pf.set_sci_parm(mode.sci_mod_parm);
        pf.set_stat_parm(mode.stat_mod_parm);
        post = id + ".fit";
        pf.write_blks(parmfile,pre,post);

        sci_sum.update(sci_parm_chain,sci_func_chain,stat_parm_chain,
          stat_func_chain,mode);
        sci_sum.write(pathname + prefix); 
      }

      if (tag == 60) {

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        assess_stat_parm_chain.resize(dim[1],dim[2]);
        buf = assess_stat_parm_chain.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".as_st_parm." + id;
        vecwrite(filename.c_str(),assess_stat_parm_chain);
          
        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        assess_stat_func_chain.resize(dim[1],dim[2]);
        buf = assess_stat_func_chain.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".as_st_func." + id;
        vecwrite(filename.c_str(),assess_stat_func_chain);

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        assess_stat_prior_chain.resize(dim[1],dim[2]);
        buf = assess_stat_prior_chain.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".as_st_prior." + id;
        vecwrite(filename.c_str(),assess_stat_prior_chain);

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        assess_stat_logl_chain.resize(dim[1],dim[2]);
        buf = assess_stat_logl_chain.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".as_st_logl." + id;
        vecwrite(filename.c_str(),assess_stat_logl_chain);

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        assess_prior_chain.resize(dim[1],dim[2]);
        buf = assess_prior_chain.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".as_prior." + id;
        vecwrite(filename.c_str(),assess_prior_chain);

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        reject.resize(dim[1],dim[2]);
        buf = reject.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        filename = pathname + prefix + ".as_rej." + id;
        vecwrite(filename.c_str(),reject);

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        realmat packed_mode(dim[1],dim[2]);
        buf = packed_mode.get_x();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        stat_val mode = unpack_stat_val(packed_mode);
        assess_sum.update(assess_stat_parm_chain, assess_stat_func_chain, mode);
        assess_sum.write(pathname + prefix); 
      }

      if (tag == 90) {
      count++;
      cout << "Process " << src << " has finished\n";
      }
    }

  }
  else {

    realmat sci_parm_start = pf.get_sci_parm();
    realmat stat_parm_start = pf.get_stat_parm();
    realmat stat_parm_start_save = stat_parm_start;

    if (est_blk.run_sci_chain) {

      sci_mcmc schain(sci_mod,sci_mod_prop,stat_mod,stat_mod_prop,ap,&data,im);
      schain.set_draw_from_posterior(!est_blk.sci_draw_from_prior);
      schain.set_len_chain(sci_blk.len_sci_chain);
      schain.set_len_sub_chain(sci_blk.len_sub_chain);
      schain.set_num_sub_chains(sci_blk.num_sub_chains);
      schain.set_num_polish_iter(sci_blk.num_polish_iter);
      schain.set_polish_toler(sci_blk.polish_toler);
      schain.set_analytic_mle(sci_blk.analytic_mle);

      INT_32BIT seed = sci_blk.sci_seed + 100*my_rank;

      for (INTEGER ifile = 0; ifile<=sci_blk.num_sci_files; ++ifile) {

        reject = schain.draw(seed, sci_parm_start, stat_parm_start,
              sci_parm_chain, sci_func_chain,
              sci_prior_chain,
              sci_sub_parm_chain, sci_sub_logl_chain, 
              sci_sub_rej_chain,
              stat_parm_chain, stat_func_chain,
              stat_prior_chain, stat_logl_chain,
              stat_sub_logl_chain);

        int tag = 50;
        int dest = 0;
        REAL* buf = 0;

        MPI_Send(&ifile,1,MPI_INT,dest,tag,MPI_COMM_WORLD);

        dim[0] = sci_parm_chain.size();
        dim[1] = sci_parm_chain.get_rows();
        dim[2] = sci_parm_chain.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = sci_parm_chain.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        dim[0] = sci_func_chain.size();
        dim[1] = sci_func_chain.get_rows();
        dim[2] = sci_func_chain.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = sci_func_chain.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        dim[0] = sci_prior_chain.size();
        dim[1] = sci_prior_chain.get_rows();
        dim[2] = sci_prior_chain.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = sci_prior_chain.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        dim[0] = stat_parm_chain.size();
        dim[1] = stat_parm_chain.get_rows();
        dim[2] = stat_parm_chain.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = stat_parm_chain.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        dim[0] = stat_func_chain.size();
        dim[1] = stat_func_chain.get_rows();
        dim[2] = stat_func_chain.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = stat_func_chain.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        dim[0] = stat_prior_chain.size();
        dim[1] = stat_prior_chain.get_rows();
        dim[2] = stat_prior_chain.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = stat_prior_chain.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        dim[0] = stat_logl_chain.size();
        dim[1] = stat_logl_chain.get_rows();
        dim[2] = stat_logl_chain.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = stat_logl_chain.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        dim[0] = reject.size();
        dim[1] = reject.get_rows();
        dim[2] = reject.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = reject.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        implied_map_buffer send_buf(im); 
        imbuf = send_buf;
        dim[0] = imbuf.size();
        dim[1] = imbuf.get_rows();
        dim[2] = imbuf.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = imbuf.get_ptr();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        MPI_Status status;

        MPI_Recv(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD,&status);
        imbuf.resize(dim[1],dim[2]);
        buf = imbuf.get_ptr();
        MPI_Recv(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD,&status);
        implied_map recv_map;
        for (INTEGER j=1; j<=imbuf.get_cols(); ++j) {
          sci_val sv = unpack_sci_val(imbuf.get_col(j));
          recv_map.set_val(sv);
        }

        im.update_map(recv_map);

        realmat sigma = ap.get_sigma();
        dim[0] = sigma.size();
        dim[1] = sigma.get_rows();
        dim[2] = sigma.get_cols();
        dim[3] = ap.get_n();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = sigma.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);
      }
    }

    if (est_blk.run_stat_chain) {

      ap.set_implied_map(im);

      stat_mcmc achain(stat_mod, stat_mod_prop, ap);
      achain.set_draw_from_posterior(!est_blk.stat_draw_from_prior);
      achain.set_data_ptr(&data);
      achain.set_use_stat_mod_prior(true);
      achain.set_use_assess_prior(true);
      achain.set_len_chain(stat_blk.len_stat_chain);

      INT_32BIT seed = stat_blk.stat_seed + 100*my_rank;

      stat_parm_start = stat_parm_start_save;

      for (INTEGER ifile = 0; ifile<=stat_blk.num_stat_files; ++ifile) {

        reject = achain.draw(seed, stat_parm_start, assess_stat_parm_chain,
            assess_stat_func_chain, assess_stat_prior_chain,
            assess_stat_logl_chain, assess_prior_chain);

        int tag = 60;
        int dest = 0;
        REAL* buf = 0;

        MPI_Send(&ifile,1,MPI_INT,dest,tag,MPI_COMM_WORLD);

        dim[0] = assess_stat_parm_chain.size();
        dim[1] = assess_stat_parm_chain.get_rows();
        dim[2] = assess_stat_parm_chain.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = assess_stat_parm_chain.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        dim[0] = assess_stat_func_chain.size();
        dim[1] = assess_stat_func_chain.get_rows();
        dim[2] = assess_stat_func_chain.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = assess_stat_func_chain.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        dim[0] = assess_stat_prior_chain.size();
        dim[1] = assess_stat_prior_chain.get_rows();
        dim[2] = assess_stat_prior_chain.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = assess_stat_prior_chain.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        dim[0] = assess_stat_logl_chain.size();
        dim[1] = assess_stat_logl_chain.get_rows();
        dim[2] = assess_stat_logl_chain.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = assess_stat_logl_chain.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        dim[0] = assess_prior_chain.size();
        dim[1] = assess_prior_chain.get_rows();
        dim[2] = assess_prior_chain.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = assess_prior_chain.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        dim[0] = reject.size();
        dim[1] = reject.get_rows();
        dim[2] = reject.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = reject.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);

        realmat packed_mode = pack_stat_val(achain.get_mode());
        dim[0] = packed_mode.size();
        dim[1] = packed_mode.get_rows();
        dim[2] = packed_mode.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        buf = packed_mode.get_x();
        MPI_Send(buf,dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);
      }
    }

    int tag = 90;
    int dest = 0;
    int i = 0;
    MPI_Send(&i,1,MPI_INT,dest,tag,MPI_COMM_WORLD);
  }

  if (my_rank == 0) delete detail_ptr;

  MPI_Finalize();

  previous_error = set_lib_error_handler(previous_error);
  previous_warn = set_lib_warn_handler(previous_warn);

  return 0;
}


