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

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

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.

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

#undef SCI_MCMC_DEBUG                    // Messages to cout and cerr
#undef STAT_MCMC_DEBUG                   // Messages to cout
#undef STICKY_CHAIN_DEBUG                // Messages to cerr via warn()
#undef BEGINNING_AND_ENDING_DEBUG        // Messages to cerr via warn()

#include "libgsm.h"
using namespace std;  
using namespace scl;
using namespace libgsm;

INTEGER libgsm::sci_mcmc::make_sci_val(INT_32BIT& sub_chain_seed, 
  realmat& sci_parm_start, realmat& stat_parm_start,
  sci_val& sv, string& errmsg, realmat& sci_sub_parm_chain, 
  realmat& sci_sub_logl_chain, realmat& sci_sub_rej_chain)
{

  if (analytic_mle) error("Error, sci_mcmc::make_sci_val, erroneous call");

  errmsg = "errmsg not set #1";

  /*
  INTEGER len_sci_parm = sci_mod.len_parm();
  INTEGER len_sci_func = sci_mod.len_func();
  INTEGER len_stat_func = stat_mod.len_func();
  */
  INTEGER len_stat_parm = stat_mod.len_parm();

  #if defined SCI_MCMC_DEBUG
    cout << "\n\n Debug: Entering sci_mcmc.make_sci_val that uses stat_mcmc\n";
  #endif

  sci_val scv(sci_parm_start);
  sv = scv;

  sci_mod.set_parm(sv.sci_mod_parm);

  /*
  Setting the parameter before checking support is intentional.
  Do not change it.  SNP will not work if you do.
  */

  sv.sci_mod_support = sci_mod.support(sv.sci_mod_parm);

  if (!sv.sci_mod_support) {
    errmsg ="sci_parm_start not in support of sci_mod";
    return 1;
  }

  realmat sci_mod_sim;
  sv.sci_mod_simulate = sci_mod.gen_sim(sci_mod_sim, sv.sci_mod_func);

  if (!sv.sci_mod_simulate) {
    errmsg = "cannot simulate sci_mod at sci_parm_start";
    return 2;
  }

  sv.sci_mod_prior = sci_mod.prior(sv.sci_mod_parm, sv.sci_mod_func);
  
  if (!sv.sci_mod_prior.positive) {
    errmsg = "sci_parm_start not in support of sci_mod_prior";
    return 3;
  }

  stat_mod.set_parm(stat_parm_start);

  /*
  Setting the parameter before checking support is intentional.
  Do not change it.  SNP will not work if you do.
  */

  if (!stat_mod.support(stat_parm_start)) {
    errmsg = "stat_parm_start not in support of stat_mod";
  }

  #if defined SCI_MCMC_DEBUG
    cout << "\n Debug: sv before running subchain\n\n" << sv << endl;
    cout << "\n Debug: stat_parm_start before\n" << stat_parm_start << endl;
  #endif

  stat_mcmc subchain(stat_mod, stat_prop, aprior); 

  subchain.set_data_ptr(&sci_mod_sim);
  subchain.set_use_stat_mod_prior(false);
  subchain.set_use_assess_prior(false);
  subchain.set_draw_from_posterior(true);
  subchain.set_len_chain(len_sub_chain);

  realmat stat_parm_subchain;
  realmat stat_func_subchain;
  realmat stat_prior_subchain;
  realmat stat_logl_subchain;
  realmat assess_prior_subchain;

  realmat statrej = subchain.draw(sub_chain_seed, stat_parm_start,
                      stat_parm_subchain, stat_func_subchain, 
                      stat_prior_subchain, stat_logl_subchain, 
                      assess_prior_subchain);

  #if defined SCI_MCMC_DEBUG
    cout << "\n Debug: statrej from subchain\n" << statrej << endl;
    cout << "\n Debug: stat_parm_start after\n" << stat_parm_start << endl;
    cout << "\n Debug: stat_parm_sv after\n" << subchain.get_mode() << endl;
    cout << "\n Debug: stat_parm_logl after\n" 
         << subchain.get_mode().stat_mod_logl << endl;
  #endif

  stat_val stat_mode = subchain.get_mode();

  if (!stat_mode.stat_mod_logl.positive) {
    errmsg = "implied map computation failed in subchain";
    return 4;
  }

  bool chain_good = false;
  for (INTEGER i=1; i<=len_stat_parm; ++i) {
    if ((statrej(i,1) < 1.0) && (statrej(i,2) > 0.0)) chain_good = true;
  }

  if (chain_good) {
    INTEGER midway = len_sub_chain >= 2 ? len_sub_chain/2 : 1;
    realmat residual 
       = stat_parm_subchain("",len_sub_chain) - stat_parm_subchain("",midway);
    aprior.increment_sigma(residual);
  }

  sv.stat_mod_parm = stat_mode.stat_mod_parm;
  sv.stat_mod_sub_logl = stat_mode.stat_mod_logl;

  stat_mod_eqns polish_eqns(stat_mod);
  nlopt minimizer(polish_eqns);

  realmat polish_stat_parm = stat_mode.stat_mod_parm;
  den_val polish_logl = stat_mode.stat_mod_logl;

  if (num_polish_iter > 0) {

    minimizer.set_iter_limit(num_polish_iter);
    minimizer.set_solution_tolerance(polish_toler);
    
    #if defined SCI_MCMC_DEBUG
      cout << "\n Debug: polish subchain output follows" << endl;
      cout << "\n Debug: polish_stat_parm = " << polish_stat_parm << endl;
      cout << "\n Debug: polish_logl = " << polish_logl << endl;
      cout << "\n Debug: polish_logl.log_den/n = " 
           << polish_logl.log_den/stat_mod.num_obs() << endl;
      minimizer.set_warning_messages(false);
      minimizer.set_check_derivatives(false);
      minimizer.set_output(true,&cout,false);
    #else
      minimizer.set_warning_messages(false);
      minimizer.set_check_derivatives(false);
      minimizer.set_output(false);
    #endif

    realmat stat_parm_stop;
    minimizer.minimize(stat_mode.stat_mod_parm, stat_parm_stop);

    #if defined SCI_MCMC_DEBUG
      cout << "\n Debug: stat_parm_stop = " << stat_parm_stop << endl;
    #endif

    stat_mod.set_parm(stat_parm_stop);

    /*
    Setting the parameter before checking support is intentional.
    Do not change it.  SNP will not work if you do.
    */

    if (stat_mod.support(stat_parm_stop)) {

      den_val dv = stat_mod.loglikelihood();
    
      #if defined SCI_MCMC_DEBUG
        cout << '\n';
        cout << " Debug: stat_mod.loglikelihood() = " << dv << endl;
        cout << " Debug: stat_mod.loglikelihood()/n = "
             << dv.log_den/stat_mod.num_obs() << endl;
      #endif

      if (dv.positive) {
        polish_stat_parm = stat_parm_stop;
        polish_logl = dv;
      }

      if (polish_logl.log_den > stat_mode.stat_mod_logl.log_den) {

        sv.stat_mod_parm = polish_stat_parm;
        sv.stat_mod_sub_logl = polish_logl;

        /* Not a great idea
        Statements like the following statement sprinkled throughout can
        cause the subchain and hence the sci chain to hang. It seems like
        a good idea to start at previous peaks, but it isn't.

        stat_parm_start = sv.stat_mod_parm;

        Not resetting stat_parm_start has the effect of using the last
        value of the previous sub chain as the start for the next.
        */

      }
    }
  }

  // The last element of sci_sub_* is overwritten when num_polish_iter > 0

  if (num_sub_chain_saved < num_sub_chain_to_save) {
    INTEGER offset = num_sub_chain_saved;
    for (INTEGER j=1; j<=len_sub_chain; ++j) {
      for (INTEGER i=1; i<=len_stat_parm; ++i) {
        sci_sub_parm_chain(i,len_sub_chain*offset+j) = stat_parm_subchain(i,j);
      }
    }
    if (num_polish_iter > 0) {
      for (INTEGER i=1; i<=len_stat_parm; ++i) {
        sci_sub_parm_chain(i,len_sub_chain*offset+len_sub_chain) 
          = polish_stat_parm[i];
      }
    }
    ++offset;
    for (INTEGER i=1; i<=len_sub_chain; ++i) {
      sci_sub_logl_chain(i,offset) = stat_logl_subchain[i];
    }
    if (num_polish_iter > 0) {
     sci_sub_logl_chain(len_sub_chain,offset) = polish_logl.log_den;
    }
    for (INTEGER i=1; i<=len_stat_parm+1; ++i) {
      sci_sub_rej_chain(i,offset) = statrej(i,1);
    }
  }

  ++num_sub_chain_saved;

  stat_mod.set_data_ptr(data_ptr);
  stat_mod.set_parm(sv.stat_mod_parm);

  /*
  Setting the parameter before checking support is intentional.
  Do not change it.  SNP will not work if you do.
  */

  sv.stat_mod_support = stat_mod.support(sv.stat_mod_parm);

  if (!sv.stat_mod_support) {
    errmsg = "stat_parm not in support of stat_mod";
    return 5;
  }

  sv.stat_mod_logl = stat_mod.loglikelihood();

  realmat stat_mod_sim;
  sv.stat_mod_simulate = stat_mod.gen_sim(stat_mod_sim, sv.stat_mod_func);

  if (!sv.stat_mod_simulate) {
    errmsg = "stat_mod simulation failed";
    return 6;
  }

  sv.stat_mod_prior = stat_mod.prior(sv.stat_mod_parm, sv.stat_mod_func);

  if (!sv.stat_mod_prior.positive) {
    errmsg = "stat sub chain mode not in support of stat_mod_prior";
    return 7;
  }

  #if defined SCI_MCMC_DEBUG
    cout << "\n Debug: sv after running subchain\n\n" << sv << endl;
    cerr << "\n Debug: sv.stat_mod_sub_logl" << sv.stat_mod_sub_logl;
    cerr << "\n Debug: sv.stat_mod_logl" << sv.stat_mod_logl << endl;
  #endif

  return 0;
}

