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

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

  string filename;
  ifstream fin;

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

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

  realmat EXX;

  var1EXX(b0, B, S, EXX);

  cout << b0 << B << S << EXX;

  INT_32BIT seed = 34059;

  realmat y(M,M);
  realmat z(M,1);
  realmat ylag(M,1,0.0);

  for (INTEGER t=1; t<=10000; ++t) {
    for (INTEGER i=1; i<=M; ++i) z[i] = unsk(seed);
    realmat y = b0 + B*ylag + R*z;
    ylag = y;
  }

  realmat sum(M+1,M+1,0.0);
  realmat x(M+1,1);

  INTEGER N = 1000000;
  for (INTEGER t=1; t<=N; ++t) {
    for (INTEGER i=1; i<=M; ++i) z[i] = unsk(seed);
    realmat y = b0 + B*ylag + R*z;
    ylag = y;
    x[1] = 1.0;
    for (INTEGER i=1; i<=M; ++i) x[1+i] = y[i];
    sum += x*T(x);
  }

  cout << sum/N;

  return 0;
}

