#include "libscl.h"

using namespace scl;
using namespace std;

int main(int argc, char** argp, char** envp)
{

  INTEGER M = 3;
  INTEGER n = 3;

  realmat mu(M,1,0.0);
  realmat sig(M,M,0.0);
  for (INTEGER i=1; i<=M; ++i) sig(i,i) = 1.0;

  realmat xvals;
  realmat probs;

  mvnquad(n,M,mu,sig,xvals,probs);

  INTEGER len = probs.size();

  realmat muhat(M,1,0.0);
  realmat sighat(M,M,0.0);
  realmat x(M,1);

  for (INTEGER i=1; i<=len; ++i) {
    for (INTEGER j=1; j<=M; ++j) x[j] = xvals(i,j);
    muhat += probs[i]*x;
    sighat += probs[i]*(x*T(x));
  }

  cout << mu << muhat;
  cout << sig << sighat - mu*T(mu);

  for (INTEGER i=1; i<=M; ++i) {
    mu[i] = i;
    for (INTEGER j=1; j<=M; ++j) {
      if (i == j) {
        sig(i,j) = i;
      }
      else {
        sig(i,j) = (i + j) / 10.0; 
      }
    }
  }

  mvnquad(n,M,mu,sig,xvals,probs);

  len = probs.size();

  for (INTEGER i=1; i<=len; ++i) {
    for (INTEGER j=1; j<=M; ++j) x[j] = xvals(i,j);
    muhat += probs[i]*x;
    sighat += probs[i]*(x*T(x));
  }

  cout << mu << muhat;
  cout << sig << sighat - mu*T(mu);

  M = 2;
  n = 5;

  mu.resize(M,1);
  sig.resize(M,M,0.0);

  mu[1] = 0.95319;
  mu[2] = 24.503;

  sig(1,1) = pow(0.054691,2);
  sig(2,2) = pow(28.243,2);

  mvnquad(n,M,mu,sig,xvals,probs);

  cout << xvals << probs;

 
  return 0;
}
