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

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.

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

#include <cerrno>
#include "libsmm.h"
#include "emm.h"

using namespace libcpp;
using namespace libsmm;
using namespace emm;
using namespace std;

emm::sv_usrmod::sv_usrmod
  (const realmat& dat, INTEGER len_mod_parm, INTEGER len_mod_func,
   const std::vector<std::string>& mod_pfvec,
   const std::vector<std::string>& mod_alvec, 
   std::ostream& detail)
: data(dat), lrho(6), lstats(8), 
  variable_seed(740726), fixed_seed(770116)
{
  vector<string>::const_iterator usr_ptr = mod_alvec.begin();
  ++usr_ptr;
  slen = atoi((usr_ptr++)->substr(0,12).c_str());
  spin = atoi((usr_ptr++)->substr(0,12).c_str());

  if (lrho != len_mod_parm) {
    error("Error, usrmod, constructor, len_mod_parm is set wrong in parmfile");
  }

  if (lstats != len_mod_func) {
    error("Error, usrmod, constructor, len_mod_parm is set wrong in parmfile");
  }

}

bool emm::sv_usrmod::sv_simulate
       (INT_32BIT& seed, realmat& sim, realmat& stats, realmat& latent)
{
  if (sim.get_rows() != 1 && sim.get_cols() != slen) {
    sim.resize(1,slen);
  }
  if (latent.get_rows() != 1 && latent.get_cols() != slen) {
    latent.resize(1,slen);
  }

  REAL a0 = rho[1];
  REAL a1 = rho[2];
  REAL b0 = rho[3];
  REAL b1 = rho[4];
  REAL s  = rho[5];
  REAL r  = rho[6];

  REAL rr2 = sqrt(1.0 - pow(r,2));

  REAL vlag = 0.0;
  REAL ylag = 0.0;
  errno = 0;  

  for (INTEGER t=1; t<=spin; ++t) {
    REAL z1 = unsk(seed);
    REAL z2 = unsk(seed);
    REAL u1 = z1;
    REAL u2 = s*(r*z1 + rr2*z2);
    REAL v = b0 + b1*(vlag - b0) + u2;
    REAL y = a0 + a1*(ylag - a0) + u1*exp(v); // or u1*exp(vlag)
    vlag = v;                                 // see User's Guide
    ylag = y;
  }

  if (errno == ERANGE) return false; 

  REAL ymin = REAL_MAX;
  REAL ymax = -REAL_MAX;
  REAL ymean = 0.0;
  REAL ysdev = 0.0;
  REAL vmin = REAL_MAX;
  REAL vmax = -REAL_MAX;
  REAL vmean = 0.0;
  REAL vsdev = 0.0;

  for (INTEGER t=1; t<=slen; ++t) {
    REAL z1 = unsk(seed);
    REAL z2 = unsk(seed);
    REAL u1 = z1;
    REAL u2 = s*(r*z1 + rr2*z2);
    REAL v = b0 + b1*(vlag - b0) + u2;
    REAL y = a0 + a1*(ylag - a0) + u1*exp(v); // or u1*exp(vlag)
    vlag = v;                                 // see User's Guide
    ylag = y;
    sim[t] = y;
    latent[t] = v;
    ymin = y < ymin ? y : ymin;
    ymax = y > ymax ? y : ymax;
    ymean += y;
    ysdev += pow(y,2);
    vmin = v < vmin ? v : vmin;
    vmax = v > vmax ? v : vmax;
    vmean += v;
    vsdev += pow(v,2);
  }

  if (errno == ERANGE) return false; 

  if (lstats != 8) error("Error, emmusr, wrong size for stats");
  if (stats.size() != lstats) stats.resize(lstats,1);

  ymean = ymean/REAL(slen);
  ysdev = sqrt( (ysdev - REAL(slen)*pow(ymean,2))/REAL(slen) );

  vmean = vmean/REAL(slen);
  vsdev = sqrt( (vsdev - REAL(slen)*pow(vmean,2))/REAL(slen) );

  stats[1] = ymin;
  stats[2] = ymax;
  stats[3] = ymean;
  stats[4] = ysdev;
  stats[5] = vmin;
  stats[6] = vmax;
  stats[7] = vmean;
  stats[8] = vsdev;

  return true;
}

bool emm::sv_usrmod::gen_bootstrap(vector<realmat>& bs)
{
  realmat sim, stats, latent;;

  if (!support(rho)) return false;

  if (!sv_simulate(variable_seed,sim,stats,latent)) return false;

  INTEGER n_sim = sim.get_cols();
  INTEGER n_dat = data.get_cols();

  if (n_sim < n_dat) return false;

  INTEGER len_vec = 2*lrho+1;
  INTEGER inc, len;

  if (n_sim/n_dat < len_vec) {
    len = n_sim/n_dat;
    inc = n_dat;
  }
  else {
    len = len_vec;
    inc = n_sim/len;
  }

  len_vec = bs.size();

  if (len_vec != len) {
    bs.resize(len);
  }

  INTEGER s = 0;
  for (INTEGER i=0; i<len; i++) {
    bs[i] = sim("",seq(s+1,s+n_dat));
    s+=inc;
  }

  return true;
}

bool emm::sv_usrmod::support(const realmat& parm) 
{
  if (parm[5] < 0.0) return false;
  if (fabs(parm[6]) > 1.0) return false;
  return true; 
}

den_val emm::sv_usrmod::prior(const realmat& parm, const realmat& stats) 
{
  return den_val(true, 0.0);
}