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

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

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA   

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

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

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

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

int main(int argc, char** argp, char** envp)
{
  istream* ctrl_ptr;
  if (argc == 2) {
    ctrl_ptr = new(nothrow) ifstream(argp[1]);
    string msg = "Error, gsm, " + string(argp[1]) + " open failed";
    if( (ctrl_ptr == 0) || (!*ctrl_ptr) ) error(msg);
  }
  else {
    ctrl_ptr = new(nothrow) ifstream("control.dat");
    string msg = "Error, gsm, control.dat open failed";
    if( (ctrl_ptr == 0) || (!*ctrl_ptr) ) error(msg);
  }
  istream& ctrl = *ctrl_ptr;

  string parmfile, prefix;
 
  while(ctrl >> parmfile >> prefix) {

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

    ofstream detail_ofs(filename.c_str());
    if(!detail_ofs) error("Error, gsm, detail.dat open failed");
    ostream& detail = detail_ofs;

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

    string line;

    vector<string> pfvec;
    while (getline(pf_ifs, line)) pfvec.push_back(line);

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

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

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

    if (est_blk.print) {
      detail << starbox("/First 12 observations//");
      detail << data("",seq(1,12));
      detail << starbox("/Last 12 observations//");
      detail << data("",seq(dat_blk.n-11, dat_blk.n));
      detail.flush();
    }
        
    if (est_blk.print) {
      detail << starbox("/Begin instantiation of stat_mod//");
      detail.flush();
    }
        
    vector<string> stat_pfvec;

    if (stat_blk.is_stat_parmfile) {

      string filename = stat_blk.stat_parmfile;

      ifstream stat_pf_ifs(filename.c_str());
      if (!stat_pf_ifs) error("Error, gsm, cannot open " + filename);
  
      while (getline(stat_pf_ifs, line)) stat_pfvec.push_back(line);

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

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

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

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

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

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

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

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

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

    if (est_blk.print) {
      detail << starbox("/End instantiation of stat_mod//");
      detail.flush();
    }
        
    if (est_blk.print) {
      detail << starbox("/Begin instantiation of sci_mod//");
      detail.flush();
    }
        
    vector<string> sci_pfvec;

    if (sci_blk.is_sci_parmfile) {

      string filename = sci_blk.sci_parmfile;

      ifstream sci_pf_ifs(filename.c_str());
      if (!sci_pf_ifs) error("Error, gsm, cannot open " + filename);
  
      while (getline(sci_pf_ifs, line)) sci_pfvec.push_back(line);

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

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

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

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

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

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

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

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

    sci_mod.set_parm(pf.get_sci_parm());
    
    if (est_blk.print) {
      detail << starbox("/Parameter values//");
      realmat parm;
      sci_mod.get_parm(parm);
      detail << parm << '\n';
      realmat sim, func;
      sci_mod.gen_sim(sim, func);
      detail << starbox("/First 12 sci_mod simulations//");
      detail << sim("",seq(1,12));
      detail << starbox("/Last 12 sci_mod simulations//");
      intvec idx = seq(sim.get_cols()-11, sim.get_cols());
      detail << sim("",idx);
      detail << starbox("/sci_mod func//");
      detail << func << '\n';
      detail.flush();
    }

    if (est_blk.print) {
      detail << starbox("/End instantiation of sci_mod//");
      detail.flush();
    }

    if (est_blk.print) {
      detail << starbox("/Instantiate the proposal densities//");
      detail.flush();
    }

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

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

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

    filename = prefix + ".implied_map.dat";

    if(!im.read_map(filename.c_str())) {
    
      if (est_blk.run_sci_chain) {
        warn("Warning, gsm, could not read implied map, will compute afresh");
      }
      else if (est_blk.run_stat_chain) {
        error("Error, gsm, could not read implied map, you must run sci_chain");
      }

    }

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

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

    filename = prefix + ".assess_sigma.dat";

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

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

    }
    else {

      if (est_blk.run_sci_chain) {
        warn("Warning, gsm, could not read assess sigma, will compute afresh");
      }
      else if (est_blk.run_stat_chain) {
        error("Error, gsm, could not read sigma, you must run sci_chain");
      }
    }

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

    if (est_blk.run_sci_chain) {

      if (est_blk.print) {
        detail << starbox("/Instantiate sci_mcmc chain//");
        detail.flush();
      }

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

      realmat sci_sub_parm_chain;
      realmat sci_sub_logl_chain;
      realmat sci_sub_rej_chain;

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

      realmat assess_prior_chain;

      INT_32BIT seed = sci_blk.sci_seed;

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

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

       if (est_blk.print) {
	 if (ifile == 0) detail << '\n';
         detail << "\tCumulative implied map hit rate " 
	        << im.get_hit_rate() << '\n';
         detail.flush();
        }

        char number[5];
        sprintf(number,"%03d",ifile);

        filename = prefix + ".sci_parm." + number + ".dat";
        vecwrite(filename.c_str(),sci_parm_chain);

        filename = prefix + ".sci_func." + number + ".dat";
        vecwrite(filename.c_str(),sci_func_chain);

        filename = prefix + ".sci_prior." + number + ".dat";
        vecwrite(filename.c_str(),sci_prior_chain);

	if (! sci_blk.analytic_mle) {
          filename = prefix + ".sci_sub_parm." + number + ".dat";
          vecwrite(filename.c_str(),sci_sub_parm_chain);
  
          filename = prefix + ".sci_sub_logl." + number + ".dat";
          vecwrite(filename.c_str(),sci_sub_logl_chain);
  
          filename = prefix + ".sci_sub_rej." + number + ".dat";
          vecwrite(filename.c_str(),sci_sub_rej_chain);
	}
  
        filename = prefix + ".stat_parm." + number + ".dat";
        vecwrite(filename.c_str(),stat_parm_chain);

        filename = prefix + ".stat_func." + number + ".dat";
        vecwrite(filename.c_str(),stat_func_chain);

        filename = prefix + ".stat_prior." + number + ".dat";
        vecwrite(filename.c_str(),stat_prior_chain);

        filename = prefix + ".stat_logl." + number + ".dat";
        vecwrite(filename.c_str(),stat_logl_chain);

        filename = prefix + ".stat_sub_logl." + number + ".dat";
        vecwrite(filename.c_str(),stat_sub_logl_chain);

        filename = prefix + ".sci_rej." + number + ".dat";
        vecwrite(filename.c_str(),sreject);

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

        filename = prefix + ".assess_sigma.new";
        ap.write_sigma(filename.c_str());

        sci_val mode = schain.get_mode();

        filename = prefix + ".sci_mode_sci_val." + number + ".dat";
        ofstream mode_ofs(filename.c_str());
        if (!mode_ofs) error("Error, gsm, open failed for " + filename);
        mode_ofs << mode;

        pf.set_sci_parm(sci_parm_start);
        pf.set_stat_parm(stat_parm_start);
        
        if(!pf.write_blks(parmfile, prefix, string("end"))) {
          detail.flush();
          error("Error, gsm, cannot write new gsm parmfile");
        }

        string suffix = number + string(".end");
        
        if(!pf.write_blks(parmfile, prefix, suffix)) {
          detail.flush();
          error("Error, gsm, cannot write new gsm parmfile");
        }

        sci_val overall_mode = im.get_mode();
        filename = prefix + ".sci_mode_sci_val.dat";
        ofstream overall_mode_ofs(filename.c_str());
        if (!mode_ofs) error("Error, gsm, open failed for " + filename);
        overall_mode_ofs << mode;

        pf.set_sci_parm(overall_mode.sci_mod_parm);
        pf.set_stat_parm(overall_mode.stat_mod_parm);

        if(!pf.write_blks(parmfile, prefix, string("fit"))) {
          detail.flush();
          error("Error, gsm, cannot write new gsm parmfile");
        }

        sci_sum.update(sci_parm_chain,sci_func_chain,stat_parm_chain,
          stat_func_chain,overall_mode); 
        if(!sci_sum.write(prefix)) {
          error("Error, gsm, cannot write summary");
        }
      }
    }

    if (est_blk.run_stat_chain) {

      ap.set_implied_map(im);

      if (est_blk.print) {
        detail << starbox("/Instantiate stat_mcmc assess chain//");
        detail.flush();
      }

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

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

      INT_32BIT seed = stat_blk.stat_seed;

      stat_parm_start = stat_parm_start_save;

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

        realmat areject 
          = achain.draw(seed, stat_parm_start, assess_stat_parm_chain,
            assess_stat_func_chain, assess_stat_prior_chain,
            assess_stat_logl_chain, assess_prior_chain);
                                                                             
        char number[5];
        sprintf(number,"%03d",ifile);

        filename = prefix + ".as_st_parm." + number + ".dat";
        vecwrite(filename.c_str(),assess_stat_parm_chain);

        filename = prefix + ".as_st_func." + number + ".dat";
        vecwrite(filename.c_str(),assess_stat_func_chain);

        filename = prefix + ".as_st_prior." + number + ".dat";
        vecwrite(filename.c_str(),assess_stat_prior_chain);

        filename = prefix + ".as_st_logl." + number + ".dat";
        vecwrite(filename.c_str(),assess_stat_logl_chain);

        filename = prefix + ".as_prior." + number + ".dat";
        vecwrite(filename.c_str(),assess_prior_chain);

        filename = prefix + ".as_rej." + number + ".dat";
        vecwrite(filename.c_str(),areject);

        stat_val mode = achain.get_mode();

        filename = prefix + ".as_st_mode_stat_val.dat";
        ofstream mode_ofs(filename.c_str());
        if (!mode_ofs) error("Error, gsm, open failed for " + filename);
        mode_ofs << mode;

        pf.set_stat_parm(mode.stat_mod_parm);

        if(!pf.write_blks(parmfile, prefix, string("alt"))) {
          detail.flush();
          error("Error, gsm, cannot write new gsm parmfile");
        }

        assess_sum.update(assess_stat_parm_chain, assess_stat_func_chain,mode); 
        if(!assess_sum.write(prefix)) {
          error("Error, gsm, cannot write summary");
        }
      }
    }
  }

  delete ctrl_ptr;
  return 0;
}


