stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Removing low discrepancy sequence DLL (licensing problems) 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@3004 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
sebastien 2009-10-01 08:44:47 +00:00
parent 6f55900bb9
commit e173e94839
5 changed files with 0 additions and 1932 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

59
mex/sources/qmc/build_qmc.m
View File

@ -1,59 +0,0 @@
addpath '../../../matlab';
MATLAB_PATH = matlabroot;
COMPILE_OPTIONS = '';
% mwSize, mwIndex and mwSignedIndex appeared in Matlab 7.3
if matlab_ver_less_than('7.3')
COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -DMWTYPES_NOT_DEFINED' ];
end
% Large array dims for 64 bits platforms appeared in Matlab 7.3
if (strcmpi('GLNXA64', computer) || strcmpi('PCWIN64', computer)) ...
&& ~matlab_ver_less_than('7.3')
COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -largeArrayDims' ];
end
if matlab_ver_less_than('7.5')
OUTPUT_DIR = '../../2007a';
else
OUTPUT_DIR = '../../2007b';
end
disp(' ')
if exist(OUTPUT_DIR,'dir')
disp('Delete old qmc mex file.')
delete([OUTPUT_DIR '/qmc_sequence.' mexext]);
else
whereami = pwd;
disp(['Create directory ' whereami(1:end-7) OUTPUT_DIR(4:end) '.'])
mkdir(OUTPUT_DIR);
end
disp(' ')
% Comment next line to suppress compilation debugging info
% COMPILE_OPTIONS = [ COMPILE_OPTIONS ' -v ' ];
% Set Optimization and Debug flags
CXXOPTIMFLAGS = ' CXXOPTIMFLAGS=-O3 ';
COPTIMFLAGS = ' COPTIMFLAGS=-O3 ';
CXXDEBUGFLAGS = ' CXXDEBUGFLAGS= ';
CDEBUGFLAGS = ' CDEBUGFLAGS= ';
LDOPTIMFLAGS = ' LDOPTIMFLAGS=-O3 ';
LDDEBUGFLAGS = ' LDDEBUGFLAGS= ';
COMPILE_OPTIONS = [ COMPILE_OPTIONS CDEBUGFLAGS COPTIMFLAGS CXXDEBUGFLAGS CXXOPTIMFLAGS LDDEBUGFLAGS LDOPTIMFLAGS];
FC='FC=/usr/bin/gfortran ' ;
FFLAGS= 'FFLAGS= ' ;
FLIBS = 'FLIBS=-L/usr/lib -lgfortran ' ;
FOPTIMFLAGS = 'FOPTIMFLAGS=-O ' ;
COMPILE_OPTIONS = [ COMPILE_OPTIONS FC FFLAGS FOPTIMFLAGS FLIBS ];
COMPILE_COMMAND = [ 'mex ' COMPILE_OPTIONS ' -outdir ' OUTPUT_DIR ] ;
disp('Compiling qmc')
system('gfortran -c -O3 low_discrepancy.f');
eval([ COMPILE_COMMAND ' qmc_sequence.cc low_discrepancy.o' ]) ;

50
mex/sources/qmc/doc/qmc_sequence.m
View File

@ -1,50 +0,0 @@
% Sobol sequence of quasi random numbers.
%
% This function produces a quasi monte carlo sequences.
%
% INPUT (for initialization)
% o dimension [uint32 scalar] dimension {2,3,...,1111}
% o flag [uint32 scalar] flag {0,1,2,3}
% o seed [uint32 scalar] seed (needed if flag>0)
% o transform [uint32 scalar] transform {0,1}
%
% OUTPUT (initialization)
%
% o qmc [structure]
% *** qmc.dimension [uint32 scalar]
% *** qmc.flag [uint32 scalar]
% *** qmc.transform [uint32 scalar]
% *** qmc.iteration [uint32 scalar]
% *** qmc.table_of_direction_numbers [uint32 matrix]
% *** qmc.table_common_denominator [uint32 scalar]
% *** qmc.last [double matrix]
% *** qmc.seed [uint32 scalar]
%
% INPUT ()
% o qmc [structure]
% o number_of_simulations [uint32 scalar] number_of_simulations*dimension
%
% OUTPUT ()
% o sArray [double matrix]
% o qmc [structure]
%
% ALGORITHMS
% Algorithm 659, Collected Algorithms From ACM. Published in
% Transactions on Mathematical Software, Vol. 14, n°1, P.88.
% Copyright (C) 2009 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/>.

