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

Copyright (C) 2004, 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 "snp.h"

using namespace scl;
using namespace libsnp;
using namespace snp;
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, parmfile, cannot interpret this line: \"" +s+ "\"");
      warn("Warning, parmfile, 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, parmfile, cannot interpret this line: \"" +s+ "\"");
      warn("Warning, parmfile, " +first+ " 1 substituted for \"" +s+ "\"");
      strtod(first.c_str(),endptr);
      if (**endptr) error("Error, parmfile, 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, parmfile, cannot read line as numeric: "+s);
    strtod(second.c_str(),endptr);
    if (**endptr) error("Error, parmfile, cannot read line as numeric: "+s);
    return pair<string,string>(first,second);
  }
    
  bool check_optparms(optparms op) 
  {
    if (op.pname.size() != 10) return false;
    if (op.snpver != 9.0) return false;
    if (op.toler <= 0) return false;
    if (op.task < 0 || op.task > 6) return false;
    if (op.extra < 0) return false;
    if (op.sfac <= 0) return false;
    if (op.ngrid <= 0) return false;
    return true;
  }
  
  bool check_datparms(datparms dp)
  {
    if (dp.M <= 0) return false;
    if (dp.n <= 0) return false;
    if (dp.drop <= 3 || dp.drop >= dp.n) return false;
    if (dp.cond > dp.n) return false;
    return true;
  }
  
  bool check_tranparms(tranparms tp) //needs work
  {
    if (tp.inflec <= 0) return false;
    realmat null;
    if (tp.mean == null) return false;
    if (tp.variance == null) return false;
    return true;
  }
  
}

snp::parmfile::parmfile()
: fn(snpden()), 
  afn(afunc()), ufn(ufunc()), rfn(rfunc()), 
  afn_mask(afunc()), ufn_mask(ufunc()), rfn_mask(rfunc()), 
  rv(false)
{ }

void snp::parmfile::set_optparms(const optparms& op)
{  
  if (check_optparms(op)) {
    opm = op;
  }
  else {
    error("Error, parmfile, set_optparms, bad input");
  }
}

void snp::parmfile::set_datparms(const datparms& dp)
{  
  if (check_datparms(dp)) {
    dpm = dp;
  }
  else {
    error("Error, parmfile, set_datparms, bad input");
  }
}

void snp::parmfile::set_tranparms(const tranparms& tp)
{  
  if (check_tranparms(tp)) {
    tpm = tp;
  }
  else {
    error("Error, parmfile, set_tranparms, bad input");
  }
}

void snp::parmfile::set_afunc(const libsnp::afunc& af)
{ 
  afn = af;
}

void snp::parmfile::set_ufunc(const libsnp::ufunc& uf)
{ 
  ufn = uf;
}

void snp::parmfile::set_rfunc(const libsnp::rfunc& rf)
{ 
  rfn = rf;
}

