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

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 "lpeqns.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;

  string filename;
  ifstream fin;

  filename = "../data/DCFdataL.txt";
  fin.open(filename.c_str());
  if (!fin) error("Error, cannot open " + filename);

  const INTEGER data_row = 85;  
  const INTEGER data_col = 4;

  realmat data(data_row,data_col);
  
  for (INTEGER t=1; t<=data_row; ++t) {
    REAL year, log_mrs, log_gdp_prof, del_log_prof;
    fin >> year >> log_mrs >> log_gdp_prof >> del_log_prof;
    data(t,1) = year;
    data(t,2) = log_mrs;
    data(t,3) = log_gdp_prof;
    data(t,4) = del_log_prof;
  }

  intvec vdx(2);
  vdx[1] = 2;
  vdx[2] = 4;

  if (M != vdx.size()) error("Error, bad dimensions");

  realmat b0, B, S, C, R;
  realmat theta(2,1);
  REAL beta;
  REAL gamma;

  //INTEGER rank = varcoef(data, vdx, L, b0, B, S, C);

  realmat rho;
  vecread("dcf.rho_mode.dat",rho);

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

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

  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 << '\n';
  cout << "beta = " << beta << '\n';
  cout << "gamma = " << gamma << '\n';

  //const REAL pi = 4.0*atan(1.0);
  //const REAL sqrt2 = sqrt(2.0);

  const REAL pi = M_PI;
  const REAL sqrt2 = M_SQRT2;
  
  realmat IB = -B;
  for (INTEGER i=1; i<=M; ++i) IB(i,i) += 1.0;
  realmat IBinv;
  IBinv = inv(IB);

//  realmat R=S;
//  if (factor(R)!=0) error("Error, factor failed");