void libgsm::sci_mcmc::remake_sci_val_pair(INT_32BIT& sub_chain_seed,
       sci_val& first, sci_val& second,
       realmat& sci_sub_parm_chain, realmat& sci_sub_logl_chain, 
       realmat& sci_sub_rej_chain)
{

  #if defined STICKY_CHAIN_DEBUG
    string fsl0 = fmt('e',27,16,first.stat_mod_sub_logl.log_den).get_ostr();
    string fl0 = fmt('e',27,16,first.stat_mod_logl.log_den).get_ostr();
    string ssl0 = fmt('e',27,16,second.stat_mod_sub_logl.log_den).get_ostr();
    string sl0 = fmt('e',27,16,second.stat_mod_logl.log_den).get_ostr();
    REAL dsl = second.stat_mod_sub_logl.log_den-first.stat_mod_sub_logl.log_den;
    REAL dl = second.stat_mod_logl.log_den - first.stat_mod_logl.log_den;
    string rsl = fmt('e',27,16,exp(dsl)).get_ostr();
    string rl = fmt('e',27,16,exp(dl)).get_ostr();
    string msg = "\nDebug, sci_mcmc, entering remake";
    msg += "\nfirst";
    msg += "\n  old sub_logl and logl " + fsl0 + "  " + fl0;
    warn(msg);
  #endif 

  string errmsg = "errmsg not set #2";

  sci_val scv;

  realmat scv_sci_mod_parm = first.sci_mod_parm;
  realmat scv_stat_parm_start = first.stat_mod_parm;
  INTEGER rv = make_sci_val(sub_chain_seed, 
         scv_sci_mod_parm, scv_stat_parm_start,
         scv, errmsg, sci_sub_parm_chain, sci_sub_logl_chain,
         sci_sub_rej_chain);
  if (rv) {
    string msg = "Warning, sci_mcmc, " + errmsg + " at sci_mod_parm\n";
    for (INTEGER i=1; i<=scv.sci_mod_parm.size(); ++i) {
      msg += fmt('f',15,5,scv.sci_mod_parm[i])(); 
      msg += fmt('e',27,16,scv.sci_mod_parm[i])(); 
      msg += "\n";
    }
    warn(msg);
  }
  else if (scv.stat_mod_sub_logl.log_den > first.stat_mod_sub_logl.log_den) {
    first = scv;
    imap.set_val(scv);
    #if defined STICKY_CHAIN_DEBUG
      string fsl = fmt('e',27,16,first.stat_mod_sub_logl.log_den).get_ostr();
      string fl = fmt('e',27,16,first.stat_mod_logl.log_den).get_ostr();
      string msg = "first started at first, ";
      msg += "\n  new sub_logl and logl " + fsl  + "  " + fl;
      warn(msg);
    #endif
  }
  else {
    #if defined STICKY_CHAIN_DEBUG
      string msg = "first started at first was not better";
      warn(msg);
    #endif
  }

  scv_sci_mod_parm = first.sci_mod_parm;
  scv_stat_parm_start = second.stat_mod_parm;
  rv = make_sci_val(sub_chain_seed, 
         scv_sci_mod_parm, scv_stat_parm_start,
         scv, errmsg, sci_sub_parm_chain, sci_sub_logl_chain,
         sci_sub_rej_chain);
  if (rv) {
    string msg = "Warning, sci_mcmc, " + errmsg + " at sci_mod_parm\n";
    for (INTEGER i=1; i<=scv.sci_mod_parm.size(); ++i) {
      msg += fmt('f',15,5,scv.sci_mod_parm[i])(); 
      msg += fmt('e',27,16,scv.sci_mod_parm[i])(); 
      msg += "\n";
    }
    warn(msg);
  }
  else if (scv.stat_mod_sub_logl.log_den > first.stat_mod_sub_logl.log_den) {
    first = scv;
    imap.set_val(scv);
    #if defined STICKY_CHAIN_DEBUG
      string fsl = fmt('e',27,16,first.stat_mod_sub_logl.log_den).get_ostr();
      string fl = fmt('e',27,16,first.stat_mod_logl.log_den).get_ostr();
      string msg = "first started at second, ";
      msg += "\n  new sub_logl and logl " + fsl  + "  " + fl;
      warn(msg);
    #endif
  }
  else {
    #if defined STICKY_CHAIN_DEBUG
      string msg = "first started at second was not better";
      warn(msg);
    #endif
  }

  #if defined STICKY_CHAIN_DEBUG
    msg = "second";
    msg += "\n  old sub_logl and logl " + ssl0 + "  " + sl0;
    warn(msg);
  #endif 

  scv_sci_mod_parm = second.sci_mod_parm;
  scv_stat_parm_start = second.stat_mod_parm;
  rv = make_sci_val(sub_chain_seed, 
         scv_sci_mod_parm, scv_stat_parm_start,
         scv, errmsg, sci_sub_parm_chain, sci_sub_logl_chain,
         sci_sub_rej_chain);
  if (rv) {
    string msg = "Warning, sci_mcmc, " + errmsg + " at sci_mod_parm\n";
    for (INTEGER i=1; i<=scv.sci_mod_parm.size(); ++i) {
      msg += fmt('f',15,5,scv.sci_mod_parm[i])(); 
      msg += fmt('e',27,16,scv.sci_mod_parm[i])(); 
      msg += "\n";
    }
    warn(msg);
  }
  else if (scv.stat_mod_sub_logl.log_den > second.stat_mod_sub_logl.log_den) {
    second = scv;
    imap.set_val(scv);
    #if defined STICKY_CHAIN_DEBUG
      string ssl = fmt('e',27,16,second.stat_mod_sub_logl.log_den).get_ostr();
      string sl = fmt('e',27,16,second.stat_mod_logl.log_den).get_ostr();
      string msg = "second started at second, ";
      msg += "\n  new sub_logl and logl " + ssl  + "  " + sl;
      warn(msg);
    #endif
  }
  else {
    #if defined STICKY_CHAIN_DEBUG
      string msg = "second started at second was not better";
      warn(msg);
    #endif
  }

  scv_sci_mod_parm = second.sci_mod_parm;
  scv_stat_parm_start = first.stat_mod_parm;
  rv = make_sci_val(sub_chain_seed, 
         scv_sci_mod_parm, scv_stat_parm_start,
         scv, errmsg, sci_sub_parm_chain, sci_sub_logl_chain,
         sci_sub_rej_chain);
  if (rv) {
    string msg = "Warning, sci_mcmc, " + errmsg + " at sci_mod_parm\n";
    for (INTEGER i=1; i<=scv.sci_mod_parm.size(); ++i) {
      msg += fmt('f',15,5,scv.sci_mod_parm[i])(); 
      msg += fmt('e',27,16,scv.sci_mod_parm[i])(); 
      msg += "\n";
    }
    warn(msg);
  }
  else if (scv.stat_mod_sub_logl.log_den > second.stat_mod_sub_logl.log_den) {
    second = scv;
    imap.set_val(scv);
    #if defined STICKY_CHAIN_DEBUG
      string ssl = fmt('e',27,16,second.stat_mod_sub_logl.log_den).get_ostr();
      string sl = fmt('e',27,16,second.stat_mod_logl.log_den).get_ostr();
      string msg = "second started at first, ";
      msg += "\n  new sub_logl and logl " + ssl  + "  " + sl;
      warn(msg);
    #endif
  }
  else {
    #if defined STICKY_CHAIN_DEBUG
      string msg = "second started at first was not better";
      warn(msg);
    #endif
  }

  #if defined STICKY_CHAIN_DEBUG
    msg = "ratio, exp(second-first)";
    msg += "\n old  sub_logl and logl " + rsl + "  " + rl;
    dsl = second.stat_mod_sub_logl.log_den-first.stat_mod_sub_logl.log_den;
    dl = second.stat_mod_logl.log_den - first.stat_mod_logl.log_den;
    rsl = fmt('e',27,16,exp(dsl)).get_ostr();
    rl = fmt('e',27,16,exp(dl)).get_ostr();
    msg += "\n new  sub_logl and logl " + rsl + "  " + rl;
    warn(msg);
  #endif
}