bool snp::parmfile::read_parms(const char* filename, ostream& detail) 
{
  ifstream pfstream(filename);
  if (!pfstream) return rv = false;

  vector<string> pf; 

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

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

bool snp::parmfile::set_parms(const vector<string>& pf, std::ostream& detail)
{
  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()) {
    if (kw_ptr + 12 >= pf.end()) error("Error, parmfile, " + header);
    kw_ptr += 12;
    while (kw_ptr != pf.end() && (*kw_ptr)[0] == '#') {
      history.push_back(*kw_ptr); 
      ++kw_ptr;
    }
  }
  
  header=kw.set_keyword("OPTIMIZATION DESCRIPTION");  //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) return rv = false;
  if (kw_ptr + 12 >= pf.end()) error("Error, parmfile, " + header);

  opm.pname = (++kw_ptr)->substr(0,10); 
  opm.snpver = atof((++kw_ptr)->substr(0,10).c_str());
  opm.itmax0 = atoi((++kw_ptr)->substr(0,10).c_str()); 
  opm.itmax1 = atoi((++kw_ptr)->substr(0,10).c_str()); 
  opm.toler = atof((++kw_ptr)->substr(0,10).c_str());
  opm.print = atoi((++kw_ptr)->substr(0,10).c_str()); 
  opm.task = atoi((++kw_ptr)->substr(0,10).c_str()); 
  opm.extra = atoi((++kw_ptr)->substr(0,10).c_str()); 
  opm.sfac = atof((++kw_ptr)->substr(0,10).c_str());
  opm.iseed = atoi((++kw_ptr)->substr(0,10).c_str());
  opm.ngrid = atoi((++kw_ptr)->substr(0,10).c_str());
  opm.kilse = atoi((++kw_ptr)->substr(0,10).c_str());

  if (opm.print) {
    detail << starbox("/Input parmfile//") << '\n';
    for (pf_itr itr=pf.begin(); itr!=pf.end(); ++itr) {
      detail << *itr << '\n';
    }
    detail.flush();
    detail << starbox("/Parmfile interpretation//") << '\n';
    detail << "     In the code, some notational conventions differ from "
           << "the documentation.\n     These differences are noted in the "
           << "output; e.g., M = ly or Lg = q.\n";
  }

  if (opm.print) {
    detail << '\n' << header << '\n';
    #if defined COMPILER_HAS_BOOLALPHA
      detail << boolalpha;
    #endif
    detail << "\t pname = " << opm.pname << '\n';
    detail << "\t snpver = " << fmt('f',3,1,opm.snpver) << '\n';
    detail << "\t itmax0 = " << opm.itmax0 << '\n';
    detail << "\t itmax1 = " << opm.itmax1 << '\n';
    detail << "\t toler = " << opm.toler << '\n';
    detail << "\t print = " << opm.print << '\n';
    detail << "\t task = " << opm.task << '\n';
    detail << "\t extra = " << opm.extra << '\n';
    detail << "\t sfac = " << opm.sfac << '\n';
    detail << "\t iseed = " << opm.iseed << '\n';
    detail << "\t ngrid = " << opm.ngrid << '\n';
    detail << "\t kilse = " << opm.kilse << '\n';
    detail.flush();
  }

  if (!check_optparms(opm)) {
     warn(string("Warning, bad or too few items in ")+header);
     return rv = false;
  }

  header=kw.set_keyword("DATA DESCRIPTION");  //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) return rv = false;
  if (kw_ptr + 7 >= pf.end()) error("Error, parmfile, " + header);

  dpm.M = atoi((++kw_ptr)->substr(0,10).c_str());
  dpm.n = atoi((++kw_ptr)->substr(0,10).c_str());
  dpm.drop = atoi((++kw_ptr)->substr(0,10).c_str());
  dpm.cond = atoi((++kw_ptr)->substr(0,10).c_str());
  dpm.reread = atoi((++kw_ptr)->substr(0,10).c_str());
  string ss = *(++kw_ptr);
  dpm.dsn = cutstr(ss);
  dpm.fields.resize(dpm.M,0);
  ++kw_ptr;
  INTEGER j=0;
  for (INTEGER i=1; i<=dpm.M; ++i) {
    while (isspace((*kw_ptr)[j])) ++j;
    INTEGER k=0;
    while (isdigit((*kw_ptr)[j+k])) ++k;
    dpm.fields[i] = atoi(kw_ptr->substr(j,k).c_str());
    j+=k;
  }
  if (dpm.fields[dpm.M]==0) {
    error("Error, parmfile, not enough fields this line:\n" + *kw_ptr);
  }
    
  if (dpm.cond < 0) dpm.cond = 0;
  if (dpm.cond > dpm.n) dpm.cond = dpm.n;

  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "\t M = " << dpm.M << '\n';
    detail << "\t n = " << dpm.n << '\n';
    detail << "\t drop = " << dpm.drop << '\n';
    detail << "\t cond = " << dpm.cond << '\n';
    detail << "\t reread = " << dpm.reread << '\n';
    detail << "\t dsn = " << dpm.dsn << '\n';
    detail << "\t fields: "; 
    for (INTEGER i=1; i<=dpm.fields.size(); ++i) detail <<" "<< dpm.fields[i];
    detail << '\n';
    detail.flush();
  }

  if (!check_datparms(dpm)) {
    warn(string("Warning, bad or too few items in ")+header);
    return rv = false;
  }

  header=kw.set_keyword("TRANSFORM DESCRIPTION");  //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) return rv = false;
  if (kw_ptr + 4 >= pf.end()) error("Error, parmfile, " + header);

  tpm.useold = atoi((++kw_ptr)->substr(0,10).c_str());
  tpm.diag = atoi((++kw_ptr)->substr(0,10).c_str());
  tpm.squash = atoi((++kw_ptr)->substr(0,10).c_str());
  tpm.inflec = atof((++kw_ptr)->substr(0,10).c_str());

  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "\t useold = " << tpm.useold << '\n';
    detail << "\t diag = " << tpm.diag << '\n';
    detail << "\t squash = " << tpm.squash << '\n';
    detail << "\t inflec = " << tpm.inflec << '\n';
    detail.flush();
  }

  header=kw.set_keyword("TRANSFORM START VALUES FOR mean");  //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) {
    if (tpm.useold) {
      warn("Warning, parmfile, no starts found, mean put to zero"); 
    }
    tpm.mean.resize(dpm.M,1,0.0);
  }
  else {
    tpm.mean.resize(dpm.M,1); 
    if (kw_ptr+tpm.mean.size()>=pf.end()) error("Error, parmfile, " + header);
    for (INTEGER i=1; i<=tpm.mean.size(); ++i) {
      tpm.mean[i]= atof((++kw_ptr)->c_str());
    }
    if (opm.print) {
      detail << '\n' << header << '\n';
      detail << "(Values are from constructed object)";
      detail << tpm.mean << '\n';
      detail.flush();
    }
  }

  header=kw.set_keyword("TRANSFORM START VALUES FOR variance");  //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) {
    if (tpm.useold) {
      warn("Warning, parmfile, no starts found, variance put to identity"); 
    }
    tpm.variance.resize(dpm.M,dpm.M,0.0);
    for (INTEGER i=1; i<=dpm.M; ++i) tpm.variance(i,i) = 1.0;
  }
  else {
    tpm.variance.resize(dpm.M,dpm.M); 
    if (kw_ptr+tpm.variance.size()>=pf.end()) error("Error, parmfile, "+header);
    for (INTEGER j=1; j<=dpm.M; ++j) {
      for (INTEGER i=1; i<=dpm.M; ++i) {
        tpm.variance(i,j) = atof((++kw_ptr)->c_str());
        if (tpm.diag && i!=j) tpm.variance(i,j) = 0.0;
      }
    }
  
    if (opm.print) {
      detail << '\n' << header << '\n';
      detail << "(Values are from constructed object)";
      detail << tpm.variance << '\n';
      detail.flush();
    }
  }

  if (!check_tranparms(tpm)) {
    warn(string("Warning, bad or too few items in ")+header);
    return rv = false;
  }
  header=kw.set_keyword("POLYNOMIAL START VALUE DESCRIPTION");  //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) return rv = false;
  if (kw_ptr + 8 >= pf.end()) error("Error, parmfile, " + header);

  INTEGER Kz = atoi((++kw_ptr)->substr(0,10).c_str());
  INTEGER Iz = atoi((++kw_ptr)->substr(0,10).c_str());
  REAL eps0 = atof((++kw_ptr)->substr(0,10).c_str());
  INTEGER lagP = atoi((++kw_ptr)->substr(0,10).c_str());
  INTEGER maxKz = atoi((++kw_ptr)->substr(0,10).c_str());
  INTEGER maxIz = atoi((++kw_ptr)->substr(0,10).c_str());
  INTEGER Kx = atoi((++kw_ptr)->substr(0,10).c_str());
  INTEGER Ix = atoi((++kw_ptr)->substr(0,10).c_str());

  if (lagP <= 0 ) {
    lagP = 1;
    warn("Warning, parmfile, lagP reset to 1");
  }

  snpden f_start(Kz, Iz, dpm.M);
  f_start.set_e0(eps0);
  vector<intvec> alpha = f_start.get_alpha();

  afunc af_start(alpha, maxKz, maxIz, Kx, Ix, dpm.M, lagP); 
  afunc af_mask_start(alpha, maxKz, maxIz, Kx, Ix, dpm.M, lagP); 
  vector<intvec> beta = af_start.get_beta();


  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "(Values are from constructed object)\n";
    detail << "\t Kz = " << f_start.get_Kz() << '\n';
    detail << "\t Iz = " << f_start.get_Iz() << '\n';
    detail << "\t e0 = " << f_start.get_e0() << " = eps0\n";
    detail << "\t lag = " << af_start.get_lag() << " = lagP\n";
    detail << "\t maxKz = " << af_start.get_maxKz() << '\n';
    detail << "\t maxIz = " << af_start.get_maxIz() << '\n';
    detail << "\t Kx = " << af_start.get_Kx() << '\n';
    detail << "\t Ix = " << af_start.get_Ix() << '\n';
    detail << '\n';
    detail << "\t ly = " << f_start.get_ly() << " = M\n";
    detail << "\t lx = " << af_start.get_lx() << " = M\n";
    detail << "\t la = " << f_start.get_la() << '\n';
    detail << "\t lR = " << f_start.get_lR() << '\n';
    detail << "\t nrowa0 = " << af_start.get_nrowa0() << '\n';
    detail << "\t nrowA = " << af_start.get_nrowA() << '\n';
    detail << "\t ncolA = " << af_start.get_ncolA() << '\n';
    detail << '\n';
    detail << "\t alpha:\n";
    for (INTEGER i=0; i<f_start.get_la(); i++) {
      detail << "\t\t   Row " << fmt('d',3,i+1) << "   Lambda ";
      for (INTEGER j=1; j<=f_start.get_ly(); ++j) {
        detail << alpha[i][j];
      }
      detail << '\n';
    }
    detail << '\n';
    detail << "\t Note: If you make modifications, be aware that alpha is\n";
    detail << "\t       indexed from zero in the code.\n\n";
    detail << "\t odx: " << af_start.get_odx() << '\n';
    detail << "\t Note: odx indicates rows of alpha that a0 determines.\n";
    detail << '\n';
    detail << "\t idx: " << af_start.get_idx() << '\n';
    detail << "\t Note: idx indicates rows of alpha that A determines.\n";
    detail << '\n';
    detail << "\t beta:\n";
    for (INTEGER i=0; i<af_start.get_ncolA(); i++) {
      detail << "\t\t   Row " << fmt('d',3,i+1) << "   lambda ";
      for (INTEGER j=1; j<=beta[i].size(); ++j) {
        detail << beta[i][j];
      }
      detail << '\n';
    }
    detail << '\n';
    detail << "\t Note: If you make modifications, be aware that beta is\n";
    detail << "\t       indexed from zero in the code.\n\n";
    detail.flush();
  }

  realmat mask_a0 = af_start.get_a0();  fill(mask_a0);

  header=kw.set_keyword("POLYNOMIAL START VALUES FOR a0");  //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
  
    realmat a0 = af_start.get_a0();
    if (kw_ptr+a0.size()>=pf.end()) error("Error, parmfile, " + header);
    for (INTEGER i=1; i<=a0.size(); ++i) {
      pair<string,string> ss = find_ss((++kw_ptr)->c_str());
      a0[i] = atof(ss.first.c_str());
      mask_a0[i] = atof(ss.second.c_str()) - 1.0;   //Coded this way because
    }                                               //newly added params are
    af_start.set_a0(a0);                            //set to 0 by constructor
  }
  
  if (mask_a0.size()>0) af_mask_start.set_a0(mask_a0);

  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "(Values are from constructed object)";
    detail << af_start.get_a0();
    detail << "\n(Mask, coded 0 for fixed, 1 for active)";
    detail << "\n(But, in the code sometimes coded -1 for fixed, 0 for active)";
    if (mask_a0.size()>0) ++mask_a0;      //External representaion
    detail << mask_a0 << '\n';  
    if (mask_a0.size()>0) --mask_a0;      //Restore internal representaion
    detail.flush();
  }

  realmat mask_A = af_start.get_A();  fill(mask_A);

  header=kw.set_keyword("POLYNOMIAL START VALUES FOR A");  //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
  
    realmat A = af_start.get_A();
    if (kw_ptr+A.size()>=pf.end()) error("Error, parmfile, " + header);
    for (INTEGER i=1; i<=A.size(); ++i) {
      pair<string,string> ss = find_ss((++kw_ptr)->c_str());
      A[i] = atof(ss.first.c_str());
      mask_A[i] = atof(ss.second.c_str()) - 1.0;
    }
    if (A[1] != 1.0) {
      A[1] = 1.0;
      warn("Warning, parmfile, A(1,1) reset to 1.0");
    }
    af_start.set_A(A);
  }
  
  mask_A[1] = -1.0;
  af_mask_start.set_A(mask_A);

  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "(Values are from constructed object)";
    detail << af_start.get_A();
    detail << "\n(Mask, coded 0 for fixed, 1 for active)";
    detail << "\n(But, in the code sometimes coded -1 for fixed, 0 for active)";
    detail << ++mask_A << '\n';  --mask_A;
    detail.flush();
  }

  header=kw.set_keyword("POLYNOMIAL DESCRIPTION");  //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
    if (kw_ptr + 8 >= pf.end()) error("Error, parmfile, " + header);
  
    Kz += atoi((++kw_ptr)->substr(0,10).c_str());
    Iz += atoi((++kw_ptr)->substr(0,10).c_str());
    eps0 += atof((++kw_ptr)->substr(0,10).c_str());
    lagP += atoi((++kw_ptr)->substr(0,10).c_str());
    maxKz += atoi((++kw_ptr)->substr(0,10).c_str());
    maxIz += atoi((++kw_ptr)->substr(0,10).c_str());
    Kx += atoi((++kw_ptr)->substr(0,10).c_str());
    Ix += atoi((++kw_ptr)->substr(0,10).c_str());
  
    if (lagP <= 0 ) {
      lagP = 1;
      warn("Warning, parmfile, lagP reset to 1");
    }
  } 

  snpden f(f_start, Kz, Iz, dpm.M);
  f.set_e0(eps0);
  alpha = f.get_alpha();
  
  afunc af(af_start, alpha, maxKz, maxIz, Kx, Ix, dpm.M, lagP); 
  afunc af_mask(af_mask_start, alpha, maxKz, maxIz, Kx, Ix, dpm.M, lagP); 
  beta = af.get_beta();

  fn = f;
  afn = af;
  afn_mask = af_mask;
  
  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "(Values are after application of increment or decrement)\n";
    detail << "(Values are from reconstructed object)\n";
    detail << "\t Kz = " << f.get_Kz() << '\n';
    detail << "\t Iz = " << f.get_Iz() << '\n';
    detail << "\t e0 = " << f.get_e0() << " = eps0\n";
    detail << "\t lag = " << af.get_lag() << " = lagP\n";
    detail << "\t maxKz = " << af.get_maxKz() << '\n';
    detail << "\t maxIz = " << af.get_maxIz() << '\n';
    detail << "\t Kx = " << af.get_Kx() << '\n';
    detail << "\t Ix = " << af.get_Ix() << '\n';
    detail << '\n';
    detail << "\t ly = " << f.get_ly() << " = M\n";
    detail << "\t lx = " << af.get_lx() << " = M\n";
    detail << "\t la = " << f.get_la() << '\n';
    detail << "\t lR = " << f.get_lR() << '\n';
    detail << "\t nrowa0 = " << af.get_nrowa0() << '\n';
    detail << "\t nrowA = " << af.get_nrowA() << '\n';
    detail << "\t ncolA = " << af.get_ncolA() << '\n';
    detail << '\n';
    detail << "\t alpha:\n";
    for (INTEGER i=0; i<f.get_la(); i++) {
      detail << "\t\t   Row " << fmt('d',3,i+1) << "   Lambda ";
      for (INTEGER j=1; j<=f.get_ly(); ++j) {
        detail << alpha[i][j];
      }
      detail << '\n';
    }
    detail << '\n';
    detail << "\t Note: If you make modifications, be aware that alpha is\n";
    detail << "\t       indexed from zero in the code.\n\n";
    detail << "\t odx: " << af.get_odx() << '\n';
    detail << "\t Note: odx indicates rows of alpha that a0 determines.\n";
    detail << '\n';
    detail << "\t idx: " << af.get_idx() << '\n';
    detail << "\t Note: idx indicates rows of alpha that A determines.\n";
    detail << '\n';
    detail << "\t beta:\n";
    for (INTEGER i=0; i<af.get_ncolA(); i++) {
      detail << "\t\t   Row " << fmt('d',3,i+1) << "   Lambda ";
      for (INTEGER j=1; j<=beta[i].size(); ++j) {
        detail << beta[i][j];
      }
      detail << '\n';
    }
    detail << '\n';
    detail << "\t Note: If you make modifications, be aware that beta is\n";
    detail << "\t       indexed from zero in the code.\n\n";
    detail.flush();
  }
  
  if (opm.print) {
    detail << '\n' << "Revised starting values for a0" << '\n';
    detail << "(Values are from reconstructed object)";
    detail << af.get_a0();
    detail << "\n(Mask, coded 0 for fixed, 1 for active)";
    detail << "\n(But, in the code sometimes coded -1 for fixed, 0 for active)";
    mask_a0 = af_mask.get_a0();
    if (mask_a0.size()>0) ++mask_a0;
    detail << mask_a0 << '\n';  
    if (mask_a0.size()>0) --mask_a0;
    detail.flush();
  }
  
  if (opm.print) {
    detail << '\n' << "Revised starting values for A" << '\n';
    detail << "(Values are from reconstructed object)";
    detail << af.get_A();
    detail << "\n(Mask, coded 0 for fixed, 1 for active)";
    detail << "\n(But, in the code sometimes coded -1 for fixed, 0 for active)";
    mask_A = af_mask.get_A();
    mask_A[1] = -1.0;
    af_mask.set_A(mask_A);
    detail << ++mask_A << '\n';  --mask_A;
    detail.flush();
  }

  header=kw.set_keyword("MEAN FUNCTION START VALUE DESCRIPTION"); //required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) return rv = false;
  if (kw_ptr + 2 >= pf.end()) error("Error, parmfile, " + header);

  INTEGER lagu = atoi((++kw_ptr)->substr(0,10).c_str());
  bool icept = atoi((++kw_ptr)->substr(0,10).c_str());

  INTEGER lagu_final = lagu;
  
  header=kw.set_keyword("MEAN FUNCTION DESCRIPTION"); //will be reread below
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
    if (kw_ptr + 2 >= pf.end()) error("Error, parmfile, " + header);
    lagu_final += atoi((++kw_ptr)->substr(0,10).c_str());
  }

  INTEGER ilx = lagu_final > 0 ? dpm.M : 0; 
  
  ufunc uf_start(dpm.M, ilx, lagu, icept);
  ufunc uf_mask_start(dpm.M, ilx, lagu, icept);

  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "(Values are from constructed object)\n";
    detail << "\t lag = " << uf_start.get_lag() << " = Lu\n";
    detail << "\t intercept = " << uf_start.is_intercept() << '\n';
    detail << '\n';
    detail << "\t regression = " << uf_start.is_regression() << '\n';
    detail << "\t ly = " << uf_start.get_ly() << " = M\n";
    detail << "\t lx = " << uf_start.get_lx() << " = M if Lu > 0 below\n";
    detail.flush();
  }

  realmat mask_b0 = uf_start.get_b0();  fill(mask_b0);

  header=kw.set_keyword("MEAN FUNCTION START VALUES FOR b0"); //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
  
    realmat b0 = uf_start.get_b0();
    if (kw_ptr+b0.size()>=pf.end()) error("Error, parmfile, " + header);
    for (INTEGER i=1; i<=b0.size(); ++i) {
      pair<string,string> ss = find_ss((++kw_ptr)->c_str());
      b0[i] = atof(ss.first.c_str());
      mask_b0[i] = atof(ss.second.c_str()) - 1.0;
    }
    uf_start.set_b0(b0);
  }

  if (mask_b0.size()>0) uf_mask_start.set_b0(mask_b0);

  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "(Values are from constructed object)";
    detail << uf_start.get_b0();
    detail << "\n(Mask, coded 0 for fixed, 1 for active)";
    detail << "\n(But, in the code sometimes coded -1 for fixed, 0 for active)";
    mask_b0 = uf_mask_start.get_b0();
    if (mask_b0.size()>0) ++mask_b0;
    detail << mask_b0 << '\n';
    if (mask_b0.size()>0) --mask_b0;
    detail.flush();
  }
  
  realmat mask_B = uf_start.get_B();  fill(mask_B);

  header=kw.set_keyword("MEAN FUNCTION START VALUES FOR B"); //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
  
    realmat B = uf_start.get_B();
    if (kw_ptr+B.size()>=pf.end()) error("Error, parmfile, " + header);
    for (INTEGER i=1; i<=B.size(); ++i) {
      pair<string,string> ss = find_ss((++kw_ptr)->c_str());
      B[i] = atof(ss.first.c_str());
      mask_B[i] = atof(ss.second.c_str()) - 1.0;
    }
    uf_start.set_B(B);
  }

  if (mask_B.size()>0) uf_mask_start.set_B(mask_B);
  
  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "(Values are from constructed object)";
    detail << uf_start.get_B();
    detail << "\n(Mask, coded 0 for fixed, 1 for active)";
    detail << "\n(But, in the code sometimes coded -1 for fixed, 0 for active)";
    mask_B = uf_mask_start.get_B();
    if (mask_B.size()>0) ++mask_B;
    detail << mask_B << '\n';
    if (mask_B.size()>0) --mask_B;
    detail.flush();
  }

  header=kw.set_keyword("MEAN FUNCTION DESCRIPTION"); //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
    if (kw_ptr + 2 >= pf.end()) error("Error, parmfile, " + header);
  
    lagu += atoi((++kw_ptr)->substr(0,10).c_str());
    int icept_int = int(icept);
    icept_int += atoi((++kw_ptr)->substr(0,10).c_str());
    if (icept_int == 0) {icept = false;} else {icept = true;}
  }

  ufunc uf(uf_start, dpm.M, ilx, lagu, icept);
  ufunc uf_mask(uf_mask_start, dpm.M, ilx, lagu, icept);

  ufn = uf;
  ufn_mask = uf_mask;
  
  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "(Values are after application of increment or decrement)\n";
    detail << "(Values are from constructed object)\n";
    detail << "\t lag = " << uf.get_lag() << " = Lu\n";
    detail << "\t intercept = " << uf.is_intercept() << '\n';
    detail << '\n';
    detail << "\t regression = " << uf.is_regression() << '\n';
    detail << "\t ly = " << uf.get_ly() << " = M\n";
    detail << "\t lx = " << uf.get_lx() << " = M if Lu > 0\n";
    detail.flush();
  }
  
  if (opm.print) {
    detail << '\n' << "Revised starting values for b0" << '\n';
    detail << "(Values are from reconstructed object)";
    detail << uf.get_b0();
    detail << "\n(Mask, coded 0 for fixed, 1 for active)";
    detail << "\n(But, in the code sometimes coded -1 for fixed, 0 for active)";
    mask_b0 = uf_mask.get_b0();
    if (mask_b0.size()>0) ++mask_b0;
    detail << mask_b0 << '\n';
    if (mask_b0.size()>0) --mask_b0;
    detail.flush();
  }
    
  if (opm.print) {
    detail << '\n' << "Revised starting values for B" << '\n';
    detail << "(Values are from reconstructed object)";
    detail << uf.get_B();
    detail << "\n(Mask, coded 0 for fixed, 1 for active)";
    detail << "\n(But, in the code sometimes coded -1 for fixed, 0 for active)";
    mask_B = uf_mask.get_B();
    if (mask_B.size()>0) ++mask_B;
    detail << mask_B << '\n';
    if (mask_B.size()>0) --mask_B;
    detail.flush();
  }

  header=kw.set_keyword("VARIANCE FUNCTION START VALUE DESCRIPTION");//required
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr == pf.end()) return rv = false;
  if (kw_ptr + 8 >= pf.end()) error("Error, parmfile, " + header);

  INTEGER q = atoi((++kw_ptr)->substr(0,10).c_str());
  char Qtype = last_char((++kw_ptr)->substr(0,10));
  INTEGER p = atoi((++kw_ptr)->substr(0,10).c_str());
  char Ptype = last_char((++kw_ptr)->substr(0,10));
  INTEGER v = atoi((++kw_ptr)->substr(0,10).c_str());
  char Vtype = last_char((++kw_ptr)->substr(0,10));
  INTEGER w = atoi((++kw_ptr)->substr(0,10).c_str());
  char Wtype = last_char((++kw_ptr)->substr(0,10));

  rfunc rf_start(dpm.M, dpm.M, tpm.squash, tpm.inflec,
    p, Ptype, q, Qtype, v, Vtype, w, Wtype);
  rfunc rf_mask_start(dpm.M, dpm.M, tpm.squash, tpm.inflec,
    p, Ptype, q, Qtype, v, Vtype, w, Wtype);

  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "(Values are from constructed object)\n";
    detail << "\t p = " << rf_start.get_p() << " = Lr\n";
    detail << "\t Ptype = " << rf_start.get_Ptype() << '\n';
    detail << "\t q = " << rf_start.get_q() << " = Lg\n";
    detail << "\t Qtype = " << rf_start.get_Qtype() << '\n';
    detail << "\t v = " << rf_start.get_v() << " = Lv\n";
    detail << "\t Vtype = " << rf_start.get_Vtype() << '\n';
    detail << "\t w = " << rf_start.get_w() << " = Lw\n";
    detail << "\t Wtype = " << rf_start.get_Wtype() << '\n';
    detail << '\n';
    detail << "\t ly = " << rf_start.get_ly() << " = M\n";
    detail << "\t lx = " << rf_start.get_lx() << " = M\n";
    detail << "\t lR = " << rf_start.get_lR() << '\n';
    detail << "\t rowsP = " << rf_start.get_rowsP() << '\n';
    detail << "\t colsP = " << rf_start.get_colsP() << '\n';
    detail << "\t rowsQ = " << rf_start.get_rowsQ() << '\n';
    detail << "\t colsQ = " << rf_start.get_colsQ() << '\n';
    detail << "\t rowsV = " << rf_start.get_rowsV() << '\n';
    detail << "\t colsV = " << rf_start.get_colsV() << '\n';
    detail << "\t rowsW = " << rf_start.get_rowsW() << '\n';
    detail << "\t colsW = " << rf_start.get_colsW() << '\n';
    detail << "\t lRparms = " << rf_start.get_lRparms() << '\n';
    detail.flush();
  }
  
  realmat mask_Rparms = rf_start.get_Rparms();  fill(mask_Rparms);

  header=kw.set_keyword("VARIANCE FUNCTION START VALUES FOR Rparms");//optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
  
    realmat Rparms = rf_start.get_Rparms();
    if (kw_ptr+Rparms.size()>=pf.end()) error("Error, parmfile, " + header);
    for (INTEGER i=1; i<=Rparms.size(); ++i) {
      pair<string,string> ss = find_ss((++kw_ptr)->c_str());
      Rparms[i] = atof(ss.first.c_str());
      mask_Rparms[i] = atof(ss.second.c_str()) - 1.0;
    }
    rf_start.set_Rparms(Rparms);
  }
  
  rf_mask_start.set_Rparms(mask_Rparms);

  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "(Values are from constructed object)\n\n";
    detail << "\t Rparms: " << rf_start.get_Rparms();
    detail << "\n(Mask, coded 0 for fixed, 1 for active)";
    detail << "\n(But, in the code sometimes coded -1 for fixed, 0 for active)";
    detail << "\n\n";
    detail << "\t Rparms mask: " << ++mask_Rparms << '\n'; --mask_Rparms;
    detail << '\n';
    detail << "(Rparms is the vec of the following matrices)\n\n";
    detail << "\t R0: " << rf_start.get_R0() << '\n';
    detail << "\t P: " << rf_start.get_P() << '\n';
    detail << "\t Q: " << rf_start.get_Q() << '\n';
    detail << "\t V: " << rf_start.get_V() << '\n';
    detail << "\t W: " << rf_start.get_W() << '\n';
    detail.flush();
  }

  header=kw.set_keyword("VARIANCE FUNCTION DESCRIPTION"); //optional
  kw_ptr = find_if(pf.begin(), pf.end(), kw);
  if (kw_ptr != pf.end()) {
  if (kw_ptr + 8 >= pf.end()) error("Error, parmfile, " + header);
  
    q += atoi((++kw_ptr)->substr(0,10).c_str());
    Qtype = last_char((++kw_ptr)->substr(0,10));
    p += atoi((++kw_ptr)->substr(0,10).c_str());
    Ptype = last_char((++kw_ptr)->substr(0,10));
    v += atoi((++kw_ptr)->substr(0,10).c_str());
    Vtype = last_char((++kw_ptr)->substr(0,10));
    w += atoi((++kw_ptr)->substr(0,10).c_str());
    Wtype = last_char((++kw_ptr)->substr(0,10));
  }
  
  rfunc rf(rf_start, dpm.M, dpm.M, tpm.squash, tpm.inflec, 
    p, Ptype, q, Qtype, v, Vtype, w, Wtype);
  rfunc rf_mask(rf_mask_start, dpm.M, dpm.M, tpm.squash, tpm.inflec, 
    p, Ptype, q, Qtype, v, Vtype, w, Wtype);

  rfn = rf;
  rfn_mask = rf_mask;
  
  if (opm.print) {
    detail << '\n' << header << '\n';
    detail << "(Values are after application of increment or decrement)\n";
    detail << "(Values are from reconstructed object)\n";
    detail << "\t p = " << rf.get_p() << " = Lr\n";
    detail << "\t Ptype = " << rf.get_Ptype() << '\n';
    detail << "\t q = " << rf.get_q() << " = Lg\n";
    detail << "\t Qtype = " << rf.get_Qtype() << '\n';
    detail << "\t v = " << rf.get_v() << " = Lv\n";
    detail << "\t Vtype = " << rf.get_Vtype() << '\n';
    detail << "\t w = " << rf.get_w() << " = Lw\n";
    detail << "\t Wtype = " << rf.get_Wtype() << '\n';
    detail << '\n';
    detail << "\t ly = " << rf.get_ly() << " = M\n";
    detail << "\t lx = " << rf.get_lx() << " = M\n";
    detail << "\t lR = " << rf.get_lR() << '\n';
    detail << "\t rowsP = " << rf.get_rowsP() << '\n';
    detail << "\t colsP = " << rf.get_colsP() << '\n';
    detail << "\t rowsQ = " << rf.get_rowsQ() << '\n';
    detail << "\t colsQ = " << rf.get_colsQ() << '\n';
    detail << "\t rowsV = " << rf.get_rowsV() << '\n';
    detail << "\t colsV = " << rf.get_colsV() << '\n';
    detail << "\t rowsW = " << rf.get_rowsW() << '\n';
    detail << "\t colsW = " << rf.get_colsW() << '\n';
    detail << "\t lRparms = " << rf.get_lRparms() << '\n';
    detail.flush();
  }
  
  if (opm.print) {
    detail << '\n' << "Revised starting values for Rparms" << '\n';
    detail << "(Values are from reconstructed object)";
    detail << rf.get_Rparms();
    detail << "\n(Mask, coded 0 for fixed, 1 for active)";
    detail << "\n(But, in the code sometimes coded -1 for fixed, 0 for active)";
    mask_Rparms = rf_mask.get_Rparms();
    detail << ++mask_Rparms << '\n'; --mask_Rparms;
    detail << '\n';
    detail << "(Rparms is the vec of the following matrices)\n\n";
    detail << "\t R0: " << rf.get_R0() << '\n';
    detail << "\t P: " << rf.get_P() << '\n';
    detail << "\t Q: " << rf.get_Q() << '\n';
    detail << "\t V: " << rf.get_V() << '\n';
    detail << "\t W: " << rf.get_W() << '\n';
    detail.flush();
  }

  return rv = true;
}

