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

Copyright (C) 2005, 2006, 2008, 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.

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

#define  COMPILER_HAS_BOOLALPHA

#include "libsnp.h"
#include "libgsm.h"
#include "snp_base.h"
#include "gsm.h"

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


/*
None of the objects that need to be constructed on a slave node on a 
parallel machine can read files.  Thus, set_blks does not read files
even though that would be the more logical way to code it.
*/

namespace {

  /*
  char last_char(string s) 
  {
    char c = ' ';
    for (string::iterator i=s.begin(); i!=s.end(); ++i) {
      if (!isspace(*i)) c=*i;
    }
    return c;
  }
  */

  pair<string,string> find_ss(string s)
  {
    string zero = "0"; 
    string one = "1";
    char val = 'x'; 
    char* end = &val;
    char** endptr=&end;
    int lim = s.size();
    int b=0;
    while(b<lim && isspace(s[b])) ++b; 
    if (b==lim) {
      warn("Warning, gsmblocks, cannot interpret this line: \"" +s+ "\"");
      warn("Warning, gsmblocks, 0 1 substituted for \"" +s+ "\"");
      return pair<string,string>(zero,one);
    }
    int l=0;
    while(b+l<lim && !isspace(s[b+l])) ++l;
    string first=s.substr(b,l);
    b=b+l;
    while(b<lim && isspace(s[b])) ++b; 
    if (b==lim) { 
      warn("Warning, gsmblocks, cannot interpret this line: \"" +s+ "\"");
      warn("Warning, gsmblocks, " +first+ " 1 substituted for \"" +s+ "\"");
      strtod(first.c_str(),endptr);
      if (**endptr) error("Error, gsmblocks, cannot read this as numeric: "+s);
      return pair<string,string>(first,one);
    }
    l=0;
    while(b+l<lim && !isspace(s[b+l])) ++l;
    string second=s.substr(b,l);
    strtod(first.c_str(),endptr);
    if (**endptr) error("Error, gsmblocks, cannot read line as numeric: "+s);
    strtod(second.c_str(),endptr);
    if (**endptr) error("Error, gsmblocks, cannot read line as numeric: "+s);
    return pair<string,string>(first,second);
  }

  bool check_estblock(const estblock& eb) 
  {
    if (eb.pname.size() != 12) return false;
    if (eb.gsmver != 1.0 && eb.gsmver != 1.1 && eb.gsmver != 1.2 
        && eb.gsmver != 1.3 && eb.gsmver != 1.4 && eb.gsmver != 1.5
       ) return false;
    if (eb.kappa < 0.0) return false;
    if (eb.sci_draw_from_prior && !eb.run_sci_chain) return false;
    if (eb.stat_draw_from_prior && !eb.run_stat_chain) return false;
    return true;
  }
  
  bool check_datblock(const datblock& db) 
  {
    if (db.M <= 0) return false;
    if (db.n <= 0) return false;
    if (db.fields == intvec()) return false;
    return true;
  }

  bool check_statblock(const statblock& stb) 
  {
    if (stb.len_stat_parm <= 0) return false;
    if (stb.len_stat_func <= 0) return false;
    if (stb.len_stat_chain <= 0) return false;
    if (stb.num_stat_files < 0) return false;
    if (stb.stat_sclfac <= 0.0) return false;
    return true;
  }
  
  bool check_sciblock(sciblock& scb) 
  {
    if (scb.len_sci_parm <= 0) return false;
    if (scb.len_sci_func <= 0) return false;
    if (scb.len_sci_chain <= 0) return false;
    if (scb.num_sci_files < 0) return false;
    if (scb.analytic_mle != 1) {
      if (scb.len_sub_chain <= 0) return false;
      if (scb.num_sub_chains < 0) return false;
      /*
      if (scb.len_sub_chain < 50) {
        scb.len_sub_chain = 50;
        warn("Warning, check_sciblock, len_sci_sub_chain reset to 50");
      }
      */
      if (scb.num_polish_iter < 0) return false;
      if (scb.polish_toler <= 0) return false;
    }
    if (scb.sci_sclfac <= 0.0) return false;
    if (scb.sci_incfac <= 0.0) return false;
    const REAL x = 5.0;
    if ((x*scb.sci_incfac)/scb.sci_incfac != x) {
      warn("Warning, check_sciblock, incfac not a (fractional) power of two");
    }
    return true;
  }
  
}

void gsm::gsmblocks::set_stat_block(const statblock& s)
{  
  if (check_statblock(s)) {
    stat_blk = s;
  }
  else {
    error("Error, gsmblocks, set_stat_block, bad input");
  }
}

void gsm::gsmblocks::set_sci_block(const sciblock& s)
{  
  sciblock ss = s; 
  if (check_sciblock(ss)) {
    sci_blk = ss;
  }
  else {
    error("Error, gsmblocks, set_sci_block, bad input");
  }
}


void gsm::gsmblocks::set_stat_parm(const realmat& p)
{  
  if (p.get_rows() != stat_parm.get_rows() || p.get_cols() != 1) {
    error("Error, gsmblks, set_stat_parm, bad input");
  }
  else {
    stat_parm = p;
  }
}

void gsm::gsmblocks::set_sci_parm(const realmat& p)
{  
  if (p.get_rows() != sci_parm.get_rows() || p.get_cols() != 1) {
    error("Error, gsmblks, set_sci_parm, bad input");
  }
  else {
    sci_parm = p;
  }
}

bool gsm::gsmblocks::read_blks(const char* filename, ostream& detail)
{
  ifstream pfstream(filename);
  if (!pfstream) error("Error, gsmblocks, cannot open " + string(filename));

  vector<string> pf; 

  string line;
  while (getline(pfstream, line)) pf.push_back(line);

  return this->set_blks(pf, detail);
}

