Added new mex file for computing Quasi Monte-Carlo sequences (Sobol), texinfo...

Added new mex file for computing Quasi Monte-Carlo sequences (Sobol), texinfo headers and unitary tests are available in <DYNARE_PATH>/matlab/qmc_sequence.m.

matlab/dynare_config.m

@@ -246,5 +246,15 @@ else
 end
 if verbose
     disp([ message 'k-order solution simulation.' ])
+end
+
+% Test if qmc_sequence DLL is present
+if exist('qmc_sequence', 'file') == 3
+    message = '[mex] ';
+else
+    message = '[no]  ';
+end
+if verbose
+    disp([ message 'Quasi Monte-Carlo sequence (Sobol).' ])
     disp(' ')
 end

matlab/qmc_sequence.m

+%@info:
+%! @deftypefn {Mex File} {[@var{a}, @var{s}] =} qmc_sequence (@var{d}, @var{s}, @var{t}, @var{n}, @var{lu})
+%! @anchor{qmc_sequence}
+%! @sp 1
+%! Computes quasi Monte-Carlo sequence (Sobol numbers).
+%! @sp 2
+%! @strong{Inputs}
+%! @sp 1
+%! @table @ @var
+%! @item d
+%! Integer scalar, dimension.
+%! @item s
+%! Integer scalar (int32 or int64), seed.
+%! @item t
+%! Integer scalar, sequence type:
+%!  @sp 1
+%!  @table @ @samp
+%!  @item @var{t}=0
+%!  Uniform numbers in a n-dimensional (unit by default) hypercube 
+%!  @item @var{t}=1
+%!  Gaussian numbers
+%!  @item @var{t}=2
+%!  Uniform numbers on a n-dimensional (unit by default) hypersphere
+%!  @end table
+%! @item n
+%! Integer scalar, number of elements in the sequence.
+%! @item lu
+%! Optional argument, the interpretation depends on its size:
+%!  @sp 1
+%!  @table @ @samp
+%!  @item @var{d}x2 array of doubles
+%!  Lower and upper bounds of the hypercube (default is 0-1 in all dimensions). @var{t} must be equal to zero.
+%!  @item @var{d}x@var{d} array of doubles
+%!  Covariance matrix of the Gaussian variates (default is the identity matrix). @var{t} must be equal to one.
+%!  @item scalar double
+%!  Radius of the hypershere (default is one). @var{t} must be equal to two.
+%!  @end table
+%! @end table
+%! @sp 2
+%! @strong{Outputs}
+%! @sp 1
+%! @table @ @var
+%! @item a
+%! @var{n}x@var{d} matrix of doubles, the Sobol sequence.
+%! @item s
+%! Integer scalar, current value of the seed.
+%! @end table
+%! @sp 2
+%! @strong{This function is called by:}
+%! @sp 2
+%! @strong{This function calls:}
+%! @sp 2
+%! @end deftypefn
+%@eod:
+
+%@test:1
+%$ t = ones(3,1);
+%$
+%$ d = 2;
+%$ n = 100;
+%$ s = int64(0);
+%$
+%$ try
+%$   [draws, S] = qmc_sequence(d,s,0,n);
+%$ catch
+%$   t(1) = 0;
+%$ end
+%$
+%$ try
+%$   [draws, S] = qmc_sequence(d,s,1,n);
+%$ catch
+%$   t(2) = 0;
+%$ end
+%$
+%$ try
+%$   [draws, S] = qmc_sequence(d,s,2,n);
+%$ catch
+%$   t(3) = 0;
+%$ end
+%$
+%$ T = all(1);
+%@eof:1
\ No newline at end of file

mex/build/matlab/Makefile.am

@@ -2,7 +2,7 @@ ACLOCAL_AMFLAGS = -I ../../../m4
 
 # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
 if DO_SOMETHING
-SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ estimation block_kalman_filter
+SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ estimation block_kalman_filter sobol
 
 if HAVE_GSL
 SUBDIRS += ms_sbvar

mex/build/matlab/configure.ac

@@ -143,6 +143,7 @@ AC_CONFIG_FILES([Makefile
                  libslicot/Makefile
                  kalman_steady_state/Makefile
                  ms_sbvar/Makefile
-                 block_kalman_filter/Makefile])
+                 block_kalman_filter/Makefile
+	         sobol/Makefile])
 
 AC_OUTPUT