1354
mex/sources/qmc/low_discrepancy.f
View File

File diff suppressed because it is too large Load Diff

405
mex/sources/qmc/qmc_sequence.cc
View File

@ -1,405 +0,0 @@
/*
* Copyright (C) 2009 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 "matrix.h"
#include "mex.h"
#if defined(__linux__)
#define SOBOL sobol_
#define INITSOBOL initsobol_
#else
#define SOBOL sobol
#define INITSOBOL initsobol
#endif
#define maxbit 30
extern "C"{
int SOBOL(double*, unsigned int*, unsigned int*, double*, unsigned int*, unsigned int*, unsigned int*, unsigned int*);
int INITSOBOL(unsigned int*, double*, unsigned int*, unsigned int*, unsigned int*, unsigned int*, unsigned int*);
}
void sobol_init(unsigned int* dimension, unsigned int* flag, unsigned int* seed, unsigned int* iteration, double* quasi, unsigned int* LL, unsigned int* SV)
{
INITSOBOL(dimension, quasi, LL, iteration, SV, flag, seed);
}
void sobol_array(unsigned int* number_of_simulations, unsigned int* dimension, unsigned int* transform,
unsigned int* LL, unsigned int* SV, double *sArray, unsigned int* iteration, double *quasi)
{
SOBOL(sArray, number_of_simulations, dimension, quasi, LL, iteration, SV, transform);
}
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
/*
** CASE 1 (nrhs==3 and nlhs==1) Initialization:
**
** prhs[0] ==> [integer scalar] dimension
** prhs[1] ==> [integer scalar] flag {0,1,2,3}
** prhs[2] ==> [integer scalar] seed
** prhs[3] ==> [integer scalar] transform {0,1}
**
** plhs[0] ==> [matlab structure] qmc
**
** qmc.dimension
** qmc.flag
** qmc.transform
** qmc.iteration
** qmc.table_of_direction_numbers
** qmc.table_common_denominator
** qmc.last
** qmc.seed
**
** CASE 2 (nrhs==2 and nlhs=2)
**
** prhs[0] ==> [matlab structure] qmc
** prhs[1] ==> [integer scalar] sequence_size
**
** plhs[0] ==> [matrix of doubles] sArray (sequence_size*dimension)
** plhs[1] ==> [matlab structure] qmc
*/
// Get and Check input and output:
if ( !(nrhs == 4 | nrhs == 2) )
{
mexErrMsgTxt("Four or two input arguments are required!");
}
if ( !(nlhs == 1 | nlhs == 2) )
{
mexErrMsgTxt("The number of output arguments should be two or one!");
}
if (nrhs==4)// CASE 1.
{
if (nlhs>1)
{
mexErrMsgTxt("The number of output arguments should be one!");
}
// Test of the first input argument (type):
if ( !( mxIsNumeric(prhs[0]) & mxIsClass(prhs[0],"uint32") ) )
{
mexPrintf("\t First input (dimension) has to be an integer [int32]. \n");
mexErrMsgTxt("\t Fatal error.");
}
// Test of the second input argument (type):
if ( !( mxIsNumeric(prhs[1]) & mxIsClass(prhs[1],"uint32") ) )
{
mexPrintf("\t First input (flag) has to be an integer [int32]. \n");
mexErrMsgTxt("\t Fatal error.");
}
// Test of the third input argument (type):
if ( !( mxIsNumeric(prhs[1]) & mxIsClass(prhs[1],"uint32") ) )
{
mexPrintf("\t First input (seed) has to be an integer [int32]. \n");
mexErrMsgTxt("\t Fatal error.");
}
// Test of the fourth input argument (type):
if ( !( mxIsNumeric(prhs[1]) & mxIsClass(prhs[1],"uint32") ) )
{
mexPrintf("\t First input (transform) has to be an integer [int32]. \n");
mexErrMsgTxt("\t Fatal error.");
}
/*
** Get the input variables.
*/
unsigned int dimension = (unsigned int) mxGetScalar(prhs[0]);
if (dimension>1111 | dimension < 2)
{
mexErrMsgTxt("dimension has to be between 2 and 1111!");
}
mxArray* DIMENSION = mxDuplicateArray(prhs[0]);
unsigned int flag = (unsigned int) mxGetScalar(prhs[1]);
if ( !( flag==0 | flag==1 | flag==2 | flag==3 ) )
{
mexErrMsgTxt("flag has to be equal to 0, 1, 2 or 3!");
}
mxArray* FLAG = mxDuplicateArray(prhs[1]);
unsigned int seed = ( unsigned int) mxGetScalar( prhs[2] );
unsigned int transform = (unsigned int) mxGetScalar(prhs[3]);
if ( !( transform==0 | transform==1) )
{
mexErrMsgTxt("transform has to be equal to 0 or 1!");
}
mxArray* TRANSFORM = mxDuplicateArray(prhs[3]);
/*
** Declaration of the variables returned by INITSOBOL.
*/
mxArray* t0;
unsigned int* LL = (unsigned int*) mxMalloc(sizeof(mxUINT32_CLASS));
t0 = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
LL = (unsigned int*) mxGetData(t0);
mxArray* t1;
unsigned int* SV = (unsigned int*) mxMalloc(dimension*maxbit*sizeof(mxUINT32_CLASS));
t1 = mxCreateNumericMatrix(dimension, maxbit, mxUINT32_CLASS, mxREAL);
SV = (unsigned int*) mxGetData(t1);
mxArray* t2;
unsigned int* iteration = (unsigned int*) mxMalloc((sizeof(mxUINT32_CLASS)));
t2 = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS, mxREAL);
iteration = (unsigned int*) mxGetData(t2);
mxArray* t3;
double* quasi = (double*) mxMalloc((dimension*sizeof(double)));
t3 = mxCreateNumericMatrix(dimension,1,mxDOUBLE_CLASS,mxREAL);
quasi = (double*) mxGetData(t3);
mxArray* t4;
unsigned int* seedout = (unsigned int*) mxMalloc(sizeof(mxUINT32_CLASS));
t4 = mxCreateNumericMatrix(1, 1, mxUINT32_CLASS,mxREAL);
seedout = (unsigned int*) mxGetData(t4);
/*
** Call to INITSOBOL subroutine.
*/
sobol_init(&dimension, &flag, &seed, iteration, quasi, LL, SV);
mxSetData(t0,LL);
mxSetData(t1,SV);
mxSetData(t2,iteration);
mxSetPr(t3,quasi);
*seedout = seed;
mxSetData(t4,seedout);
/*
** Now I build the returned matlab structure.
*/
char *fieldnames[8]; //This will hold field names.
fieldnames[0] = (char*) mxMalloc(sizeof("dimension"));
memcpy(fieldnames[0],"dimension",sizeof("dimension"));
fieldnames[1] = (char*) mxMalloc(sizeof("flag"));
memcpy(fieldnames[1],"flag",sizeof("flag"));
fieldnames[2] = (char*) mxMalloc(sizeof("transform"));
memcpy(fieldnames[2],"transform",sizeof("transform"));
fieldnames[3] = (char*) mxMalloc(sizeof("iteration"));
memcpy(fieldnames[3],"iteration",sizeof("iteration"));
fieldnames[4] = (char*) mxMalloc(sizeof("table_of_direction_numbers"));
memcpy(fieldnames[4],"table_of_direction_numbers",sizeof("table_of_direction_numbers"));
fieldnames[5] = (char*) mxMalloc(sizeof("table_common_denominator"));
memcpy(fieldnames[5],"table_common_denominator",sizeof("table_common_denominator"));
fieldnames[6] = (char*) mxMalloc(sizeof("last"));
memcpy(fieldnames[6],"last",sizeof("last"));
fieldnames[7] = (char*) mxMalloc(sizeof("seed"));
memcpy(fieldnames[7],"seed",sizeof("seed"));
plhs[0] = mxCreateStructMatrix(1,1,8,(const char**)fieldnames);
mxFree( fieldnames[0] );
mxFree( fieldnames[1] );
mxFree( fieldnames[2] );
mxFree( fieldnames[3] );
mxFree( fieldnames[4] );
mxFree( fieldnames[5] );
mxFree( fieldnames[6] );
mxFree( fieldnames[7] );
mxSetField(plhs[0], 0, "dimension", DIMENSION);
mxSetField(plhs[0], 0, "flag", FLAG);
mxSetField(plhs[0], 0, "transform", TRANSFORM);
mxSetField(plhs[0], 0, "iteration", t2);
mxSetField(plhs[0], 0, "table_of_direction_numbers", t1);
mxSetField(plhs[0], 0, "table_common_denominator", t0);
mxSetField(plhs[0], 0, "last", t3);
mxSetField(plhs[0], 0, "seed", t4);
}
if (nrhs==2)// CASE 2.
{
if (nlhs!=2)
{
mexErrMsgTxt("The number of output arguments has to be two!");
}
if ( !mxIsStruct(prhs[0]) )
{
mexErrMsgTxt("The first argument must be a matlab structure!");
}
else
{
if (mxGetNumberOfFields(prhs[0]) != 8)
{
mexErrMsgTxt("The matlab structure must have eight fields!");
}
}
/*
** Get the fields of the matlab structure.
*/
plhs[1] = mxDuplicateArray(prhs[0]);
mxArray* t0;// dimension
t0 = mxGetField( plhs[1], 0, "dimension");
if (t0 == NULL)
{
mexPrintf("\t Field ""dimension"" is empty. \n");
mexErrMsgTxt("\t Fatal error.");
}
if (!mxIsClass(t0,"uint32"))
{
mexPrintf("\t Field ""dimension"" is not uint32. \n");
mexErrMsgTxt("\t Fatal error.");
}
unsigned int* dimension = (unsigned int*) mxMalloc(sizeof( mxUINT32_CLASS ));
dimension = (unsigned int*) mxGetData( t0 );
// mxArray* t1;// flag
// t1 = mxGetField(plhs[1], 0, "flag");
// if (t1 == NULL)
// {
// mexPrintf("\t Field ""flag"" is empty. \n");
// mexErrMsgTxt("\t Fatal error.");
// }
// if (!mxIsClass(t1,"uint32"))
// {
// mexPrintf("\t Field ""flag"" is not uint32. \n");
// mexErrMsgTxt("\t Fatal error.");
// }
// unsigned int* flag = (unsigned int*) mxMalloc(sizeof( mxUINT32_CLASS ));
// flag = (unsigned int*) mxGetData( t1 );
mxArray* t2;// transform
t2 = mxGetField(plhs[1], 0, "transform");
if (t2 == NULL)
{
mexPrintf("\t Field ""transform"" is empty. \n");
mexErrMsgTxt("\t Fatal error.");
}
if (!mxIsClass(t2,"uint32"))
{
mexPrintf("\t Field ""transform"" is not uint32. \n");
mexErrMsgTxt("\t Fatal error.");
}
unsigned int* transform = (unsigned int*) mxMalloc(sizeof( mxUINT32_CLASS ));
transform = (unsigned int*) mxGetData( t2 );
mxArray* t3;// iteration
t3 = mxGetField(plhs[1], 0, "iteration");
if (t3 == NULL)
{
mexPrintf("\t Field ""iteration"" is empty. \n");
mexErrMsgTxt("\t Fatal error.");
}
if (!mxIsClass(t3,"uint32"))
{
mexPrintf("\t Field ""iteration"" is not uint32. \n");
mexErrMsgTxt("\t Fatal error.");
}
unsigned int* iteration = (unsigned int*) mxMalloc(sizeof( mxUINT32_CLASS ));
iteration = (unsigned int*) mxGetData( t3 );
mxArray* t4;// table_of_direction_numbers
t4 = mxGetField(plhs[1], 0, "table_of_direction_numbers");
if (t4 == NULL)
{
mexPrintf("\t Field ""table_of_direction_numbers"" is empty. \n");
mexErrMsgTxt("\t Fatal error.");
}
if (!mxIsClass(t4,"uint32"))
{
mexPrintf("\t Field ""table_of_direction_numbers"" is not uint32. \n");
mexErrMsgTxt("\t Fatal error.");
}
unsigned int* SV = (unsigned int*) mxMalloc(*dimension*maxbit*sizeof( mxUINT32_CLASS ));
SV = (unsigned int*) mxGetData( t4 );
mxArray* t5;//table_common_denominator
t5 = mxGetField(plhs[1], 0, "table_common_denominator");
if (t5 == NULL)
{
mexPrintf("\t Field ""table_common_denominator"" is empty. \n");
mexErrMsgTxt("\t Fatal error.");
}
if (!mxIsClass(t5,"uint32"))
{
mexPrintf("\t Field ""table_common_denominator"" is not uint32. \n");
mexErrMsgTxt("\t Fatal error.");
}
unsigned int* LL = (unsigned int*) mxMalloc(sizeof( mxUINT32_CLASS ));
LL = (unsigned int*) mxGetData( t5 );
mxArray* t6;// last
t6 = mxGetField(plhs[1], 0, "last");
if (t6 == NULL)
{
mexPrintf("\t Field ""last"" is empty. \n");
mexErrMsgTxt("\t Fatal error.");
}
if (!mxIsClass(t6,"double"))
{
mexPrintf("\t Field ""last"" is not double. \n");
mexErrMsgTxt("\t Fatal error.");
}
double* last = (double*) mxMalloc(*dimension*sizeof(double));
last = (double*) mxGetData( t6 );
// mxArray* t7;// seed
// t7 = mxGetField(prhs[0], 0, "seed");
// if (t7 == NULL)
// {
// mexPrintf("\t Field ""seed"" is empty. \n");
// mexErrMsgTxt("\t Fatal error.");
// }
// if (!mxIsClass(t7,"uint32"))
// {
// mexPrintf("\t Field ""seed"" is not uint32. \n");
// mexErrMsgTxt("\t Fatal error.");
// }
// unsigned int* seed = (unsigned int*) mxMalloc(sizeof( mxUINT32_CLASS ));
// seed = (unsigned int*) mxGetData( t7 );
/*
** Get the second input.
*/
unsigned int number_of_simulations = (unsigned int) mxGetScalar(prhs[1]);
/*
** Initialization of the first output argument.
*/
plhs[0] = mxCreateNumericMatrix(number_of_simulations, *dimension, mxDOUBLE_CLASS, mxREAL);
double* sArray = (double*)mxMalloc(*dimension*number_of_simulations*sizeof(double));
mxSetPr(plhs[0],sArray);
/*
** Call to SOBOL subroutine.
*/
sobol_array(&number_of_simulations, dimension, transform, LL, SV, sArray, iteration, last);
mxSetData(t6,last);
mxSetField(plhs[1], 0, "last", t6);
mxSetData(t3,iteration);
mxSetField(plhs[1], 0, "iteration", t3);
}
}

64
mex/sources/qmc/tests/test_qmc_sequence.m
View File

@ -1,64 +0,0 @@
function test_qmc_sequence()
% Copyright (C) 2009 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/>.
transform = uint32(0);
dimension = uint32(2);
flag = uint32(3);
seed = uint32(1477);
count = 0;
for i=1:20
disp('Step 1...')
qmc_init = qmc_sequence(dimension,flag,seed,transform);
disp('Done!')
disp('Step 2...')
number_of_simulations = uint32(20);
[Q, qmc_new] = qmc_sequence(qmc_init,number_of_simulations);
disp('Done!')
number_of_simulations = uint32(10);
disp('Step 3...')
[Q1, qmc_new1] = qmc_sequence(qmc_init,number_of_simulations);
disp('Done!')
disp('Step 4...')
[Q2, qmc_new2] = qmc_sequence(qmc_new1,number_of_simulations);
disp('Done!')
Q3 = Q - [Q1;Q2] ;
test = max(max(abs(Q3)));
if test<2*eps
count = count+1;
end
if i==5 | i== 10
clear mex;
close all;
end
end
disp('qmc_sequence:: No error!')
clear all