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

Copyright (C) 2014

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 "libmle.h"
#include "mle.h"
#include "mleusr.h"

using namespace scl;
using namespace libmle;
using namespace mle;
using namespace std;

mle::binomial_usrmod::binomial_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), rho(1,1,0.5), lrho(1), lstats(1)
{
  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_func is set wrong in parmfile");
  }

  vector<string>::const_iterator mod_alvec_ptr = mod_alvec.begin();
  if (*mod_alvec_ptr != "#begin additional lines") {
    cerr << "Line read " << *mod_alvec_ptr << '\n';
    error("Error, binomial_usrmod, bad mod_alvec");
  }

  prior_mean = atof((++mod_alvec_ptr)->substr(0,12).c_str());
  prior_sdev = atof((++mod_alvec_ptr)->substr(0,12).c_str());

  if (*(++mod_alvec_ptr) != "#end additional lines") {
    cerr << "Line read " << *mod_alvec_ptr << '\n';
    error("Error, binomial_usrmod, bad mod_alvec");
  }
}

den_val mle::binomial_usrmod::likelihood() 
{
  INTEGER n = data.get_cols();

  REAL sum = 0.0;
  for (INTEGER i=1; i<=n; ++i) {
    sum += data[i]*log(rho[1]) + (1.0-data[i])*log(1.0-rho[1]);
  }

  if (IsFinite(sum)) {
    return den_val(true,sum);
  }
  else {
    return den_val(false,-REAL_MAX);;
  }
}

bool mle::binomial_usrmod::support(const realmat& parm) 
{
  if (parm[1] <= 0.0) return false;
  if (parm[1] >= 1.0) return false;
  return true; 
}

bool mle::binomial_usrmod::get_stats(scl::realmat& stats) 
{ 
  stats.resize(1,1,0.0);
  return true;
}     

den_val mle::binomial_usrmod::prior(const realmat& parm, const realmat& stats)
{
  const REAL minus_log_root_two_pi = -9.1893853320467278e-01;

  REAL z = (parm[1] - prior_mean)/prior_sdev;
  REAL e = minus_log_root_two_pi - log(prior_sdev) - 0.5*pow(z,2);

  return den_val(true,e);
}