mex/build/matlab/sobol/Makefile.am

+include ../mex.am
+include ../../qmc_sequence.am

mex/build/octave/Makefile.am

@@ -2,7 +2,7 @@ ACLOCAL_AMFLAGS = -I ../../../m4
 
 # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
 if DO_SOMETHING
-SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ qzcomplex ordschur block_kalman_filter
+SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ qzcomplex ordschur block_kalman_filter sobol
 
 if HAVE_GSL
 SUBDIRS += ms_sbvar

mex/build/octave/configure.ac

@@ -130,6 +130,7 @@ AC_CONFIG_FILES([Makefile
                  libslicot/Makefile
                  kalman_steady_state/Makefile
                  ms_sbvar/Makefile
-                 block_kalman_filter/Makefile])
+                 block_kalman_filter/Makefile
+		 sobol/Makefile])
 
 AC_OUTPUT

mex/build/octave/sobol/Makefile.am

+include ../mex.am
+include ../../qmc_sequence.am

mex/build/qmc_sequence.am

+vpath %.cc $(top_srcdir)/../../sources/sobol
+
+noinst_PROGRAMS = qmc_sequence
+
+nodist_qmc_sequence_SOURCES = qmc_sequence.cc

mex/sources/Makefile.am

@@ -13,7 +13,8 @@ EXTRA_DIST = \
 	libslicot \
 	kalman_steady_state \
 	ms-sbvar \
-	block_kalman_filter
+	block_kalman_filter \
+	sobol
 
 clean-local:
 	rm -rf `find mex/sources -name *.o`

mex/sources/sobol/gaussian.hh

