#ifndef __FILE_LIBSMM_H_SEEN__
#define __FILE_LIBSMM_H_SEEN__ 

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

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

A. Ronald Gallant

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 "libsmm_base.h"
#include <queue>

namespace libsmm {

class smm_mcmc : public mcmc_base {
private:
  proposal_base& T;
  cachemgr_base& L;
  usrmod_base& U;
  INTEGER simulation_size;
  INTEGER stride;
  bool draw_from_posterior;
  REAL temp;
  scl::realmat mode;
  REAL high;
public:
  smm_mcmc(proposal_base& T_fn, cachemgr_base& L_fn, usrmod_base& U_mod) 
    : T(T_fn), L(L_fn), U(U_mod), simulation_size(1), stride(1),
      draw_from_posterior(true), temp(1.0), mode(), high(-REAL_MAX) { } 
  scl::realmat draw(INT_32BIT& seed, scl::realmat& rho_start,
    scl::realmat& rho_sim, scl::realmat& stats_sim, 
    scl::realmat& pi_sim);
  void set_simulation_size(INTEGER n);
  void set_stride(INTEGER k);
  void set_draw_from_posterior(bool from_posterior);
  void set_temp(REAL temperature);
  void set_mode(const scl::realmat& new_mode, REAL new_high);
  REAL get_temp() { return temp; }
  scl::realmat get_mode() { return mode; }
  REAL get_high() { return high; }
};


class grid_group_move : public proposal_base {
private:
  typedef std::map<scl::intvec,REAL,scl::intvec_cmp> pdf_type;
  typedef std::map<scl::intvec,REAL,scl::intvec_cmp>::iterator pdf_itr;
  struct prop_aux {        //The elements of increment must be (fractional)
    REAL prob;             //powers of two for a grid_group_move proposal to
    pdf_type pdf;          //work correctly. This, in turn, requires that
  };                       //the elements of ginc in the prop_def argument of
  scl::realmat increment;  //the constructor be (fractional) powers of two.
  prop_def pd;
  INTEGER  lparm;
  INTEGER  lpd;
  std::vector<prop_aux> pa;
  bool get_grid(const scl::realmat& parm, scl::realmat& parm_real, 
    scl::intvec& parm_int); 
  void make_pdf (INTEGER d, const scl::realmat& ginc, 
    const scl::realmat& Vinv, const scl::intvec& lo, 
    const scl::intvec& hi, scl::intvec& index, pdf_type& pdf);
public:
  grid_group_move(const prop_def& def);
  scl::den_val operator()(const scl::realmat& parm_old, 
    const scl::realmat& parm_new);
  void draw(INT_32BIT& seed, const scl::realmat& parm_old, 
    scl::realmat& parm_new); 
  INTEGER len_rho() { return lparm; }
  bool transition_is_symmetric() { return true; }
  std::ostream& write_proposal(std::ostream& os);
};


class group_move : public proposal_base {
private:
  struct prop_aux {
    REAL    prob;
    scl::realmat Rmat;
    scl::realmat Vinv;
    REAL    scale;
  };
  prop_def pd;
  INTEGER  lrho;
  INTEGER  lpd;
  std::vector<prop_aux> pa;
public:
  group_move(prop_def def);
  scl::den_val operator()(const scl::realmat& rho_old, 
    const scl::realmat& rho_new);
  void draw(INT_32BIT& seed, const scl::realmat& rho_old, 
    scl::realmat& rho_new); 
  INTEGER len_rho() { return lrho; }
  bool transition_is_symmetric() { return true; }
};

class conditional_move : public proposal_base {
private:
  struct possible_move {
    INTEGER y;
    scl::intvec  x;
    REAL    b0;
    scl::realmat b;
    REAL    s;
  };
  INTEGER  lrho;
  INTEGER  lpm;
  std::vector<possible_move> pm;
public:
  conditional_move(prop_def def, std::ostream& detail, bool print);
  scl::den_val operator()(const scl::realmat& rho_old, 
    const scl::realmat& rho_new);
  void draw(INT_32BIT& seed, const scl::realmat& rho_old, 
    scl::realmat& rho_new); 
  INTEGER len_rho() { return lrho; }
};

class cachemgr : public cachemgr_base{
private:
  usrmod_base& model;
  objfun_base& obj;
  std::queue< scl::realmat,std::list<scl::realmat> >::size_type cache_size;
  INTEGER hit;
  INTEGER miss;
  std::queue< scl::realmat,std::list<scl::realmat> > cache;
  struct val {
    REAL v;
    scl::realmat stats;
    val() : v(0), stats() {}
    val(REAL rv, scl::realmat st) : v(rv), stats(st) {}
  };
  std::map<scl::realmat,val,scl::realmat_cmp> past_values;
  REAL evaluate(const scl::realmat& rho, scl::realmat& stats);
public:
  cachemgr(usrmod_base& m, objfun_base& o, INTEGER csz) 
  : model(m), obj(o), cache_size(csz), hit(0), miss(0)
  { 
    cache_size = cache_size < 10 ? 10 : cache_size;
  }
  scl::den_val operator()(const scl::realmat& rho, scl::realmat& stats); 
  bool read_cache(const char* filename);
  bool write_cache(const char* filename) const;
  bool set_cache(const scl::realmat& cache_new);
  bool set_cache(const scl::realmat& cache_new, INTEGER limit);
  void get_cache(scl::realmat& cache_sav) const; 
  INTEGER len_rho() { return model.len_rho(); }
  INTEGER len_stats() { return model.len_stats(); }
  REAL cache_hit_rate() const { return REAL(hit)/REAL(hit + miss); }
};

class bstrap_asymptotics : public asymptotics_base {
private:
  const scl::realmat& data;
  INTEGER n_dat;
  usrmod_base& model;
  objfun_base& objfun;
  mcmc_base& mcmc;
  INTEGER p;
  scl::realmat rho_sum;
  scl::realmat rho_ctr;
  scl::realmat rho_sse;
  scl::realmat mean;
  scl::realmat mode;
  scl::realmat sscp;
  scl::realmat foc;
  INTEGER n_tot;
  scl::realmat I_mat;
  INTEGER I_reps;
  REAL high; // temperature does not affect, prior does ifdef ADJUST_FOR_PRIOR
  struct deriv_item {
    objfun_base* obj;
    scl::realmat der_obj;
  };
  std::vector<deriv_item> deriv_vec;
  struct info_item {
    scl::realmat score;
    REAL lo;
    REAL hi;
  };
  std::vector<info_item> info_vec;
public:
  bstrap_asymptotics(const scl::realmat& dat, usrmod_base& mod, 
    objfun_base& obfn, mcmc_base& mc) 
    : data(dat), n_dat(dat.get_cols()), model(mod), objfun(obfn), mcmc(mc),
      p(mod.len_rho()), rho_sum(p,1,0.0), rho_ctr(p,1,0.0), rho_sse(p,p,0.0), 
      mean(p,1,0.0), mode(p,1,0.0), sscp(p,p,0.0), foc(p,1,0.0), n_tot(0),
      I_mat(p,p,0.0), I_reps(0), high(-REAL_MAX) 
    { };
  bool set_asymptotics(const scl::realmat& sim); 
    // Aggregates if called repeatedly.
  void get_asymptotics(scl::realmat& rho_hat, scl::realmat& V_hat, 
    INTEGER& n);
  void get_asymptotics(scl::realmat& rho_mean, scl::realmat& rho_mode, 
    REAL& post_high, scl::realmat& I, scl::realmat& invJ, 
    scl::realmat& foc_hat, INTEGER& reps);
};

class minimal_asymptotics : public asymptotics_base {
private:
  const scl::realmat& data;
  INTEGER n_dat;
  mcmc_base& mcmc;
  INTEGER p;
  scl::realmat rho_sum;
  scl::realmat rho_ctr;
  scl::realmat rho_sse;
  scl::realmat mean;
  scl::realmat mode;
  scl::realmat sscp;
  scl::realmat foc;
  INTEGER n_tot;
  REAL high; // temperature does not affect, prior does ifdef ADJUST_FOR_PRIOR
public:
  minimal_asymptotics(const scl::realmat& dat, usrmod_base& mod, 
    mcmc_base& mc) 
    : data(dat), n_dat(dat.get_cols()), mcmc(mc), p(mod.len_rho()), 
      rho_sum(p,1,0.0), rho_ctr(p,1,0.0), rho_sse(p,p,0.0), mean(p,1,0.0), 
      mode(p,1,0.0), sscp(p,p,0.0), foc(p,1,0.0), n_tot(0), high(-REAL_MAX)
    { };
  bool set_asymptotics(const scl::realmat& sim); 
    // Aggregates if called repeatedly.
  void get_asymptotics(scl::realmat& rho_hat, scl::realmat& V_hat, 
    INTEGER& n);
  void get_asymptotics(scl::realmat& rho_mean, scl::realmat& rho_mode, 
    REAL& post_high, scl::realmat& I, scl::realmat& invJ, 
    scl::realmat& foc_hat, INTEGER& reps);
};

}

#endif