bool 
snp::parmfile::write_parms(const char* filename, string ctrl, snpll& ll) const
{
  ofstream ps(filename);
  if (!ps) return false; 

  ps << "PARMFILE HISTORY (optional)\n"
     << "#\n"
     << "# This parmfile was written by SNP Version " 
     << fmt('f',3,1,opm.snpver) << ' '
     << "using the following line from\n" 
     << "# control.dat, which was read as " 
     << "char*, char*, float, float, int, int\n"
     << "# -------------------------------------------"
     << "----------------------------------\n"
     << "#      input_file          output_file         fnew      fold" 
     << "   nstart   jseed\n"
     << "# -------------------- -------------------- --------- ---------"
     << " ------ --------\n"
     << "# " << ctrl << '\n'
     << "# -------------------------------------------"
     << "----------------------------------\n"
     << "# If fnew is negative, only the polynomial part of the model "
     << "is perturbed.\n"
     << "# Similarly for fold.\n"
     << "#\n";
  for (vector<string>::const_iterator i=history.begin();i!=history.end();++i) {
    ps << *i << '\n';
  }

  ps << "OPTIMIZATION DESCRIPTION (required)\n"
     << opm.pname << "     "
     << "Project name, pname, char*\n"
     << fmt('f', 10, 1, opm.snpver) << "     "
     << "SNP version, defines format of this file, snpver, float\n"
     << fmt('d', 10, opm.itmax0) << "     "
     << "Maximum number of primary iterations, itmax0, int\n"
     << fmt('d', 10, opm.itmax1) << "     "
     << "Maximum number of secondary iterations, itmax1, int\n"
     << fmt('e', 10, 2, opm.toler) << "     "
     << "Convergence tolerance, toler, float\n"
     << fmt('d', 10, opm.print) << "     "
     << "Write detailed output if print=1, int\n"
     << fmt('d', 10, opm.task) << "     "
     << "task, 0 fit, 1 res, 2 mu, 3 sig, 4 plt, 5 sim, 6 usr, int\n"
     << fmt('d', 10, opm.extra) << "     "
     << "Increase simulation length by extra, int\n"
     << fmt('e', 10, 2, opm.sfac) << "     "
     << "Scale factor for plots, sfac, float\n"
     << fmt('d', 10, opm.iseed) << "     "
     << "Seed for simulations, iseed, int\n"
     << fmt('d', 10, opm.ngrid) << "     "
     << "Number of plot grid points, ngrid, int\n"
     << fmt('d', 10, opm.kilse) << "     "
     << "Statistics not computed if kilse=1, int\n";

  ps << "DATA DESCRIPTION (required)\n"
     << fmt('d', 10, dpm.M) << "     "
     << "Dimension of the time series, M, int\n"
     << fmt('d', 10, dpm.n) << "     "
     << "Number of observations, n, int\n"
     << fmt('d', 10, dpm.drop) << "     "
     << "Provision for initial lags, must have 3<drop<n, int\n"
     << fmt('d', 10, dpm.cond) << "     "
     << "Condition set for plt is mean if cond=0, it-th obs if it, int\n"
     << fmt('d', 10, dpm.reread) << "     "
     << "Reread, do not use data from prior fit, if reread=1, int\n"
     << dpm.dsn << ' ';
  for (INTEGER i=dpm.dsn.size(); i<14; ++i) ps << ' ';
  ps << "File name, any length, no embedded blanks, dsn, string\n";
  for (INTEGER i=1; i<=dpm.M; ++i) ps << dpm.fields[i] << ' ';
  for (INTEGER i=2*dpm.M; i<15; ++i) ps << ' ';
  ps << "Read these white space separated fields, fields, intvec\n";

  ps << "TRANSFORM DESCRIPTION (required)\n"
     << fmt('d', 10, tpm.useold) << "     "
     << "Normalize using start values if useold=1 else compute, int\n"
     << fmt('d', 10, tpm.diag) << "     "
     << "Make variance matrix diagonal if diag=1, int\n"
     << fmt('d', 10, tpm.squash) << "     "
     << "Spline transform x if squash=1, logistic if squash=2, int\n"
     << fmt('e', 10, 2, tpm.inflec) << "     "
     << "Inflection point of transform in normalized x, inflec, float\n";

  ps << "POLYNOMIAL START VALUE DESCRIPTION (required)\n"
     << fmt('d', 10, fn.get_Kz()) << "     "
     << "Degree of the polynomial in z, Kz, int\n"
     << fmt('d', 10, fn.get_Iz()) << "     "
     << "Degree of interactions in z, Iz, int\n"
     << fmt('e', 10, 2, fn.get_e0()) << "     "
     << "Zero or positive to get positive SNP for EMM, eps0, float\n"
     << fmt('d', 10, afn.get_lag()) << "     "
     << "Lags in polynomial part, Lp, int\n"
     << fmt('d', 10, afn.get_maxKz()) << "     "
     << "Max degree of z polynomial that depends on x, maxKz, int\n"
     << fmt('d', 10, afn.get_maxIz()) << "     "
     << "Max interaction of z polynomial that depends on x, maxIz, int\n"
     << fmt('d', 10, afn.get_Kx()) << "     "
     << "Degree of the polynomial in x, Kx, int\n"
     << fmt('d', 10, afn.get_Ix()) << "     "
     << "Degree of the interactions x, Ix, int\n";

  ps << "MEAN FUNCTION START VALUE DESCRIPTION (required)\n"
     << fmt('d', 10, ufn.get_lag()) << "     "
     << "Lags in VAR part, Lu, int\n"
     << fmt('d', 10, ufn.is_intercept()) << "     "
     << "Intercept if icept=1, int\n";

  ps << "VARIANCE FUNCTION START VALUE DESCRIPTION (required)\n"
     << fmt('d', 10, rfn.get_q()) << "     "
     << "Lags in GARCH (autoregressive) part, may be zero, Lg, int\n"
     << "         " << rfn.get_Qtype() << "     "
     << "Coded 's','d','f' for scalar, diagonal, full, Qtype, char\n"
     << fmt('d', 10, rfn.get_p()) << "     "
     << "Lags in ARCH (moving average) part, may be zero, Lr, int\n"
     << "         " << rfn.get_Ptype() << "     "
     << "Coded 's','d','f' for scalar, diagonal, full, Ptype, char\n"
     << fmt('d', 10, rfn.get_v()) << "     "
     << "Lags in leverage effect of GARCH, may be zero, Lv, int\n"
     << "         " << rfn.get_Vtype() << "     "
     << "Coded 's','d','f' for scalar, diagonal, full, Vtype, char\n"
     << fmt('d', 10, rfn.get_w()) << "     "
     << "Lags in additive level effect, may be zero, Lw, int\n"
     << "         " << rfn.get_Wtype() << "     "
     << "Coded 's','d','f' for scalar, diagonal, full, Wtype, char\n";

  ps << "POLYNOMIAL DESCRIPTION (optional)\n"
     << "         0     Increment or decrement to Kz, int\n"
     << "         0     Increment or decrement to Iz, int\n"
     << "  0.00e+00     Increment or decrement to eps0, float\n"
     << "         0     Increment or decrement to Lp, int\n"
     << "         0     Increment or decrement to maxKz, int\n"
     << "         0     Increment or decrement to maxIz, int\n"
     << "         0     Increment or decrement to Kx, int\n"
     << "         0     Increment or decrement to Ix, int\n";

  ps << "MEAN FUNCTION DESCRIPTION (optional)\n"
     << "         0     Increment or decrement to Lu, int\n"
     << "         0     Increment or decrement to icept, int\n";

  ps << "VARIANCE FUNCTION DESCRIPTION (optional)\n"
     << "         0     Increment or decrement to GARCH lag Lg, int\n"
     << "         " << rfn.get_Qtype() << "     "
     << "Coded 's','d','f' for scalar, diagonal, full, Qtype, char\n"
     << "         0     Increment or decrement to ARCH lag Lr, int\n"
     << "         " << rfn.get_Ptype() << "     "
     << "Coded 's','d','f' for scalar, diagonal, full, Ptype, char\n"
     << "         0     Increment or decrement to leverage effect lag Lv, int\n"
     << "         " << rfn.get_Vtype() << "     "
     << "Coded 's','d','f' for scalar, diagonal, full, Vtype, char\n"
     << "         0     Increment or decrement to level effect lag Lw, int\n"
     << "         " << rfn.get_Wtype() << "     "
     << "Coded 's','d','f' for scalar, diagonal, full, Wtype, char\n";

  realmat a0 = afn.get_a0();
  realmat mask_a0 = afn_mask.get_a0();
  if (a0.size() > 0) {
    ps << "POLYNOMIAL START VALUES FOR a0 (optional)\n";
    ++mask_a0;
    for (INTEGER i=1; i<=a0.size(); ++i) {
      ps << fmt('e',26,17,a0[i]) << fmt('d',5,INTEGER(mask_a0[i])) << '\n';
    }
  }

  realmat A = afn.get_A();
  realmat mask_A = afn_mask.get_A();
  ++mask_A;
  mask_A[1] = 0;
  if (A.size() > 0) {
    ps << "POLYNOMIAL START VALUES FOR A (optional)\n";
    for (INTEGER i=1; i<=A.size(); ++i) {
      ps << fmt('e',26,17,A[i]) << fmt('d',5,INTEGER(mask_A[i])) << '\n';
    }
  }

  realmat b0 = ufn.get_b0();
  realmat mask_b0 = ufn_mask.get_b0();
  if (b0.size() > 0) {
    ps << "MEAN FUNCTION START VALUES FOR b0 (optional)\n";
    ++mask_b0;
    for (INTEGER i=1; i<=b0.size(); ++i) {
      ps << fmt('e',26,17,b0[i]) << fmt('d',5,INTEGER(mask_b0[i])) << '\n';
    }
  }

  realmat B = ufn.get_B();
  realmat mask_B = ufn_mask.get_B();
  if (B.size() > 0) {
    ps << "MEAN FUNCTION START VALUES FOR B (optional)\n";
    ++mask_B;
    for (INTEGER i=1; i<=B.size(); ++i) {
      ps << fmt('e',26,17,B[i]) << fmt('d',5,INTEGER(mask_B[i])) << '\n';
    }
  }

  rfunc rf = rfn;  //Because write_parms is const and get_Rparms() is not.
  rfunc rf_mask = rfn_mask; 
  realmat Rparms = rf.get_Rparms();
  realmat mask_Rparms = rf_mask.get_Rparms();
  ++mask_Rparms;
  ps << "VARIANCE FUNCTION START VALUES FOR Rparms (optional)\n";
  for (INTEGER i=1; i<=Rparms.size(); ++i) {
      ps << fmt('e',26,17,Rparms[i]) 
          << fmt('d',5,INTEGER(mask_Rparms[i])) << '\n';
  }

  ps << "TRANSFORM START VALUES FOR mean (optional)\n";
  for (INTEGER i=1; i<=tpm.mean.size(); ++i) {
    ps << fmt('e',26,17,tpm.mean[i]) << '\n';
  }

  ps << "TRANSFORM START VALUES FOR variance (optional)\n";
  for (INTEGER i=1; i<=tpm.variance.size(); ++i) {
    ps << fmt('e',26,17,tpm.variance[i]) << '\n';
  }

  if (!opm.kilse) {
    ps << "SUMMARY STATISTICS (optional)\n";
    ll.write_stats(ps);
  }

  return ps.good();
}