+/*  Generates gaussian random deviates from uniform random deviates.
+** 
+**  Pseudo code of the algorithm is given at http://home.online.no/~pjacklam/notes/invnorm  
+**
+**  Copiright (C)  2010 Dynare Team
+**
+** This file is part of Dynare.
+**
+** Dynare 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 3 of the License, or
+** (at your option) any later version.
+**
+** Dynare 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 Dynare.  If not, see <http://www.gnu.org/licenses/>.  
+*/
+
+#include <cstdlib>
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <ctime>
+
+#include <dynblas.h>
+
+using namespace std;
+
+#define lb  .02425
+#define ub  .97575
+
+#ifdef USE_OMP
+# include <omp.h>
+#endif
+#define DEBUG_OMP 0
+
+template<typename T> T icdf( const T uniform )
+/*
+**  This function invert the gaussian cumulative distribution function.
+**
+*/
+{ 
+  static T A[6] = 
+    {
+      -3.969683028665376e+01,
+       2.209460984245205e+02,
+      -2.759285104469687e+02,
+       1.383577518672690e+02,
+      -3.066479806614716e+01,
+       2.506628277459239e+00
+    };
+  static T B[5] = 
+    {
+      -5.447609879822406e+01,
+       1.615858368580409e+02,
+      -1.556989798598866e+02,
+       6.680131188771972e+01,
+      -1.328068155288572e+01
+    };
+  static T C[6] = 
+    {
+      -7.784894002430293e-03,
+      -3.223964580411365e-01,
+      -2.400758277161838e+00,
+      -2.549732539343734e+00,
+       4.374664141464968e+00,
+       2.938163982698783e+00
+    };
+  static T D[4] = 
+    {
+      7.784695709041462e-03,
+      3.224671290700398e-01,
+      2.445134137142996e+00,
+      3.754408661907416e+00
+    };
+  T gaussian = (T)0.0;
+  if ( (0<uniform)  && (uniform<lb) )
+    {
+      T tmp;
+      tmp = sqrt(-2*log(uniform));
+      gaussian = (((((C[0]*tmp+C[1])*tmp+C[2])*tmp+C[3])*tmp+C[4])*tmp+C[5])/((((D[0]*tmp+D[1])*tmp+D[2])*tmp+D[3])*tmp+1);
+    }
+  else
+    {
+      if ((lb <= uniform) && (uniform <= ub))
+        {
+          T tmp, TMP;
+          tmp = uniform - .5;
+          TMP = tmp*tmp;
+          gaussian = (((((A[0]*TMP+A[1])*TMP+A[2])*TMP+A[3])*TMP+A[4])*TMP+A[5])*tmp/(((((B[0]*TMP+B[1])*TMP+B[2])*TMP+B[3])*TMP+B[4])*TMP+1);
+        }
+      else
+        {
+          if ((ub < uniform) && (uniform < 1))
+            {
+              T tmp;
+              tmp = sqrt(-2*log(1-uniform));
+              gaussian = -(((((C[0]*tmp+C[1])*tmp+C[2])*tmp+C[3])*tmp+C[4])*tmp+C[5])/((((D[0]*tmp+D[1])*tmp+D[2])*tmp+D[3])*tmp+1);
+            }
+        }
+    }
+  if ( (0<uniform) && (uniform<1) )
+    {
+      T tmp, tmp_;
+      tmp = .5*erfc(-gaussian/sqrt(2.0))-uniform;
+      tmp_ = tmp*sqrt(2*M_PI)*exp(.5*gaussian*gaussian);
+      gaussian = gaussian - tmp_/(1+.5*gaussian*tmp_);
+    }
+ return(gaussian);
+}
+
+template<typename T> void icdfm( const int n, T *U)
+{ 
+  #if USE_OMP
+  #pragma omp parallel for num_threads(omp_get_num_threads())
+  #endif
+  for(int i=0; i<n; i++)
+    {
+      U[i] = icdf(U[i]);
+    }
+  return;
+}
+
+template<typename T> void icdfmSigma( const int d, const int n, T *U, const double *LowerCholSigma)
+{ 
+  double one = 1.0;
+  icdfm(n*d, U);
+  double tmp[n*d];
+  dgemm("N","N",(blas_int*)d,(blas_int*)n,(blas_int*)d,&one,LowerCholSigma,(blas_int*)d,U,(blas_int*)d,&one,tmp,(blas_int*)d);
+  memcpy(U,tmp,d*n*sizeof(double));
+  return;
+}
+
+template<typename T> void usphere( const int d, const int n, T *U)
+{ 
+  icdfm(n*d, U);
+  #if USE_OMP
+  #pragma omp parallel for num_threads(omp_get_num_threads())
+  #endif
+  for (int j=0; j<n; j++)// sequence index.
+    {
+      int k = j*d;
+      double norm = 0.0;
+      for(int i=0; i<d; i++)// dimension index.
+        {
+          norm = norm + U[k+i]*U[k+i];
+        }
+      norm = sqrt(norm);
+      for(int i=0; i<d; i++)// dimension index.
+        {
+          U[k+i] = U[k+i]/norm;
+        }
+    }
+  return;
+}
+
+template<typename T> void usphereRadius( const int d, const int n, double radius, T *U)
+{ 
+  icdfm(n*d, U);
+  #if USE_OMP
+  #pragma omp parallel for num_threads(omp_get_num_threads())
+  #endif
+  for (int j=0; j<n; j++)// sequence index.
+    {
+      int k = j*d;
+      double norm = 0.0;
+      for(int i=0; i<d; i++)// dimension index.
+        {
+          norm = norm + U[k+i]*U[k+i];
+        }
+      norm = sqrt(norm);
+      for(int i=0; i<d; i++)// dimension index.
+        {
+          U[k+i] = radius*U[k+i]/norm;
+        }
+    }
+  return;
+}

mex/sources/sobol/initialize_v_array.hh

+template<typename T> int initialize_v_array (int dim_max, int log_max, T **v)
+/*
+** This function initializes the v array used in the sobol routine. 
+**
+**  Original files downloaded from http://people.sc.fsu.edu/~burkardt/cpp_src/sobol/ (version 17-Feb-2009 09:46)
+**
+**  Copyright (C) -2009 John Burkardt
+**  Copiright (C)  2009 Stéphane Adjemian
+**
+*/
+{
+  #include "initialize_v_array.inc"
+  return 1;
+}

mex/sources/sobol/qmc_sequence.cc

