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

Copyright (C) 2015.

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

using namespace scl;
using namespace std;

namespace mle {

  sims_usrmod::sims_usrmod
    (const realmat& dat, INTEGER len_mod_parm, INTEGER len_mod_func,
     const vector<string>& mod_pfvec, const vector<string>& mod_alvec,
     ostream& detail)
  : data(dat), theta(n_parms,1), prior_mean(n_parms,1), prior_sdev(n_parms,1) 
  {
    detail << starbox("/sims_usrmod constructor now controls print//");
    detail << '\n';

    if (len_mod_parm != n_parms) error("Error, sims_usrmod, bad len_mod_parm");
    if (len_mod_func != n_funcs) error("Error, sims_usrmod, bad len_mod_func");

    if (mod_alvec.size() != 5) {
      error("Error, sims_usrmod, prior_scale or mnt_eq not in parmfile");
    }
    prior_mean[1] = atof(mod_alvec[1].substr(0,12).c_str());
    prior_sdev[1] = atof(mod_alvec[2].substr(0,12).c_str());
    mnt_eq = atoi(mod_alvec[3].substr(0,12).c_str());

    detail << boolalpha << '\n';
    detail << "\t prior_mean  = " << prior_mean[1] << '\n';
    detail << "\t prior_sdev  = " << prior_sdev[1] << '\n';
    detail << "\t mnt_eq = " << mnt_eq << '\n';
    detail << '\n';

    detail << starbox("/mle now controls print//");
    detail.flush();
  }

  bool sims_usrmod::get_stats(realmat& stats)
  {
    stats.resize(n_funcs,1,1.0);
    return true;
  }

  den_val sims_usrmod::likelihood()
  {
    //const REAL minus_log_root_two_pi = -9.1893853320467278e-01;

    INTEGER n = data.ncol();
    realmat gam(n,1);
    realmat score(n,1);
    realmat logp(n,1);
    for (INTEGER i=1; i<=n; ++i) {
      REAL x = data[i];
      gam[i] = pow(x,theta[1]);
      logp[i] = log(theta[1]) + (theta[1]-1.0)*log(x) - gam[i];
      score[i] = 1.0/theta[1] + log(x) - log(x)*gam[i];
    }

    REAL log_likelihood = 0.0;
    REAL mean;
    REAL var;
    REAL sum;
    REAL adj;

    switch (mnt_eq) {

      case 1 : 
        for (INTEGER i=1; i<=n; ++i) {
          log_likelihood += logp[i];
        }
        break;

      case 2 : 
        sum = 0.0;
        for (INTEGER i=1; i<=n; ++i) {
          sum += gam[i];
        }
        log_likelihood = REAL(n-1)*log(sum) - sum;
        break;

      case 3 :
        mean = 0.0;
        for (INTEGER i=1; i<=n; ++i) {
          mean += gam[i] - 1.0;
        }
        mean /= REAL(n);
        var = 0.0;
        for (INTEGER i=1; i<=n; ++i) {
          var += pow(gam[i] - 1.0 - mean, 2);
        }
        var /= REAL(n-1);
        log_likelihood = -0.5*REAL(n)*pow(mean,2)/var;
        break;

      case 4 :
        mean = 0.0;
        for (INTEGER i=1; i<=n; ++i) {
          mean += score[i];
        }
        mean /= REAL(n);
        var = 0.0;
        for (INTEGER i=1; i<=n; ++i) {
          var += pow(score[i] - mean, 2);
        }
        var /= REAL(n-1);
        log_likelihood = -0.5*REAL(n)*pow(mean,2)/var;
        break;

      case 5 : 
        sum = 0.0;
        for (INTEGER i=1; i<=n; ++i) {
          sum += gam[i];
        }
        adj = log(theta[1]) + log(sum)*(1.0-1.0/theta[1]);
        log_likelihood = adj + REAL(n-1)*log(sum) - sum;
        break;

      default : error("Error, sims, invalid mnt_eq");
        break;
    }

    return denval(true, log_likelihood);
  }

  bool sims_usrmod::support(const realmat& theta)
  {
    return (theta[1] > 1.0);
  }

  den_val sims_usrmod::prior(const realmat& theta, const realmat& stats)
  {
    const REAL minus_log_root_two_pi = -9.1893853320467278e-01;
    REAL z = (theta[1] - prior_mean[1])/prior_sdev[1];
    REAL e = minus_log_root_two_pi - log(prior_sdev[1]) - 0.5*pow(z,2);
    return den_val(true,e);
  }

  void sims_usrmod::write_usrvar(const char* filename)
  {
    /*
    string stem = filename;

    string::size_type suffix = stem.find(string("dat"));
    if (suffix != std::string::npos) {
      stem = stem.substr(0, suffix);
    }
    
    ofstream fout;
    fout.open(filename);
    if (!fout) error("Error, usrvar, could not open " + string(filename));

    fout.clear(); fout.close();
    */
  }

}
