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

Copyright (C) 2004, 2005, 2006, 2007.

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 <sstream>
#include "libsnp.h"
#include "libsmm.h"
#include "snpusr.h"
#include "snp.h"
#include "emm.h"

using namespace scl;
using namespace libsnp;
using namespace libsmm;
using namespace snp;
using namespace emm;
using namespace std;


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, emmparms, cannot interpret this line: \"" +s+ "\"");
      warn("Warning, emmparms, 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, emmparms, cannot interpret this line: \"" +s+ "\"");
      warn("Warning, emmparms, " +first+ " 1 substituted for \"" +s+ "\"");
      strtod(first.c_str(),endptr);
      if (**endptr) error("Error, emmparms, 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, emmparms, cannot read line as numeric: "+s);
    strtod(second.c_str(),endptr);
    if (**endptr) error("Error, emmparms, cannot read line as numeric: "+s);
    return pair<string,string>(first,second);
  }

  bool check_estblock(estblock eb) 
  {
    if (eb.pname.size() != 12) return false;
    if (eb.emmver != 2.6) return false;
    if (eb.objtype < 0 || eb.objtype > 3) return false;
    if (eb.proptype < 0 || eb.proptype > 1) return false;
    if (eb.lchain < 0) return false;
    if (eb.nfile < 0) return false;
    if (eb.incfac <= 0) return false;
    const REAL x = 5.0;
    if ((x*eb.incfac)/eb.incfac != x) {
      warn("Warning, check_estparms, incfac not a (fractional) power of two");
    }
    if (eb.sclfac <= 0) return false;
    if (eb.temperature <= 0) return false;
    if (eb.stride < 1) return false;
    if (eb.max_cache_size < 0) 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_modblock(const modblock& mb) 
  {
    if (mb.len_mod_parm <= 0) return false;
    if (mb.len_mod_func <= 0) return false;
    return true;
  }
}

void emm::emmparms::set_rho(const realmat& r)
{  
  if (rho.get_rows() != r.get_rows() || rho.get_cols() != r.get_cols() ) {
    rho = r;
  }
  else {
    error("Error, emmparms, set_rho, bad input");
  }
}

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

  vector<string> pf; 

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

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