realmat libgsm::sci_mcmc::draw(INT_32BIT& seed,
          realmat& sci_parm_start, realmat& stat_parm_start,
          realmat& sci_parm_chain, realmat& sci_func_chain, 
          realmat& sci_prior_chain, 
          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)
{

  #if defined SCI_MCMC_DEBUG
    cout << "\n\n Debug: Entering sci_mcmc.draw with analytic_mle = " 
         << boolalpha << analytic_mle << '\n';
  #endif

  INT_32BIT sub_chain_seed = seed;
  string errmsg = "errmsg not set #3";
  num_sub_chain_saved = 0;

  INTEGER len_sci_parm = sci_mod.len_parm();
  INTEGER len_sci_func = sci_mod.len_func();
  INTEGER len_stat_parm = stat_mod.len_parm();
  INTEGER len_stat_func = stat_mod.len_func();

  INTEGER sci_sub_total = len_sub_chain*num_sub_chain_to_save;
  sci_sub_parm_chain.resize(len_stat_parm,sci_sub_total,0.0);
  sci_sub_logl_chain.resize(len_sub_chain,num_sub_chain_to_save,0.0);
  sci_sub_rej_chain.resize(len_stat_parm+1,num_sub_chain_to_save,0.0);

  if (sci_parm_start.get_rows() != len_sci_parm) {
    error("Error, sci_mcmc, bad sci_parm_start or sci_mod");
  }

  if (sci_parm_start.get_cols() != 1) {
    error("Error, sci_mcmc, sci_parm_start has more than one column");
  }

  if (len_sci_func <= 0) {
    string err = "Error, sci_mcmc, len_sci_func cannot be zero.\n";
    err += "Recode sci_mod to return realmat func(1,1,-REAL_MAX).";
    error(err); 
  }

  if (stat_parm_start.get_rows() != len_stat_parm) {
    error("Error, sci_mcmc, bad stat_parm_start or stat_mod");
  }
 
  if (stat_parm_start.get_cols() != 1) {
    error("Error, sci_mcmc, stat_parm_start has more than one column");
  }

  if (len_stat_func <= 0) {
    string err = "Error, sci_mcmc, stat_mod.len_func cannot be zero.\n";
    err += "Recode stat_mod to return realmat func(1,1,-REAL_MAX).";
    error(err);
  }

  INTEGER rv;

  if (analytic_mle) {
    rv = make_sci_val(sci_parm_start, stat_parm_start, mode, errmsg); 
    if (mode.stat_mod_parm.get_rows() != len_stat_parm) {
      error("Error, sci_mcmc, stat_mod.mle() appears to be coded wrong");
    }
  }
  else {
    rv = make_sci_val(sub_chain_seed, 
           sci_parm_start, stat_parm_start,
           mode, errmsg, sci_sub_parm_chain, sci_sub_logl_chain,
           sci_sub_rej_chain);
  
    if (num_sub_chains > 0) {
      for (INTEGER run = 1; run <= num_sub_chains; ++run) {
        realmat mode_sci_mod_parm = mode.sci_mod_parm;
        if (rv == 0) stat_parm_start = mode.stat_mod_parm;
        rv = make_sci_val(sub_chain_seed, 
               mode_sci_mod_parm, stat_parm_start,
               mode, errmsg, sci_sub_parm_chain, sci_sub_logl_chain,
               sci_sub_rej_chain);
      }
    }
  }

  if (rv) error("Error, sci_mcmc, " + errmsg);

  sci_val chk = mode;
  
  if (imap.empty()) {
    goto map_good;
  }
  else {
    chk = imap.get_mode();
  }

  if (mode.sci_mod_parm.get_rows() != chk.sci_mod_parm.get_rows()) {
    warn("Warning, sci_mcmc, bad implied map was cleared, bad len_sci_parm");
    imap.clear();
    goto map_good;
  }

  if (mode.stat_mod_parm.get_rows() != chk.stat_mod_parm.get_rows()) {
    warn("Warning, sci_mcmc, bad implied map was cleared, bad len_stat_parm");
    imap.clear();
    goto map_good;
  }

  stat_mod.set_data_ptr(data_ptr);
  stat_mod.set_parm(chk.stat_mod_parm);

  /*
  Setting the parameter before checking support is intentional.
  Do not change it.  SNP will not work if you do.
  */

  if (!stat_mod.support(chk.stat_mod_parm)) {
    warn("Warning, sci_mcmc, bad implied map was cleared, bad stat_support");
    imap.clear();
    goto map_good;
  }
  
  if (stat_mod.loglikelihood() != chk.stat_mod_logl) {
    warn("Warning, sci_mcmc, bad implied map was cleared, bad stat_logl");
    imap.clear();
    goto map_good;
  }

  map_good: 
  
  map_val mv;
  if (imap.get_val(mode.sci_mod_parm, mv)) {
    sci_val scv(mode.sci_mod_parm, mv);
    if (scv.sci_mod_support && scv.sci_mod_simulate 
      && scv.sci_mod_prior.positive && scv.stat_mod_sub_logl.positive
      && (scv.stat_mod_sub_logl.log_den > mode.stat_mod_sub_logl.log_den)
      && scv.stat_mod_logl.positive && scv.stat_mod_support
      && scv.stat_mod_simulate && scv.stat_mod_prior.positive) {
      mode = scv;
    }
  }

  imap.set_val(mode);

  den_val pi_mode = mode.sci_mod_prior;
  pi_mode += mode.stat_mod_prior;

  if (draw_from_posterior) {
    pi_mode += mode.stat_mod_logl;
  }

  sci_val val_old = mode;
  den_val pi_old = pi_mode;

  sci_parm_chain.resize(len_sci_parm,len_chain); 
  sci_func_chain.resize(len_sci_func,len_chain);
  sci_prior_chain.resize(1,len_chain);
  stat_parm_chain.resize(len_stat_parm,len_chain); 
  stat_func_chain.resize(len_stat_func,len_chain);
  stat_prior_chain.resize(1,len_chain); 
  stat_logl_chain.resize(1,len_chain);
  stat_sub_logl_chain.resize(1,len_chain);

  realmat reject(len_sci_parm+1,4,0.0);

  sci_val val_new = val_old;
  den_val pi_new = pi_old;
 
  #if defined SCI_MCMC_DEBUG
    INTEGER loop_count = 0;
    INTEGER max_loop_count = 10;
  #endif

  /* Not a great idea
  Statements like the following statement sprinkled throughout can 
  cause the subchain and hence the sci chain to hang. It seems like 
  a good idea to start at previous peaks, but it isn't.

    stat_parm_start = mode.stat_mod_parm;

  Not resetting stat_parm_start has the effect of using the last
  value of the previous sub chain as the start for the next.
  */

  INTEGER sticky_chain_count = 0;

  realmat sci_mod_sim;
  stat_mcmc subchain(stat_mod, stat_prop, aprior);
  subchain.set_data_ptr(&sci_mod_sim);
  subchain.set_use_stat_mod_prior(false);
  subchain.set_use_assess_prior(false);
  subchain.set_draw_from_posterior(true);
  subchain.set_len_chain(len_sub_chain);
  realmat stat_parm_subchain;
  realmat stat_func_subchain;
  realmat stat_prior_subchain;
  realmat stat_logl_subchain;
  realmat assess_prior_subchain;
  realmat stat_mod_sim;
  stat_mod_eqns polish_eqns(stat_mod);
  nlopt minimizer(polish_eqns);
  minimizer.set_warning_messages(false);
  minimizer.set_check_derivatives(false);
  minimizer.set_output(false);

  sci_val val_prn = val_old;

  #if defined BEGINNING_AND_ENDING_DEBUG
    string sb = val_prn.annotated_sci_val();
    warn("\nDebug, sci_mcmc, sci_val at beginning of chain\n" + sb);
  #endif

  for (INTEGER t=1; t <= len_chain; ++t) {
   
    #if defined SCI_MCMC_DEBUG
      ++loop_count;
      if (loop_count == len_chain || loop_count < max_loop_count) {
        cout << "\n Debug: val_old within chain\n\n" << val_old << endl;
        cout << "\n Debug: pi_old within chain\n\n\t " << pi_old << endl;
      }
    #endif

    sci_prop.draw(seed, val_old.sci_mod_parm, val_new.sci_mod_parm);

    #if defined SCI_MCMC_DEBUG
      ++loop_count;
      if (loop_count == len_chain || loop_count < max_loop_count) {
        cout << "\n Debug: val_old within chain\n\n" << val_old << '\n';
        cout << "\n Debug: pi_old within chain\n\n\t " << pi_old << '\n';
        cout << "\n Debug: val_old.sci_mod_parm, val_new.sci_mod_parm, diff";
        cout << val_old.sci_mod_parm << val_new.sci_mod_parm 
             << val_old.sci_mod_parm - val_new.sci_mod_parm << endl;
      }
    #endif

    /*
    Either there is a sci_val for the proposed sci_mod_parm in the
    implied map or it has to be constructed and put there.  One or 
    the other is the next task.
    */

    map_val mv;
    if (imap.get_val(val_new.sci_mod_parm, mv)) {
      sci_val scv(val_new.sci_mod_parm, mv);
      val_new = scv;
    }
    else {
      realmat val_new_sci_mod_parm = val_new.sci_mod_parm;
      if (analytic_mle) {
        rv = make_sci_val(val_new_sci_mod_parm,stat_parm_start,val_new,errmsg);
      }
      else {
        rv = make_sci_val(sub_chain_seed, 
               val_new_sci_mod_parm, stat_parm_start,
               val_new, errmsg, sci_sub_parm_chain, sci_sub_logl_chain,
               sci_sub_rej_chain);
  
        if (num_sub_chains > 0) {
          for (INTEGER run = 1; run <= num_sub_chains; ++run) {
            val_new_sci_mod_parm = val_new.sci_mod_parm;
            if (rv == 0) stat_parm_start = val_new.stat_mod_parm;
            rv = make_sci_val(sub_chain_seed, 
                   val_new_sci_mod_parm, stat_parm_start,
                   val_new, errmsg, sci_sub_parm_chain, sci_sub_logl_chain,
                   sci_sub_rej_chain);
          }
        }
      }

      if (rv) {
        string msg = "Warning, sci_mcmc, " + errmsg + " at sci_mod_parm\n";
        for (INTEGER i=1; i<=val_new.sci_mod_parm.size(); ++i) {
          msg += fmt('f',15,5,val_new.sci_mod_parm[i])(); 
          msg += fmt('e',27,16,val_new.sci_mod_parm[i])(); 
          msg += "\n";
        }
        warn(msg);
      }
      imap.set_val(val_new);
    }

    bool move = false;

    bool move_possible = val_new.sci_mod_support
                         && val_new.sci_mod_simulate
                         && val_new.sci_mod_prior.positive
                         && val_new.stat_mod_logl.positive
                         && val_new.stat_mod_support
                         && val_new.stat_mod_simulate
                         && val_new.stat_mod_prior.positive;

    #if defined SCI_MCMC_DEBUG
      if (loop_count == len_chain || loop_count < max_loop_count) {
        cout << "\n Debug: val_new within chain\n\n" << val_new << endl;
        
        pi_new = val_new.sci_mod_prior;
        pi_new += val_new.stat_mod_prior;;

        if (draw_from_posterior) {
          pi_new += val_new.stat_mod_logl;
          cout << "\n Debug: draw is for posterior" << endl;
        }
        else {
          cout << "\n Debug: pi is for prior only" << endl;
        }
        cout << " Debug: pi_new within chain\n\n\t " << pi_new << endl;
      }
    #endif

    if (move_possible) {

      bool stuck = (num_sub_chains > 0)||(sticky_chain_count > 2*len_sci_parm);

      if ( stuck && (!analytic_mle) ) {
        
        #if defined STICKY_CHAIN_DEBUG
          warn("\nNow at t = " + fmt('d',6,t).get_ostr());
        #endif

        remake_sci_val_pair(sub_chain_seed, val_old, val_new,
          sci_sub_parm_chain, sci_sub_logl_chain, sci_sub_rej_chain);
      
        pi_old = val_old.sci_mod_prior;
        pi_old += val_old.stat_mod_prior;

        if (draw_from_posterior) {
          pi_old += val_old.stat_mod_logl;
        }

        if (pi_old.log_den > pi_mode.log_den) {
          mode = val_old;
          pi_mode = pi_old;
        }
      }  

      pi_new = val_new.sci_mod_prior;
      pi_new += val_new.stat_mod_prior;

      if (draw_from_posterior) {
        pi_new += val_new.stat_mod_logl;
      }

      if (pi_new.log_den > pi_mode.log_den) {
        mode = val_new;
        pi_mode = pi_new;
      }

      den_val top = pi_new;
      den_val bot = pi_old;
      if (!sci_prop.transition_is_symmetric()) {
        top += sci_prop(val_new.sci_mod_parm, val_old.sci_mod_parm);
        bot += sci_prop(val_old.sci_mod_parm, val_new.sci_mod_parm);
      }

      REAL diff = top.log_den - bot.log_den;
      REAL r = (0.0 < diff) ? 1.0 : exp(diff);
      REAL u = scl::ran(seed);
      if (u <= r) move = true;

      #if defined STICKY_CHAIN_DEBUG
        REAL pdiff = 0.0;
        pdiff += val_new.sci_mod_prior.log_den + val_new.stat_mod_prior.log_den;
        pdiff -= val_old.sci_mod_prior.log_den + val_old.stat_mod_prior.log_den;
        REAL ldiff = 0.0;
        ldiff += val_new.stat_mod_logl.log_den - val_old.stat_mod_logl.log_den;
        string strpdiff = fmt('e',27,16,pdiff).get_ostr();
        string strldiff = fmt('e',27,16,ldiff).get_ostr();
        string strdiff = fmt('e',27,16,diff).get_ostr();
        string strr = fmt('e',27,16,r).get_ostr();
        string stru = fmt('e',27,16,u).get_ostr();
        warn("accept/reject computation (assumes symmetric proposal density)");
        warn(" pdiff" + strpdiff);
        warn(" ldiff" + strldiff);
        warn(" diff " + strdiff);
        warn(" r    " + strr);
        warn(" u    " + stru);
        if (move) warn(" move      true"); else warn(" move      false");
      #endif

      #if defined SCI_MCMC_DEBUG
        if (loop_count == len_chain || loop_count < max_loop_count) {
          cout<<" Debug: diff, r, u = "<<diff<<", "<< r <<", "<< u << endl;
        }
      #endif

    } 

    if (move) {
       sticky_chain_count = 0;

      #if defined SCI_MCMC_DEBUG
        if (loop_count == len_chain || loop_count < max_loop_count) {
          cout << "\n Debug: Move accepted" << endl;
        }
      #endif

      val_prn = val_new;

      for (INTEGER i=1; i<=len_sci_parm; ++i) {
        sci_parm_chain(i,t) = val_new.sci_mod_parm[i];
      }
      for (INTEGER i=1; i<=len_sci_func; ++i) {
        sci_func_chain(i,t) = val_new.sci_mod_func[i];
      }
      sci_prior_chain[t] = val_new.sci_mod_prior.log_den;

      for (INTEGER i=1; i<=len_stat_parm; ++i) {
        stat_parm_chain(i,t) = val_new.stat_mod_parm[i];
      }
      for (INTEGER i=1; i<=len_stat_func; ++i) {
        stat_func_chain(i,t) = val_new.stat_mod_func[i];
      }
      stat_prior_chain[t] = val_new.stat_mod_prior.log_den;
      
      stat_logl_chain[t] = val_new.stat_mod_logl.log_den;

      stat_sub_logl_chain[t] = val_new.stat_mod_sub_logl.log_den;

      reject(len_sci_parm+1,4) += 1;
      for (INTEGER i=1; i<=len_sci_parm; ++i) {
        if (val_old.sci_mod_parm[i] != val_new.sci_mod_parm[i]) {
          reject(i,4) += 1;
        }
      }

      val_old = val_new;
      pi_old = pi_new;

      /* Not a great idea
      stat_parm_start = val_new.stat_mod_parm;
      */
    }
    else {
      ++sticky_chain_count;
        
      #if defined SCI_MCMC_DEBUG
        if (loop_count == len_chain || loop_count < max_loop_count) {
          cout << "\n Debug: Move rejected" << endl;
        }
      #endif

      val_prn = val_old;

      for (INTEGER i=1; i<=len_sci_parm; ++i) {
        sci_parm_chain(i,t) = val_old.sci_mod_parm[i];
      }
      for (INTEGER i=1; i<=len_sci_func; ++i) {
        sci_func_chain(i,t) = val_old.sci_mod_func[i];
      }
      sci_prior_chain[t] = val_old.sci_mod_prior.log_den;

      for (INTEGER i=1; i<=len_stat_parm; ++i) {
        stat_parm_chain(i,t) = val_old.stat_mod_parm[i];
      }
      for (INTEGER i=1; i<=len_stat_func; ++i) {
        stat_func_chain(i,t) = val_old.stat_mod_func[i];
      }
      stat_prior_chain[t] = val_old.stat_mod_prior.log_den;

      stat_logl_chain[t] = val_old.stat_mod_logl.log_den;

      stat_sub_logl_chain[t] = val_old.stat_mod_sub_logl.log_den;

      reject(len_sci_parm+1,4) += 1;
      reject(len_sci_parm+1,3) += 1;
      for (INTEGER i=1; i<=len_sci_parm; ++i) {
        if (val_old.sci_mod_parm[i] != val_new.sci_mod_parm[i]) {
          reject(i,4) += 1;
          reject(i,3) += 1;
        }
      }

      /* Not a great idea
      stat_parm_start = val_old.stat_mod_parm;
      */
    }

    #if defined SCI_MCMC_DEBUG
       if (loop_count == len_chain || loop_count < max_loop_count) {
         cout << "\n Debug: reject within chain\n\n\t " << reject << endl;
       }
    #endif
  }

  for (INTEGER i=1; i<=len_sci_parm; ++i) {
    sci_parm_start[i] = sci_parm_chain(i,len_chain);
  }

  /* Not a good idea
  for (INTEGER i=1; i<=len_stat_parm; ++i) {
    stat_parm_start[i] = stat_parm_chain(i,len_chain);
  }
  */

  for (INTEGER i=1; i<=len_sci_parm+1; ++i) {
    REAL bot = reject(i,4);
    if (bot>0.0) reject(i,1) = reject(i,3)/bot;
    bot = reject(len_sci_parm+1,4);
    if (bot>0.0) reject(i,2) = reject(i,4)/bot;
  }

  #if defined BEGINNING_AND_ENDING_DEBUG
    string se = val_prn.annotated_sci_val();
    warn("\nDebug, sci_mcmc, sci_val at end of chain\n" + se);
  #endif

  return reject;
}