+/*
+ * Copyright (C) 2010 Dynare Team
+ *
+ * This file is part of Dynare (can be used outside Dynare).
+ *
+ * Dynare 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 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare 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 Dynare.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <string.h>
+#include <dynmex.h>
+
+#include "sobol.hh"
+#include "gaussian.hh"
+
+void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+  /* 
+  ** INPUTS:
+  ** prhs[0]    [integer]    scalar, dimension.
+  ** prhs[1]    [integer]    scalar, seed.
+  ** prhs[2]    [integer]    scalar, sequence type:
+  **                                   0 ==> uniform,
+  **                                   1 ==> gaussian,
+  **                                   2 ==> uniform on an hypershere.
+  ** prhs[3]    [integer]    scalar, sequence size.
+  ** prhs[4]    [double]     dimension*2 array, lower and upper bounds of the hypercube (default is 0-1 in all dimensions) if prhs[2]==0,
+  **                         dimension*dimension array, covariance of the multivariate gaussian distribution of prhs[2]==1 (default is the identity matrix),
+  **                         scalar, radius of the hypershere if prhs[2]==2 (default is one).
+  **
+  ** OUTPUTS:
+  ** plhs[0]    [double]     sequence_size*dimension array, the Sobol sequence.
+  ** plhs[1]    [integer]    scalar, seed. 
+  ** 
+  */
+  /*
+  ** Check the number of input and output arguments.
+  */
+  if ( !( (nrhs==5) | (nrhs==4) | (nrhs==3) )  )
+    {
+      mexErrMsgTxt("Five, four or three input arguments are required!");
+    }
+  if ( !(nlhs == 2) )
+    {
+      mexErrMsgTxt("The number of output arguments has to be two");
+    }
+  /*
+  ** Test the first input argument and assign it to dimension.
+  */
+  if (  !( mxIsNumeric(prhs[0]) ) )
+    {
+      mexPrintf("\t First input (dimension) has to be a positive integer. \n");
+      mexErrMsgTxt("\t Fatal error.");
+    }
+  int dimension = ( int ) mxGetScalar(prhs[0]);
+  /*
+  ** Test the second input argument and assign it to seed.
+  */
+  if ( !( mxIsNumeric(prhs[1]) && mxIsClass(prhs[1],"int64") ) )
+    {
+      mexPrintf("\t Second input (seed) has to be an integer [int64]. \n");
+      mexErrMsgTxt("\t Fatal error.");
+    }
+  int64_T seed = ( int64_T ) mxGetScalar( prhs[1] );
+  /*
+  ** Test the third input argument and assign it to type (kind of QMC sequence).
+  */
+  int error_flag_3 = 0;
+  if (  !(mxIsNumeric(prhs[2])) )
+    {
+      error_flag_3 = 1;
+    }
+  int type = ( int ) mxGetScalar(prhs[2]);
+  if ( !(type==0 || type==1 || type==2) )
+    {
+      error_flag_3 = 1;
+    }
+  if (error_flag_3==1)
+    {
+      mexPrintf("\t Third input (type of QMC sequence) has to be an integer equal to 0, 1 or 2: \n");
+      mexPrintf("\t   0 ==> Points uniformly distributed in an hypercube. \n");
+      mexPrintf("\t   1 ==> Points normally distributed in R^n. \n");
+      mexPrintf("\t   2 ==> Points uniformly distributed on an hypershere. \n");
+      mexErrMsgTxt("\t Fatal error.");
+    }
+  /*
+  ** Test dimension>=2 when type==2
+  */
+  if ( (type==2) && (dimension<2) )
+    {
+      mexPrintf("\t First input (dimension) has to be greater than 1 for a uniform QMC on an hypershere.\n");
+      mexErrMsgTxt("\t Fatal error.");
+    }
+  /*
+  ** Test the optional fourth input argument and assign it to sequence_size.
+  */
+  if  ( ( nrhs>3 )  &&  !mxIsNumeric(prhs[3])  ) 
+    {
+      mexPrintf("\t Fourth input (qmc sequence size) has to be a positive integer. \n");
+      mexErrMsgTxt("\t Fatal error.");
+    }
+  int sequence_size;
+  if ( nrhs>3)
+    {
+      sequence_size = ( int ) mxGetScalar( prhs[3] );
+    }
+  else
+    {
+      sequence_size = 1;
+    }
+  /*
+  ** Test the optional fifth input argument and assign it to lower_and_upper_bounds.
+  */
+  if  (  ( nrhs>4 )  &&  (type==0) && ( !( mxGetN(prhs[4])==2) )  )// Sequence of uniformly distributed numbers in an hypercube.
+    {
+      mexPrintf("\t The fifth input argument must be an array with two columns. \n");
+      mexErrMsgTxt("\t Fatal error.");
+    }
+  if  ( (nrhs>4)   &&  (type==0) &&  ( !( (int)mxGetM(prhs[4])==dimension) ) )
+    {
+      mexPrintf("\t The fourth input argument must be an array with a number of lines equal to dimension (first input argument). \n");
+      mexErrMsgTxt("\t Fatal error.");
+    }
+  if  ( ( nrhs>4 )  &&  (type==1) && ( !( ((int)mxGetN(prhs[4])==dimension) && ((int)mxGetM(prhs[4])==dimension) ) ) )// Sequence of normally distributed numbers.
+    {
+      mexPrintf("\t The fifth input argument must be a squared matrix (whose dimension is given by the first input argument). \n");
+      mexErrMsgTxt("\t Fatal error.");
+    }
+  if  ( ( nrhs>4 )  &&  (type==2) && ( !( (mxGetN(prhs[4])==1) && (mxGetM(prhs[4])==1) ) ) )// Sequence of uniformly distributed numbers on an hypershere.
+    {
+      mexPrintf("\t The fifth input argument must be a positive scalar. \n");
+      mexErrMsgTxt("\t Fatal error.");
+    }
+  double *lower_bounds, *upper_bounds;
+  int unit_hypercube_flag = 1;
+  if ( (type==0) && (nrhs>4) )
+    {
+      lower_bounds = (double *) mxCalloc(dimension,sizeof(double));
+      upper_bounds = (double *) mxCalloc(dimension,sizeof(double));
+      double *tmp;
+      tmp = (double *) mxCalloc(dimension*2,sizeof(double));
+      memcpy(tmp,mxGetPr(prhs[4]),dimension*2*sizeof(double));
+      lower_bounds = &tmp[0];
+      upper_bounds = &tmp[dimension];
+      unit_hypercube_flag = 0;
+    }
+  double *cholcov;
+  int identity_covariance_matrix = 1;
+  if ( (type==1) && (nrhs>4) )
+    {
+      cholcov = (double *) mxCalloc(dimension*dimension,sizeof(double));
+      double *tmp;
+      tmp = (double *) mxCalloc(dimension*dimension,sizeof(double));
+      memcpy(tmp,mxGetPr(prhs[4]),dimension*2*sizeof(double));
+      cholcov = &tmp[0];
+      identity_covariance_matrix = 0;
+    }
+  double radius = 1.0;
+  int unit_radius = 1;
+  if ( (type==2) && (nrhs>4) )
+    {
+      double *tmp;
+      tmp = (double *) mxCalloc(1,sizeof(double));
+      memcpy(tmp,mxGetPr(prhs[4]),dimension*2*sizeof(double));
+      radius = tmp[0];
+      unit_radius = 0;
+    }
+  /*
+  ** Initialize outputs of the mex file.
+  */
+  double *qmc_draws;
+  plhs[0] = mxCreateDoubleMatrix(dimension,sequence_size,mxREAL);
+  qmc_draws = mxGetPr(plhs[0]);
+  int64_T *seed_out;
+  plhs[1] = mxCreateNumericMatrix(1, 1, mxINT64_CLASS, mxREAL);
+  seed_out = (int64_T *) mxGetData(plhs[1]);
+  if (type==0)// Uniform QMC sequence in an hypercube.
+    {
+      if (sequence_size==1)
+        {
+          next_sobol ( dimension, &seed, qmc_draws );
+          *seed_out = seed;
+        }
+      else
+        {
+          *seed_out = sobol_block( dimension, sequence_size, seed, qmc_draws);
+        }
+      if (unit_hypercube_flag==0)
+        {
+          expand_unit_hypercube( dimension, sequence_size, qmc_draws, lower_bounds, upper_bounds);
+        }
+      return;
+    }
+  if (type==1)// Normal QMC sequance in R^n.
+    {
+      if (sequence_size==1)
+        {
+          next_sobol ( dimension, &seed, qmc_draws );
+          *seed_out = seed;
+        }
+      else
+        {
+          *seed_out = sobol_block( dimension, sequence_size, seed, qmc_draws);
+        }
+      if (identity_covariance_matrix==1)
+        {
+          icdfm(dimension*sequence_size, qmc_draws);
+        }
+      else
+        {
+          icdfmSigma(dimension,sequence_size, qmc_draws, cholcov);
+        }
+      return;
+    }
+  if (type==2)// Uniform QMC sequence on an hypershere.
+    {
+      if (sequence_size==1)
+        {
+          next_sobol ( dimension, &seed, qmc_draws );
+          *seed_out = seed;
+        }
+      else
+        {
+          *seed_out = sobol_block( dimension, sequence_size, seed, qmc_draws);
+        }
+      if (unit_r