bool emm::emmparms::set_parms(const vector<string>& pf, ostream& detail)
{

  history.erase(history.begin(),history.end());
  grouptxt.erase(grouptxt.begin(),grouptxt.end());
  parmrhs.erase(parmrhs.begin(),parmrhs.end());
  scalerhs.erase(scalerhs.begin(),scalerhs.end());
  incrhs.erase(incrhs.begin(),incrhs.end());
  usrparms.erase(usrparms.begin(),usrparms.end());
  mod_alvec.erase(mod_alvec.begin(),mod_alvec.end());
  obj_alvec.erase(obj_alvec.begin(),obj_alvec.end());
  
  typedef vector<string>::const_iterator pf_itr;

  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, emmparms, " + 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()) return false;
  if (kw_ptr + 12 >= pf.end()) error("Error, emmparms, " + header);

  est_blk.pname = (++kw_ptr)->substr(0,12); 
  est_blk.emmver = atof((++kw_ptr)->substr(0,12).c_str());
  est_blk.objtype = atoi((++kw_ptr)->substr(0,12).c_str()); 
  est_blk.proptype = atoi((++kw_ptr)->substr(0,12).c_str()); 
  est_blk.print = atoi((++kw_ptr)->substr(0,12).c_str()); 
  est_blk.seed = atoi((++kw_ptr)->substr(0,12).c_str()); 
  est_blk.lchain = atoi((++kw_ptr)->substr(0,12).c_str()); 
  est_blk.nfile = atoi((++kw_ptr)->substr(0,12).c_str()); 
  if (est_blk.emmver == 2.1) {
    if (est_blk.objtype == 1) error("Error, emmparms, EDF objfun removed");
    est_blk.incfac = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.sclfac = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.temperature = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.kilse = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.stride = 1;
    est_blk.draw_from_prior = false;
    est_blk.max_cache_size = 50000;
    est_blk.emmver = 2.6;
  }
  else if (est_blk.emmver == 2.2) {
    if (kw_ptr + 13 >= pf.end()) error("Error, emmparms, " + header);
    if (est_blk.objtype == 1) error("Error, emmparms, EDF objfun removed");
    est_blk.incfac = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.sclfac = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.temperature = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.kilse = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.stride = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.draw_from_prior = false;
    est_blk.max_cache_size = 50000;
    est_blk.emmver = 2.6;
  }
  else if (est_blk.emmver == 2.3) {
    if (kw_ptr + 13 >= pf.end()) error("Error, emmparms, " + header);
    est_blk.incfac = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.sclfac = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.temperature = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.kilse = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.stride = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.draw_from_prior = false;
    est_blk.max_cache_size = 50000;
    est_blk.emmver = 2.6;
  }
  else if (est_blk.emmver == 2.4) {
    if (kw_ptr + 14 >= pf.end()) error("Error, emmparms, " + header);
    est_blk.incfac = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.sclfac = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.temperature = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.kilse = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.stride = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.draw_from_prior = false;
    est_blk.max_cache_size = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.emmver = 2.6;
  }
  else if (est_blk.emmver == 2.5) {
    if (kw_ptr + 14 >= pf.end()) error("Error, emmparms, " + header);
    est_blk.incfac = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.sclfac = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.temperature = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.kilse = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.stride = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.draw_from_prior = false;
    est_blk.max_cache_size = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.emmver = 2.6;
  }
  else {
    if (kw_ptr + 14 >= pf.end()) error("Error, emmparms, " + header);
    est_blk.sclfac = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.incfac = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.temperature = atof((++kw_ptr)->substr(0,12).c_str());
    est_blk.kilse = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.stride = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.draw_from_prior = atoi((++kw_ptr)->substr(0,12).c_str());
    est_blk.max_cache_size = atoi((++kw_ptr)->substr(0,12).c_str());
  }

  if (est_blk.print) {
    detail << starbox("/Input EMM parmfile//") << '\n';
    for (pf_itr itr=pf.begin(); itr!=pf.end(); ++itr) {
      detail << *itr << '\n';
    }
    detail.flush();
    detail << starbox("/EMM 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 emmver = " << fmt('f',3,1,est_blk.emmver) << '\n';
    detail << "\t objtype = " << est_blk.objtype << '\n';
    detail << "\t proptype = " << est_blk.proptype << '\n';
    detail << "\t print = " << est_blk.print << '\n';
    detail << "\t seed = " << est_blk.seed << '\n';
    detail << "\t lchain = " << est_blk.lchain << '\n';
    detail << "\t nfile = " << est_blk.nfile << '\n';
    detail << "\t sclfac = " << est_blk.sclfac << '\n';
    detail << "\t incfac = " << est_blk.incfac 
           << " (should have exact binary representation)\n";
    detail << "\t temperature = " << est_blk.temperature << '\n';
    detail << "\t kilse = " << est_blk.kilse << '\n';
    detail << "\t stride = " << est_blk.stride << '\n';
    detail << "\t draw_from_prior = " << est_blk.draw_from_prior << '\n';
    detail << "\t max_cache_size = " << est_blk.max_cache_size << '\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, emmparms, " + 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, emmparms, 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("MODEL 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, emmparms, " + header);

  mod_blk.len_mod_parm = atoi((++kw_ptr)->substr(0,12).c_str());
  mod_blk.len_mod_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_mod_parm = " << mod_blk.len_mod_parm << '\n';
    detail << "\t len_mod_func = " << mod_blk.len_mod_func << '\n';
    detail.flush();
  }

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

  header=kw.set_keyword("MODEL 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, emmparms, " + header);
  string line = *(++kw_ptr);
  mod_blk.mod_parmfile = cutstr(line);
  pf_itr start = ++kw_ptr;  // No idea why djgpp compiler insists on these
  pf_itr stop = pf.end();   // explicit definitions for find_if below.
  if (mod_blk.mod_parmfile != string("__none__")) {
    mod_blk.is_mod_parmfile = true;
  }
  else {
    mod_blk.is_mod_parmfile = false;
  }
  keyword endkw("#end additional lines");
  pf_itr end_ptr = find_if(start, stop, endkw); 
  if (end_ptr==stop || end_ptr<start) {
    error("Error, emmparms, " + header);
  }
  ++end_ptr;
  for (kw_ptr=start; kw_ptr != end_ptr; ++kw_ptr) {
    mod_alvec.push_back(*kw_ptr); 
  }

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

  rho.resize(mod_blk.len_mod_parm,1);
  rho_mask.resize(rho.size());
  
  INTEGER nfixed = 0;

  header=kw.set_keyword("PARAMETER START VALUES");   //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr + rho.size() >= pf.end()) error("Error, emmparms, " + header);
  for (INTEGER i=1; i<=rho.size(); ++i) {
    pair<string,string> ss = find_ss((++kw_ptr)->c_str());
    rho[i] = atof(ss.first.c_str());
    rho_mask[i] = atoi(ss.second.c_str());
    if (rho_mask[i] == 0) ++nfixed;
    string line = *kw_ptr;
    cutstr(line); cutstr(line);
    parmrhs.push_back(line);
  }
  
  INTEGER nactive = rho.size() - nfixed;

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

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    detail << "\n\t rho = ";
    detail << rho;
    detail << "\n\t rho_mask = ";
    detail << rho_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.flush();
  }

  rho_increment.resize(rho.size(),1);
  
  header=kw.set_keyword("PARAMETER INCREMENTS");   //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
    est_blk.is_inc_block = true;
    if (est_blk.max_cache_size <= 100) {
      warn("Warning, max_cache_size must be > 100 if PARAMETER INCREMENTS set");
      warn("max_cache_size reset to 50000");
      est_blk.max_cache_size = 50000;
    }
    if (kw_ptr + rho.size() >= pf.end()) error("Error, emmparms, " + header);
    for (INTEGER i=1; i<=rho.size(); ++i) {
      rho_increment[i] = atof((++kw_ptr)->c_str());
      rho_increment[i] *= est_blk.incfac;
      if (rho_increment[i]<=0) error("Error, emmparms, " + header);
      string line = *kw_ptr;
      cutstr(line);
      incrhs.push_back(line);
    }
  }
  else {
    est_blk.is_inc_block = false;
    for (INTEGER i=1; i<=rho.size(); ++i) {
      rho_increment[i] = 9.53674316406250000e-07;
    }
    warn("Warning, no PARAMETER INCREMENTS blk, parameters not put on a grid");
    warn("Warning, no PARAMETER INCREMENTS blk, cache is ineffective");
  }

  if (est_blk.print) {
    if (est_blk.is_inc_block) {
      detail << '\n' << header << '\n';
      detail << "(after multiplication by incfac)\n";
      detail << "(output parmfile this increment and incfac=1.0)\n";
      detail << "(should have exact binary representation)\n";
      detail << "\n\t rho_increment = ";
      detail << rho_increment;
      detail.flush();
    }
    else {
      detail << '\n' << header << " not found." << '\n';
      detail << "Parameters are not put on a grid." << '\n';
      detail << "Cache is ineffective." << '\n';
      detail.flush();
    }
  }

  prop_scale.resize(rho.size(),1);
  
  header=kw.set_keyword("PROPOSAL SCALING");   //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr + rho.size() >= pf.end()) error("Error, emmparms, " + header);
  for (INTEGER i=1; i<=rho.size(); ++i) {
    prop_scale[i] = atof((++kw_ptr)->c_str());
    prop_scale[i] *= est_blk.sclfac;
    if (prop_scale[i]<=0) error("Error, emmparms, " + header);
    string line = *kw_ptr;
    cutstr(line);
    scalerhs.push_back(line);
  }

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    detail << "(after multiplication by sclfac)\n";
    detail << "(output parmfile this scale and sclfac=1.0)\n";
    detail << "\n\t prop_scale = ";
    detail << prop_scale;
    detail.flush();
  }

  header=kw.set_keyword("OBJFUN 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, emmparms, " + header);
  line = *(++kw_ptr);
  obj_blk.obj_parmfile = cutstr(line);
  start = ++kw_ptr;  // No idea why djgpp compiler insists on these
  stop = pf.end();   // explicit definitions for find_if below.
  if (obj_blk.obj_parmfile != string("__none__")) {
    obj_blk.is_obj_parmfile = true;
  }
  else {
    obj_blk.is_obj_parmfile = false;
  }
  endkw.set_keyword("#end additional lines");
  end_ptr = find_if(start, stop, endkw); 
  if (end_ptr==stop || end_ptr<start) {
    error("Error, emmparms, " + header);
  }
  ++end_ptr;
  for (kw_ptr=start; kw_ptr != end_ptr; ++kw_ptr) {
    obj_alvec.push_back(*kw_ptr); 
  }

  if (est_blk.print) {
    detail << '\n' << header << '\n';
    #if defined COMPILER_HAS_BOOLALPHA
      detail << boolalpha;
    #endif
    detail << "\t is_obj_parmfile = " << obj_blk.is_obj_parmfile << '\n';
    if (obj_blk.is_obj_parmfile) {
      detail << "\t obj_parmfile = " << obj_blk.obj_parmfile << '\n';
    }
    vector<string>::const_iterator itr;
    for (itr=obj_alvec.begin(); itr!=obj_alvec.end(); ++itr ) {
      detail << *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] = rho_increment[fixed[i]];
      ufixed[i] = rho[fixed[i]];
      Vfixed(i,i) = pow(prop_scale[fixed[i]],2);
    }
    prop_groups.push_back(prop_group(0.0,fixed,ifixed,ufixed,Vfixed));
  }

  header=kw.set_keyword("PROPOSAL GROUPING");  //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
    INTEGER count = 0;
    while (count < rho.size()) {
      INTEGER top = rho.size();
      INTEGER r=1;
      realmat gmat(rho.size()+1,rho.size()+1,0.0);
      while (r <= top) {
        ++kw_ptr;
        grouptxt.push_back(*kw_ptr);
        if (kw_ptr == pf.end()) error("Error, emmparms, " + 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 <= rho.size()) {
          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, emmparms, " + 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] = rho_increment[gvec[i]];
          mean[i] = rho[gvec[i]];
        }
        for (INTEGER j=1; j<=gvec.size(); ++j) {
          for (INTEGER i=1; i<=gvec.size(); ++i) {
            Vmat(i,j) = Vmat(i,j)*prop_scale[gvec[i]]*prop_scale[gvec[j]];
          }
        }
        prop_groups.push_back(prop_group(freq, gvec, ginc, mean, Vmat));
      }
    }
  }
  else {
    for (INTEGER i=1; i<=nactive; ++i) {
      intvec gvec(1,active[i]);
      REAL var = pow(prop_scale[active[i]],2);
      REAL mu = rho[active[i]];
      REAL inc = rho_increment[active[i]];
      prop_groups.push_back(prop_group(1.0, gvec, realmat(1,1,inc), 
        realmat(1,1,mu), realmat(1,1,var)));
    }
  }
  if (est_blk.print) {
    detail << starbox("/Prop groups//") << '\n';
    detail << "\t\t Constructed from "<<header<<", if present; \n";
    detail << "\t\t If not, move one at a time grouping is used.\n";
    detail << "\t\t If rho has fixed elements, they are in Group 0\n";
    detail << "\t\t and not changed during the MCMC iterations.\n";
    detail << "\t\t "<<header<<" correlations have been converted\n";
    detail << "\t\t to variances using the prop_scale vector above.\n";
    detail << "\t\t Frequencies are relative, not absolute.\n\n";
    prop_def::const_iterator itr;
    INTEGER i=0;
    for (itr=prop_groups.begin(); itr!=prop_groups.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';
    }
    detail.flush();
  }

  return true;
}