realmat libgsm::stat_mcmc::draw(INT_32BIT& seed, realmat& stat_parm_start,
          realmat& stat_parm_chain, realmat& stat_func_chain, 
          realmat& stat_prior_chain, realmat& stat_logl_chain,
          realmat& assess_prior_chain) 
{
  #if defined STAT_MCMC_DEBUG
    cout << "\n\n Debug: Entering stat_mcmc.draw\n";
  #endif

  INTEGER len_parm = stat_mod.len_parm();
  INTEGER len_func = stat_mod.len_func();

  if (stat_parm_start.get_rows() != len_parm) {
    error("Error, stat_mcmc, bad stat_parm_start or stat_mod");
  }
  if (stat_parm_start.get_cols() != 1) {
    error("Error, stat_mcmc, stat_parm_start has more than one column");
  }

  if (len_func <= 0) {
    string err = "Error, stat_mcmc, len_func cannot be zero.\n";
    err += "Recode stat_mod to return realmat func(1,1,-REAL_MAX).";
    error(err); 
  }

  mode.stat_mod_parm = stat_parm_start;

  stat_mod.set_data_ptr(data_ptr);
  stat_mod.set_parm(mode.stat_mod_parm);

  /*
  Setting the parameter before checking support is intentional.
  Do not change it.  SNP will not work if you do.
  */

  mode.stat_mod_support = stat_mod.support(mode.stat_mod_parm);

  if (!mode.stat_mod_support) {
    error("Error, stat_mcmc, stat_parm not in support of stat_mod");
  }

  mode.stat_mod_logl = stat_mod.loglikelihood();

  if (use_stat_mod_prior) {
    realmat stat_mod_sim;
    mode.stat_mod_simulate = stat_mod.gen_sim(stat_mod_sim, mode.stat_mod_func);
  
    if (!mode.stat_mod_simulate) {
      error("Error, stat_mcmc, stat_mod simulation failed");
    }
  
    mode.stat_mod_prior=stat_mod.prior(mode.stat_mod_parm, mode.stat_mod_func);
  
    if (!mode.stat_mod_prior.positive) {
      error("Error, gsm_mcmc, map_val not in support of stat_mod_prior");
    }
  } 
  else {
    mode.stat_mod_simulate = true;
    mode.stat_mod_func.resize(len_func,1,-REAL_MAX);
    mode.stat_mod_prior.positive = true;
    mode.stat_mod_prior.log_den = 0.0;
  }

  if (use_assess_prior) {
    mode.assess_prior = aprior(mode.stat_mod_parm);
  }
  else {
    mode.assess_prior.positive = true;
    mode.assess_prior.log_den = 0.0;
  }

  den_val pi_mode = mode.stat_mod_prior;
  pi_mode += mode.assess_prior;

  if (draw_from_posterior) {
    pi_mode += mode.stat_mod_logl;
  }
  
  stat_parm_chain.resize(len_parm,len_chain);
  stat_func_chain.resize(len_func,len_chain);
  stat_prior_chain.resize(1,len_chain);
  stat_logl_chain.resize(1,len_chain);
  assess_prior_chain.resize(1,len_chain);

  realmat reject(len_parm+1,4,0.0);

  stat_val val_old = mode;
  den_val pi_old = pi_mode;

  stat_val val_new = val_old;
  den_val pi_new = pi_old;

  #if defined STAT_MCMC_DEBUG
    cout << "\n\t N = " << data_ptr->get_cols() << '\n';
    cout << "\n Debug: first six observations = "<<(*data_ptr)("","1:6")<< endl;
    intvec ivec = seq(data_ptr->get_cols()-5,data_ptr->get_cols());
    cout << "\n Debug: last six observation = " << (*data_ptr)("",ivec) << '\n';
    cout << "\n Debug: val_old prior to running chain\n\n" << val_old << endl;
    cout << "\n Debug: pi_old prior to running chain\n\n\t " << pi_old << endl;
    cout << "\n Debug: pi_old parts " 
         << "\n\t\t logl            " << mode.stat_mod_logl 
         << "\n\t\t stat_mod_prior  " << mode.stat_mod_prior
         << "\n\t\t assess_prior    " << mode.assess_prior 
         << endl;
    INTEGER loop_count = 0;
    INTEGER max_loop_count = 10;
  #endif

  for (INTEGER t=1; t <= len_chain; t++) {

    #if defined STAT_MCMC_DEBUG
      ++loop_count;
      if (loop_count == len_chain || loop_count < max_loop_count) {
        cout << "\n Debug: val_old within chain\n\n" << val_old << endl;
        cout << "\n Debug: pi_old within chain\n\n\t " << pi_old << endl;
      }
    #endif

    stat_prop.draw(seed, val_old.stat_mod_parm, val_new.stat_mod_parm);

    stat_mod.set_parm(val_new.stat_mod_parm);

    /*
    Setting the parameter before checking support is intentional.
    Do not change it.  SNP will not work if you do.
    */

    val_new.stat_mod_support = stat_mod.support(val_new.stat_mod_parm);

    bool move = false;

    if (val_new.stat_mod_support) {
      
      val_new.stat_mod_logl = stat_mod.loglikelihood();

      if (use_stat_mod_prior) {
        realmat stat_mod_sim;
        val_new.stat_mod_simulate 
          = stat_mod.gen_sim(stat_mod_sim, val_new.stat_mod_func);
        if (val_new.stat_mod_simulate) {
          val_new.stat_mod_prior
            = stat_mod.prior(val_new.stat_mod_parm, val_new.stat_mod_func);
        }
        else {
          val_new.stat_mod_prior.positive = false;
          val_new.stat_mod_prior.log_den = -REAL_MAX;
        }
      }
      else {
        val_new.stat_mod_simulate = true;
        val_new.stat_mod_func.resize(len_func,1,-REAL_MAX);
        val_new.stat_mod_prior.positive = true;
        val_new.stat_mod_prior.log_den = 0.0;
      }

      if (use_assess_prior) {
        val_new.assess_prior = aprior(val_new.stat_mod_parm);
      }
      else {
        val_new.assess_prior.positive = true;
        val_new.assess_prior.log_den = 0.0;
      }

      pi_new = val_new.stat_mod_prior;
      pi_new += val_new.assess_prior;

      if (draw_from_posterior) {
        pi_new += val_new.stat_mod_logl;
      }

      #if defined STAT_MCMC_DEBUG
        if (loop_count == len_chain || loop_count < max_loop_count) {
          cout << "\n Debug: val_new within chain\n\n" << val_new << endl;
          cout << "\n Debug: pi_new within chain\n\n\t " << pi_new << endl;
          cout << "\n Debug: pi_new parts " 
               << "\n\t\t logl            " << val_new.stat_mod_logl 
               << "\n\t\t stat_mod_prior  " << val_new.stat_mod_prior
               << "\n\t\t assess_prior    " << val_new.assess_prior 
               << endl;
        }
      #endif

      if (pi_new.positive) {
      
        if (pi_new.log_den > pi_mode.log_den) {
          mode = val_new;
          pi_mode = pi_new;
        }

        den_val top = pi_new;
        den_val bot = pi_old;
        if (!stat_prop.transition_is_symmetric()) {
          top += stat_prop(val_new.stat_mod_parm, val_old.stat_mod_parm);
          bot += stat_prop(val_old.stat_mod_parm, val_new.stat_mod_parm);
        }

        REAL diff = top.log_den - bot.log_den;
        REAL r = (0.0 < diff) ? 1.0 : exp(diff);
        REAL u = scl::ran(seed);
        if (u <= r) move = true;

        #if defined STAT_MCMC_DEBUG
          if (loop_count == len_chain || loop_count < max_loop_count) {
            cout<<" Debug: diff, r, u = "<<diff<<", "<< r <<", "<< u << endl;
          }
        #endif
      }
    }

    if (move) {
      #if defined STAT_MCMC_DEBUG
        if (loop_count == len_chain || loop_count < max_loop_count) {
          cout << "\n Debug: Move accepted" << endl;
        }
      #endif

      for (INTEGER i=1; i<=len_parm; ++i) {
        stat_parm_chain(i,t) = val_new.stat_mod_parm[i];
      }
      for (INTEGER i=1; i<=len_func; ++i) {
        stat_func_chain(i,t) = val_new.stat_mod_func[i];
      }
      stat_prior_chain[t] = val_new.stat_mod_prior.log_den;
      
      stat_logl_chain[t] = val_new.stat_mod_logl.log_den;
      
      assess_prior_chain[t] = val_new.assess_prior.log_den;

      reject(len_parm+1,4) += 1;
      for (INTEGER i=1; i<=len_parm; ++i) {
        if (val_old.stat_mod_parm[i] != val_new.stat_mod_parm[i]) {
          reject(i,4) += 1;
        }
      }

      val_old = val_new;
      pi_old = pi_new;
    }
    else {
      #if defined STAT_MCMC_DEBUG
        if (loop_count == len_chain || loop_count < max_loop_count) {
          cout << "\n Debug: Move rejected" << endl;
        }
      #endif

      for (INTEGER i=1; i<=len_parm; ++i) {
        stat_parm_chain(i,t) = val_old.stat_mod_parm[i];
      }
      for (INTEGER i=1; i<=len_func; ++i) {
        stat_func_chain(i,t) = val_old.stat_mod_func[i];
      }

      stat_prior_chain[t] = val_old.stat_mod_prior.log_den;

      stat_logl_chain[t] = val_old.stat_mod_logl.log_den;

      assess_prior_chain[t] = val_old.assess_prior.log_den;

      reject(len_parm+1,4) += 1;
      reject(len_parm+1,3) += 1;
      for (INTEGER i=1; i<=len_parm; ++i) {
        if (val_old.stat_mod_parm[i] != val_new.stat_mod_parm[i]) {
          reject(i,4) += 1;
          reject(i,3) += 1;
        }
      }
   }

   #if defined STAT_MCMC_DEBUG
      if (loop_count == len_chain || loop_count < max_loop_count) {
        cout << "\n Debug: reject within chain\n\n\t " << reject << endl;
      }
   #endif

  }

  for (INTEGER i=1; i<=len_parm; ++i) {
    stat_parm_start[i] = stat_parm_chain(i,len_chain);
  }

  for (INTEGER i=1; i<=len_parm+1; ++i) {
    REAL bot = reject(i,4);
    if (bot>0.0) reject(i,1) = reject(i,3)/bot;
    bot = reject(len_parm+1,4);
    if (bot>0.0) reject(i,2) = reject(i,4)/bot;
  }
  return reject;
}