bool gsm::gsmblocks::set_blks(const vector<string>& pf, ostream& detail)
{
  history.erase(history.begin(),history.end());
  stat_parm_rhs.erase(stat_parm_rhs.begin(),stat_parm_rhs.end());
  stat_scale_rhs.erase(stat_scale_rhs.begin(),stat_scale_rhs.end());
  stat_groups_txt.erase(stat_groups_txt.begin(),stat_groups_txt.end());
  sci_parm_rhs.erase(sci_parm_rhs.begin(),sci_parm_rhs.end());
  sci_scale_rhs.erase(sci_scale_rhs.begin(),sci_scale_rhs.end());
  sci_incr_rhs.erase(sci_incr_rhs.begin(),sci_incr_rhs.end());
  sci_groups_txt.erase(sci_groups_txt.begin(),sci_groups_txt.end());
  stat_prop_def.erase(stat_prop_def.begin(),stat_prop_def.end());
  stat_alvec.erase(stat_alvec.begin(),stat_alvec.end());
  sci_prop_def.erase(sci_prop_def.begin(),sci_prop_def.end());
  sci_alvec.erase(sci_alvec.begin(),sci_alvec.end());
  
  typedef vector<string>::const_iterator pf_itr;

  snp::keyword kw;
  string header;

  pf_itr kw_ptr;
  
  header=kw.set_keyword("PARMFILE HISTORY");  //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
    kw_ptr += 8;    // Number of lines output history
    if (kw_ptr > pf.end()) error("Error, gsmblocks, " + header);
    while (kw_ptr != pf.end() && (*kw_ptr)[0] == '#') {
      history.push_back(*kw_ptr); 
      ++kw_ptr;
    }
  }
  
  header=kw.set_keyword("ESTIMATION DESCRIPTION");  //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  if (kw_ptr + 6 >= pf.end()) error("Error, gsmblocks, " + header);

  est_blk.pname = (++kw_ptr)->substr(0,12); 
  est_blk.gsmver = atof((++kw_ptr)->substr(0,12).c_str());
  est_blk.print = atoi((++kw_ptr)->substr(0,12).c_str()); 
  if (est_blk.gsmver >= 1.2) {
    est_blk.sci_draw_from_prior = atoi((++kw_ptr)->substr(0,12).c_str()); 
  }
  else {
    est_blk.sci_draw_from_prior = false;
  }
  est_blk.run_sci_chain = atoi((++kw_ptr)->substr(0,12).c_str()); 
  if (est_blk.gsmver >= 1.3) {
    est_blk.stat_draw_from_prior = atoi((++kw_ptr)->substr(0,12).c_str()); 
  }
  else {
    est_blk.stat_draw_from_prior = false;
  }
  est_blk.run_stat_chain = atoi((++kw_ptr)->substr(0,12).c_str()); 
  est_blk.kappa = atof((++kw_ptr)->substr(0,12).c_str());

  if (est_blk.print) {
    detail << starbox("/Input GSM parmfile//") << '\n';
    for (pf_itr itr=pf.begin(); itr!=pf.end(); ++itr) {
      detail << *itr << '\n';
    }
    detail.flush();
    detail << starbox("/GSM parmfile interpretation//") << '\n';
  }

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    #if defined COMPILER_HAS_BOOLALPHA
      detail << boolalpha;
    #endif
    detail << "\t pname = " << est_blk.pname << '\n';
    detail << "\t gsmver = " << fmt('f',3,1,est_blk.gsmver) << '\n';
    detail << "\t print = " << est_blk.print << '\n';
    detail << "\t sci_draw_from_prior = " << est_blk.sci_draw_from_prior
           << " (run_sci_chain must be true if true)\n";
    detail << "\t run_sci_chain = " << est_blk.run_sci_chain 
           << " (implied map from previous run required if false)\n";
    detail << "\t stat_draw_from_prior = " << est_blk.stat_draw_from_prior
           << " (run_stat_chain must be true if true)\n";
    detail << "\t run_stat_chain = " << est_blk.run_stat_chain 
           << " (this is the chain that uses the assess_prior)\n";
    detail << "\t kappa = " << est_blk.kappa << '\n';
    detail.flush();
  }

  if (!check_estblock(est_blk)) error("Error, bad or too few items in "+header);

  header=kw.set_keyword("DATA DESCRIPTION");  //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  if (kw_ptr + 4 >= pf.end()) error("Error, parmfile, " + header);

  dat_blk.M = atoi((++kw_ptr)->substr(0,12).c_str());
  dat_blk.n = atoi((++kw_ptr)->substr(0,12).c_str());
  string ss = *(++kw_ptr);
  dat_blk.dsn = cutstr(ss);
  dat_blk.fields.resize(dat_blk.M,0);
  ++kw_ptr;
  dat_blk.fieldline = *kw_ptr;
  INTEGER j=0;
  INTEGER i=1;
  while (i<=dat_blk.M) {
    while (isspace((*kw_ptr)[j])) ++j;
    INTEGER k=0;
    while (isdigit((*kw_ptr)[j+k])) ++k;
    if ((*kw_ptr)[j+k] != ':') {
      dat_blk.fields[i++] = atoi(kw_ptr->substr(j,k).c_str());
    }
    else {
      INTEGER start = atoi(kw_ptr->substr(j,k).c_str());
      j+=k+1;
      k=0;
      while (isdigit((*kw_ptr)[j+k])) ++k;
      INTEGER stop = atoi(kw_ptr->substr(j,k).c_str());
      for (INTEGER ii=start; ii<=stop; ++ii) {
        if (i <= dat_blk.M) dat_blk.fields[i++] = ii;
      }
    }
    j+=k;
  }
  if (dat_blk.fields[dat_blk.M]==0) {
    error("Error, parmfile, not enough fields this line:\n" + *kw_ptr);
  }
    
  if (est_blk.print) {
    detail << '\n' << header << '\n';
    detail << "\t M = " << dat_blk.M << '\n';
    detail << "\t n = " << dat_blk.n << '\n';
    detail << "\t dsn = " << dat_blk.dsn << '\n';
    detail << "\t fieldline: " << dat_blk.fieldline.substr(0,54) << '\n';
    detail << "\t fields: "; 
    for (INTEGER i=1; i<=dat_blk.fields.size(); ++i) {
      detail <<" "<< dat_blk.fields[i];
    }
    detail << '\n';
    detail.flush();
  }

  if (!check_datblock(dat_blk)) error("Error, bad or too few items in "+header);

  header=kw.set_keyword("STAT_MOD DESCRIPTION");  //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  if (kw_ptr + 2 >= pf.end()) error("Error, gsmblocks, " + header);

  stat_blk.len_stat_parm = atoi((++kw_ptr)->substr(0,12).c_str());
  stat_blk.len_stat_func = atoi((++kw_ptr)->substr(0,12).c_str());

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    #if defined COMPILER_HAS_BOOLALPHA
      detail << boolalpha;
    #endif
    detail << "\t len_stat_parm = " << stat_blk.len_stat_parm << '\n';
    detail << "\t len_stat_func = " << stat_blk.len_stat_func << '\n';
    detail.flush();
  }

  string h1 = header;

  header=kw.set_keyword("STAT_MOD MCMC DESCRIPTION");  //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  if (kw_ptr + 4 >= pf.end()) error("Error, gsmblocks, " + header);

  stat_blk.stat_seed = atoi((++kw_ptr)->substr(0,12).c_str());
  stat_blk.len_stat_chain = atoi((++kw_ptr)->substr(0,12).c_str());
  stat_blk.num_stat_files = atoi((++kw_ptr)->substr(0,12).c_str());
  stat_blk.stat_sclfac = atof((++kw_ptr)->substr(0,12).c_str());

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    #if defined COMPILER_HAS_BOOLALPHA
      detail << boolalpha;
    #endif
    detail << "\t stat_seed = " << stat_blk.stat_seed << '\n';
    detail << "\t len_stat_chain = " << stat_blk.len_stat_chain << '\n';
    detail << "\t num_stat_files = " << stat_blk.num_stat_files 
           << " (in addition to the first)\n";
    detail << "\t stat_sclfac = " << stat_blk.stat_sclfac << '\n';
    detail.flush();
  }

  string h2 = header;

  if (!check_statblock(stat_blk)) {
    error("Error, bad or too few items in\n" + h1 + " or " + h2);
  }

  header=kw.set_keyword("STAT_MOD PARMFILE");  //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  if (kw_ptr + 3 >= pf.end()) error("Error, gsmblocks, " + header);
  string line = *(++kw_ptr);
  stat_blk.stat_parmfile = cutstr(line);
  pf_itr beg_ptr = ++kw_ptr;
  if (stat_blk.stat_parmfile != string("__none__")) {
    stat_blk.is_stat_parmfile = true;
  }
  else {
    stat_blk.is_stat_parmfile = false;
  }
  snp::keyword endkw("#end additional lines");
  pf_itr end_ptr = find_if(beg_ptr, pf.end(), endkw);
  if (end_ptr==pf.end() || end_ptr<beg_ptr) {
    error("Error, gsmblocks, " + header);
  }
  ++end_ptr;
  for (kw_ptr=beg_ptr; kw_ptr != end_ptr; ++kw_ptr) {
    stat_alvec.push_back(*kw_ptr); 
  }

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    #if defined COMPILER_HAS_BOOLALPHA
      detail << boolalpha;
    #endif
    detail << "\t is_stat_parmfile = " << stat_blk.is_stat_parmfile << '\n';
    if (stat_blk.is_stat_parmfile) {
      detail << "\t stat_parmfile = " << stat_blk.stat_parmfile << '\n';
    }
    vector<string>::const_iterator itr;
    for (itr=stat_alvec.begin(); itr!=stat_alvec.end(); ++itr ) {
      detail << *itr << '\n';
    }
    detail.flush();
  }

  INTEGER len_stat_parm = stat_blk.len_stat_parm;
  stat_parm.resize(len_stat_parm,1);
  stat_mask.resize(len_stat_parm,1);
  
  INTEGER nfixed = 0;

  header=kw.set_keyword("STAT_MOD PARAMETER START VALUES");   //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  if (kw_ptr+len_stat_parm >= pf.end()) error("Error, gsmblocks, " + header);
  for (INTEGER i=1; i<=len_stat_parm; ++i) {
    pair<string,string> ss = find_ss((++kw_ptr)->c_str());
    stat_parm[i] = atof(ss.first.c_str());
    stat_mask[i] = atoi(ss.second.c_str());
    if (stat_mask[i] == 0) ++nfixed;
    string line = *kw_ptr;
    cutstr(line); cutstr(line);
    stat_parm_rhs.push_back(line);
  }
  
  INTEGER nactive = len_stat_parm - nfixed;

  intvec fixed;
  intvec active = seq(1,len_stat_parm);
  if (nfixed>0) {
    fixed.resize(nfixed);
    active.resize(nactive);
    INTEGER j=0;
    INTEGER k=0;
    for (INTEGER i=1; i<=len_stat_parm; ++i) {
      if (stat_mask[i] == 0) {
        fixed[++j]=i;
      }
      else {
        active[++k]=i;
      }
    }
  }

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    detail << "\n\t stat_parm = ";
    detail << stat_parm;
    detail << "\n\t stat_mask = ";
    detail << stat_mask << '\n';
    detail << "\t Number fixed = " << nfixed << '\n';
    detail << "\t Number active = " << nactive << '\n';
    detail << "\n\t fixed = ";
    detail << fixed << '\n';
    detail << "\n\t active = ";
    detail << active << '\n';
    detail << "\n\t stat_parm_rhs = \n\n";
    vector<string>::const_iterator itr;
    for (itr=stat_parm_rhs.begin(); itr!=stat_parm_rhs.end(); ++itr) {
      detail << "\t " << *itr << '\n';
    }
    detail.flush();
  }

  stat_scale.resize(len_stat_parm,1);
  
  header=kw.set_keyword("STAT_MOD PROPOSAL SCALING");   //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  if (kw_ptr + len_stat_parm >= pf.end()) error("Error, gsmblocks, " + header);
  for (INTEGER i=1; i<=len_stat_parm; ++i) {
    string line = *(++kw_ptr);
    stat_scale[i] = atof(cutstr(line).c_str());
    stat_scale[i] *= stat_blk.stat_sclfac;
    if (stat_scale[i]<=0) error("Error, gsmblocks, " + header);
    stat_scale_rhs.push_back(line);
  }

  stat_blk.stat_sclfac = 1.0;

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    detail << "(after multiplication by stat_sclfac)\n";
    detail << "(output parmfile this scale and sclfac=1.0)\n";
    detail << "\n\t stat_scale = ";
    detail << stat_scale;
    detail << "\n\t stat_scale_rhs = \n\n";
    vector<string>::const_iterator itr;
    for (itr=stat_scale_rhs.begin(); itr!=stat_scale_rhs.end(); ++itr) {
      detail << "\t " << *itr << '\n';
    }
    detail.flush();
  }

  if (nfixed > 0) {
    realmat ifixed(nfixed,1);
    realmat ufixed(nfixed,1);
    realmat Vfixed(nfixed,nfixed,0.0);
    for (INTEGER i=1; i<=nfixed; ++i) {
      ifixed[i] = 1.0;
      ufixed[i] = stat_parm[fixed[i]];
      Vfixed(i,i) = pow(stat_scale[fixed[i]],2);
    }
    stat_prop_def.push_back(prop_group(0.0,fixed,ifixed,ufixed,Vfixed));
  }
  
  header=kw.set_keyword("STAT_MOD PROPOSAL GROUPING");  //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
    INTEGER count = 0;
    while (count < len_stat_parm) {
      INTEGER top = len_stat_parm;
      INTEGER r=1;
      realmat gmat(len_stat_parm+1,len_stat_parm+1,0.0);
      while (r <= top) {
        ++kw_ptr;
        stat_groups_txt.push_back(*kw_ptr);
        if (kw_ptr == pf.end()) error("Error, gsmblocks, " + header); 
        INTEGER j=0;
        INTEGER c=0;
        INTEGER lim = kw_ptr->size();
        while(lim > 0 && isspace((*kw_ptr)[lim-1])) --lim;
        while (j < lim && c <= len_stat_parm) {
          while (j < lim && isspace((*kw_ptr)[j])) ++j;
          INTEGER k=0;
          while (j+k < lim && !isspace((*kw_ptr)[j+k])) ++k;
          gmat(r,++c) = atof(kw_ptr->substr(j,k).c_str());
          j+=k;
        }
        if (r != 1) ++count;
        if (r == 1) top = c;
        ++r;
      }
      intvec idx = seq(1,top);
      gmat = gmat(idx,idx);
      for (INTEGER i=1; i<=top; ++i) {
        if (gmat(1,i) != gmat(i,1)) error("Error, gsmblocks, " + header);
      }
      for (INTEGER i=2; i<=top; ++i) {
        for (INTEGER j=1; j<=fixed.size(); ++j) {
          if (INTEGER(gmat(1,idx[i])) == fixed[j]) idx[i] = -idx[i];
        }
      }
      realmat gmat_active = gmat(idx,idx);
      if (gmat_active.size() > 1) {
        idx[1] = -idx[1];
        REAL freq = gmat(1,1);
        realmat rvec = gmat(1,idx);
        realmat Vmat = gmat(idx,idx);
        intvec gvec(rvec.size());
        realmat ginc(rvec.size(),1);
        realmat mean(rvec.size(),1);
        for (INTEGER i=1; i<=gvec.size(); ++i) {
          gvec[i] = INTEGER(rvec[i]);
          ginc[i] = 1.0;
          mean[i] = stat_parm[gvec[i]];
        }
        for (INTEGER j=1; j<=gvec.size(); ++j) {
          for (INTEGER i=1; i<=gvec.size(); ++i) {
            Vmat(i,j) = Vmat(i,j)*stat_scale[gvec[i]]*stat_scale[gvec[j]];
          }
        }
        stat_prop_def.push_back(prop_group(freq, gvec, ginc, mean, Vmat));
      }
    }
  }
  else {
    for (INTEGER i=1; i<=nactive; ++i) {
      REAL freq = 1.0;
      intvec gvec(1,active[i]);
      realmat ginc(1,1, 1.0);
      realmat mean(1,1, stat_parm[active[i]]);
      realmat Vmat(1,1, pow(stat_scale[active[i]],2));
      stat_prop_def.push_back(prop_group(freq, gvec, ginc, mean, Vmat));
    }
  }

  if (est_blk.print) {
    detail << starbox("/Prop groups//") << '\n';
    detail << "\t Constructed from "<<header<<", if present; \n";
    detail << "\t If not, move one at a time grouping is used.\n";
    detail << "\t If parm has fixed elements, they are in Group 0\n";
    detail << "\t and not changed during the MCMC iterations.\n";
    detail << "\t "<<header<<" correlations have been converted\n";
    detail << "\t to variances using the prop_scale vector above.\n";
    detail << "\t Frequencies are relative, not absolute.\n\n";
    prop_def::const_iterator itr;
    INTEGER i=0;
    for (itr=stat_prop_def.begin(); itr!=stat_prop_def.end(); ++itr) {
      detail << "     Group " << i++ << "\n\n"; 
      detail << "\t freq = " << itr->freq << '\n';
      detail << "\t gvec = " << itr->gvec << '\n';
      detail << "\t ginc = " << itr->ginc << '\n';
      detail << "\t mean = " << itr->mean << '\n';
      detail << "\t Vmat = " << itr->Vmat << '\n';
    }
  }

  header=kw.set_keyword("SCI_MOD DESCRIPTION");  //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  if (kw_ptr + 2 >= pf.end()) error("Error, gsmblocks, " + header);

  sci_blk.len_sci_parm = atoi((++kw_ptr)->substr(0,12).c_str());
  sci_blk.len_sci_func = atoi((++kw_ptr)->substr(0,12).c_str());

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    #if defined COMPILER_HAS_BOOLALPHA
      detail << boolalpha;
    #endif
    detail << "\t len_sci_parm = " << sci_blk.len_sci_parm << '\n';
    detail << "\t len_sci_func = " << sci_blk.len_sci_func << '\n';
    detail.flush();
  }

  h1 = header;

  header=kw.set_keyword("SCI_MOD MCMC DESCRIPTION");  //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  
  if (est_blk.gsmver == 1.0) {
    if (kw_ptr + 6 >= pf.end()) error("Error, gsmblocks, " + header);
    sci_blk.sci_seed = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.len_sci_chain = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.num_sci_files = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.analytic_mle = 0;
    sci_blk.len_sub_chain = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.num_sub_chains = 0;
    sci_blk.num_polish_iter = 0;
    sci_blk.polish_toler = 1.0e-10;
    sci_blk.sci_sclfac = atof((++kw_ptr)->substr(0,12).c_str());
    sci_blk.sci_incfac = atof((++kw_ptr)->substr(0,12).c_str());
  }
  else if (est_blk.gsmver <= 1.3) {
    if (kw_ptr + 8 >= pf.end()) error("Error, gsmblocks, " + header);
    sci_blk.sci_seed = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.len_sci_chain = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.num_sci_files = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.analytic_mle = 0;
    sci_blk.len_sub_chain = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.num_sub_chains = 0;
    sci_blk.num_polish_iter = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.polish_toler = atof((++kw_ptr)->substr(0,12).c_str());
    sci_blk.sci_sclfac = atof((++kw_ptr)->substr(0,12).c_str());
    sci_blk.sci_incfac = atof((++kw_ptr)->substr(0,12).c_str());
  }
  else if (est_blk.gsmver <= 1.4) {
    if (kw_ptr + 9 >= pf.end()) error("Error, gsmblocks, " + header);
    sci_blk.sci_seed = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.len_sci_chain = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.num_sci_files = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.analytic_mle = 0;
    sci_blk.len_sub_chain = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.num_sub_chains = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.num_polish_iter = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.polish_toler = atof((++kw_ptr)->substr(0,12).c_str());
    sci_blk.sci_sclfac = atof((++kw_ptr)->substr(0,12).c_str());
    sci_blk.sci_incfac = atof((++kw_ptr)->substr(0,12).c_str());
  }
  else {
    if (kw_ptr + 10 >= pf.end()) error("Error, gsmblocks, " + header);
    sci_blk.sci_seed = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.len_sci_chain = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.num_sci_files = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.analytic_mle = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.len_sub_chain = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.num_sub_chains = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.num_polish_iter = atoi((++kw_ptr)->substr(0,12).c_str());
    sci_blk.polish_toler = atof((++kw_ptr)->substr(0,12).c_str());
    sci_blk.sci_sclfac = atof((++kw_ptr)->substr(0,12).c_str());
    sci_blk.sci_incfac = atof((++kw_ptr)->substr(0,12).c_str());
  }

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    #if defined COMPILER_HAS_BOOLALPHA
      detail << boolalpha;
    #endif
    detail << "\t sci_seed = " << sci_blk.sci_seed << '\n';
    detail << "\t len_sci_chain = " << sci_blk.len_sci_chain << '\n';
    detail << "\t num_sci_files = " << sci_blk.num_sci_files 
           << " (in addition to the first)\n";
    detail << "\t analytic_mle = " << sci_blk.analytic_mle << '\n';
    detail << "\t len_sub_chain = " << sci_blk.len_sub_chain << '\n';
    detail << "\t num_sub_chains = " << sci_blk.num_sub_chains << '\n';
    detail << "\t num_polish_iter = " << sci_blk.num_polish_iter << '\n';
    detail << "\t polish_toler = " << sci_blk.polish_toler << '\n';
    detail << "\t sci_sclfac = " << sci_blk.sci_sclfac << '\n';
    detail << "\t sci_incfac = " << sci_blk.sci_incfac << '\n';
    detail.flush();
  }

  h2 = header;

  if (!check_sciblock(sci_blk)) {
    error("Error, bad or too few items in\n" + h1 + " or " + h2);
  }

  header=kw.set_keyword("SCI_MOD PARMFILE");  //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  if (kw_ptr + 3 >= pf.end()) error("Error, gsmblocks, " + header);
  line = *(++kw_ptr);
  sci_blk.sci_parmfile = cutstr(line);
  beg_ptr = ++kw_ptr;
  if (sci_blk.sci_parmfile != string("__none__")) {
    sci_blk.is_sci_parmfile = true;
  }
  else {
    sci_blk.is_sci_parmfile = false;
  }
  endkw("#end additional lines");
  end_ptr = find_if(beg_ptr, pf.end(), endkw);
  if (end_ptr==pf.end() || end_ptr<beg_ptr) {
    error("Error, gsmblocks, " + header);
  }
  ++end_ptr;
  for (kw_ptr=beg_ptr; kw_ptr != end_ptr; ++kw_ptr) {
    sci_alvec.push_back(*kw_ptr); 
  }

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    #if defined COMPILER_HAS_BOOLALPHA
      detail << boolalpha;
    #endif
    detail << "\t is_sci_parmfile = " << sci_blk.is_sci_parmfile << '\n';
    if (sci_blk.is_sci_parmfile) {
      detail << "\t sci_parmfile = " << sci_blk.sci_parmfile << '\n';
    }
    vector<string>::const_iterator itr;
    for (itr=sci_alvec.begin(); itr!=sci_alvec.end(); ++itr ) {
      detail << *itr << '\n';
    }
    detail.flush();
  }

  INTEGER len_sci_parm = sci_blk.len_sci_parm;
  sci_parm.resize(len_sci_parm,1);
  sci_mask.resize(len_sci_parm,1);
  
  nfixed = 0;

  header=kw.set_keyword("SCI_MOD PARAMETER START VALUES");   //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  if (kw_ptr+len_sci_parm >= pf.end()) error("Error, gsmblocks, " + header);
  for (INTEGER i=1; i<=len_sci_parm; ++i) {
    pair<string,string> ss = find_ss((++kw_ptr)->c_str());
    sci_parm[i] = atof(ss.first.c_str());
    sci_mask[i] = atoi(ss.second.c_str());
    if (sci_mask[i] == 0) ++nfixed;
    string line = *kw_ptr;
    cutstr(line); cutstr(line);
    sci_parm_rhs.push_back(line);
  }
  
  nactive = len_sci_parm - nfixed;

  active = seq(1,len_sci_parm);
  if (nfixed>0) {
    fixed.resize(nfixed);
    active.resize(nactive);
    INTEGER j=0;
    INTEGER k=0;
    for (INTEGER i=1; i<=len_sci_parm; ++i) {
      if (sci_mask[i] == 0) {
        fixed[++j]=i;
      }
      else {
        active[++k]=i;
      }
    }
  }

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    detail << "\n\t sci_parm = ";
    detail << sci_parm;
    detail << "\n\t sci_mask = ";
    detail << sci_mask << '\n';
    detail << "\t Number fixed = " << nfixed << '\n';
    detail << "\t Number active = " << nactive << '\n';
    detail << "\n\t fixed = ";
    detail << fixed << '\n';
    detail << "\n\t active = ";
    detail << active << '\n';
    detail << "\n\t sci_parm_rhs = \n\n";
    vector<string>::const_iterator itr;
    for (itr=sci_parm_rhs.begin(); itr!=sci_parm_rhs.end(); ++itr) {
      detail << "\t " << *itr << '\n';
    }
    detail.flush();
  }

  sci_scale.resize(len_sci_parm,1);

  header=kw.set_keyword("SCI_MOD PROPOSAL SCALING");   //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  if (kw_ptr + len_sci_parm >= pf.end()) error("Error, gsmblocks, " + header);
  for (INTEGER i=1; i<=len_sci_parm; ++i) {
    string line = *(++kw_ptr);
    sci_scale[i] = atof(cutstr(line).c_str());
    sci_scale[i] *= sci_blk.sci_sclfac;
    if (sci_scale[i]<=0) error("Error, gsmblocks, " + header);
    sci_scale_rhs.push_back(line);
  }

  sci_blk.sci_sclfac = 1.0;

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    detail << "(after multiplication by sci_sclfac)\n";
    detail << "(output parmfile this scale and sclfac=1.0)\n";
    detail << "\n\t sci_scale = ";
    detail << sci_scale;
    detail << "\n\t sci_scale_rhs = \n\n";
    vector<string>::const_iterator itr;
    for (itr=sci_scale_rhs.begin(); itr!=sci_scale_rhs.end(); ++itr) {
      detail << "\t " << *itr << '\n';
    }
    detail.flush();
  }

  sci_incr.resize(len_sci_parm,1);
  
  header=kw.set_keyword("SCI_PARAMETER INCREMENTS");   //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) error("Error, " + header + " block not found");
  if (kw_ptr + len_sci_parm >= pf.end()) error("Error, gsmblocks, " + header);
  for (INTEGER i=1; i<=len_sci_parm; ++i) {
    string line = *(++kw_ptr);
    sci_incr[i] = atof(cutstr(line).c_str());
    sci_incr[i] *= sci_blk.sci_incfac;
    if (sci_incr[i]<=0) error("Error, gsmblocks, " + header);
    sci_incr_rhs.push_back(line);
  }

  sci_blk.sci_incfac = 1.0;
  
  if (est_blk.print) {
    detail << '\n' << header << '\n';
    detail << "(after multiplication by sci_sclfac)\n";
    detail << "(output parmfile this incr and incfac=1.0)\n";
    detail << "\n\t sci_incr = ";
    detail << sci_incr;
    detail << "\n\t sci_incr_rhs = \n\n";
    vector<string>::const_iterator itr;
    for (itr=sci_incr_rhs.begin(); itr!=sci_incr_rhs.end(); ++itr) {
      detail << "\t " << *itr << '\n';
    }
    detail.flush();
  }

  if (est_blk.print) {
    string msg = "\nImplied maximum move of grid_group_move proposal\n";
    msg += "formula is INTEGER(2.0*sci_scale[i]/sci_incr[i])\n\n";
    detail << starbox(msg.c_str(),'\n');
    intvec out(sci_scale.get_rows());
    for (INTEGER i=1; i<=sci_scale.get_rows(); ++i) {
      out[i] =INTEGER(2.0*sci_scale[i]/sci_incr[i]);
    }
    detail << out;
  }

  if (nfixed > 0) {
    realmat ifixed(nfixed,1);
    realmat ufixed(nfixed,1);
    realmat Vfixed(nfixed,nfixed,0.0);
    for (INTEGER i=1; i<=nfixed; ++i) {
      ifixed[i] = sci_incr[fixed[i]];
      ufixed[i] = sci_parm[fixed[i]];
      Vfixed(i,i) = pow(sci_scale[fixed[i]],2);
    }
    sci_prop_def.push_back(prop_group(0.0,fixed,ifixed,ufixed,Vfixed));
  }
  
  header=kw.set_keyword("SCI_MOD PROPOSAL GROUPING");  //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
    INTEGER count = 0;
    while (count < len_sci_parm) {
      INTEGER top = len_sci_parm;
      INTEGER r=1;
      realmat gmat(len_sci_parm+1,len_sci_parm+1,0.0);
      while (r <= top) {
        ++kw_ptr;
        sci_groups_txt.push_back(*kw_ptr);
        if (kw_ptr == pf.end()) error("Error, gsmblocks, " + header); 
        INTEGER j=0;
        INTEGER c=0;
        INTEGER lim = kw_ptr->size();
        while(lim > 0 && isspace((*kw_ptr)[lim-1])) --lim;
        while (j < lim && c <= len_sci_parm) {
          while (j < lim && isspace((*kw_ptr)[j])) ++j;
          INTEGER k=0;
          while (j+k < lim && !isspace((*kw_ptr)[j+k])) ++k;
          gmat(r,++c) = atof(kw_ptr->substr(j,k).c_str());
          j+=k;
        }
        if (r != 1) ++count;
        if (r == 1) top = c;
        ++r;
      }
      intvec idx = seq(1,top);
      gmat = gmat(idx,idx);
      for (INTEGER i=1; i<=top; ++i) {
        if (gmat(1,i) != gmat(i,1)) error("Error, gsmblocks, " + header);
      }
      for (INTEGER i=2; i<=top; ++i) {
        for (INTEGER j=1; j<=fixed.size(); ++j) {
          if (INTEGER(gmat(1,idx[i])) == fixed[j]) idx[i] = -idx[i];
        }
      }
      realmat gmat_active = gmat(idx,idx);
      if (gmat_active.size() > 1) {
        idx[1] = -idx[1];
        REAL freq = gmat(1,1);
        realmat rvec = gmat(1,idx);
        realmat Vmat = gmat(idx,idx);
        intvec gvec(rvec.size());
        realmat ginc(rvec.size(),1);
        realmat mean(rvec.size(),1);
        for (INTEGER i=1; i<=gvec.size(); ++i) {
          gvec[i] = INTEGER(rvec[i]);
          ginc[i] = sci_incr[gvec[i]];
          mean[i] = sci_parm[gvec[i]];
        }
        for (INTEGER j=1; j<=gvec.size(); ++j) {
          for (INTEGER i=1; i<=gvec.size(); ++i) {
            Vmat(i,j) = Vmat(i,j)*sci_scale[gvec[i]]*sci_scale[gvec[j]];
          }
        }
        sci_prop_def.push_back(prop_group(freq, gvec, ginc, mean, Vmat));
      }
    }
  }
  else {
    for (INTEGER i=1; i<=nactive; ++i) {
      REAL freq = 1.0;
      intvec gvec(1,active[i]);
      realmat ginc(1,1, sci_incr[active[i]]);
      realmat mean(1,1, sci_parm[active[i]]);
      realmat Vmat(1,1, pow(sci_scale[active[i]],2));
      sci_prop_def.push_back(prop_group(freq, gvec, ginc, mean, Vmat));
    }
  }

  if (est_blk.print) {
    detail << starbox("/Prop groups//") << '\n';
    detail << "\t Constructed from "<<header<<", if present; \n";
    detail << "\t If not, move one at a time grouping is used.\n";
    detail << "\t If parm has fixed elements, they are in Group 0\n";
    detail << "\t and not changed during the MCMC iterations.\n";
    detail << "\t "<<header<<" correlations have been converted\n";
    detail << "\t to variances using the prop_scale vector above.\n";
    detail << "\t Frequencies are relative, not absolute.\n\n";
    prop_def::const_iterator itr;
    INTEGER i=0;
    for (itr=sci_prop_def.begin(); itr!=sci_prop_def.end(); ++itr) {
      detail << "     Group " << i++ << "\n\n"; 
      detail << "\t freq = " << itr->freq << '\n';
      detail << "\t gvec = " << itr->gvec << '\n';
      detail << "\t ginc = " << itr->ginc << '\n';
      detail << "\t mean = " << itr->mean << '\n';
      detail << "\t Vmat = " << itr->Vmat << '\n';
    }
  }
  return true;
}


