/*!\file: random
 * \brief random number generating functions
 */ 

/*Headers*/
/*{{{*/
#include <stdio.h>
#include <sys/types.h>
#include <math.h>
#include <float.h>    /*  DBL_EPSILON  */
#include <chrono>
#include <cstdarg>
#include <iostream>
#include <random>

#include "../Matrix/matrix.h"
#include "../Exceptions/exceptions.h"
#include "../MemOps/MemOps.h"
#include "../io/io.h"
/*}}}*/

void univariateNormal(IssmDouble* prand, IssmDouble mean, IssmDouble sdev) { /*{{{*/
	/*univariateNormal	generates a random value follwoing Normal distribution*/
	unsigned seed = std::chrono::steady_clock::now().time_since_epoch().count(); //random seed using time_since_epoch
   std::default_random_engine generator(seed); //generator of random numbers
   std::normal_distribution<double> normdistri(mean,sdev); //Normal probability distribution
	double tfunc;
	*prand = normdistri(generator);
} /*}}}*/
void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix) { /*{{{*/
   IssmDouble* sampleStandardNormal     = xNew<IssmDouble>(dim);
   IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
   IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
   for(int ii{0};ii<dim;ii++){univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);}
   CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
   for(int ii{0};ii<dim;ii++){ //matrix by vector multiplication
      double sum{0.0};
      for(int jj{0};jj<dim;jj++){sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];} //entry-by-entry multiplication along matrix row
      sampleMultivariateNormal[ii] = mean+sum; //assign value
   }
   *prand = sampleMultivariateNormal;
   xDelete<IssmDouble>(sampleStandardNormal);
   xDelete<IssmDouble>(Lchol);
} /*}}}*/
void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix) { /*{{{*/
	IssmDouble* sampleStandardNormal     = xNew<IssmDouble>(dim);
	IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
	IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
	for(int ii{0};ii<dim;ii++){univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);}
	CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
	for(int ii{0};ii<dim;ii++){ //matrix by vector multiplication
		double sum{0.0};
      for(int jj{0};jj<dim;jj++){sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];} //entry-by-entry multiplication along matrix row
      sampleMultivariateNormal[ii] = mean[ii]+sum; //assign value
	}
	*prand = sampleMultivariateNormal;
	xDelete<IssmDouble>(sampleStandardNormal);
	xDelete<IssmDouble>(Lchol);
} /*}}}*/