void libgsm::sci_mcmc::set_len_chain(INTEGER N)
{
  if (N > 0) {
    len_chain = N;
  }
  else {
    error("Error, sci_mcmc, set_len_chain, N must be positive");
  }
}

void libgsm::sci_mcmc::set_len_sub_chain(INTEGER N)
{
  if (N > 0) {
    len_sub_chain = N;
  }
  else {
    error("Error, sci_mcmc, set_len_sub_chain, N must be positive");
  }
}

void libgsm::sci_mcmc::set_num_sub_chains(INTEGER N)
{
  if (N >= 0) {
    num_sub_chains = N;
  }
  else {
    error("Error, sci_mcmc, set_num_sub_chains, N cannot be negative");
  }
}

void libgsm::sci_mcmc::set_num_polish_iter(INTEGER N)
{
  if (N >= 0) {
    num_polish_iter = N;
  }
  else {
    error("Error, sci_mcmc, set_num_polish_iter, N must be zero or positive");
  }
}

void libgsm::sci_mcmc::set_polish_toler(REAL tol)
{
  if (tol > 0.0) {
    polish_toler = tol;
  }
  else {
    error("Error, sci_mcmc, polish_toler must be postitive");
  }
}

void libgsm::sci_mcmc::set_draw_from_posterior(bool from_posterior)
{
  draw_from_posterior = from_posterior;
}

