...

Commits (103)
 ... ... @@ -104,7 +104,6 @@ doc/internals/ltxpng /dynare++/kord/tests /dynare++/kord/tests.exe /dynare++/kord/out.txt /dynare++/sylv/testing/*.mm /dynare++/sylv/testing/tests /dynare++/sylv/testing/tests.exe /dynare++/parser/cc/*_ll.cc ... ... @@ -117,6 +116,10 @@ doc/internals/ltxpng /dynare++/src/dynglob_tab.hh /dynare++/tl/testing/tests /dynare++/tl/testing/tests.exe /dynare++/tests/*.jnl /dynare++/tests/*.m /dynare++/tests/*.mat /dynare++/tests/*.dump !/dynare++/extern/R/Makefile # Windows ... ...
 ... ... @@ -107,6 +107,10 @@ test_dynare++: artifacts: paths: - dynare++/kord/out.txt - dynare++/tests/*.jnl - dynare++/tests/*.m - dynare++/tests/*.mat - dynare++/tests/*.dump deploy_manual_unstable: stage: deploy ... ...
 ... ... @@ -9,7 +9,7 @@ Described on the homepage: Most users should use the precompiled package available for your OS, also Most users should use the precompiled package available for their OS, also available via the Dynare homepage: . # Contributions ... ... @@ -38,13 +38,19 @@ Note that if you obtain the source code via git, you will need to install more t The first section of this page gives general instructions, which apply to all platforms. Then some specific platforms are discussed. **Note:** Here, when we refer to 32-bit or 64-bit, we refer to the type of MATLAB installation, not the type of Windows installation. It is perfectly possible to run a 32-bit MATLAB on a 64-bit Windows: in that case, instructions for Windows 32-bit should be followed. To determine the type of your MATLAB installation, type: **Note:** Here, when we refer to 32-bit or 64-bit, we refer to the type of MATLAB or Octave installation, not the type of operating system installation. For example, it is perfectly possible to run a 32-bit MATLAB on a 64-bit Windows: in that case, instructions for Windows 32-bit should be followed. To determine the type of your MATLAB/Octave installation, type: matlab >> computer  at the MATLAB prompt: if it returns PCWIN64, GLNX64 or MACI64, then you have a 64-bit MATLAB; if it returns PCWIN, MACI or GLNX, then you have a 32-bit MATLAB. at the MATLAB/Octave prompt. Under MATLAB, if it returns PCWIN64, GLNX64 or MACI64, then it is a 64-bit MATLAB; if it returns PCWIN, MACI or GLNX, then it is a 32-bit MATLAB. Under Octave, if it returns a string that begins with x86_64, it is a 64-bit Octave; if the strings begins with i686, it is a 32-bit Octave. **Contents** ... ... @@ -102,28 +108,39 @@ Simply launch the configure script from a terminal:  If you have MATLAB, you need to indicate both the MATLAB location and version. For example, on GNU/Linux:  ./configure --with-matlab=/usr/local/MATLAB/R2013a MATLAB_VERSION=8.1 ./configure --with-matlab=/usr/local/MATLAB/R2019a MATLAB_VERSION=9.6  Note that the MATLAB version can also be specified via the MATLAB family product release (R2009a, R2008b, ...). Note that the MATLAB version can also be specified via the MATLAB family product release (R2019a, R2018b, …). Alternatively, you can disable the compilation of MEX files for MATLAB with the --disable-matlab flag, and MEX files for Octave with --disable-octave. You may need to specify additional options to the configure script, see the platform specific instructions below. You may need to specify additional options to the configure script, see the output of the --help option, and also the platform specific instructions below. Note that if you don't want to compile the C/C++ programs with debugging information, you can specify the CFLAGS and CXXFLAGS variables to the configure script, such as:  ./configure CFLAGS="-O3" CXXFLAGS="-O3"  To remove debugging information for MATLAB MEX functions, the analagous call would be: To remove debugging information for MATLAB MEX functions, the analogous call would be:  ./configure MATLAB_MEX_CFLAGS="-O3" MATLAB_MEX_CXXFLAGS="-O3"  If you want to give a try to the parallelized versions of some mex files (A_times_B_kronecker_C and sparse_hessian_times_B_kronecker_C used to get the reduced form of the second order approximation of the model) you can add the --enable-openmp flag, for instance:  ./configure --with-matlab=/usr/local/MATLAB/R2013a MATLAB_VERSION=8.1 --enable-openmp If the configuration goes well, the script will tell you which components are correctly configured and will be built. Note that it is possible that some MEX files cannot be compiled, due to missing build dependencies. If you find no way of installing the missing dependencies, a workaround can be to give up on compiling these MEX files and rather use slower implementations (in the MATLAB/Octave language) that are available under the matlab/missing/mex/ subdirectories. For example, if you fail to compile the gensylv MEX, you can type the following at the MATLAB/Octave prompt before running Dynare: matlab addpath /matlab/missing/mex/gensylv  If the configuration goes well, the script will tell you which components are correctly configured and will be built. (where you need to replace  with the full path to your Dynare copy). ### Building ... ... @@ -151,7 +168,7 @@ The Git source comes with unit tests (in the MATLAB functions) and integration t make check  In the tests subfolder. If Dynare has been compiled against MATLAB and Octave, the tests will be run with MATLAB and Octave. Depending on your PC, this can take several hours. It is possible to run the tests only with MATLAB: the performance of your machine, this can take several hours. It is possible to run the tests only with MATLAB:  make check-matlab  ... ... @@ -230,7 +247,8 @@ apt install build-essential gfortran liboctave-dev libboost-graph-dev libboost-f ## Windows - Install [MSYS2](http://www.msys2.org) (pick the 64-bit version) - Install [MSYS2](http://www.msys2.org) (pick the 64-bit version, unless you have a 32-bit Windows, in which case see below) - Run a MSYS MinGW 64-bit shell - Update the system:  ... ...
 ... ... @@ -136,6 +136,7 @@ AC_CONFIG_FILES([Makefile dynare++/integ/testing/Makefile dynare++/kord/Makefile dynare++/src/Makefile dynare++/tests/Makefile mex/sources/Makefile ]) ... ...
 ... ... @@ -2,7 +2,7 @@ EXTRA_DIST = source \ utils/dynare_dom.py \ utils/dynare_lex.py SRC = $(wildcard src/source/*.rst) SRC =$(wildcard source/*.rst) html-local: build/html/index.html ... ...
 ... ... @@ -44,7 +44,7 @@ description, please refer to the comments inside the files themselves. Multi-country RBC model with time to build, presented in *Backus, Kehoe and Kydland (1992)*. The file shows how to use Dynare’s macro-processor. macro processor. agtrend.mod ... ...
 ... ... @@ -7,7 +7,7 @@ Installation and configuration Software requirements ===================== Packaged versions of Dynare are available for Windows XP/Vista/7/8/10, Packaged versions of Dynare are available for Windows 7/8/10, Debian GNU/Linux _, Ubuntu_ and macOS 10.8 or later. Dynare should work on other systems, but some compilation steps are necessary in that case. ... ... @@ -197,8 +197,7 @@ Homebrew, type:: If you don’t want to type this command every time you run Octave, you can put it in a file called .octaverc in your home directory (under Windows this will generally be c:\Documents and Settings\USERNAME\ while under macOS it is (under Windows this will generally be c:\Users\USERNAME while under macOS it is /Users/USERNAME/). This file is run by Octave at every startup. ... ...
 ... ... @@ -141,22 +141,22 @@ by the dynare command. .. option:: savemacro[=FILENAME] Instructs dynare to save the intermediary file which is obtained after macro-processing (see :ref:macro-proc-lang); obtained after macro processing (see :ref:macro-proc-lang); the saved output will go in the file specified, or if no file is specified in FILENAME-macroexp.mod .. option:: onlymacro Instructs the preprocessor to only perform the macro-processing step, and stop just after. Mainly useful for debugging purposes or for using the macro-processor macro processing step, and stop just after. Useful for debugging purposes or for using the macro processor independently of the rest of Dynare toolbox. .. option:: nolinemacro Instructs the macro-preprocessor to omit line numbering Instructs the macro preprocessor to omit line numbering information in the intermediary .mod file created after the macro-processing step. Useful in conjunction with the macro processing step. Useful in conjunction with :opt:savemacro  when one wants that to reuse the intermediary .mod file, without having it cluttered by line numbering directives. ... ... @@ -321,7 +321,7 @@ by the dynare command. .. option:: -I<> Defines a path to search for files to be included by the macroprocessor (using the @#include command). Multiple macro processor (using the @#include command). Multiple -I flags can be passed on the command line. The paths will be searched in the order that the -I flags are passed and the first matching file will be used. The flags passed here ... ...
 ... ... @@ -17,9 +17,7 @@ computations. On Linux and macOS, the default location of the configuration file is $HOME/.dynare, while on Windows it is %APPDATA%\dynare.ini (typically C:\Documents and Settings\USERNAME\Application Data\dynare.ini under Windows XP, or C:\Users\USERNAME\AppData\dynare.ini under Windows Vista/7/8). You (typically c:\Users\USERNAME\AppData\dynare.ini). You can specify a non standard location using the conffile option of the dynare command (see :ref:dyn-invoc). ... ... This diff is collapsed.  ... ... @@ -74,7 +74,7 @@ below. Basic operations can be performed on dates: **horzcat operator ([,])** Concatenates dates objects without removing repetitions. For instance [1950Q1, 1950Q2] is a a dates object with two instance [1950Q1, 1950Q2] is a dates object with two elements (1950Q1 and 1950Q2). **vertcat operator ([;])** ... ... @@ -132,7 +132,7 @@ he would extract some elements from a vector in Matlab/Octave. Let a the two last elements of a (by instantiating the dates object [1950Q4, 1951Q1]). Remark Dynare substitutes any occurrence of dates in the .mod file Remark: Dynare substitutes any occurrence of dates in the .mod file into an instantiation of the dates class regardless of the context. For instance, d = 1950Q1 will be translated as d = dates('1950Q1');. This automatic substitution can lead to a crash if ... ... @@ -776,7 +776,7 @@ The dseries class |br| The Matlab/Octave dseries class handles time series data. As any Matlab/Octave statements, this class can be used in a Dynare’s mod file. A dseries object has eight members: Dynare’s mod file. A dseries object has six members: :arg name: A nobs*1 cell of strings or a nobs*p character array, the names of the variables. ... ... @@ -784,6 +784,8 @@ The dseries class array, the tex names of the variables. :arg dates dates: An object with nobs elements, the dates of the sample. :arg double data: A nobs by vobs array, the data. :arg ops: The history of operations on the variables. :arg tags: The user-defined tags on the variables. data, name, tex are private members. The following constructors are available: ... ... @@ -791,7 +793,7 @@ The dseries class .. construct:: dseries () dseries (INITIAL_DATE) |br| Instantiates an empty dseries object, with, if |br| Instantiates an empty dseries object with, if defined, an initial date given by the single element dates object *INITIAL_DATE.* ... ... @@ -967,7 +969,7 @@ The dseries class deterministic_trend = .1*transpose(1:200); x = zeros(200,1); for i=2:200 x(i) = .75*x(i-1) + e(i); x(i) = .75*x(i-1) + u(i); end y = x + stochastic_trend + deterministic_trend; ... ... @@ -1075,7 +1077,7 @@ The dseries class 6Y | 64 7Y | 128 >> ts3 = ts1.cumsum(dates('3Y')); >> ts3 = ts1.cumprod(dates('3Y')); >> ts3 ts3 is a dseries object: ... ... @@ -1089,7 +1091,7 @@ The dseries class 6Y | 8 7Y | 16 >> ts4 = ts1.cumsum(dates('3Y'),dseries(pi)); >> ts4 = ts1.cumprod(dates('3Y'),dseries(pi)); >> ts4 ts4 is a dseries object: ... ... @@ -1351,7 +1353,7 @@ The dseries class deterministic_trend = .1*transpose(1:200); x = zeros(200,1); for i=2:200 x(i) = .75*x(i-1) + e(i); x(i) = .75*x(i-1) + u(i); end y = x + stochastic_trend + deterministic_trend; ... ...  SUBDIRS = utils/cc sylv parser/cc tl doc integ kord src EXTRA_DIST = tests SUBDIRS = utils/cc sylv parser/cc tl doc integ kord src tests  EXTRA_DIST = dynare++-ramsey.tex dynare++-tutorial.tex changelog-old.html EXTRA_DIST = \ dynare++-ramsey.tex \ dynare++-tutorial.tex \ sylvester.tex \ tl.tex \ changelog-old.html \ changelog-sylv-old.html if HAVE_PDFLATEX pdf-local: dynare++-ramsey.pdf dynare++-tutorial.pdf pdf-local: dynare++-ramsey.pdf dynare++-tutorial.pdf sylvester.pdf tl.pdf endif %.pdf: %.tex ... ...  //$Id: precalc_quadrature.dat 431 2005-08-16 15:41:01Z kamenik $// Copyright 2005, Ondra Kamenik /* * Copyright © 2005 Ondra Kamenik * Copyright © 2019 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 . */ // The file contains one dimensional quadrature points and weights for // a few quadratures. The format of data is clear. There is a class // OneDPrecalcQuadrature which implements an interface OneDQuadrature // using the data of this format. // The file contains one dimensional quadrature points and weights for a few // quadratures. The format of data is clear. There is a class // OneDPrecalcQuadrature which implements an interface OneDQuadrature using the // data of this format. // Gauss-Hermite quadrature; prefix gh // number of levels // Number of levels static const int gh_num_levels = 26; // number of points in each level // Number of points in each level static const int gh_num_points[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 30, 32, 40, 50, 60, 64 }; // weights, starting with the first level // Weights, starting with the first level static const double gh_weights[] = { // weights 1 = sqrt(pi) // weights 1 = √π 1.77245385090551588191942755656782537698745727539062, // weights 2 0.886226925452758013649083741671e+00, ... ... @@ -533,7 +550,7 @@ static const double gh_weights[] = { 0.553570653584e-48 }; // points, starting with the first level // Points, starting with the first level static const double gh_points[] = { // points 1 0.0, ... ... @@ -1051,16 +1068,16 @@ static const double gh_points[] = { // Gauss-Legendre quadrature; prefix gl // number of levels // Number of levels static const int gl_num_levels = 22; // number of points in each level // Number of points in each level static const int gl_num_points[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 32, 64 }; // weights, starting with the first level // Weights, starting with the first level static const double gl_weights[] = { // weight 1 2.0e+00, ... ... @@ -1392,7 +1409,7 @@ static const double gl_weights[] = { 0.178328072169643294729607914497e-02 }; // points, starting with the first level // Points, starting with the first level static const double gl_points[] = { // points 1 0.0e+00, ... ...  // Copyright 2005, Ondra Kamenik /* * Copyright © 2005 Ondra Kamenik * Copyright © 2019 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 . */ #include "product.hh" #include "symmetry.hh" ... ... @@ -6,7 +24,7 @@ #include #include /* This constructs a product iterator corresponding to index$(j0,0\ldots,0)$. */ /* This constructs a product iterator corresponding to index (j0,0,…,0). */ prodpit::prodpit(const ProductQuadrature &q, int j0, int l) : prodq(q), level(l), npoints(q.uquad.numPoints(l)), ... ... @@ -33,7 +51,6 @@ prodpit::operator==(const prodpit &ppit) const prodpit & prodpit::operator++() { // todo: throw if |prodq==NULL| or |jseq==NULL| or |sig==NULL| or |end_flag==true| int i = prodq.dimen()-1; jseq[i]++; while (i >= 0 && jseq[i] == npoints) ... ... @@ -59,8 +76,6 @@ prodpit::operator++() void prodpit::setPointAndWeight() { // todo: raise if |prodq==NULL| or |jseq==NULL| or |sig==NULL| or // |p==NULL| or |end_flag==true| w = 1.0; for (int i = 0; i < prodq.dimen(); i++) { ... ... @@ -89,26 +104,26 @@ prodpit::print() const ProductQuadrature::ProductQuadrature(int d, const OneDQuadrature &uq) : QuadratureImpl(d), uquad(uq) { // todo: check |d>=1| // TODO: check d≥1 } /* This calls |prodpit| constructor to return an iterator which points approximatelly at |ti|-th portion out of |tn| portions. First we find /* This calls prodpit constructor to return an iterator which points approximatelly at ‘ti’-th portion out of ‘tn’ portions. First we find out how many points are in the level, and then construct an interator$(j0,0,\ldots,0)$where$j0=$|ti*npoints/tn|. */ (j0,0,…,0) where j0=ti·npoints/tn. */ prodpit ProductQuadrature::begin(int ti, int tn, int l) const { // todo: raise is |l  // Copyright 2005, Ondra Kamenik /* * Copyright © 2005 Ondra Kamenik * Copyright © 2019 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 . */ // Product quadrature. /* This file defines a product multidimensional quadrature. If$Q_k$denotes the one dimensional quadrature, then the product quadrature$Q$of$k$level and dimension$d$takes the form $$Qf=\sum_{i_1=1}^{n_k}\ldots\sum_{i_d=1}^{n^k}w_{i_1}\cdot\ldots\cdot w_{i_d} f(x_{i_1},\ldots,x_{i_d})$$ which can be written in terms of the one dimensional quadrature$Q_k$as $$Qf=(Q_k\otimes\ldots\otimes Q_k)f$$ /* This file defines a product multidimensional quadrature. If Qₖ$ denotes the one dimensional quadrature, then the product quadrature Q of k level and dimension d takes the form: Here we define the product quadrature iterator |prodpit| and plug it into |QuadratureImpl| to obtains |ProductQuadrature|. */ nₖ nₖ Qf = ∑ … ∑ w_i₁·…·w_{i_d} f(x_i₁,…,x_{i_d}) i₁=1 i_d=1 which can be written in terms of the one dimensional quadrature Qₖ as Qf=(Qₖ⊗…⊗Qₖ)f Here we define the product quadrature iterator prodpit and plug it into QuadratureImpl to obtains ProductQuadrature. */ #ifndef PRODUCT_H #define PRODUCT_H ... ... @@ -20,18 +42,18 @@ #include "vector_function.hh" #include "quadrature.hh" /* This defines a product point iterator. We have to maintain the following: a pointer to product quadrature in order to know the dimension and the underlying one dimensional quadrature, then level, number of points in the level, integer sequence of indices, signal, the coordinates of the point and the weight. /* This defines a product point iterator. We have to maintain the following: a pointer to product quadrature in order to know the dimension and the underlying one dimensional quadrature, then level, number of points in the level, integer sequence of indices, signal, the coordinates of the point and the weight. The point indices, signal, and point coordinates are implmented as pointers in order to allow for empty constructor. The point indices, signal, and point coordinates are implmented as pointers in order to allow for empty constructor. The constructor |prodpit(const ProductQuadrature& q, int j0, int l)| constructs an iterator pointing to $(j0,0,\ldots,0)$, which is used by |begin| dictated by |QuadratureImpl|. */ The constructor prodpit(const ProductQuadrature& q, int j0, int l) constructs an iterator pointing to (j0,0,…,0), which is used by begin() dictated by QuadratureImpl. */ class ProductQuadrature; ... ... @@ -79,13 +101,12 @@ protected: void setPointAndWeight(); }; /* The product quadrature is just |QuadratureImpl| with the product iterator plugged in. The object is constructed by just giving the underlying one dimensional quadrature, and the dimension. The only extra method is |designLevelForEvals| which for the given maximum number of evaluations (and dimension and underlying quadrature from the object) returns a maximum level yeilding number of evaluations less than the given number. */ /* The product quadrature is just QuadratureImpl with the product iterator plugged in. The object is constructed by just giving the underlying one dimensional quadrature, and the dimension. The only extra method is designLevelForEvals() which for the given maximum number of evaluations (and dimension and underlying quadrature from the object) returns a maximum level yeilding number of evaluations less than the given number. */ class ProductQuadrature : public QuadratureImpl { ... ...
 // Copyright 2005, Ondra Kamenik /* * Copyright © 2005 Ondra Kamenik * Copyright © 2019 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 . */ #include "quasi_mcarlo.hh" ... ... @@ -7,9 +25,9 @@ #include #include /* Here in the constructor, we have to calculate a maximum length of |coeff| array for a given |base| and given maximum |maxn|. After allocation, we calculate the coefficients. */ /* Here in the constructor, we have to calculate a maximum length of ‘coeff’ array for a given ‘base’ and given maximum ‘maxn’. After allocation, we calculate the coefficients. */ RadicalInverse::RadicalInverse(int n, int b, int mxn) : num(n), base(b), maxn(mxn), ... ... @@ -26,23 +44,26 @@ RadicalInverse::RadicalInverse(int n, int b, int mxn) while (nr > 0); } /* This evaluates the radical inverse. If there was no permutation, we have to calculate $${c_0\over b}+{c_1\over b^2}+\ldots+{c_j\over b^{j+1}}$$ which is evaluated as $$\left(\ldots\left(\left({c_j\over b}\cdot{1\over b}+{c_{j-1}\over b}\right) \cdot{1\over b}+{c_{j-2}\over b}\right) \ldots\right)\cdot{1\over b}+{c_0\over b}$$ Now with permutation $\pi$, we have $$\left(\ldots\left(\left({\pi(c_j)\over b}\cdot{1\over b}+ {\pi(c_{j-1})\over b}\right)\cdot{1\over b}+ {\pi(c_{j-2})\over b}\right) \ldots\right)\cdot{1\over b}+{\pi(c_0)\over b}$$ */ /* This evaluates the radical inverse. If there was no permutation, we have to calculate: c₀ c₁ cⱼ ── + ── + … + ──── b b² bʲ⁺¹ which is evaluated as: ⎛ ⎛⎛cⱼ 1 cⱼ₋₁⎞ 1 cⱼ₋₂⎞ ⎞ 1 c₀ ⎢…⎢⎢──·─ + ────⎥·─ + ────⎥…⎥·─ + ── ⎝ ⎝⎝ b b b ⎠ b b ⎠ ⎠ b b Now with permutation π, we have: ⎛ ⎛⎛π(cⱼ) 1 π(cⱼ₋₁)⎞ 1 π(cⱼ₋₂)⎞ ⎞ 1 π(c₀) ⎢…⎢⎢─────·─ + ───────⎥·─ + ───────⎥…⎥·─ + ───── ⎝ ⎝⎝ b b b ⎠ b b ⎠ ⎠ b b */ double RadicalInverse::eval(const PermutationScheme &p) const ... ... @@ -61,7 +82,7 @@ RadicalInverse::eval(const PermutationScheme &p) const void RadicalInverse::increase() { // todo: raise if |num+1 > maxn| // TODO: raise if num+1 > maxn num++; int i = 0; coeff[i]++; ... ... @@ -82,8 +103,8 @@ RadicalInverse::print() const coeff.print(); } /* Here we have the first 170 primes. This means that we are not able to integrate dimensions greater than 170. */ /* Here we have the first 170 primes. This means that we are not able to integrate dimensions greater than 170. */ std::array HaltonSequence::primes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, ... ... @@ -105,21 +126,21 @@ std::array HaltonSequence::primes = { 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013 }; /* This takes first |dim| primes and constructs |dim| radical inverses and calls |eval|. */ /* This takes first ‘dim’ primes and constructs ‘dim’ radical inverses and calls eval(). */ HaltonSequence::HaltonSequence(int n, int mxn, int dim, const PermutationScheme &p) : num(n), maxn(mxn), per(p), pt(dim) { // todo: raise if |dim > num_primes| // todo: raise if |n > mxn| // TODO: raise if dim > num_primes // TODO: raise if n > mxn for (int i = 0; i < dim; i++) ri.emplace_back(num, primes[i], maxn); eval(); } /* This calls |RadicalInverse::increase| for all radical inverses and calls |eval|. */ /* This calls RadicalInverse::increase() for all radical inverses and calls eval(). */ void HaltonSequence::increase() ... ... @@ -131,7 +152,7 @@ HaltonSequence::increase() eval(); } /* This sets point |pt| to radical inverse evaluations in each dimension. */ /* This sets point ‘pt’ to radical inverse evaluations in each dimension. */ void HaltonSequence::eval() { ... ... @@ -169,7 +190,6 @@ qmcpit::operator==(const qmcpit &qpit) const qmcpit & qmcpit::operator++() { // todo: raise if |halton == null || qmcq == NULL| halton.increase(); return *this; } ... ... @@ -180,14 +200,12 @@ qmcpit::weight() const return 1.0/spec.level(); } /* Clear from code. */ int WarnockPerScheme::permute(int i, int base, int c) const { return (c+i) % base; } /* Clear from code. */ int ReversePerScheme::permute(int i, int base, int c) const { ... ...
 // Copyright 2005, Ondra Kamenik /* * Copyright © 2005 Ondra Kamenik * Copyright © 2019 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 . */ // Quasi Monte Carlo quadrature. /* This defines quasi Monte Carlo quadratures for cube and for a function multiplied by normal density. The quadrature for a cube is named |QMCarloCubeQuadrature| and integrates: $$\int_{\langle 0,1\rangle^n}f(x){\rm d}x$$ The quadrature for a function of normally distributed parameters is named |QMCarloNormalQuadrature| and integrates: $${1\over\sqrt{(2\pi)^n}}\int_{(-\infty,\infty)^n}f(x)e^{-{1\over 2}x^Tx}{\rm d}x$$ For a cube we define |qmcpit| as iterator of |QMCarloCubeQuadrature|, and for the normal density multiplied function we define |qmcnpit| as iterator of |QMCarloNormalQuadrature|. The quasi Monte Carlo method generates low discrepancy points with equal weights. The one dimensional low discrepancy sequences are generated by |RadicalInverse| class, the sequences are combined for higher dimensions by |HaltonSequence| class. The Halton sequence can use a permutation scheme; |PermutattionScheme| is an abstract class for all permutaton schemes. We have three implementations: |WarnockPerScheme|, |ReversePerScheme|, and |IdentityPerScheme|. */ QMCarloCubeQuadrature and integrates: ∫ f(x)dx [0,1]ⁿ The quadrature for a function of normally distributed parameters is named QMCarloNormalQuadrature and integrates: 1 ──────── ∫ f(x)e^{−½xᵀx}dx √{(2π)ⁿ} [−∞,+∞]ⁿ For a cube we define qmcpit as iterator of QMCarloCubeQuadrature, and for the normal density multiplied function we define qmcnpit as iterator of QMCarloNormalQuadrature. The quasi Monte Carlo method generates low discrepancy points with equal weights. The one dimensional low discrepancy sequences are generated by RadicalInverse class, the sequences are combined for higher dimensions by HaltonSequence class. The Halton sequence can use a permutation scheme; PermutattionScheme is an abstract class for all permutaton schemes. We have three implementations: WarnockPerScheme, ReversePerScheme, and IdentityPerScheme. */ #ifndef QUASI_MCARLO_H #define QUASI_MCARLO_H ... ... @@ -32,9 +56,9 @@ #include /* This abstract class declares |permute| method which permutes coefficient |c| having index of |i| fro the base |base| and returns the permuted coefficient which must be in $0,\ldots,base-1$. */ /* This abstract class declares permute() method which permutes coefficient ‘c’ having index of ‘i’ fro the base ‘base’ and returns the permuted coefficient which must be in 0,…,base−1. */ class PermutationScheme { ... ... @@ -44,15 +68,14 @@ public: virtual int permute(int i, int base, int c) const = 0; }; /* This class represents an integer number |num| as $c_0+c_1b+c_2b^2+\ldots+c_jb^j$, where $b$ is |base| and $c_0,\ldots,c_j$ is stored in |coeff|. The size of |IntSequence| coeff does not grow with growing |num|, but is fixed from the very beginning and is set according to supplied maximum |maxn|. /* This class represents an integer number ‘num’ as c₀+c₁b+c₂b²+…+cⱼbʲ, where b is ‘base’ and c₀,…,cⱼ is stored in ‘coeff’. The size of IntSequence coeff does not grow with growing ‘num’, but is fixed from the very beginning and is set according to supplied maximum ‘maxn’. The basic method is |eval| which evaluates the |RadicalInverse| with a given permutation scheme and returns the point, and |increase| which increases |num| and recalculates the coefficients. */ The basic method is eval() which evaluates the RadicalInverse with a given permutation scheme and returns the point, and increase() which increases ‘num’ and recalculates the coefficients. */ class RadicalInverse { ... ... @@ -70,11 +93,11 @@ public: void print() const; }; /* This is a vector of |RadicalInverse|s, each |RadicalInverse| has a different prime as its base. The static members |primes| and |num_primes| define a precalculated array of primes. The |increase| method of the class increases indices in all |RadicalInverse|s and sets point |pt| to contain the points in each dimension. */ /* This is a vector of RadicalInverses, each RadicalInverse has a different prime as its base. The static members ‘primes’ and ‘num_primes’ define a precalculated array of primes. The increase() method of the class increases indices in all RadicalInverse’s and sets point ‘pt’ to contain the points in each dimension. */ class HaltonSequence { ... ... @@ -106,10 +129,9 @@ protected: void eval(); }; /* This is a specification of quasi Monte Carlo quadrature. It consists of dimension |dim|, number of points (or level) |lev|, and the permutation scheme. This class is common to all quasi Monte Carlo classes. */ /* This is a specification of quasi Monte Carlo quadrature. It consists of dimension ‘dim’, number of points (or level) ‘lev’, and the permutation scheme. This class is common to all quasi Monte Carlo classes. */ class QMCSpecification { ... ... @@ -140,10 +162,10 @@ public: } }; /* This is an iterator for quasi Monte Carlo over a cube |QMCarloCubeQuadrature|. The iterator maintains |HaltonSequence| of the same dimension as given by the specification. An iterator can be constructed from a given number |n|, or by a copy constructor. */ /* This is an iterator for quasi Monte Carlo over a cube QMCarloCubeQuadrature. The iterator maintains HaltonSequence of the same dimension as given by the specification. An iterator can be constructed from a given number ‘n’, or by a copy constructor. */ class qmcpit { ... ... @@ -181,10 +203,9 @@ public: } }; /* This is an easy declaration of quasi Monte Carlo quadrature for a cube. Everything important has been done in its iterator |qmcpit|, so we only inherit from general |Quadrature| and reimplement |begin| and |numEvals|. */ /* This is an easy declaration of quasi Monte Carlo quadrature for a cube. Everything important has been done in its iterator qmcpit, so we only inherit from general Quadrature and reimplement begin() and numEvals(). */ class QMCarloCubeQuadrature : public QuadratureImpl, public QMCSpecification { ... ...
 // Copyright 2005, Ondra Kamenik /* * Copyright © 2005 Ondra Kamenik * Copyright © 2019 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 . */ #include "smolyak.hh" #include "symmetry.hh" ... ... @@ -6,9 +24,9 @@ #include #include /* This constructs a beginning of |isum| summand in |smolq|. We must be careful here, since |isum| can be past-the-end, so no reference to vectors in |smolq| by |isum| must be done in this case. */ /* This constructs a beginning of ‘isum’ summand in ‘smolq’. We must be careful here, since ‘isum’ can be past-the-end, so no reference to vectors in ‘smolq’ by ‘isum’ must be done in this case. */ smolpit::smolpit(const SmolyakQuadrature &q, unsigned int isum) : smolq(q), isummand(isum), ... ... @@ -26,14 +44,13 @@ smolpit::operator==(const smolpit &spit) const return &smolq == &spit.smolq && isummand == spit.isummand && jseq == spit.jseq; } /* We first try to increase index within the current summand. If we are at maximum, we go to a subsequent summand. Note that in this case all indices in |jseq| will be zero, so no change is needed. */ /* We first try to increase index within the current summand. If we are at maximum, we go to a subsequent summand. Note that in this case all indices in ‘jseq’ will be zero, so no change is needed. */ smolpit & smolpit::operator++() { // todo: throw if |smolq==NULL| or |jseq==NULL| or |sig==NULL| const IntSequence &levpts = smolq.levpoints[isummand]; int i = smolq.dimen()-1; jseq[i]++; ... ... @@ -55,14 +72,13 @@ smolpit::operator++() return *this; } /* Here we set the point coordinates according to |jseq| and |isummand|. Also the weight is set here. */ /* Here we set the point coordinates according to ‘jseq’ and ‘isummand’. Also the weight is set here. */ void smolpit::setPointAndWeight() { // todo: raise if |smolq==NULL| or |jseq==NULL| or |sig==NULL| or // |p==NULL| or |isummand>=smolq.numSummands()| // todo: raise if isummand ≥ smolq.numSummands() int l = smolq.level; int d = smolq.dimen(); int sumk = (smolq.levels[isummand]).sum(); ... ... @@ -96,21 +112,19 @@ smolpit::print() const std::cout.flags(ff); } /* Here is the constructor of |SmolyakQuadrature|. We have to setup |levels|, |levpoints| and |cumevals|. We have to go through all $d$-dimensional sequences $k$, such that $l\leq \vert k\vert\leq l+d-1$ and all $k_i$ are positive integers. This is equivalent to going through all $k$ such that $l-d\leq\vert k\vert\leq l-1$ and all $k_i$ are non-negative integers. This is equivalent to going through $d+1$ dimensional sequences $(k,x)$ such that $\vert(k,x)\vert =l-1$ and $x=0,\ldots,d-1$. The resulting sequence of positive integers is obtained by adding $1$ to all $k_i$. */ /* Here is the constructor of SmolyakQuadrature. We have to setup ‘levels’, ‘levpoints’ and ‘cumevals’. We have to go through all d-dimensional sequences k, such that l≤|k|≤l+d−1 and all kᵢ are positive integers. This is equivalent to going through all k such that l−d≤|k|≤l−1 and all kᵢ are non-negative integers. This is equivalent to going through d+1 dimensional sequences (k,x) such that |(k,x)|=l−1 and x=0,…,d−1. The resulting sequence of positive integers is obtained by adding 1 to all kᵢ. */ SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature &uq) : QuadratureImpl(d), level(l), uquad(uq) { // todo: check |l>1|, |l>=d| // todo: check |l>=uquad.miLevel()|, |l<=uquad.maxLevel()| // TODO: check l>1, l≥d // TODO: check l≥uquad.miLevel(), l≤uquad.maxLevel() int cum = 0; for (const auto &si : SymmetrySet(l-1, d+1)) { ... ... @@ -129,10 +143,10 @@ SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature &uq) } } /* Here we return a number of evalutions of the quadrature for the given level. If the given level is the current one, we simply return the maximum cumulative number of evaluations. Otherwise we call costly |calcNumEvaluations| method. */ /* Here we return a number of evalutions of the quadrature for the given level. If the given level is the current one, we simply return the maximum cumulative number of evaluations. Otherwise we call costly calcNumEvaluations() method. */ int SmolyakQuadrature::numEvals(int l) const ... ... @@ -143,14 +157,14 @@ SmolyakQuadrature::numEvals(int l) const return cumevals[numSummands()-1]; } /* This divides all the evaluations to |tn| approximately equal groups, and returns the beginning of the specified group |ti|. The granularity of divisions are summands as listed by |levels|. */ /* This divides all the evaluations to ‘tn’ approximately equal groups, and returns the beginning of the specified group ‘ti’. The granularity of divisions are summands as listed by ‘levels’. */ smolpit SmolyakQuadrature::begin(int ti, int tn, int l) const { // todo: raise is |level!=l| // TODO: raise is level≠l if (ti == tn) return smolpit(*this, numSummands()); ... ... @@ -162,9 +176,9 @@ SmolyakQuadrature::begin(int ti, int tn, int l) const return smolpit(*this, isum); } /* This is the same in a structure as |@<|SmolyakQuadrature| constructor@>|. We have to go through all summands and calculate a number of evaluations in each summand. */ /* This is the same in a structure as SmolyakQuadrature constructor. We have to go through all summands and calculate a number of evaluations in each summand. */ int SmolyakQuadrature::calcNumEvaluations(int lev) const ... ... @@ -185,8 +199,8 @@ SmolyakQuadrature::calcNumEvaluations(int lev) const return cum; } /* This returns a maximum level such that the number of evaluations is less than the given number. */ /* This returns a maximum level such that the number of evaluations is less than the given number. */ void SmolyakQuadrature::designLevelForEvals(int max_evals, int &lev, int &evals) const ... ...