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

Copyright (C) 2018

A. Ronald Gallant
Post Office Box 659
Raleigh 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 "lprior.h"
#include "slpeqns.h"

using namespace scl;
using namespace std;

int main(int argc, char** argp, char** envp)
{
  const INTEGER L = 1;
  const INTEGER M = 2;
  const INTEGER mrs_pos = 1;
  const INTEGER cf_pos = 2;
  const INTEGER years = 30;
  const INTEGER n_stats = 10;
  const INTEGER reps = 1000;
  const REAL lambda = 1.0;
  const REAL weight = 1.0e-2;
  //const REAL weight = 0.0;

  string filename;
  ifstream fin;

  realmat rho_hat, var_hat;
  
  if(vecread("dcf.rho_mode.dat", rho_hat)==0) error("Error, vecread failed");
  if(vecread("dcf.V_hat_hess.dat", var_hat)==0) error("Error, vecread failed");

  realmat R_hat = var_hat;
  if(factor(R_hat)!=0) error("Error, factor failed");

  realmat rho = rho_hat;

  realmat b0, B, S, R, theta;
  frac_rho(M, L, rho, b0, B, S, R, theta);

  realmat E;
  INTEGER ier;
  REAL maxlam = eigen(B,E,ier);

  cout << '\n';
  //cout << "rank = " << rank << '\n';
  cout << "ier = " << ier << '\n';
  cout << "maxlam = " << maxlam << '\n';
  cout << '\n';
  cout << "b0 = " << b0 << '\n';
  cout << "B = " << B << '\n';
  cout << "R = " << R << '\n';
  cout << "S = " << S << '\n';
  cout << "theta = " << theta << '\n';

  REAL beta = theta[1];
  REAL gamma = theta[2];

  intvec idx(n_stats);
  realmat stats(n_stats,1);
  realmat target(n_stats,1);

  target[ 1] = log(beta)/gamma;
  target[ 2] = 1.0/gamma;
  target[ 3] = 0.01;
  target[ 4] = 0.02;
  target[ 5] = 0.03;
  target[ 6] = 0.02;
  target[ 7] = 0.07;
  target[ 8] = 0.14;
  target[ 9] = 0.07;
  target[10] = 0.04;

  idx[ 1] = 0;
  idx[ 2] = 0;
  idx[ 3] = 1;
  idx[ 4] = 1;
  idx[ 5] = 1;
  idx[ 6] = 0;
  idx[ 7] = 1;
  idx[ 8] = 1;
  idx[ 9] = 1;
  idx[10] = 0;

  denval lp=lprior(M,L,mrs_pos,cf_pos,years,rho,lambda,target,idx,stats);

  cout << '\n';
  cout << "stat        idx    value    target" << '\n';
  cout << "lcg_mean     "  << idx[ 1]
    << fmt('f',10,5,stats[ 1]) << fmt('f',10,5,target[ 1]) << '\n';
  cout << "lcg_sdev     "  << idx[ 2]
    << fmt('f',10,5,stats[ 2]) << fmt('f',10,5,target[ 2]) << '\n';
  cout << "yld_mean_01  "  << idx[ 3]
    << fmt('f',10,5,stats[ 3]) << fmt('f',10,5,target[ 3]) << '\n';
  cout << "yld_sdev_01  "  << idx[ 4]
    << fmt('f',10,5,stats[ 4]) << fmt('f',10,5,target[ 4]) << '\n';
  cout << "yld_mean_yr  "  << idx[ 5]
    << fmt('f',10,5,stats[ 5]) << fmt('f',10,5,target[ 5]) << '\n';
  cout << "yld_sdev_yr  "  << idx[ 6]
    << fmt('f',10,5,stats[ 6]) << fmt('f',10,5,target[ 6]) << '\n';
  cout << "stk_mean_01  "  << idx[ 7]
    << fmt('f',10,5,stats[ 7]) << fmt('f',10,5,target[ 7]) << '\n';
  cout << "stk_sdev_01  "  << idx[ 8]
    << fmt('f',10,5,stats[ 8]) << fmt('f',10,5,target[ 8]) << '\n';
  cout << "stk_mean_yr  "  << idx[ 9]
    << fmt('f',10,5,stats[ 9]) << fmt('f',10,5,target[ 9]) << '\n';
  cout << "stk_sdev_yr  "  << idx[10]
    << fmt('f',10,5,stats[10]) << fmt('f',10,5,target[10]) << '\n';

  cout << boolalpha;
  cout << "lp = (" << lp.positive << ", " << lp.log_den << ')' << '\n';

  INTEGER p = rho_hat.size();
  INTEGER pless2 = p - 2;

  INT_32BIT seed = 8903735;
  realmat z_hat(p,1);

  REAL max_logprior = lp.log_den;
  realmat max_rho = rho;
  realmat max_stats = stats;
  denval max_lp = lp;
  INTEGER failcount = 0;

  for (INTEGER r=1; r<=reps; ++r) {

    for (INTEGER i=1; i<=p; ++i) z_hat[i] = unsk(seed);

    rho = rho_hat + R_hat*z_hat;

    rho[p-1] = beta;
    rho[p] = gamma;

    lp = lprior(M,L,mrs_pos,cf_pos,years,rho,lambda,target,idx,stats);

    if (lp.positive) {
      if (lp.log_den > max_logprior) {
        max_logprior = lp.log_den;
        max_rho = rho;
        max_lp = lp;
        max_stats = stats;
      }
    }
    else {
      ++failcount;
    }
  }

  rho = max_rho;
  stats = max_stats;
  realmat parms = rho;
  lp = max_lp;

  cout << '\n';
  cout << "stat        idx    value    target" << '\n';
  cout << "lcg_mean     "  << idx[ 1]
    << fmt('f',10,5,stats[ 1]) << fmt('f',10,5,target[ 1]) << '\n';
  cout << "lcg_sdev     "  << idx[ 2]
    << fmt('f',10,5,stats[ 2]) << fmt('f',10,5,target[ 2]) << '\n';
  cout << "yld_mean_01  "  << idx[ 3]
    << fmt('f',10,5,stats[ 3]) << fmt('f',10,5,target[ 3]) << '\n';
  cout << "yld_sdev_01  "  << idx[ 4]
    << fmt('f',10,5,stats[ 4]) << fmt('f',10,5,target[ 4]) << '\n';
  cout << "yld_mean_yr  "  << idx[ 5]
    << fmt('f',10,5,stats[ 5]) << fmt('f',10,5,target[ 5]) << '\n';
  cout << "yld_sdev_yr  "  << idx[ 6]
    << fmt('f',10,5,stats[ 6]) << fmt('f',10,5,target[ 6]) << '\n';
  cout << "stk_mean_01  "  << idx[ 7]
    << fmt('f',10,5,stats[ 7]) << fmt('f',10,5,target[ 7]) << '\n';
  cout << "stk_sdev_01  "  << idx[ 8]
    << fmt('f',10,5,stats[ 8]) << fmt('f',10,5,target[ 8]) << '\n';
  cout << "stk_mean_yr  "  << idx[ 9]
    << fmt('f',10,5,stats[ 9]) << fmt('f',10,5,target[ 9]) << '\n';
  cout << "stk_sdev_yr  "  << idx[10]
    << fmt('f',10,5,stats[10]) << fmt('f',10,5,target[10]) << '\n';
  cout << "lp = (" << lp.positive << ", " << lp.log_den << ')' << '\n';
  cout << "failcount = " << failcount << '\n';

  slpeqns lpe(M,L,mrs_pos,cf_pos,years,lambda,weight,theta,target,idx);

  nlopt lpopt(lpe);
  //lpopt.set_lower_bound(0.0);
  //lpopt.set_output(true,&cout,true);
  lpopt.set_output(true);
  lpopt.set_warning_messages(true);
  lpopt.set_iter_limit(100);
  //lpopt.set_solution_tolerance(0.0001);

  realmat x_start(pless2,1); 
  for (INTEGER i=1; i<=pless2; ++i) x_start[i] = parms[i]; 
  realmat x_stop;

  //cout << '\n';
  //cout << "x_start = " << x_start << '\n';

  lpopt.minimize(x_start,x_stop);

  cout << '\n';
  cout << "x_stop = " << x_stop << '\n';
  cout << "theta = " << theta << '\n';

  for (INTEGER i=1; i<=pless2; ++i) parms[i] = x_stop[i]; 

  lp = lprior(M,L,mrs_pos,cf_pos,years,parms,lambda,target,idx,stats);

  cout << '\n';
  cout << "lp = (" << lp.positive << ", " << lp.log_den << ')' << '\n';

  cout << '\n';
  cout << "stat        idx    value    target" << '\n';
  cout << "lcg_mean     "  << idx[ 1]
    << fmt('f',10,5,stats[ 1]) << fmt('f',10,5,target[ 1]) << '\n';
  cout << "lcg_sdev     "  << idx[ 2]
    << fmt('f',10,5,stats[ 2]) << fmt('f',10,5,target[ 2]) << '\n';
  cout << "yld_mean_01  "  << idx[ 3]
    << fmt('f',10,5,stats[ 3]) << fmt('f',10,5,target[ 3]) << '\n';
  cout << "yld_sdev_01  "  << idx[ 4]
    << fmt('f',10,5,stats[ 4]) << fmt('f',10,5,target[ 4]) << '\n';
  cout << "yld_mean_yr  "  << idx[ 5]
    << fmt('f',10,5,stats[ 5]) << fmt('f',10,5,target[ 5]) << '\n';
  cout << "yld_sdev_yr  "  << idx[ 6]
    << fmt('f',10,5,stats[ 6]) << fmt('f',10,5,target[ 6]) << '\n';
  cout << "stk_mean_01  "  << idx[ 7]
    << fmt('f',10,5,stats[ 7]) << fmt('f',10,5,target[ 7]) << '\n';
  cout << "stk_sdev_01  "  << idx[ 8]
    << fmt('f',10,5,stats[ 8]) << fmt('f',10,5,target[ 8]) << '\n';
  cout << "stk_mean_yr  "  << idx[ 9]
    << fmt('f',10,5,stats[ 9]) << fmt('f',10,5,target[ 9]) << '\n';
  cout << "stk_sdev_yr  "  << idx[10]
    << fmt('f',10,5,stats[10]) << fmt('f',10,5,target[10]) << '\n';

  realmat f,F;
  bool success = lpe.get_f(x_stop,f);

  cout << '\n';
  cout << "beginning maxlam = " << maxlam << '\n';
  frac_tau(M, L, x_stop, b0, B, S, R);
  maxlam = eigen(B,E,ier);
  cout << "ending maxlam    = " << maxlam << '\n';

  cout << '\n';
  cout << "success = " << success << '\n';
  cout << "f = " << f << '\n';
  success = lpe.get_F(x_stop,f,F);
  cout << "success = " << success << '\n';
  cout << "F = " << F << '\n';

  if(vecwrite("dcf.rho_mode.old",rho_hat)==0) error("Error, vecwrite failed");
  if(vecwrite("dcf.V_hat_hess.old",var_hat)==0) error("Error, vecwrite failed");

  rho_hat = parms;
  realmat H_matrix = lpopt.get_H_matrix();
  fill(var_hat);
  
  for (INTEGER j=1; j<=pless2; ++j) {
    for (INTEGER i=1; i<=pless2; ++i) {
      var_hat(i,j) = H_matrix(i,j);
    }
  }
  var_hat(p-1,p-1) = 1.0;
  var_hat(p,p) = 1.0;

  if(vecwrite("dcf.rho_mode.new",rho_hat)==0) error("Error, vecwrite failed");
  if(vecwrite("dcf.V_hat_hess.new",var_hat)==0) error("Error, vecwrite failed");

  return 0;
}