bool emm::emmparms::write_parms (string parmfile, string prefix,
                 INT_32BIT seed, const realmat& rho_last, 
		 const realmat& rho_mode, const realmat& sig) const
{
  stringstream ps;
  string filename;
  ofstream fout;

  ps << "PARMFILE HISTORY (optional)\n"
     << "#\n"
     << "# This parmfile was written by EMM Version " 
     << fmt('f',3,1,est_blk.emmver) << ' '
     << "using the following line from\n" 
     << "# control.dat, which was read as " 
     << "char*, char*\n"
     << "# -----------------------------------------"
     << "----------------------------------\n";
//ps << "#      input_parmfile       output_prefix\n"   // Change kw_ptr
//   << "#      --------------       -------------\n";  // line if altered
  ps << '#';
  ps << ' ';
  for (INTEGER i=1; i<=INTEGER(19-parmfile.size()); ++i) ps << ' ';
  ps << parmfile;
  ps << ' ';
  for (INTEGER i=1; i<=INTEGER(19-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, est_blk.emmver) <<  "   "
     << "EMM version, defines format of this file, emmver, float\n"
     << fmt('d', 12, est_blk.objtype) <<  "   "
     << "Objfun type, 0 EMM, 1 MLE, 2 usr, objtype, int\n";

  string reghead = ps.str();
  string althead = ps.str();
  string endhead = ps.str();
  ps.str("");

  ps << fmt('d', 12, 1) <<  "   "
     << "Proposal type, 0 group_move, 1 cond_move, 2 usr, proptype, int\n";

  althead += ps.str();
  ps.str("");
    
  ps << fmt('d', 12, 0) <<  "   "
     << "Proposal type, 0 group_move, 1 cond_move, 2 usr, proptype, int\n";

  reghead += ps.str();
  endhead += ps.str();
  ps.str("");

  ps << fmt('d', 12, est_blk.print) <<  "   "
     << "Write detailed output if print=1, int\n";

  reghead += ps.str();
  althead += ps.str();
  endhead += ps.str();
  ps.str("");

  ps << fmt('d', 12, est_blk.seed) <<  "   "
     << "Seed for MCMC simulations, iseed, int\n";

  reghead += ps.str();
  althead += ps.str();
  ps.str("");

  ps << fmt('d', 12, seed) <<  "   "
     << "Seed for MCMC simulations, iseed, int\n";

  endhead += ps.str();
  ps.str("");

  ps << fmt('d', 12, est_blk.lchain) <<  "   " 
     << "Number of MCMC simulations per file, lchain, int\n"
     << fmt('d', 12, est_blk.nfile) <<  "   "
     << "Number of MCMC simulation files beyond the first, nfile, int\n"
     << fmt('f', 12, 1, 1.0) << "   "
     << "Rescale proposal scaling by this value, sclfac, float\n"
     << fmt('f', 12, 1, 1.0) << "   "
     << "Rescale parameter increments by this value, incfac, float\n";
  if (est_blk.temperature == 1.0) {
    ps << fmt('f', 12, 1, est_blk.temperature) << "   ";
  } 
  else {
    ps << fmt('g', 12, 5, est_blk.temperature) << "   ";
  }
  ps << "Rescale objfun by this value, temperature, float\n"
     << fmt('d', 12, est_blk.kilse) << "   "
     << "Sandwich variance not computed if kilse=1, int\n"
     << fmt('d', 12, est_blk.stride) << "   "
     << "The stride used to write MCMC simulations, stride, int\n"
     << fmt('d', 12, est_blk.draw_from_prior) << "   "
     << "Draw from prior if draw_from_prior=1, int\n"
     << fmt('d', 12, est_blk.max_cache_size) << "   "
     << "Max cache size, max_cache_size, int\n";

  ps << "DATA DESCRIPTION (required) "
     << "(mod and obj constructors see realmat data(M,n))\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 << "MODEL DESCRIPTION (required)\n"
     << fmt('d', 12, mod_blk.len_mod_parm) <<  "   "
     << "Number of modal parameters, len_mod_parm, int\n"
     << fmt('d', 12, mod_blk.len_mod_func) <<  "   "
     << "Number of model functionals, len_mod_func, int\n";

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

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

  reghead += ps.str();
  althead += ps.str();
  endhead += ps.str();
  ps.str("");
  
  INTEGER lrho = mod_blk.len_mod_parm;

  ps << "PARAMETER START VALUES (required)\n";
  itr = parmrhs.begin();
  for (INTEGER i=1; i<=lrho; ++i) {
    ps << fmt('e',26,17,rho_mode[i]) << fmt('d',5, rho_mask[i]) 
       << (itr++)->substr(0,49) << '\n';
  }
  
  string mode = ps.str();
  ps.str("");

  ps << "PARAMETER START VALUES (required)\n";
  itr = parmrhs.begin();
  for (INTEGER i=1; i<=lrho; ++i) {
    ps << fmt('e',26,17,rho_last[i]) << fmt('d',5, rho_mask[i]) 
       << (itr++)->substr(0,49) << '\n';
  }

  string last = ps.str();
  ps.str("");

  if (lrho <=20) {
    realmat corr(lrho,lrho,0.0);
    realmat rho_scale(lrho,1);
    for (INTEGER i=1; i<=lrho; ++i) {
      rho_scale[i] = sqrt(sig(i,i));
      for (INTEGER j=1; j<=lrho; ++j) {
        if (sig(i,i) > 0 && sig(j,j) > 0) {
          corr(i,j) = sig(i,j)/sqrt(sig(i,i)*sig(j,j));
        }
      }
    }

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

    if (est_blk.is_inc_block) {
      ps << "PARAMETER INCREMENTS (optional) " 
         << "(fractional powers of two recommended)\n";
      itr = incrhs.begin();
      for (INTEGER i=1; i<=lrho; ++i) {
        ps << fmt('e',26,17,rho_increment[i]) << (itr++)->substr(0,49) << '\n';
      }
    }

    ps << "PROPOSAL GROUPING (optional) (frequencies are relative)\n";
    ps << " 1.0 ";
    for (INTEGER j=1; j<=lrho; ++j) ps << j << " ";
    ps << '\n';
    for (INTEGER i=1; i<=lrho; ++i) {
      ps << "  " << i << ' ';
      for (INTEGER j=1; j<=lrho; ++j) {
        ps << fmt('e',26,17,corr(i,j)) << ' ';  
      }
      ps << '\n'; 
    }

    string alttail = ps.str();
    ps.str("");

    filename = prefix + ".parmfile.alt";

    fout.open(filename.c_str());
    if (!fout) return false; 

    fout << althead << mode << alttail;

    fout.clear();
    fout.close();
  }
  
  ps << "PROPOSAL SCALING (required)" << '\n';
  itr = scalerhs.begin();
  for (INTEGER i=1; i<=lrho; ++i) {
    ps << fmt('e',26,17,prop_scale[i]) << (itr++)->substr(0,49) << '\n';
  }

  if (est_blk.is_inc_block) {
    ps << "PARAMETER INCREMENTS (optional) " 
       << "(fractional powers of two recommended)\n";
    itr = incrhs.begin();
    for (INTEGER i=1; i<=lrho; ++i) {
      ps << fmt('e',26,17,rho_increment[i]) << (itr++)->substr(0,49) << '\n';
    }
  }

  if (grouptxt.size() > 0) {
    ps << "PROPOSAL GROUPING (optional) (frequencies are relative)\n";
    vector<string>::const_iterator grpitr;
    for (grpitr=grouptxt.begin(); grpitr!=grouptxt.end(); ++grpitr){ 
      ps << *grpitr << '\n';
    }
  
  }

  string regtail = ps.str();
  ps.str("");

  filename = prefix + ".parmfile.fit";

  fout.open(filename.c_str());
  if (!fout) return false; 

  fout << reghead << mode << regtail;

  fout.clear();
  fout.close();

  filename = prefix + ".parmfile.end";

  fout.open(filename.c_str());
  if (!fout) return false; 

  fout << endhead << last << regtail;

  fout.clear();
  fout.close();

  return true;
}