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

Copyright (C) 2005, 2006, 2008.

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"
using namespace std;  
using namespace scl;
using namespace libgsm;

bool libgsm::implied_map::read_map(const char* filename)
{
  ifstream is(filename);
  if (!is.good()) return false;
  sci_val sv;
  while (is >> sv) this->set_val(sv);
  return true;
}

bool libgsm::implied_map::write_map(const char* filename) const
{
  ofstream os(filename);
  if (!os.good()) return os.good();
  for (im_c_itr itr = past_values.begin(); itr != past_values.end(); ++itr) {
    sci_val sv(itr->first, itr->second);
    os << sv;
  }
  return os.good();
}

sci_val libgsm::implied_map::get_mode() const
{
  sci_val mode;
  REAL pi_mode = -REAL_MAX;
  bool first_time = true;
  for (im_c_itr itr = past_values.begin(); itr != past_values.end(); ++itr) {
    sci_val sv(itr->first, itr->second);
    bool good = sv.sci_mod_support && sv.sci_mod_simulate 
                && sv.sci_mod_prior.positive && sv.stat_mod_logl.positive
                && sv.stat_mod_support && sv.stat_mod_simulate
                && sv.stat_mod_prior.positive;
    if (good) {
      REAL pi = sv.sci_mod_prior.log_den;
      pi += sv.stat_mod_prior.log_den;
      if (draw_from_posterior) {
        pi += sv.stat_mod_logl.log_den;
      }
      if (first_time) {
        mode = sv;
        pi_mode = pi;
        first_time = false;
      }
      else if (pi > pi_mode) {
         mode = sv;
         pi_mode = pi;
      }
    }
  }
  return mode;
}

bool libgsm::implied_map::get_val(const realmat& theta, map_val& mv) const
{
  im_c_itr itr = past_values.find(theta);
  if ( itr == past_values.end() ) { ++miss ; return false; }
  ++hit;
  mv = itr->second;
  return true;
}

void libgsm::implied_map::set_val(const realmat& theta, const map_val& mv)
{
  past_values[theta] = mv;
}

void libgsm::implied_map::set_val(const sci_val& sv)
{
  map_val mv(sv);
  past_values[sv.sci_mod_parm] = mv;
}

void libgsm::implied_map::set_draw_from_posterior(bool from_posterior)
{
  draw_from_posterior = from_posterior;
}

void libgsm::implied_map::update_map(const implied_map& im)
{
  if (im.empty()) return;

  if (this->empty()) {
    implied_map::im_c_itr itr;
    for (itr=im.get_begin(); itr!=im.get_end(); ++itr) {
      this->set_val(itr->first, itr->second);
    }
  }
  else {
    implied_map::im_c_itr itr;
    for (itr=im.get_begin(); itr!=im.get_end(); ++itr) {
      realmat theta = itr->first;
      map_val mv_new = itr->second;
      map_val mv_old;
      if (this->get_val(theta,mv_old)) {
        if (mv_new.stat_mod_logl.positive) {
          if (mv_old.stat_mod_logl.positive) {
            if (mv_old.stat_mod_parm.nrow() != mv_old.stat_mod_parm.nrow()) {
              error("Error, implied_map, update, size mismatch");
            }
            if (mv_old.sci_mod_func.nrow() != mv_old.sci_mod_func.nrow()) {
              error("Error, implied_map, update, size mismatch");
            }
            if (mv_new.stat_mod_sub_logl.log_den 
                 > mv_old.stat_mod_sub_logl.log_den) {
              this->set_val(theta, mv_new);
            }
          }
          else {
            this->set_val(theta, mv_new);
          }
        }
      }
      else {
        this->set_val(theta, mv_new);
      }
    }
  }
}

void libgsm::implied_map_buffer::resize(INTEGER r, INTEGER c)
{
  buffer.resize(r,c); 
  bufcol.resize(r,1);
}

const realmat& libgsm::implied_map_buffer::get_col(INTEGER j)
{
  INTEGER rows = buffer.get_rows();
  for (INTEGER i=1; i<=rows; ++i) bufcol[i] = buffer(i,j);
  return bufcol;
}

libgsm::implied_map_buffer::implied_map_buffer(const implied_map& im)
{
  if (im.empty()) {
    buffer.resize(0,0);
    bufcol.resize(0,0);
  }
  else {
    implied_map::im_c_itr im_ptr = im.get_begin();
    sci_val sv(im_ptr->first, im_ptr->second);

    realmat rm = pack_sci_val(sv);

    INTEGER cols = im.size();

    vector<realmat> buf_vec;
    buf_vec.reserve(cols);

    INTEGER maxrows = 0;
    for (im_ptr = im.get_begin(); im_ptr != im.get_end(); ++im_ptr) {
      sci_val sv(im_ptr->first, im_ptr->second);
      realmat rm = pack_sci_val(sv); 
      maxrows = maxrows < rm.get_rows() ? rm.get_rows() : maxrows;
      buf_vec.push_back(rm);
    }

    INTEGER rows = maxrows;
    buffer.resize(rows,cols);
    bufcol.resize(rows,1);

    INTEGER k = 0;
    vector<realmat>::const_iterator itr;
    for (itr = buf_vec.begin(); itr != buf_vec.end(); ++itr) {
      realmat rm = *itr;
      ++k;
      INTEGER r = rm.get_rows();
      for (INTEGER i=1; i<=r; ++i) {
        buffer(i,k) = rm[i];
      }
    }
  }
}