bool gsm::gsmblocks::write_blks
  (string parmfile, string prefix, string suffix) const
{
  string filename = prefix + ".parmfile." + suffix;
  ofstream ps(filename.c_str());
  if (!ps) return false; 

  REAL gsmver = 1.5;

  ps << "PARMFILE HISTORY (optional)\n"
     << "#\n"
     << "# This parmfile was written by GSM Version " 
     << fmt('f',3,1,gsmver) << ' '
     << "using the following line from\n" 
     << "# control.dat, which was read as " 
     << "char*, char*\n"
     << "# -----------------------------------------"
     << "----------------------------------\n";
  ps << '#';
  for (INTEGER i=1; i<=INTEGER(20-parmfile.size()); ++i) ps << ' ';
  ps << parmfile;
  for (INTEGER i=1; i<=INTEGER(20-prefix.size()); ++i) ps << ' ';
  ps << prefix << '\n';
  ps << "# -----------------------------------------"
     << "----------------------------------\n"
     << "#\n";
  for (vector<string>::const_iterator i=history.begin();i!=history.end();++i) {
    ps << *i << '\n';
  }

  ps << "ESTIMATION DESCRIPTION (required)\n"
     << est_blk.pname << "   "
     << "Project name, pname, char*\n"
     << fmt('f', 12, 1, gsmver) <<  "   "
     << "GSM version, defines format of this file, gsmver, float\n" 
     << fmt('d', 12, est_blk.print) <<  "   "
     << "Write detailed output if print=1, int\n"
     << fmt('d', 12, est_blk.sci_draw_from_prior) <<  "   "
     << "Prior draws in sci_mod chain if sci_draw_from_prior=1, int\n"
     << fmt('d', 12, est_blk.run_sci_chain) <<  "   "
     << "Run sci_mod chain if run_sci_chain=1, int\n"
     << fmt('d', 12, est_blk.stat_draw_from_prior) <<  "   "
     << "Prior draws in stat_mod chain if stat_draw_from_prior=1, int\n"
     << fmt('d', 12, est_blk.run_stat_chain) <<  "   " 
     << "Run stat_mod chain (i.e. assess_chain) if run_stat_chain=1, int\n"
     << fmt('f', 12, 1, est_blk.kappa) << "   "
     << "Value of kappa for assess prior, kappa, float\n";

  ps << "DATA DESCRIPTION (required)\n"
     << fmt('d', 12, dat_blk.M) << "   "
     << "Dimension of the data, M, int\n"
     << fmt('d', 12, dat_blk.n) << "   "
     << "Number of observations, n, int\n";
  ps << dat_blk.dsn;
  if (dat_blk.dsn.size() <= 12) {
    for (INTEGER i=dat_blk.dsn.size(); i<15; ++i) ps << ' ';
    ps << "File name, any length, no embedded blanks, dsn, string\n";
  }
  else {
    ps << '\n';
  }
  if (dat_blk.M <= 6) {
    for (INTEGER i=1; i<=dat_blk.M; ++i) ps << dat_blk.fields[i] << ' ';
    for (INTEGER i=2*dat_blk.M; i<15; ++i) ps << ' ';
    ps << "Read these white space separated fields, fields, intvec\n";
  }
  else {
    ps << dat_blk.fieldline << '\n';
  }

  vector<string>::const_iterator itr;

  ps << "STAT_MOD DESCRIPTION (required)\n"
     << fmt('d', 12, stat_blk.len_stat_parm) <<  "   "
     << "Number of parameters, len_stat_parm, int\n"
     << fmt('d', 12, stat_blk.len_stat_func) <<  "   "
     << "Number of functionals, len_stat_func, int\n";

  ps << "STAT_MOD PARMFILE (required) " 
     << "(constructor sees as vector<string> pfvec, alvec)\n";
  ps << stat_blk.stat_parmfile;
  if (stat_blk.stat_parmfile.size() <= 12) {
    for (INTEGER i=stat_blk.stat_parmfile.size(); i<15; ++i) ps << ' ';
    ps << "File name, code __none__ if none, stat_parmfile, string\n";
  }
  else {
    ps << '\n';
  }
  for (itr=stat_alvec.begin(); itr!=stat_alvec.end(); ++itr) {
    ps << *itr << '\n';
  }

  ps << "STAT_MOD MCMC DESCRIPTION (required) (describes assess chains)\n"
     << fmt('d', 12, stat_blk.stat_seed) << "   "
     << "Seed for stat_mod MCMC simulations, stat_seed, int\n"
     << fmt('d', 12, stat_blk.len_stat_chain) << "   "
     << "Number stat_mod MCMC simulations per file, len_stat_chain, int\n"
     << fmt('d', 12, stat_blk.num_stat_files) << "   "
     << "Number of extra MCMC simulation files, num_stat_files, int\n"
     << fmt('f', 12, 1, stat_blk.stat_sclfac) << "   "
     << "Rescale proposal scaling by this value, stat_sclfac, float\n";

  ps << "STAT_MOD PARAMETER START VALUES (required)\n";
  itr = stat_parm_rhs.begin();
  for (INTEGER i=1; i<=stat_blk.len_stat_parm; ++i) {
    ps << fmt('e',26,17,stat_parm[i]) << fmt('d',5, stat_mask[i]) 
       << (itr++)->substr(0,49) << '\n';
  }

  ps << "STAT_MOD PROPOSAL SCALING (required)\n";
  itr = stat_scale_rhs.begin();
  for (INTEGER i=1; i<=stat_blk.len_stat_parm; ++i) {
    ps << fmt('e',26,17,stat_scale[i]) << (itr++)->substr(0,49) << '\n';
  }

  if (stat_groups_txt.size()>0) {
    ps << "STAT_MOD PROPOSAL GROUPING (optional) (frequencies are relative)\n";
    for (itr=stat_groups_txt.begin(); itr!=stat_groups_txt.end(); ++itr) {
      ps << *itr << '\n';
    }
  }

  ps << "SCI_MOD DESCRIPTION (required)\n"
     << fmt('d', 12, sci_blk.len_sci_parm) <<  "   "
     << "Number of parameters, len_sci_parm, int\n"
     << fmt('d', 12, sci_blk.len_sci_func) <<  "   "
     << "Number of functionals, len_sci_func, int\n";

  ps << "SCI_MOD PARMFILE (required) " 
     << "(constructor sees as vector<string> pfvec, alvec)\n";
  ps << sci_blk.sci_parmfile;
  if (sci_blk.sci_parmfile.size() <= 12) {
    for (INTEGER i=sci_blk.sci_parmfile.size(); i<15; ++i) ps << ' ';
    ps << "File name, code __none__ if none, sci_parmfile, string\n";
  }
  else {
    ps << '\n';
  }
  for (itr=sci_alvec.begin(); itr!=sci_alvec.end(); ++itr) {
    ps << *itr << '\n';
  }

  ps << "SCI_MOD MCMC DESCRIPTION (required)\n"
     << fmt('d', 12, sci_blk.sci_seed) << "   "
     << "Seed for sci_mod MCMC simulations, sci_seed, int\n"
     << fmt('d', 12, sci_blk.len_sci_chain) << "   "
     << "Number sci_mod MCMC simulations per file, len_sci_chain, int\n"
     << fmt('d', 12, sci_blk.num_sci_files) << "   "
     << "Number of extra MCMC simulation files, num_sci_files, int\n"
     << fmt('d', 12, sci_blk.analytic_mle) << "   "
     << "Use analytic expression to compute mle if analytic_mle=1, int\n"
     << fmt('d', 12, sci_blk.len_sub_chain) << "   "
     << "Length of sub chain to compute mle, len_sub_chain, int\n"
     << fmt('d', 12, sci_blk.num_sub_chains) << "   "
     << "Number of extra sub chains, num_sub_chains, int\n"
     << fmt('d', 12, sci_blk.num_polish_iter) << "   "
     << "Number of quasi-Newton iterates, num_polish_iter, int\n"
     << fmt('e',12,2,sci_blk.polish_toler) << "   "
     << "Tolerance for quasi-Newton, polish_toler, float\n"
     << fmt('f', 12, 1, sci_blk.sci_sclfac) << "   "
     << "Rescale proposal scaling by this value, sci_sclfac, float\n"
     << fmt('f', 12, 1, sci_blk.sci_incfac) << "   "
     << "Rescale parameter increments by this value, sci_incfac, float\n";

  ps << "SCI_MOD PARAMETER START VALUES (required)\n";
  itr = sci_parm_rhs.begin();
  for (INTEGER i=1; i<=sci_blk.len_sci_parm; ++i) {
    ps << fmt('e',26,17,sci_parm[i]) << fmt('d',5, sci_mask[i]) 
       << (itr++)->substr(0,49) << '\n';
  }

  ps << "SCI_MOD PROPOSAL SCALING (required)\n";
  itr = sci_scale_rhs.begin();
  for (INTEGER i=1; i<=sci_blk.len_sci_parm; ++i) {
    ps << fmt('e',26,17,sci_scale[i]) << (itr++)->substr(0,49) << '\n';
  }

  ps << "SCI_PARAMETER INCREMENTS (required) "
     << "(must be (fractional) powers of two)\n";
  itr = sci_incr_rhs.begin();
  for (INTEGER i=1; i<=sci_blk.len_sci_parm; ++i) {
    ps << fmt('e',26,17,sci_incr[i]) << (itr++)->substr(0,49) << '\n';
  }

  if (sci_groups_txt.size()>0) {
    ps << "SCI_MOD PROPOSAL GROUPING (optional) (frequencies are relative)\n";
    for (itr=sci_groups_txt.begin(); itr!=sci_groups_txt.end(); ++itr) {
      ps << *itr << '\n';
    }
  }

  return true;
}
