#ifndef __FILE_LIBGSM_H_SEEN__
#define __FILE_LIBGSM_H_SEEN__ 

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

Copyright (C) 2005, 2006, 2008, 2011, 2013.

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_base.h"

namespace libgsm {

typedef std::map<scl::intvec,REAL,scl::intvec_cmp> pdf_type;
typedef std::map<scl::intvec,REAL,scl::intvec_cmp>::iterator pdf_itr;

class grid_group_move : public proposal_base {
private:
  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);
  void draw(INT_32BIT& seed, const scl::realmat& parm_old, 
    scl::realmat& parm_new); 
  scl::den_val operator()(const scl::realmat& parm_old, 
    const scl::realmat& parm_new);
  INTEGER len_parm() { 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 lparm;
  INTEGER lpd;
  std::vector<prop_aux> pa;
public:
  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_parm() { return lparm; }
  bool transition_is_symmetric() { return true; }
};

class implied_map {
private:
  std::map<scl::realmat,map_val,scl::realmat_cmp> past_values;
  bool draw_from_posterior;
  mutable INTEGER hit;
  mutable INTEGER miss;
public:
  implied_map() : past_values(),draw_from_posterior(true),hit(0),miss(0) {}
  bool read_map(const char * filename);
  bool write_map(const char * filename) const;
  bool get_val(const scl::realmat& theta, map_val& mv) const;
  void set_val(const scl::realmat& theta, const map_val& mv);
  void set_val(const sci_val& sv);
  void set_draw_from_posterior(bool from_posterior);
  void update_map(const implied_map& im);
  void clear() { past_values.erase(past_values.begin(), past_values.end()); } 
  bool empty() const { return past_values.empty(); }
  typedef std::map<scl::realmat,map_val,scl::realmat_cmp>::const_iterator
    im_c_itr;
  im_c_itr get_begin() const { return past_values.begin(); }
  im_c_itr get_end() const { return past_values.end(); }
  typedef std::map<scl::realmat,map_val,scl::realmat_cmp>::size_type 
    im_st;
  im_st size() const { return past_values.size(); }
  sci_val get_mode() const;
  REAL get_hit_rate() const { return REAL(hit)/REAL(hit+miss); }
};  

class implied_map_buffer {   // For passing an implied_map via MPI
private:
  scl::realmat buffer;
  scl::realmat bufcol;
public:
  implied_map_buffer() : buffer(), bufcol() { }
  implied_map_buffer(const implied_map& im);
  INTEGER size() const {return buffer.size();}
  INTEGER get_rows() const {return buffer.get_rows();}
  INTEGER get_cols() const {return buffer.get_cols();}
  void resize(INTEGER r, INTEGER c);
  const scl::realmat& get_col(INTEGER j);
  const scl::realmat& get_buf() const {return buffer;}
  REAL* get_ptr() {return buffer.get_x();}
};

class assess_prior {
private:
  const INTEGER l;            // Length of stat_parm
  INTEGER n;                  // Number of times increment_sigma called
  scl::realmat sum_squares;
  scl::realmat sigma;
  REAL kappa;
  bool imap_set;
  bool kappa_set;
  bool sigma_set;
  struct man_val {
    scl::realmat stat_parm;
    scl::den_val sci_prior;
  };
  std::vector<man_val> manifold;
public:
          assess_prior(const INTEGER len_stat_parm);
  bool    read_sigma(const char * filename);
  bool    write_sigma(const char * filename) const;
  scl::den_val operator()(const scl::realmat& eta) const;
  void    increment_sigma(const scl::realmat& eta_residual);
  void    update_sigma(const scl::realmat& sig, INTEGER nobs);
  void    set_implied_map(const implied_map& im);
  void    set_kappa(const REAL kap);
  INTEGER get_n() const;
  const scl::realmat& get_sigma() const;
};

class stat_mcmc {
private:
  stat_mod_base& stat_mod;
  proposal_base& stat_prop;
  assess_prior&  aprior;
  bool     draw_from_posterior;
  INTEGER  len_chain;
  const scl::realmat* data_ptr;
  bool     use_stat_mod_prior; 
  bool     use_assess_prior;
  stat_val mode;
public:
  stat_mcmc(stat_mod_base& stm, proposal_base& stp, assess_prior& ap)
    : stat_mod(stm), stat_prop(stp), aprior(ap), 
      draw_from_posterior(true), len_chain(500), data_ptr(),
      use_stat_mod_prior(true), use_assess_prior(true), mode() { }
  scl::realmat  draw(INT_32BIT& seed, scl::realmat& stat_parm_start,
  scl::realmat& stat_parm_chain, scl::realmat& stat_func_chain, 
  scl::realmat& stat_prior_chain, scl::realmat& stat_logl_chain,
  scl::realmat& assess_prior_chain);
  void set_len_chain(INTEGER N);
  void set_use_stat_mod_prior(bool use) { use_stat_mod_prior = use; }
  void set_use_assess_prior(bool use) { use_assess_prior = use; }
  void set_data_ptr(const scl::realmat* dat_ptr) { data_ptr = dat_ptr; }
  void set_draw_from_posterior(bool from_posterior);
  stat_val get_mode() { return mode; }
};

class sci_mcmc {
private:
  sci_mod_base& sci_mod;
  proposal_base& sci_prop;
  stat_mod_base& stat_mod;
  proposal_base& stat_prop;
  assess_prior& aprior;
  const scl::realmat* data_ptr;
  implied_map& imap;
  bool    draw_from_posterior;
  INTEGER len_chain;
  INTEGER len_sub_chain;
  INTEGER num_sub_chains;
  INTEGER num_polish_iter;
  REAL    polish_toler;
  INTEGER num_sub_chain_to_save;
  INTEGER num_sub_chain_saved;
  sci_val mode;
  bool    analytic_mle;
  INTEGER make_sci_val( // Used when analytic_mle == true
    scl::realmat& sci_parm_start, scl::realmat& stat_parm_start,
    sci_val& sv, std::string& errmsg);
  INTEGER make_sci_val(INT_32BIT& sub_chain_seed,
    scl::realmat& sci_parm_start, scl::realmat& stat_parm_start,
    sci_val& sv, std::string& errmsg,
    scl::realmat& sci_sub_parm_chain, scl::realmat& sci_sub_logl_chain,
    scl::realmat& sci_sub_rej_chain);
  void remake_sci_val_pair(INT_32BIT& sub_chain_seed,
    sci_val& first, sci_val& second,
    scl::realmat& sci_sub_parm_chain, scl::realmat& sci_sub_logl_chain,
    scl::realmat& sci_sub_rej_chain);
public:
  sci_mcmc(sci_mod_base& scm, proposal_base& scp, 
    stat_mod_base& stm, proposal_base& stp, assess_prior& ap,
    const scl::realmat* dat_ptr, implied_map& im) 
    : sci_mod(scm), sci_prop(scp), 
      stat_mod(stm), stat_prop(stp), aprior(ap), 
      data_ptr(dat_ptr), imap(im), 
      draw_from_posterior(true),
      len_chain(1000), len_sub_chain(100), num_sub_chains(0), 
      num_polish_iter(0), polish_toler(1.0e-10),
      num_sub_chain_to_save(25), num_sub_chain_saved(0),
      mode(), analytic_mle(false) { }  
  scl::realmat draw(INT_32BIT& seed, 
    scl::realmat& sci_parm_start, scl::realmat& stat_parm_start,
    scl::realmat& sci_parm_chain, scl::realmat& sci_func_chain, 
    scl::realmat& sci_prior_chain, 
    scl::realmat& sci_sub_parm_chain, scl::realmat& sci_sub_logl_chain, 
    scl::realmat& sci_sub_rej_chain,
    scl::realmat& stat_parm_chain, scl::realmat& stat_func_chain, 
    scl::realmat& stat_prior_chain, scl::realmat& stat_logl_chain,
    scl::realmat& stat_sub_logl_chain);
  void set_len_chain(INTEGER N);
  void set_len_sub_chain(INTEGER N);
  void set_num_sub_chains(INTEGER N);
  void set_num_polish_iter(INTEGER N);
  void set_polish_toler(REAL tol);
  void set_draw_from_posterior(bool from_posterior);
  void set_analytic_mle(INTEGER am);
  sci_val get_mode() { return mode; }
};
}

#endif

