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

Copyright (C) 2018

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.

------------------------------------------------------------------------------

Function      mvnquad - Compute quadrature abcissae and weights for mvn

Syntax        #include "libscl.h"

              void mvnquad(INTEGER n, INTEGER M, const realmat& mu, 
	            const realmat& sig, realmat& xvals, realmat& probs);


Prototype in  libscl.h

Description   n is the number of univariate quad points mu x is an M by 1 
              vector containing the mean, and sig is an M by M positive
              definite symmetric variance-covariance matrix. xvals is an
	      n^M by M matrix of abcissae and probs is an n^M by 1 vector 
	      of corresponding weights.

Remarks       If factorization sig by scl::cholesky fails or sig is less
              than full rank, then scl::error is called.  If scl::error
              is not modified by the calling program, scl::error will
              terminate execution.  

Return value  none

Functions     Library: sqrt
called        libscl: cholesky, hquad

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


#include "libscl.h"
using namespace scl;
using namespace std;

namespace {

 // Computes all multi indexes of type (i_1, i_2, ..., i_d) where each
 // element i_l is one of the integers 1, 2, ..., n;

 // Recall that indexing midx starts with 0

  void multi(INTEGER n, INTEGER d, INTEGER k, intvec& jv, vector<intvec>& midx)
  {
    if (k == d) { // external call must have k == d
      INTEGER m = 1;
      for (INTEGER i=1; i<=d; ++i) m *= n;
      midx.reserve(m);
      jv.resize(d);
    }

    for (INTEGER j=1; j<=n; ++j) {
      jv[k] = j;
      if (k == 1) {
        midx.push_back(jv);
      }
      else {
        multi(n, d, k-1, jv, midx);
      }
    }
  }  

  void multi(INTEGER n, INTEGER d, vector<intvec>& midx)
  {
    intvec jvec;
    multi(n, d, d, jvec, midx);
  }

}

void mvnquad(INTEGER n, INTEGER M, const realmat& mu, const realmat& sig,
  realmat& xvals, realmat& probs)
{

 if (mu.size() != M) error("Error, mvnquad, mu dimension wrong");

 if (sig.nrow() != M || sig.ncol() != M) {
   error("Error, mvnquad, sig dimensions wrong");
 }

 const REAL tolerance = 16.0*EPS;

 realmat R;
 INTEGER rank = cholesky(sig,R,tolerance);

  if (rank != M) error("Error, mvnquad, sig is less than full rank");

  realmat x, w;

  INTEGER ier = hquad(n,x,w); 
  
  if (ier != 0) error("Error, mvnquad, hquad failed");

  const REAL root2 = M_SQRT2;
  const REAL rootpi = sqrt(M_PI);

  realmat abscissae(n,1);
  realmat weights(n,1);

  for (INTEGER i=1; i<=n; ++i) {
    abscissae[i] = root2*x[i];
    weights[i] = w[i]/rootpi;
  }

  vector<intvec> midx;
  multi(n,M,midx);

  INTEGER len = midx.size();

  xvals.resize(len,M);
  probs.resize(len,1);

  INTEGER i = 0;
  for(vector<intvec>::const_iterator m=midx.begin(); m!=midx.end(); ++m) {
    ++i;
    REAL p = 1.0;
    for (INTEGER j=1; j<=M; ++j) {
      xvals(i,j) = abscissae[(*m)[j]];
      p *= weights[(*m)[j]];
    }
    probs[i] = p;
  }

  if (i != len) error("Error, mvnquad, should never happen");

  xvals = xvals*R;

  for (INTEGER i=1; i<=len; ++i) {
    for (INTEGER j=1; j<=M; ++j) xvals(i,j) += mu[j];
  }

}