//  S = R*T(R);

  cout << '\n';
  //cout << "R = " << R << '\n';
  cout << "IB = " << IB << '\n';
  cout << "IBinv = " << IBinv << '\n';

  realmat y0(M,1);
  for (INTEGER i=1; i<=M; ++i) y0[i] = data(1,vdx[i]);
  realmat ecf,pvcf,pv1,cecf,dcf,yld;

  /*
  dcfyld(b0,B,S,y0,mrs_pos,cf_pos,years,ecf,pvcf,pv1,cecf,dcf,yld); 

  cout << '\n';
  cout << "years = " << years << '\n';
  cout << '\n';
  cout << "y0 = " << y0 << '\n';
  cout << "ecf = " << ecf << '\n';
  cout << "pvcf = " << pvcf << '\n';
  cout << "pv1 = " << pv1 << '\n';
  cout << "cecf = " << cecf << '\n';
  cout << "dcf = " << dcf << '\n';
  cout << "yld = " << yld << '\n';
  */

  const INTEGER quad_n = 5;

  realmat x, w;

  if (hquad(quad_n,x,w)!=0) error("Error, hquad failed");

  cout << '\n';
  cout << "hquad x = " << x << '\n';
  cout << "hquad w = " << w << '\n';

  realmat z(M,1);

  REAL lcg_mean = 0.0;
  REAL lcg_sdev = 0.0;

  REAL yld_mean_01 = 0.0;
  REAL yld_sdev_01 = 0.0;
  REAL yld_mean_yr = 0.0;
  REAL yld_sdev_yr = 0.0;

  REAL stk_mean_01 = 0.0;
  REAL stk_sdev_01 = 0.0;
  REAL stk_mean_yr = 0.0;
  REAL stk_sdev_yr = 0.0;

  for (INTEGER i=1; i<=quad_n; ++i) {
    for (INTEGER j=1; j<=quad_n; ++j) {

      z[1] = sqrt2*x[i];
      z[2] = sqrt2*x[j];

      y0 = IBinv*(b0 + R*z);

      REAL weight = w[i]*w[j]/pi;

      dcfyld(b0,B,S,y0,mrs_pos,cf_pos,years,ecf,pvcf,pv1,cecf,dcf,yld);

      lcg_mean += weight*(log(beta)-z[1])/gamma;
      lcg_sdev += weight*pow((log(beta)-z[1])/gamma,2);

      yld_mean_01 += weight*yld[1];
      yld_mean_yr += weight*yld[years];
      yld_sdev_01 += weight*pow(yld[1],2);
      yld_sdev_yr += weight*pow(yld[years],2);

      REAL ret01 = log(dcf[years]+ecf[1]) -log(dcf[years]);
      REAL retyr = ( log(dcf[years]+cecf[years]) -log(dcf[years]))/years;

      stk_mean_01 += weight*ret01;
      stk_mean_yr += weight*retyr;
      stk_sdev_01 += weight*pow(ret01,2);
      stk_sdev_yr += weight*pow(retyr,2);
    }
  }

  lcg_sdev = sqrt(lcg_sdev - pow(lcg_mean,2));

  yld_sdev_01 = sqrt(yld_sdev_01 - pow(yld_mean_01,2));
  yld_sdev_yr = sqrt(yld_sdev_yr - pow(yld_mean_yr,2));

  stk_sdev_01 = sqrt(stk_sdev_01 - pow(stk_mean_01,2));
  stk_sdev_yr = sqrt(stk_sdev_yr - pow(stk_mean_yr,2));

  realmat stats(10,1);
  realmat target(10,1);

  stats[ 1] = lcg_mean;
  stats[ 2] = lcg_sdev;
  stats[ 3] = yld_mean_01;
  stats[ 4] = yld_sdev_01;
  stats[ 5] = yld_mean_yr;
  stats[ 6] = yld_sdev_yr;
  stats[ 7] = stk_mean_01;
  stats[ 8] = stk_sdev_01;
  stats[ 9] = stk_mean_yr;
  stats[10] = stk_sdev_yr;

  target[ 1] = 0.00;
  target[ 2] = 0.10;
  target[ 3] = 0.01;
  target[ 4] = 0.02;
  target[ 5] = 0.02;
  target[ 6] = 0.05;
  target[ 7] = 0.07;
  target[ 8] = 0.14;
  target[ 9] = 0.05;
  target[10] = 0.05;

  const INTEGER lidx = 10;
  intvec idx(lidx);
  idx[ 1] = 1;
  idx[ 2] = 1;
  idx[ 3] = 1;
  idx[ 4] = 1;
  idx[ 5] = 0;
  idx[ 6] = 0;
  idx[ 7] = 1;
  idx[ 8] = 1;
  idx[ 9] = 0;
  idx[10] = 0;

  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';

  REAL ssd = 0.0;
  for (INTEGER i=1; i<=lidx; ++i) {
    if (idx[i] == 1) ssd += pow(stats[i]-target[i],2);
  }

  const REAL minus_log_root_two_pi = -log(sqrt(2.0*pi));


  REAL lambda = 1.0;
  //REAL logprior = -lambda*ssd;
  REAL logprior = minus_log_root_two_pi + 0.5*log(lambda) - 0.5*lambda*ssd;

  REAL plambda = exp(-lambda*ssd);

  cout << '\n';
  cout << "ssd = " << ssd << '\n';
  cout << "lambda = " << lambda << '\n';
  cout << "plambda = " << plambda << '\n';
  cout << "logprior = " << logprior << '\n';

  realmat parms(M+M*M+(M*M+M)/2+2,1);

  parms[ 1] = b0[1];
  parms[ 2] = b0[2];
  parms[ 3] = B[1];
  parms[ 4] = B[2];
  parms[ 5] = B[3];
  parms[ 6] = B[4];
  parms[ 7] = R[1];
  parms[ 8] = R[3];
  parms[ 9] = R[4];
  parms[10] = beta;
  parms[11] = gamma;

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

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

  lpeqns lpe(M,L,mrs_pos,cf_pos,years,lambda,target,idx);

  realmat f,F;
  bool success = lpe.get_f(parms,f);
  cout << '\n';
  cout << "success = " << success << '\n';
  cout << "f = " << f << '\n';
  success = lpe.get_F(parms,f,F);
  cout << "success = " << success << '\n';
  cout << "F = " << F << '\n';

  nlopt lpopt(lpe);
  lpopt.set_output(true);
  lpopt.set_warning_messages(true);
  lpopt.set_iter_limit(100);
  //lpopt.set_solution_tolerance(0.0001);

  realmat x_start = parms;
  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';

  parms = x_stop;

  b0[1] = parms[1];
  b0[2] = parms[2];
  B[1]  = parms[3];
  B[2]  = parms[4];
  B[3]  = parms[5]; 
  B[4]  = parms[6];
  R[1]  = parms[7];
  R[2]  = 0.0;
  R[3]  = parms[8];
  R[4]  = parms[9];
  theta[1] = beta = parms[10];
  theta[2] = gamma = parms[11];

  S = R*T(R);

  lcg_mean = 0.0;
  lcg_sdev = 0.0;

  yld_mean_01 = 0.0;
  yld_sdev_01 = 0.0;
  yld_mean_yr = 0.0;
  yld_sdev_yr = 0.0;

  stk_mean_01 = 0.0;
  stk_sdev_01 = 0.0;
  stk_mean_yr = 0.0;
  stk_sdev_yr = 0.0;


  for (INTEGER i=1; i<=quad_n; ++i) {
    for (INTEGER j=1; j<=quad_n; ++j) {

      z[1] = sqrt2*x[i];
      z[2] = sqrt2*x[j];

      y0 = IBinv*(b0 + R*z);

      REAL weight = w[i]*w[j]/pi;

      dcfyld(b0,B,S,y0,mrs_pos,cf_pos,years,ecf,pvcf,pv1,cecf,dcf,yld);

      lcg_mean += weight*(log(beta)- z[1])/gamma;
      lcg_sdev += weight*pow((log(beta)- z[1])/gamma,2);

      yld_mean_01 += weight*yld[1];
      yld_mean_yr += weight*yld[years];
      yld_sdev_01 += weight*pow(yld[1],2);
      yld_sdev_yr += weight*pow(yld[years],2);

      REAL ret01 = log(dcf[years]+ecf[1]) -log(dcf[years]);
      REAL retyr = ( log(dcf[years]+cecf[years]) -log(dcf[years]))/years;

      stk_mean_01 += weight*ret01;
      stk_mean_yr += weight*retyr;
      stk_sdev_01 += weight*pow(ret01,2);
      stk_sdev_yr += weight*pow(retyr,2);
    }
  }

  lcg_sdev = sqrt(lcg_sdev - pow(lcg_mean,2));

  yld_sdev_01 = sqrt(yld_sdev_01 - pow(yld_mean_01,2));
  yld_sdev_yr = sqrt(yld_sdev_yr - pow(yld_mean_yr,2));

  stk_sdev_01 = sqrt(stk_sdev_01 - pow(stk_mean_01,2));
  stk_sdev_yr = sqrt(stk_sdev_yr - pow(stk_mean_yr,2));

  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';

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

  return 0;
}