void libgsm::sci_mcmc::set_analytic_mle(INTEGER am)
{
  analytic_mle = (am == 1 ? true : false);
}

void libgsm::stat_mcmc::set_len_chain(INTEGER N)
{
  if (N > 0) {
    len_chain = N;
  }
  else {
    error("Error, stat_mcmc, set_len_chain, N must be positive");
  }
}

void libgsm::stat_mcmc::set_draw_from_posterior(bool from_posterior)
{
  draw_from_posterior = from_posterior;
}

INTEGER libgsm::sci_mcmc::make_sci_val( // Used when analytic_mle == true 
  realmat& sci_parm_start, realmat& stat_parm_start,
  sci_val& sv, string& errmsg)
{

  if (!analytic_mle) error("Error, sci_mcmc::make_sci_val, erroneous call");

  errmsg = "errmsg not set #4";

  /*
  INTEGER len_sci_parm = sci_mod.len_parm();
  INTEGER len_sci_func = sci_mod.len_func();
  INTEGER len_stat_parm = stat_mod.len_parm();
  INTEGER len_stat_func = stat_mod.len_func();
  */

  #if defined SCI_MCMC_DEBUG
    cout << "\n\n Debug: Entering the sci_mcmc.make_sci_val that calls mle\n";
  #endif

  sci_val scv(sci_parm_start);
  sv = scv;

  sci_mod.set_parm(sv.sci_mod_parm);

  /*
  Setting the parameter before checking support is intentional.
  Do not change it.  SNP will not work if you do.
  */

  sv.sci_mod_support = sci_mod.support(sv.sci_mod_parm);

  if (!sv.sci_mod_support) {
    errmsg ="sci_parm_start not in support of sci_mod";
    return 1;
  }

  realmat sci_mod_sim;
  sv.sci_mod_simulate = sci_mod.gen_sim(sci_mod_sim, sv.sci_mod_func);

  if (!sv.sci_mod_simulate) {
    errmsg = "cannot simulate sci_mod at sci_parm_start";
    return 2;
  }

  sv.sci_mod_prior = sci_mod.prior(sv.sci_mod_parm, sv.sci_mod_func);
  
  if (!sv.sci_mod_prior.positive) {
    errmsg = "sci_parm_start not in support of sci_mod_prior";
    return 3;
  }

  stat_mod.set_parm(stat_parm_start);
  stat_mod.set_data_ptr(&sci_mod_sim);

  /*
  Setting the parameter before checking support is intentional.
  Do not change it.  SNP will not work if you do.
  */

  if (!stat_mod.support(stat_parm_start)) {
    error("Error, make_sci_val, stat_parm_start not in support of stat_mod");
  }

  #if defined SCI_MCMC_DEBUG
    cout << "\n Debug: sv before calling stat_mod::mle\n\n" << sv << endl;
    cout << "\n Debug: stat_parm_start before calling stat_mod::mle\n"
         << stat_parm_start << endl;
  #endif

  realmat stat_parm_mle;
  realmat stat_parm_V;
  INTEGER nobs = stat_mod.mle(stat_parm_mle, stat_parm_V);

  stat_mod.set_parm(stat_parm_mle);

  /*
  Setting the parameter before checking support is intentional.
  Do not change it.  SNP will not work if you do.
  */

  if (!stat_mod.support(stat_parm_mle)) {
    errmsg = "mle computation failed, support condition violated on sim data";
    return 4;
  }

  #if defined SCI_MCMC_DEBUG
    cout << "\n Debug: stat_parm after calling mle \n" << stat_parm_mle << endl;
  #endif

  sv.stat_mod_parm = stat_parm_mle;
  sv.stat_mod_sub_logl = stat_mod.loglikelihood();

  if (!sv.stat_mod_sub_logl.positive) {
    errmsg = "mle computation failed, stat_mod_logl not positive on sim data";
    return 4;
  }

  aprior.update_sigma(stat_parm_V,nobs);

  #if defined SCI_MCMC_DEBUG
    cout << "\n Debug: update_sigma called with nobs = " << nobs 
         << " and stat_parm_V = " << stat_parm_V << endl;
  #endif

  stat_mod.set_data_ptr(data_ptr);
  stat_mod.set_parm(sv.stat_mod_parm);

  /*
  Setting the parameter before checking support is intentional.
  Do not change it.  SNP will not work if you do.
  */

  sv.stat_mod_support = stat_mod.support(sv.stat_mod_parm);

  if (!sv.stat_mod_support) {
    errmsg = "stat_parm not in support of stat_mod on actual data";
    return 5;
  }

  sv.stat_mod_logl = stat_mod.loglikelihood();

  realmat stat_mod_sim;
  sv.stat_mod_simulate = stat_mod.gen_sim(stat_mod_sim, sv.stat_mod_func);

  if (!sv.stat_mod_simulate) {
    errmsg = "stat_mod simulation failed on actual data";
    return 6;
  }

  sv.stat_mod_prior = stat_mod.prior(sv.stat_mod_parm, sv.stat_mod_func);

  if (!sv.stat_mod_prior.positive) {
    errmsg = "stat mle not in support of stat_mod_prior on actual data";
    return 7;
  }

  #if defined SCI_MCMC_DEBUG
    cout << "\n Debug: sv returned by sci_mcmc.make_sci_val\n\n" << sv << endl;
    cerr << "\n Debug: sv.stat_mod_sub_logl" << sv.stat_mod_sub_logl;
    cerr << "\n Debug: sv.stat_mod_logl" << sv.stat_mod_logl << endl;
  #endif

  return 0;
}

