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

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 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 factorizing 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 scl {

  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];
    }
  
  }
}
