...
 
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: <https://www.dynare.org/>
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: <https://www.dynare.org/download/>.
# 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 <DYNARE_ROOT>/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 `<DYNARE_ROOT>` 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 <http://www.debian.org/>`_, `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 <savemacro[=FILENAME]>` 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<<path>>
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 <http://www.gnu.org/licenses/>.
*/
// 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 <http://www.gnu.org/licenses/>.
*/
#include "product.hh"
#include "symmetry.hh"
......@@ -6,7 +24,7 @@
#include <iostream>
#include <iomanip>
/* 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<prodpit>(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<dimen()|
// todo: check |l<=uquad.numLevels()|
// TODO: raise if l<dimen()
// TODO: check l ≤ uquad.numLevels()
int npoints = uquad.numPoints(l);
return prodpit(*this, ti*npoints/tn, l);
}
/* This just starts at the first level and goes to a higher level as
long as a number of evaluations (which is $n_k^d$ for $k$ being the
level) is less than the given number of evaluations. */
/* This just starts at the first level and goes to a higher level as long as a
number of evaluations (which is nₖᵈ for k being the level) is less than the
given number of evaluations. */
void
ProductQuadrature::designLevelForEvals(int max_evals, int &lev, int &evals) const
......@@ -125,5 +140,4 @@ ProductQuadrature::designLevelForEvals(int max_evals, int &lev, int &evals) cons
while (lev < uquad.numLevels()-2 && evals < max_evals);
lev--;
evals = last_evals;
}
// 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 <http://www.gnu.org/licenses/>.
*/
// 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<prodpit>
{
......
// 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 <http://www.gnu.org/licenses/>.
*/
#include "quadrature.hh"
#include "precalc_quadrature.hh"
......
This diff is collapsed.
// 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 <http://www.gnu.org/licenses/>.
*/
#include "quasi_mcarlo.hh"
......@@ -7,9 +25,9 @@
#include <iomanip>
#include <array>
/* 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<int, 170> HaltonSequence::primes = {
2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
......@@ -105,21 +126,21 @@ std::array<int, 170> 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 <http://www.gnu.org/licenses/>.
*/
// 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 <vector>
/* 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<qmcpit>, 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 <http://www.gnu.org/licenses/>.
*/
#include "smolyak.hh"
#include "symmetry.hh"
......@@ -6,9 +24,9 @@
#include <iostream>
#include <iomanip>
/* 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<smolpit>(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
......
// 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 <http://www.gnu.org/licenses/>.