diff --git a/.gitignore b/.gitignore
index a72fda3648983226ce1fb86436845d1c6ae2c624..6714f985fadbc92a2c8cb4be29e1816adfe9007b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -120,9 +120,6 @@ mex/build/matlab/run_m2html.m
/mex/octave/
# Dynare++
-/dynare++/integ/cc/*.cpp
-/dynare++/integ/cc/*.h
-/dynare++/integ/cc/main.tex
/dynare++/integ/src/quadrature-points
/dynare++/integ/src/quadrature-points.exe
/dynare++/integ/testing/tests
diff --git a/dynare++/integ/cc/Makefile.am b/dynare++/integ/cc/Makefile.am
index ef0102ef783a64e65349b53fd6621800ad910b3f..3db1f12b7645a4187145a6d85e708c6dca5cc664 100644
--- a/dynare++/integ/cc/Makefile.am
+++ b/dynare++/integ/cc/Makefile.am
@@ -1,54 +1,16 @@
-CWEBSRC = \
- product.cweb \
- quadrature.cweb \
- quasi_mcarlo.cweb \
- smolyak.cweb \
- vector_function.cweb \
- product.hweb \
- quadrature.hweb \
- quasi_mcarlo.hweb \
- smolyak.hweb \
- vector_function.hweb
-
-GENERATED_FILES = \
- product.cpp \
- quadrature.cpp \
- quasi_mcarlo.cpp \
- smolyak.cpp \
- vector_function.cpp \
- product.h \
- quadrature.h \
- quasi_mcarlo.h \
- smolyak.h \
- vector_function.h
-
noinst_LIBRARIES = libinteg.a
-libinteg_a_SOURCES = $(CWEBSRC) $(GENERATED_FILES) precalc_quadrature.dat
+libinteg_a_SOURCES = \
+ quadrature.cc \
+ quadrature.hh \
+ quasi_mcarlo.cc \
+ quasi_mcarlo.hh \
+ product.cc \
+ product.hh \
+ smolyak.cc \
+ smolyak.hh \
+ vector_function.cc \
+ vector_function.hh \
+ precalc_quadrature.dat
libinteg_a_CPPFLAGS = -I../../sylv/cc -I../../tl/cc -I$(top_srcdir)/mex/sources
libinteg_a_CXXFLAGS = $(AM_CXXFLAGS) $(PTHREAD_CFLAGS)
-
-BUILT_SOURCES = $(GENERATED_FILES)
-
-EXTRA_DIST = main.web dummy.ch
-
-%.cpp: %.cweb dummy.ch
- $(CTANGLE) -bhp $< dummy.ch $@
-
-%.h: %.hweb dummy.ch
- $(CTANGLE) -bhp $< dummy.ch $@
-
-if HAVE_CWEAVE
-if HAVE_PDFTEX
-if HAVE_EPLAIN
-pdf-local: integ.pdf
-
-integ.pdf: main.web $(CWEBSRC)
- $(CWEAVE) -bhp main.web
- $(PDFTEX) main
- mv main.pdf integ.pdf
-endif
-endif
-endif
-
-CLEANFILES = integ.pdf main.idx main.log main.scn main.tex main.toc
diff --git a/dynare++/integ/cc/dummy.ch b/dynare++/integ/cc/dummy.ch
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/dynare++/integ/cc/main.web b/dynare++/integ/cc/main.web
deleted file mode 100644
index 416d25417b308d43b512372bbf6c70c06d5318ca..0000000000000000000000000000000000000000
--- a/dynare++/integ/cc/main.web
+++ /dev/null
@@ -1,48 +0,0 @@
-@q $Id: main.web 2333 2009-01-14 10:32:55Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@q cwebmac.tex defines its own \ifpdf, which is incompatible with the @>
-@q \ifpdf defined by eplain, so undefine it @>
-\let\ifpdf\relax
-\input eplain
-
-@q now define \ifpdf to be always false: PDF macros of cwebmac are buggy @>
-\newif\ifpdf
-\iffalse\fi
-
-\def\title{{\mainfont Numerical Integration Module}}
-
-
-@i ../../c++lib.w
-@s Vector int
-@s ConstVector int
-@s IntSequence int
-@s GeneralMatrix int
-@s THREAD int
-@s THREAD_GROUP int
-@s SYNCHRO int
-
-\titletrue
-\null\vfill
-\centerline{\titlefont Numerical Integration Module}
-\vfill\vfill
-Copyright \copyright\ 2005 by Ondra Kamenik
-
-\penalty-10000
-
-@i vector_function.hweb
-@i vector_function.cweb
-
-@i quadrature.hweb
-@i quadrature.cweb
-
-@i product.hweb
-@i product.cweb
-
-@i smolyak.hweb
-@i smolyak.cweb
-
-@i quasi_mcarlo.hweb
-@i quasi_mcarlo.cweb
-
-
diff --git a/dynare++/integ/cc/product.cc b/dynare++/integ/cc/product.cc
new file mode 100644
index 0000000000000000000000000000000000000000..08a83a9a1ef0084b575c44e7ecef5bd383d278de
--- /dev/null
+++ b/dynare++/integ/cc/product.cc
@@ -0,0 +1,191 @@
+// Copyright 2005, Ondra Kamenik
+
+#include "product.hh"
+#include "symmetry.h"
+
+prodpit::prodpit()
+ : prodq(NULL), level(0), npoints(0), jseq(NULL),
+ end_flag(true), sig(NULL), p(NULL)
+{
+}
+
+/* This constructs a product iterator corresponding to index $(j0,0\ldots,0)$. */
+
+prodpit::prodpit(const ProductQuadrature &q, int j0, int l)
+ : prodq(&q), level(l), npoints(q.uquad.numPoints(l)), jseq(new IntSequence(q.dimen(), 0)),
+ end_flag(false), sig(new ParameterSignal(q.dimen())), p(new Vector(q.dimen()))
+{
+ if (j0 < npoints)
+ {
+ (*jseq)[0] = j0;
+ setPointAndWeight();
+ }
+ else
+ {
+ end_flag = true;
+ }
+}
+
+prodpit::prodpit(const prodpit &ppit)
+ : prodq(ppit.prodq), level(ppit.level), npoints(ppit.npoints),
+ end_flag(ppit.end_flag), w(ppit.w)
+{
+ if (ppit.jseq)
+ jseq = new IntSequence(*(ppit.jseq));
+ else
+ jseq = NULL;
+ if (ppit.sig)
+ sig = new ParameterSignal(*(ppit.sig));
+ else
+ sig = NULL;
+ if (ppit.p)
+ p = new Vector(*(ppit.p));
+ else
+ p = NULL;
+}
+
+prodpit::~prodpit()
+{
+ if (jseq)
+ delete jseq;
+ if (sig)
+ delete sig;
+ if (p)
+ delete p;
+}
+
+bool
+prodpit::operator==(const prodpit &ppit) const
+{
+ bool ret = true;
+ ret = ret & prodq == ppit.prodq;
+ ret = ret & end_flag == ppit.end_flag;
+ ret = ret & ((jseq == NULL && ppit.jseq == NULL)
+ || (jseq != NULL && ppit.jseq != NULL && *jseq == *(ppit.jseq)));
+ return ret;
+}
+
+const prodpit &
+prodpit::operator=(const prodpit &ppit)
+{
+ prodq = ppit.prodq;
+ end_flag = ppit.end_flag;
+ w = ppit.w;
+
+ if (jseq)
+ delete jseq;
+ if (sig)
+ delete sig;
+ if (p)
+ delete p;
+
+ if (ppit.jseq)
+ jseq = new IntSequence(*(ppit.jseq));
+ else
+ jseq = NULL;
+ if (ppit.sig)
+ sig = new ParameterSignal(*(ppit.sig));
+ else
+ sig = NULL;
+ if (ppit.p)
+ p = new Vector(*(ppit.p));
+ else
+ p = NULL;
+
+ return *this;
+}
+
+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)
+ {
+ (*jseq)[i] = 0;
+ i--;
+ if (i >= 0)
+ (*jseq)[i]++;
+ }
+ sig->signalAfter(std::max(i, 0));
+
+ if (i == -1)
+ end_flag = true;
+
+ if (!end_flag)
+ setPointAndWeight();
+
+ return *this;
+}
+
+/* This calculates the weight and sets point coordinates from the indices. */
+
+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++)
+ {
+ (*p)[i] = (prodq->uquad).point(level, (*jseq)[i]);
+ w *= (prodq->uquad).weight(level, (*jseq)[i]);
+ }
+}
+
+/* Debug print. */
+
+void
+prodpit::print() const
+{
+ printf("j=[");
+ for (int i = 0; i < prodq->dimen(); i++)
+ printf("%2d ", (*jseq)[i]);
+ printf("] %+4.3f*(", w);
+ for (int i = 0; i < prodq->dimen()-1; i++)
+ printf("%+4.3f ", (*p)[i]);
+ printf("%+4.3f)\n", (*p)[prodq->dimen()-1]);
+}
+
+ProductQuadrature::ProductQuadrature(int d, const OneDQuadrature &uq)
+ : QuadratureImpl<prodpit>(d), uquad(uq)
+{
+ // 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
+ out how many points are in the level, and then construct an interator
+ $(j0,0,\ldots,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()|
+ 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. */
+
+void
+ProductQuadrature::designLevelForEvals(int max_evals, int &lev, int &evals) const
+{
+ int last_evals;
+ evals = 1;
+ lev = 1;
+ do
+ {
+ lev++;
+ last_evals = evals;
+ evals = numEvals(lev);
+ }
+ while (lev < uquad.numLevels()-2 && evals < max_evals);
+ lev--;
+ evals = last_evals;
+
+}
diff --git a/dynare++/integ/cc/product.cweb b/dynare++/integ/cc/product.cweb
deleted file mode 100644
index 39a6e846bf0f95e40ab4de8fa32c1dc88c5365ca..0000000000000000000000000000000000000000
--- a/dynare++/integ/cc/product.cweb
+++ /dev/null
@@ -1,213 +0,0 @@
-@q $Id: product.cweb 431 2005-08-16 15:41:01Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@ This is {\tt product.cpp} file.
-
-@c
-#include "product.h"
-#include "symmetry.h"
-
-@<|prodpit| empty constructor@>;
-@<|prodpit| regular constructor@>;
-@<|prodpit| copy constructor@>;
-@<|prodpit| destructor@>;
-@<|prodpit::operator==| code@>;
-@<|prodpit::operator=| code@>;
-@<|prodpit::operator++| code@>;
-@<|prodpit::setPointAndWeight| code@>;
-@<|prodpit::print| code@>;
-@<|ProductQuadrature| constructor@>;
-@<|ProductQuadrature::begin| code@>;
-@<|ProductQuadrature::designLevelForEvals| code@>;
-
-@
-@<|prodpit| empty constructor@>=
-prodpit::prodpit()
- : prodq(NULL), level(0), npoints(0), jseq(NULL),
- end_flag(true), sig(NULL), p(NULL)
-{
-}
-
-@ This constructs a product iterator corresponding to index $(j0,0\ldots,0)$.
-@<|prodpit| regular constructor@>=
-prodpit::prodpit(const ProductQuadrature& q, int j0, int l)
- : prodq(&q), level(l), npoints(q.uquad.numPoints(l)), jseq(new IntSequence(q.dimen(), 0)),
- end_flag(false), sig(new ParameterSignal(q.dimen())), p(new Vector(q.dimen()))
-{
- if (j0 < npoints) {
- (*jseq)[0] = j0;
- setPointAndWeight();
- } else {
- end_flag = true;
- }
-}
-
-@ Copy constructor, clear.
-@<|prodpit| copy constructor@>=
-prodpit::prodpit(const prodpit& ppit)
- : prodq(ppit.prodq), level(ppit.level), npoints(ppit.npoints),
- end_flag(ppit.end_flag), w(ppit.w)
-{
- if (ppit.jseq)
- jseq = new IntSequence(*(ppit.jseq));
- else
- jseq = NULL;
- if (ppit.sig)
- sig = new ParameterSignal(*(ppit.sig));
- else
- sig = NULL;
- if (ppit.p)
- p = new Vector(*(ppit.p));
- else
- p = NULL;
-}
-
-@
-@<|prodpit| destructor@>=
-prodpit::~prodpit()
-{
- if (jseq)
- delete jseq;
- if (sig)
- delete sig;
- if (p)
- delete p;
-}
-
-@
-@<|prodpit::operator==| code@>=
-bool prodpit::operator==(const prodpit& ppit) const
-{
- bool ret = true;
- ret = ret & prodq == ppit.prodq;
- ret = ret & end_flag == ppit.end_flag;
- ret = ret & ((jseq==NULL && ppit.jseq==NULL) ||
- (jseq!=NULL && ppit.jseq!=NULL && *jseq == *(ppit.jseq)));
- return ret;
-}
-
-@
-@<|prodpit::operator=| code@>=
-const prodpit& prodpit::operator=(const prodpit& ppit)
-{
- prodq = ppit.prodq;
- end_flag = ppit.end_flag;
- w = ppit.w;
-
- if (jseq)
- delete jseq;
- if (sig)
- delete sig;
- if (p)
- delete p;
-
- if (ppit.jseq)
- jseq = new IntSequence(*(ppit.jseq));
- else
- jseq = NULL;
- if (ppit.sig)
- sig = new ParameterSignal(*(ppit.sig));
- else
- sig = NULL;
- if (ppit.p)
- p = new Vector(*(ppit.p));
- else
- p = NULL;
-
- return *this;
-}
-
-@
-@<|prodpit::operator++| code@>=
-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) {
- (*jseq)[i] = 0;
- i--;
- if (i >= 0)
- (*jseq)[i]++;
- }
- sig->signalAfter(std::max(i,0));
-
- if (i == -1)
- end_flag = true;
-
- if (! end_flag)
- setPointAndWeight();
-
- return *this;
-}
-
-
-@ This calculates the weight and sets point coordinates from the indices.
-@<|prodpit::setPointAndWeight| code@>=
-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++) {
- (*p)[i] = (prodq->uquad).point(level, (*jseq)[i]);
- w* = (prodq->uquad).weight(level, (*jseq)[i]);
- }
-}
-
-@ Debug print.
-@<|prodpit::print| code@>=
-void prodpit::print() const
-{
- printf("j=[");
- for (int i = 0; i < prodq->dimen(); i++)
- printf("%2d ", (*jseq)[i]);
- printf("] %+4.3f*(",w);
- for (int i = 0; i < prodq->dimen()-1; i++)
- printf("%+4.3f ", (*p)[i]);
- printf("%+4.3f)\n",(*p)[prodq->dimen()-1]);
-}
-
-@
-@<|ProductQuadrature| constructor@>=
-ProductQuadrature::ProductQuadrature(int d, const OneDQuadrature& uq)
- : QuadratureImpl<prodpit>(d), uquad(uq)
-{
- // 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
-out how many points are in the level, and then construct an interator
-$(j0,0,\ldots,0)$ where $j0=$|ti*npoints/tn|.
-
-@<|ProductQuadrature::begin| code@>=
-prodpit ProductQuadrature::begin(int ti, int tn, int l) const
-{
- // todo: raise is |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.
-
-@<|ProductQuadrature::designLevelForEvals| code@>=
-void ProductQuadrature::designLevelForEvals(int max_evals, int& lev, int& evals) const
-{
- int last_evals;
- evals = 1;
- lev = 1;
- do {
- lev++;
- last_evals = evals;
- evals = numEvals(lev);
- } while (lev < uquad.numLevels()-2 && evals < max_evals);
- lev--;
- evals = last_evals;
-
-}
-
-@ End of {\tt product.cpp} file
diff --git a/dynare++/integ/cc/product.hh b/dynare++/integ/cc/product.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3a3960424ad6e021bc308df5b85b5993b9a8cea9
--- /dev/null
+++ b/dynare++/integ/cc/product.hh
@@ -0,0 +1,112 @@
+// Copyright 2005, Ondra Kamenik
+
+// 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$$
+
+ Here we define the product quadrature iterator |prodpit| and plug it
+ into |QuadratureImpl| to obtains |ProductQuadrature|. */
+
+#ifndef PRODUCT_H
+#define PRODUCT_H
+
+#include "int_sequence.h"
+#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.
+
+ 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|. */
+
+class ProductQuadrature;
+
+class prodpit
+{
+protected:
+ const ProductQuadrature *prodq;
+ int level;
+ int npoints;
+ IntSequence *jseq;
+ bool end_flag;
+ ParameterSignal *sig;
+ Vector *p;
+ double w;
+public:
+ prodpit();
+ prodpit(const ProductQuadrature &q, int j0, int l);
+ prodpit(const prodpit &ppit);
+ ~prodpit();
+ bool operator==(const prodpit &ppit) const;
+ bool
+ operator!=(const prodpit &ppit) const
+ {
+ return !operator==(ppit);
+ }
+ const prodpit &operator=(const prodpit &spit);
+ prodpit &operator++();
+ const ParameterSignal &
+ signal() const
+ {
+ return *sig;
+ }
+ const Vector &
+ point() const
+ {
+ return *p;
+ }
+ double
+ weight() const
+ {
+ return w;
+ }
+ void print() const;
+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. */
+
+class ProductQuadrature : public QuadratureImpl<prodpit>
+{
+ friend class prodpit;
+ const OneDQuadrature &uquad;
+public:
+ ProductQuadrature(int d, const OneDQuadrature &uq);
+ virtual ~ProductQuadrature()
+ {
+ }
+ int
+ numEvals(int l) const
+ {
+ int res = 1;
+ for (int i = 0; i < dimen(); i++)
+ res *= uquad.numPoints(l);
+ return res;
+ }
+ void designLevelForEvals(int max_eval, int &lev, int &evals) const;
+protected:
+ prodpit begin(int ti, int tn, int level) const;
+};
+
+#endif
diff --git a/dynare++/integ/cc/product.hweb b/dynare++/integ/cc/product.hweb
deleted file mode 100644
index 59483a76ebb8dc34744af6e3700266dd8a418ad4..0000000000000000000000000000000000000000
--- a/dynare++/integ/cc/product.hweb
+++ /dev/null
@@ -1,107 +0,0 @@
-@q $Id: product.hweb 431 2005-08-16 15:41:01Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@*2 Product quadrature. This is {\tt product.h} file
-
-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$$
-
-Here we define the product quadrature iterator |prodpit| and plug it
-into |QuadratureImpl| to obtains |ProductQuadrature|.
-
-@s prodpit int
-@s ProductQuadrature int
-
-@c
-#ifndef PRODUCT_H
-#define PRODUCT_H
-
-#include "int_sequence.h"
-#include "vector_function.h"
-#include "quadrature.h"
-
-@<|prodpit| class declaration@>;
-@<|ProductQuadrature| class declaration@>;
-
-#endif
-
-@ 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 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|.
-
-@<|prodpit| class declaration@>=
-class ProductQuadrature;
-
-class prodpit {
-protected:@;
- const ProductQuadrature* prodq;
- int level;
- int npoints;
- IntSequence* jseq;
- bool end_flag;
- ParameterSignal* sig;
- Vector* p;
- double w;
-public:@;
- prodpit();
- prodpit(const ProductQuadrature& q, int j0, int l);
- prodpit(const prodpit& ppit);
- ~prodpit();
- bool operator==(const prodpit& ppit) const;
- bool operator!=(const prodpit& ppit) const
- {@+ return ! operator==(ppit);@+}
- const prodpit& operator=(const prodpit& spit);
- prodpit& operator++();
- const ParameterSignal& signal() const
- {@+ return *sig;@+}
- const Vector& point() const
- {@+ return *p;@+}
- double weight() const
- {@+ return w;@+}
- void print() const;
-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.
-
-@<|ProductQuadrature| class declaration@>=
-class ProductQuadrature : public QuadratureImpl<prodpit> {
- friend class prodpit;
- const OneDQuadrature& uquad;
-public:@;
- ProductQuadrature(int d, const OneDQuadrature& uq);
- virtual ~ProductQuadrature()@+ {}
- int numEvals(int l) const
- {
- int res = 1;
- for (int i = 0; i < dimen(); i++)
- res *= uquad.numPoints(l);
- return res;
- }
- void designLevelForEvals(int max_eval, int& lev, int& evals) const;
-protected:@;
- prodpit begin(int ti, int tn, int level) const;
-};
-
-@ End of {\tt product.h} file
diff --git a/dynare++/integ/cc/quadrature.cc b/dynare++/integ/cc/quadrature.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8a13cec556c1e7f121bb9a17fde8e1bad618547f
--- /dev/null
+++ b/dynare++/integ/cc/quadrature.cc
@@ -0,0 +1,53 @@
+// Copyright 2005, Ondra Kamenik
+
+#include "quadrature.hh"
+#include "precalc_quadrature.dat"
+
+#include <cmath>
+
+void
+OneDPrecalcQuadrature::calcOffsets()
+{
+ offsets[0] = 0;
+ for (int i = 1; i < num_levels; i++)
+ offsets[i] = offsets[i-1] + num_points[i-1];
+}
+
+GaussHermite::GaussHermite()
+ : OneDPrecalcQuadrature(gh_num_levels, gh_num_points, gh_weights, gh_points)
+{
+}
+
+GaussLegendre::GaussLegendre()
+ : OneDPrecalcQuadrature(gl_num_levels, gl_num_points, gl_weights, gl_points)
+{
+}
+
+/* Here we transform a draw from univariate $\langle 0,1\rangle$ to the
+ draw from Gaussina $N(0,1)$. This is done by a table lookup, the table
+ is given by |normal_icdf_step|, |normal_icfd_data|, |normal_icdf_num|,
+ and a number |normal_icdf_end|. In order to avoid wrong tails for lookups close
+ to zero or one, we rescale input |x| by $(1-2*(1-end))=2*end-1$. */
+
+double
+NormalICDF::get(double x)
+{
+ double xx = (2*normal_icdf_end-1)*std::abs(x-0.5);
+ int i = (int) floor(xx/normal_icdf_step);
+ double xx1 = normal_icdf_step*i;
+ double yy1 = normal_icdf_data[i];
+ double y;
+ if (i < normal_icdf_num-1)
+ {
+ double yy2 = normal_icdf_data[i+1];
+ y = yy1 + (yy2-yy1)*(xx-xx1)/normal_icdf_step;
+ }
+ else // this should never happen
+ {
+ y = yy1;
+ }
+ if (x > 0.5)
+ return y;
+ else
+ return -y;
+}
diff --git a/dynare++/integ/cc/quadrature.cweb b/dynare++/integ/cc/quadrature.cweb
deleted file mode 100644
index bf78fc2afedbaf5341b2e5f7c514dc63ea6be888..0000000000000000000000000000000000000000
--- a/dynare++/integ/cc/quadrature.cweb
+++ /dev/null
@@ -1,63 +0,0 @@
-@q $Id: quadrature.cweb 431 2005-08-16 15:41:01Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@ This is {\tt quadrature.cpp} file.
-
-@c
-#include "quadrature.h"
-#include "precalc_quadrature.dat"
-
-#include <cmath>
-
-@<|OneDPrecalcQuadrature::calcOffsets| code@>;
-@<|GaussHermite| constructor code@>;
-@<|GaussLegendre| constructor code@>;
-@<|NormalICDF| get code@>;
-
-@
-@<|OneDPrecalcQuadrature::calcOffsets| code@>=
-void OneDPrecalcQuadrature::calcOffsets()
-{
- offsets[0] = 0;
- for (int i = 1; i < num_levels; i++)
- offsets[i] = offsets[i-1] + num_points[i-1];
-}
-
-@
-@<|GaussHermite| constructor code@>=
-GaussHermite::GaussHermite()
- : OneDPrecalcQuadrature(gh_num_levels, gh_num_points, gh_weights, gh_points)@+ {}
-
-@
-@<|GaussLegendre| constructor code@>=
-GaussLegendre::GaussLegendre()
- : OneDPrecalcQuadrature(gl_num_levels, gl_num_points, gl_weights, gl_points)@+ {}
-
-@ Here we transform a draw from univariate $\langle 0,1\rangle$ to the
-draw from Gaussina $N(0,1)$. This is done by a table lookup, the table
-is given by |normal_icdf_step|, |normal_icfd_data|, |normal_icdf_num|,
-and a number |normal_icdf_end|. In order to avoid wrong tails for lookups close
-to zero or one, we rescale input |x| by $(1-2*(1-end))=2*end-1$.
-
-@<|NormalICDF| get code@>=
-double NormalICDF::get(double x)
-{
- double xx = (2*normal_icdf_end-1)*std::abs(x-0.5);
- int i = (int)floor(xx/normal_icdf_step);
- double xx1 = normal_icdf_step*i;
- double yy1 = normal_icdf_data[i];
- double y;
- if (i < normal_icdf_num-1) {
- double yy2 = normal_icdf_data[i+1];
- y = yy1 + (yy2-yy1)*(xx-xx1)/normal_icdf_step;
- } else { // this should never happen
- y = yy1;
- }
- if (x > 0.5)
- return y;
- else
- return -y;
-}
-
-
-@ End of {\tt quadrature.cpp} file
diff --git a/dynare++/integ/cc/quadrature.hh b/dynare++/integ/cc/quadrature.hh
new file mode 100644
index 0000000000000000000000000000000000000000..73ee9d9f4214c2f5b17da1412cc9d8308347970e
--- /dev/null
+++ b/dynare++/integ/cc/quadrature.hh
@@ -0,0 +1,325 @@
+// Copyright 2005, Ondra Kamenik
+
+// Quadrature.
+
+/* This file defines an interface for one dimensional (non-nested) quadrature
+ |OneDQuadrature|, and a parent for all multi-dimensional
+ quadratures. This parent class |Quadrature| presents a general concept of
+ quadrature, this is
+ $$\int f(x){\rm d}x \approx\sum_{i=1}^N w_ix_i$$
+ The class |Quadrature| just declares this concept. The concept is
+ implemented by class |QuadratureImpl| which paralelizes the
+ summation. All implementations therefore wishing to use the parallel
+ implementation should inherit from |QuadratureImpl| and integration is
+ done.
+
+ The integration concept relies on a point iterator, which goes through
+ all $x_i$ and $w_i$ for $i=1,\ldots,N$. All the iterators must be able
+ to go through only a portion of the set $i=1,\ldots,N$. This enables
+ us to implement paralelism, for two threads for example, one iterator
+ goes from the beginning to the (approximately) half, and the other
+ goes from the half to the end.
+
+ Besides this concept of the general quadrature, this file defines also
+ one dimensional quadrature, which is basically a scheme of points and
+ weights for different levels. The class |OneDQuadrature| is a parent
+ of all such objects, the classes |GaussHermite| and |GaussLegendre|
+ are specific implementations for Gauss--Hermite and Gauss--Legendre
+ quadratures resp. */
+
+#ifndef QUADRATURE_H
+#define QUADRATURE_H
+
+#include <cstdlib>
+#include "vector_function.hh"
+#include "int_sequence.h"
+#include "sthread.h"
+
+/* This pure virtual class represents a concept of one-dimensional
+ (non-nested) quadrature. So, one dimensional quadrature must return
+ number of levels, number of points in a given level, and then a point
+ and a weight in a given level and given order. */
+
+class OneDQuadrature
+{
+public:
+ virtual ~OneDQuadrature()
+ {
+ }
+ virtual int numLevels() const = 0;
+ virtual int numPoints(int level) const = 0;
+ virtual double point(int level, int i) const = 0;
+ virtual double weight(int lelel, int i) const = 0;
+};
+
+/* This is a general concept of multidimensional quadrature. at this
+ general level, we maintain only a dimension, and declare virtual
+ functions for integration. The function take two forms; first takes a
+ constant |VectorFunction| as an argument, creates locally
+ |VectorFunctionSet| and do calculation, second one takes as an
+ argument |VectorFunctionSet|.
+
+ Part of the interface is a method returning a number of evaluations
+ for a specific level. Note two things: this assumes that the number of
+ evaluations is known apriori and thus it is not applicable for
+ adaptive quadratures, second for Monte Carlo type of quadrature, the
+ level is a number of evaluations. */
+
+class Quadrature
+{
+protected:
+ int dim;
+public:
+ Quadrature(int d) : dim(d)
+ {
+ }
+ virtual ~Quadrature()
+ {
+ }
+ int
+ dimen() const
+ {
+ return dim;
+ }
+ virtual void integrate(const VectorFunction &func, int level,
+ int tn, Vector &out) const = 0;
+ virtual void integrate(VectorFunctionSet &fs, int level, Vector &out) const = 0;
+ virtual int numEvals(int level) const = 0;
+};
+
+/* This is just an integration worker, which works over a given
+ |QuadratureImpl|. It also needs the function, level, a specification
+ of the subgroup of points, and output vector.
+
+ See |@<|QuadratureImpl| class declaration@>| for details. */
+
+template <typename _Tpit>
+class QuadratureImpl;
+
+template <typename _Tpit>
+class IntegrationWorker : public THREAD
+{
+ const QuadratureImpl<_Tpit> &quad;
+ VectorFunction &func;
+ int level;
+ int ti;
+ int tn;
+ Vector &outvec;
+public:
+ IntegrationWorker(const QuadratureImpl<_Tpit> &q, VectorFunction &f, int l,
+ int tii, int tnn, Vector &out)
+ : quad(q), func(f), level(l), ti(tii), tn(tnn), outvec(out)
+ {
+ }
+
+ /* This integrates the given portion of the integral. We obtain first
+ and last iterators for the portion (|beg| and |end|). Then we iterate
+ through the portion. and finally we add the intermediate result to the
+ result |outvec|.
+
+ This method just everything up as it is coming. This might be imply
+ large numerical errors, perhaps in future I will implement something
+ smarter. */
+
+ void
+ operator()()
+ {
+ _Tpit beg = quad.begin(ti, tn, level);
+ _Tpit end = quad.begin(ti+1, tn, level);
+ Vector tmpall(outvec.length());
+ tmpall.zeros();
+ Vector tmp(outvec.length());
+
+ // note that since beg came from begin, it has empty signal
+ // and first evaluation gets no signal
+ for (_Tpit run = beg; run != end; ++run)
+ {
+ func.eval(run.point(), run.signal(), tmp);
+ tmpall.add(run.weight(), tmp);
+ }
+
+ {
+ SYNCHRO syn(&outvec, "IntegrationWorker");
+ outvec.add(1.0, tmpall);
+ }
+ }
+};
+
+/* This is the class which implements the integration. The class is
+ templated by the iterator type. We declare a method |begin| returning
+ an iterator to the beginnning of the |ti|-th portion out of total |tn|
+ portions for a given level.
+
+ In addition, we define a method which saves all the points to a given
+ file. Only for debugging purposes. */
+
+template <typename _Tpit>
+class QuadratureImpl : public Quadrature
+{
+ friend class IntegrationWorker<_Tpit>;
+public:
+ QuadratureImpl(int d) : Quadrature(d)
+ {
+ }
+
+ /* Just fill a thread group with workes and run it. */
+ void
+ integrate(VectorFunctionSet &fs, int level, Vector &out) const
+ {
+ // todo: out.length()==func.outdim()
+ // todo: dim == func.indim()
+ out.zeros();
+ THREAD_GROUP gr;
+ for (int ti = 0; ti < fs.getNum(); ti++)
+ {
+ gr.insert(new IntegrationWorker<_Tpit>(*this, fs.getFunc(ti),
+ level, ti, fs.getNum(), out));
+ }
+ gr.run();
+ }
+ void
+ integrate(const VectorFunction &func,
+ int level, int tn, Vector &out) const
+ {
+ VectorFunctionSet fs(func, tn);
+ integrate(fs, level, out);
+ }
+
+ /* Just for debugging. */
+ void
+ savePoints(const char *fname, int level) const
+ {
+ FILE *fd;
+ if (NULL == (fd = fopen(fname, "w")))
+ {
+ // todo: raise
+ fprintf(stderr, "Cannot open file %s for writing.\n", fname);
+ exit(1);
+ }
+ _Tpit beg = begin(0, 1, level);
+ _Tpit end = begin(1, 1, level);
+ for (_Tpit run = beg; run != end; ++run)
+ {
+ fprintf(fd, "%16.12g", run.weight());
+ for (int i = 0; i < dimen(); i++)
+ fprintf(fd, "\t%16.12g", run.point()[i]);
+ fprintf(fd, "\n");
+ }
+ fclose(fd);
+ }
+
+ _Tpit
+ start(int level) const
+ {
+ return begin(0, 1, level);
+ }
+ _Tpit
+ end(int level) const
+ {
+ return begin(1, 1, level);
+ }
+protected:
+ virtual _Tpit begin(int ti, int tn, int level) const = 0;
+};
+
+/* This is only an interface to a precalculated data in file {\tt
+ precalc\_quadrature.dat} which is basically C coded static data. It
+ implements |OneDQuadrature|. The data file is supposed to define the
+ following data: number of levels, array of number of points at each
+ level, an array of weights and array of points. The both latter array
+ store data level by level. An offset for a specific level is stored in
+ |offsets| integer sequence.
+
+ The implementing subclasses just fill the necessary data from the
+ file, the rest is calculated here. */
+
+class OneDPrecalcQuadrature : public OneDQuadrature
+{
+ int num_levels;
+ const int *num_points;
+ const double *weights;
+ const double *points;
+ IntSequence offsets;
+public:
+ OneDPrecalcQuadrature(int nlevels, const int *npoints,
+ const double *wts, const double *pts)
+ : num_levels(nlevels), num_points(npoints),
+ weights(wts), points(pts), offsets(num_levels)
+ {
+ calcOffsets();
+ }
+ virtual ~OneDPrecalcQuadrature()
+ {
+ }
+ int
+ numLevels() const
+ {
+ return num_levels;
+ }
+ int
+ numPoints(int level) const
+ {
+ return num_points[level-1];
+ }
+ double
+ point(int level, int i) const
+ {
+ return points[offsets[level-1]+i];
+ }
+ double
+ weight(int level, int i) const
+ {
+ return weights[offsets[level-1]+i];
+ }
+protected:
+ void calcOffsets();
+};
+
+/* Just precalculated Gauss--Hermite quadrature. This quadrature integrates exactly integrals
+ $$\int_{-\infty}^{\infty} x^ke^{-x^2}{\rm d}x$$
+ for level $k$.
+
+ Note that if pluging this one-dimensional quadrature to product or
+ Smolyak rule in order to integrate a function $f$ through normally
+ distributed inputs, one has to wrap $f$ to
+ |GaussConverterFunction| and apply the product or Smolyak rule to the
+ new function.
+
+ Check {\tt precalc\_quadrature.dat} for available levels. */
+
+class GaussHermite : public OneDPrecalcQuadrature
+{
+public:
+ GaussHermite();
+};
+
+/* Just precalculated Gauss--Legendre quadrature. This quadrature integrates exactly integrals
+ $$\int_0^1x^k{\rm d}x$$
+ for level $k$.
+
+ Check {\tt precalc\_quadrature.dat} for available levels. */
+
+class GaussLegendre : public OneDPrecalcQuadrature
+{
+public:
+ GaussLegendre();
+};
+
+/* This is just an inverse cummulative density function of normal
+ distribution. Its only method |get| returns for a given number
+ $x\in(0,1)$ a number $y$ such that $P(z<y)=x$, where the probability
+ is taken over normal distribution $N(0,1)$.
+
+ Currently, the implementation is done by a table lookup which implies
+ that the tails had to be chopped off. This further implies that Monte
+ Carlo quadratures using this transformation to draw points from
+ multidimensional $N(0,I)$ fail to integrate with satisfactory
+ precision polynomial functions, for which the tails matter. */
+
+class NormalICDF
+{
+public:
+ static double get(double x);
+};
+
+#endif
diff --git a/dynare++/integ/cc/quadrature.hweb b/dynare++/integ/cc/quadrature.hweb
deleted file mode 100644
index 7aa22c2074d7bb3aee4e3608c428b3252b5f4be3..0000000000000000000000000000000000000000
--- a/dynare++/integ/cc/quadrature.hweb
+++ /dev/null
@@ -1,311 +0,0 @@
-@q $Id: quadrature.hweb 2269 2008-11-23 14:33:22Z michel $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@*2 Quadrature. This is {\tt quadrature.h} file
-
-This file defines an interface for one dimensional (non-nested) quadrature
-|OneDQuadrature|, and a parent for all multi-dimensional
-quadratures. This parent class |Quadrature| presents a general concept of
-quadrature, this is
-$$\int f(x){\rm d}x \approx\sum_{i=1}^N w_ix_i$$
-The class |Quadrature| just declares this concept. The concept is
-implemented by class |QuadratureImpl| which paralelizes the
-summation. All implementations therefore wishing to use the parallel
-implementation should inherit from |QuadratureImpl| and integration is
-done.
-
-The integration concept relies on a point iterator, which goes through
-all $x_i$ and $w_i$ for $i=1,\ldots,N$. All the iterators must be able
-to go through only a portion of the set $i=1,\ldots,N$. This enables
-us to implement paralelism, for two threads for example, one iterator
-goes from the beginning to the (approximately) half, and the other
-goes from the half to the end.
-
-Besides this concept of the general quadrature, this file defines also
-one dimensional quadrature, which is basically a scheme of points and
-weights for different levels. The class |OneDQuadrature| is a parent
-of all such objects, the classes |GaussHermite| and |GaussLegendre|
-are specific implementations for Gauss--Hermite and Gauss--Legendre
-quadratures resp.
-
-@s OneDQuadrature int
-@s Quadrature int
-@s IntegrationWorker int
-@s QuadratureImpl int
-@s OneDPrecalcQuadrature int
-@s GaussHermite int
-@s GaussLegendre int
-@s NormalICDF int
-@s _Tpit int
-
-@c
-#ifndef QUADRATURE_H
-#define QUADRATURE_H
-
-#include <cstdlib>
-#include "vector_function.h"
-#include "int_sequence.h"
-#include "sthread.h"
-
-@<|OneDQuadrature| class declaration@>;
-@<|Quadrature| class declaration@>;
-@<|IntegrationWorker| class declaration@>;
-@<|QuadratureImpl| class declaration@>;
-@<|OneDPrecalcQuadrature| class declaration@>;
-@<|GaussHermite| class declaration@>;
-@<|GaussLegendre| class declaration@>;
-@<|NormalICDF| class declaration@>;
-
-#endif
-
-@ This pure virtual class represents a concept of one-dimensional
-(non-nested) quadrature. So, one dimensional quadrature must return
-number of levels, number of points in a given level, and then a point
-and a weight in a given level and given order.
-
-@<|OneDQuadrature| class declaration@>=
-class OneDQuadrature {
-public:@;
- virtual ~OneDQuadrature()@+ {}
- virtual int numLevels() const =0;
- virtual int numPoints(int level) const =0;
- virtual double point(int level, int i) const =0;
- virtual double weight(int lelel, int i) const =0;
-};
-
-@ This is a general concept of multidimensional quadrature. at this
-general level, we maintain only a dimension, and declare virtual
-functions for integration. The function take two forms; first takes a
-constant |VectorFunction| as an argument, creates locally
-|VectorFunctionSet| and do calculation, second one takes as an
-argument |VectorFunctionSet|.
-
-Part of the interface is a method returning a number of evaluations
-for a specific level. Note two things: this assumes that the number of
-evaluations is known apriori and thus it is not applicable for
-adaptive quadratures, second for Monte Carlo type of quadrature, the
-level is a number of evaluations.
-
-@<|Quadrature| class declaration@>=
-class Quadrature {
-protected:@;
- int dim;
-public:@;
- Quadrature(int d) : dim(d)@+ {}
- virtual ~Quadrature()@+ {}
- int dimen() const
- {@+ return dim;@+}
- virtual void integrate(const VectorFunction& func, int level,
- int tn, Vector& out) const =0;
- virtual void integrate(VectorFunctionSet& fs, int level, Vector& out) const =0;
- virtual int numEvals(int level) const =0;
-};
-
-@ This is just an integration worker, which works over a given
-|QuadratureImpl|. It also needs the function, level, a specification
-of the subgroup of points, and output vector.
-
-See |@<|QuadratureImpl| class declaration@>| for details.
-
-@<|IntegrationWorker| class declaration@>=
-template <typename _Tpit>
-class QuadratureImpl;
-
-template <typename _Tpit>
-class IntegrationWorker : public THREAD {
- const QuadratureImpl<_Tpit>& quad;
- VectorFunction& func;
- int level;
- int ti;
- int tn;
- Vector& outvec;
-public:@;
- IntegrationWorker(const QuadratureImpl<_Tpit>& q, VectorFunction& f, int l,
- int tii, int tnn, Vector& out)
- : quad(q), func(f), level(l), ti(tii), tn(tnn), outvec(out) @+{}
- @<|IntegrationWorker::operator()()| code@>;
-};
-
-
-@ This integrates the given portion of the integral. We obtain first
-and last iterators for the portion (|beg| and |end|). Then we iterate
-through the portion. and finally we add the intermediate result to the
-result |outvec|.
-
-This method just everything up as it is coming. This might be imply
-large numerical errors, perhaps in future I will implement something
-smarter.
-
-@<|IntegrationWorker::operator()()| code@>=
-void operator()() {
- _Tpit beg = quad.begin(ti, tn, level);
- _Tpit end = quad.begin(ti+1, tn, level);
- Vector tmpall(outvec.length());
- tmpall.zeros();
- Vector tmp(outvec.length());
-
- // note that since beg came from begin, it has empty signal
- // and first evaluation gets no signal
- for (_Tpit run = beg; run != end; ++run) {
- func.eval(run.point(), run.signal(), tmp);
- tmpall.add(run.weight(), tmp);
- }
-
- {
- SYNCHRO@, syn(&outvec, "IntegrationWorker");
- outvec.add(1.0, tmpall);
- }
-}
-
-
-@ This is the class which implements the integration. The class is
-templated by the iterator type. We declare a method |begin| returning
-an iterator to the beginnning of the |ti|-th portion out of total |tn|
-portions for a given level.
-
-In addition, we define a method which saves all the points to a given
-file. Only for debugging purposes.
-
-@<|QuadratureImpl| class declaration@>=
-template <typename _Tpit>
-class QuadratureImpl : public Quadrature {
- friend class IntegrationWorker<_Tpit>;
-public:@;
- QuadratureImpl(int d) : Quadrature(d)@+ {}
- @<|QuadratureImpl::integrate| code@>;
- void integrate(const VectorFunction& func,
- int level, int tn, Vector& out) const {
- VectorFunctionSet fs(func, tn);
- integrate(fs, level, out);
- }
- @<|Quadrature::savePoints| code@>;
- _Tpit start(int level) const
- {@+ return begin(0,1,level);@+}
- _Tpit end(int level) const
- {@+ return begin(1,1,level);@+}
-protected:@;
- virtual _Tpit begin(int ti, int tn, int level) const =0;
-};
-
-@ Just fill a thread group with workes and run it.
-@<|QuadratureImpl::integrate| code@>=
-void integrate(VectorFunctionSet& fs, int level, Vector& out) const {
- // todo: out.length()==func.outdim()
- // todo: dim == func.indim()
- out.zeros();
- THREAD_GROUP@, gr;
- for (int ti = 0; ti < fs.getNum(); ti++) {
- gr.insert(new IntegrationWorker<_Tpit>(*this, fs.getFunc(ti),
- level, ti, fs.getNum(), out));
- }
- gr.run();
-}
-
-
-@ Just for debugging.
-@<|Quadrature::savePoints| code@>=
-void savePoints(const char* fname, int level) const
-{
- FILE* fd;
- if (NULL==(fd = fopen(fname,"w"))) {
- // todo: raise
- fprintf(stderr, "Cannot open file %s for writing.\n", fname);
- exit(1);
- }
- _Tpit beg = begin(0,1,level);
- _Tpit end = begin(1,1,level);
- for (_Tpit run = beg; run != end; ++run) {
- fprintf(fd, "%16.12g", run.weight());
- for (int i = 0; i < dimen(); i++)
- fprintf(fd, "\t%16.12g", run.point()[i]);
- fprintf(fd, "\n");
- }
- fclose(fd);
-}
-
-
-@ This is only an interface to a precalculated data in file {\tt
-precalc\_quadrature.dat} which is basically C coded static data. It
-implements |OneDQuadrature|. The data file is supposed to define the
-following data: number of levels, array of number of points at each
-level, an array of weights and array of points. The both latter array
-store data level by level. An offset for a specific level is stored in
-|offsets| integer sequence.
-
-The implementing subclasses just fill the necessary data from the
-file, the rest is calculated here.
-
-@<|OneDPrecalcQuadrature| class declaration@>=
-class OneDPrecalcQuadrature : public OneDQuadrature {
- int num_levels;
- const int* num_points;
- const double* weights;
- const double* points;
- IntSequence offsets;
-public:@;
- OneDPrecalcQuadrature(int nlevels, const int* npoints,
- const double* wts, const double* pts)
- : num_levels(nlevels), num_points(npoints),
- weights(wts), points(pts), offsets(num_levels)
- {@+ calcOffsets();@+}
- virtual ~OneDPrecalcQuadrature()@+ {}
- int numLevels() const
- {@+ return num_levels;@+}
- int numPoints(int level) const
- {@+ return num_points[level-1];@+}
- double point(int level, int i) const
- {@+ return points[offsets[level-1]+i];@+}
- double weight(int level, int i) const
- {@+ return weights[offsets[level-1]+i];@+}
-protected:@;
- void calcOffsets();
-};
-
-@ Just precalculated Gauss--Hermite quadrature. This quadrature integrates exactly integrals
-$$\int_{-\infty}^{\infty} x^ke^{-x^2}{\rm d}x$$
-for level $k$.
-
-Note that if pluging this one-dimensional quadrature to product or
-Smolyak rule in order to integrate a function $f$ through normally
-distributed inputs, one has to wrap $f$ to
-|GaussConverterFunction| and apply the product or Smolyak rule to the
-new function.
-
-Check {\tt precalc\_quadrature.dat} for available levels.
-
-@<|GaussHermite| class declaration@>=
-class GaussHermite : public OneDPrecalcQuadrature {
-public:@;
- GaussHermite();
-};
-
-@ Just precalculated Gauss--Legendre quadrature. This quadrature integrates exactly integrals
-$$\int_0^1x^k{\rm d}x$$
-for level $k$.
-
-Check {\tt precalc\_quadrature.dat} for available levels.
-
-@<|GaussLegendre| class declaration@>=
-class GaussLegendre : public OneDPrecalcQuadrature {
-public:@;
- GaussLegendre();
-};
-
-@ This is just an inverse cummulative density function of normal
-distribution. Its only method |get| returns for a given number
-$x\in(0,1)$ a number $y$ such that $P(z<y)=x$, where the probability
-is taken over normal distribution $N(0,1)$.
-
-Currently, the implementation is done by a table lookup which implies
-that the tails had to be chopped off. This further implies that Monte
-Carlo quadratures using this transformation to draw points from
-multidimensional $N(0,I)$ fail to integrate with satisfactory
-precision polynomial functions, for which the tails matter.
-
-@<|NormalICDF| class declaration@>=
-class NormalICDF {
-public:@;
- static double get(double x);
-};
-
-@ End of {\tt quadrature.h} file
diff --git a/dynare++/integ/cc/quasi_mcarlo.cc b/dynare++/integ/cc/quasi_mcarlo.cc
new file mode 100644
index 0000000000000000000000000000000000000000..131c398fe1d14d2ca12f47bf54f787a5ba8a0712
--- /dev/null
+++ b/dynare++/integ/cc/quasi_mcarlo.cc
@@ -0,0 +1,287 @@
+// Copyright 2005, Ondra Kamenik
+
+#include "quasi_mcarlo.hh"
+
+#include <cmath>
+
+/* 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),
+ coeff((int) (floor(log((double)maxn)/log((double)b))+2), 0)
+{
+ int nr = num;
+ j = -1;
+ do
+ {
+ j++;
+ coeff[j] = nr % base;
+ nr = nr / base;
+ }
+ 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}
+ $$ */
+
+double
+RadicalInverse::eval(const PermutationScheme &p) const
+{
+ double res = 0;
+ for (int i = j; i >= 0; i--)
+ {
+ int cper = p.permute(i, base, coeff[i]);
+ res = (cper + res)/base;
+ }
+ return res;
+}
+
+/* We just add 1 to the lowest coefficient and check for overflow with respect
+ to the base. */
+void
+RadicalInverse::increase()
+{
+ // todo: raise if |num+1 > maxn|
+ num++;
+ int i = 0;
+ coeff[i]++;
+ while (coeff[i] == base)
+ {
+ coeff[i] = 0;
+ coeff[++i]++;
+ }
+ if (i > j)
+ j = i;
+}
+
+/* Debug print. */
+void
+RadicalInverse::print() const
+{
+ printf("n=%d b=%d c=", num, base);
+ coeff.print();
+}
+
+/* Here we have the first 170 primes. This means that we are not able
+ to integrate dimensions greater than 170. */
+
+int HaltonSequence::num_primes = 170;
+int HaltonSequence::primes[] = {
+ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
+ 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
+ 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
+ 127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
+ 179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
+ 233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
+ 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
+ 353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
+ 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
+ 467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
+ 547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
+ 607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
+ 661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
+ 739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
+ 811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
+ 877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
+ 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013
+};
+
+/* 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|
+ for (int i = 0; i < dim; i++)
+ ri.push_back(RadicalInverse(num, primes[i], maxn));
+ eval();
+}
+
+const HaltonSequence &
+HaltonSequence::operator=(const HaltonSequence &hs)
+{
+ num = hs.num;
+ maxn = hs.maxn;
+ ri.clear();
+ for (unsigned int i = 0; i < hs.ri.size(); i++)
+ ri.push_back(RadicalInverse(hs.ri[i]));
+ pt = hs.pt;
+ return *this;
+}
+
+/* This calls |RadicalInverse::increase| for all radical inverses and
+ calls |eval|. */
+
+void
+HaltonSequence::increase()
+{
+ for (unsigned int i = 0; i < ri.size(); i++)
+ ri[i].increase();
+ num++;
+ if (num <= maxn)
+ eval();
+}
+
+/* This sets point |pt| to radical inverse evaluations in each dimension. */
+void
+HaltonSequence::eval()
+{
+ for (unsigned int i = 0; i < ri.size(); i++)
+ pt[i] = ri[i].eval(per);
+}
+
+/* Debug print. */
+void
+HaltonSequence::print() const
+{
+ for (unsigned int i = 0; i < ri.size(); i++)
+ ri[i].print();
+ printf("point=[ ");
+ for (unsigned int i = 0; i < ri.size(); i++)
+ printf("%7.6f ", pt[i]);
+ printf("]\n");
+}
+
+qmcpit::qmcpit()
+ : spec(NULL), halton(NULL), sig(NULL)
+{
+}
+
+qmcpit::qmcpit(const QMCSpecification &s, int n)
+ : spec(&s), halton(new HaltonSequence(n, s.level(), s.dimen(), s.getPerScheme())),
+ sig(new ParameterSignal(s.dimen()))
+{
+}
+
+qmcpit::qmcpit(const qmcpit &qpit)
+ : spec(qpit.spec), halton(NULL), sig(NULL)
+{
+ if (qpit.halton)
+ halton = new HaltonSequence(*(qpit.halton));
+ if (qpit.sig)
+ sig = new ParameterSignal(qpit.spec->dimen());
+}
+
+qmcpit::~qmcpit()
+{
+ if (halton)
+ delete halton;
+ if (sig)
+ delete sig;
+}
+
+bool
+qmcpit::operator==(const qmcpit &qpit) const
+{
+ return (spec == qpit.spec)
+ && ((halton == NULL && qpit.halton == NULL)
+ || (halton != NULL && qpit.halton != NULL && halton->getNum() == qpit.halton->getNum()));
+}
+
+const qmcpit &
+qmcpit::operator=(const qmcpit &qpit)
+{
+ spec = qpit.spec;
+ if (halton)
+ delete halton;
+ if (qpit.halton)
+ halton = new HaltonSequence(*(qpit.halton));
+ else
+ halton = NULL;
+ return *this;
+}
+
+qmcpit &
+qmcpit::operator++()
+{
+ // todo: raise if |halton == null || qmcq == NULL|
+ halton->increase();
+ return *this;
+}
+
+double
+qmcpit::weight() const
+{
+ return 1.0/spec->level();
+}
+
+qmcnpit::qmcnpit()
+ : qmcpit(), pnt(NULL)
+{
+}
+
+qmcnpit::qmcnpit(const QMCSpecification &s, int n)
+ : qmcpit(s, n), pnt(new Vector(s.dimen()))
+{
+}
+
+qmcnpit::qmcnpit(const qmcnpit &qpit)
+ : qmcpit(qpit), pnt(NULL)
+{
+ if (qpit.pnt)
+ pnt = new Vector(*(qpit.pnt));
+}
+
+qmcnpit::~qmcnpit()
+{
+ if (pnt)
+ delete pnt;
+}
+
+const qmcnpit &
+qmcnpit::operator=(const qmcnpit &qpit)
+{
+ qmcpit::operator=(qpit);
+ if (pnt)
+ delete pnt;
+ if (qpit.pnt)
+ pnt = new Vector(*(qpit.pnt));
+ else
+ pnt = NULL;
+ return *this;
+}
+
+/* Here we inccrease a point in Halton sequence ant then store images
+ of the points in |NormalICDF| function. */
+
+qmcnpit &
+qmcnpit::operator++()
+{
+ qmcpit::operator++();
+ for (int i = 0; i < halton->point().length(); i++)
+ (*pnt)[i] = NormalICDF::get(halton->point()[i]);
+ return *this;
+}
+
+/* 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
+{
+ return (base-c) % base;
+}
diff --git a/dynare++/integ/cc/quasi_mcarlo.cweb b/dynare++/integ/cc/quasi_mcarlo.cweb
deleted file mode 100644
index 55963e33bd6f61530c4c937295b37a1cb6380f07..0000000000000000000000000000000000000000
--- a/dynare++/integ/cc/quasi_mcarlo.cweb
+++ /dev/null
@@ -1,341 +0,0 @@
-@q $Id: quasi_mcarlo.cweb 431 2005-08-16 15:41:01Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@ This is {\tt quasi\_mcarlo.cpp} file.
-
-@c
-#include "quasi_mcarlo.h"
-
-#include <cmath>
-
-@<|RadicalInverse| constructor code@>;
-@<|RadicalInverse::eval| code@>;
-@<|RadicalInverse::increase| code@>;
-@<|RadicalInverse::print| code@>;
-@<|HaltonSequence| static data@>;
-@<|HaltonSequence| constructor code@>;
-@<|HaltonSequence::operator=| code@>;
-@<|HaltonSequence::increase| code@>;
-@<|HaltonSequence::eval| code@>;
-@<|HaltonSequence::print| code@>;
-@<|qmcpit| empty constructor code@>;
-@<|qmcpit| regular constructor code@>;
-@<|qmcpit| copy constructor code@>;
-@<|qmcpit| destructor@>;
-@<|qmcpit::operator==| code@>;
-@<|qmcpit::operator=| code@>;
-@<|qmcpit::operator++| code@>;
-@<|qmcpit::weight| code@>;
-@<|qmcnpit| empty constructor code@>;
-@<|qmcnpit| regular constructor code@>;
-@<|qmcnpit| copy constructor code@>;
-@<|qmcnpit| destructor@>;
-@<|qmcnpit::operator=| code@>;
-@<|qmcnpit::operator++| code@>;
-@<|WarnockPerScheme::permute| code@>;
-@<|ReversePerScheme::permute| code@>;
-
-@ 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| constructor code@>=
-RadicalInverse::RadicalInverse(int n, int b, int mxn)
- : num(n), base(b), maxn(mxn),
- coeff((int)(floor(log((double)maxn)/log((double)b))+2), 0)
-{
- int nr = num;
- j = -1;
- do {
- j++;
- coeff[j] = nr % base;
- nr = nr / base;
- } 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}
-$$
-
-
-@<|RadicalInverse::eval| code@>=
-double RadicalInverse::eval(const PermutationScheme& p) const
-{
- double res = 0;
- for (int i = j; i >= 0; i--) {
- int cper = p.permute(i, base, coeff[i]);
- res = (cper + res)/base;
- }
- return res;
-}
-
-@ We just add 1 to the lowest coefficient and check for overflow with respect to the base.
-@<|RadicalInverse::increase| code@>=
-void RadicalInverse::increase()
-{
- // todo: raise if |num+1 > maxn|
- num++;
- int i = 0;
- coeff[i]++;
- while (coeff[i] == base) {
- coeff[i] = 0;
- coeff[++i]++;
- }
- if (i > j)
- j = i;
-}
-
-@ Debug print.
-@<|RadicalInverse::print| code@>=
-void RadicalInverse::print() const
-{
- printf("n=%d b=%d c=", num, base);
- coeff.print();
-}
-
-@ Here we have the first 170 primes. This means that we are not able
-to integrate dimensions greater than 170.
-
-@<|HaltonSequence| static data@>=
-int HaltonSequence::num_primes = 170;
-int HaltonSequence::primes[] = {
- 2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
- 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
- 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
- 127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
- 179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
- 233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
- 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
- 353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
- 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
- 467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
- 547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
- 607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
- 661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
- 739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
- 811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
- 877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
- 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013
-};
-
-
-@ This takes first |dim| primes and constructs |dim| radical inverses
-and calls |eval|.
-
-@<|HaltonSequence| constructor code@>=
-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|
- for (int i = 0; i < dim; i++)
- ri.push_back(RadicalInverse(num, primes[i], maxn));
- eval();
-}
-
-@
-@<|HaltonSequence::operator=| code@>=
-const HaltonSequence& HaltonSequence::operator=(const HaltonSequence& hs)
-{
- num = hs.num;
- maxn = hs.maxn;
- ri.clear();
- for (unsigned int i = 0; i < hs.ri.size(); i++)
- ri.push_back(RadicalInverse(hs.ri[i]));
- pt = hs.pt;
- return *this;
-}
-
-
-
-@ This calls |RadicalInverse::increase| for all radical inverses and
-calls |eval|.
-
-@<|HaltonSequence::increase| code@>=
-void HaltonSequence::increase()
-{
- for (unsigned int i = 0; i < ri.size(); i++)
- ri[i].increase();
- num++;
- if (num <= maxn)
- eval();
-}
-
-@ This sets point |pt| to radical inverse evaluations in each dimension.
-@<|HaltonSequence::eval| code@>=
-void HaltonSequence::eval()
-{
- for (unsigned int i = 0; i < ri.size(); i++)
- pt[i] = ri[i].eval(per);
-}
-
-@ Debug print.
-@<|HaltonSequence::print| code@>=
-void HaltonSequence::print() const
-{
- for (unsigned int i = 0; i < ri.size(); i++)
- ri[i].print();
- printf("point=[ ");
- for (unsigned int i = 0; i < ri.size(); i++)
- printf("%7.6f ", pt[i]);
- printf("]\n");
-}
-
-@
-@<|qmcpit| empty constructor code@>=
-qmcpit::qmcpit()
- : spec(NULL), halton(NULL), sig(NULL)@+ {}
-
-@
-@<|qmcpit| regular constructor code@>=
-qmcpit::qmcpit(const QMCSpecification& s, int n)
- : spec(&s), halton(new HaltonSequence(n, s.level(), s.dimen(), s.getPerScheme())),
- sig(new ParameterSignal(s.dimen()))
-{
-}
-
-@
-@<|qmcpit| copy constructor code@>=
-qmcpit::qmcpit(const qmcpit& qpit)
- : spec(qpit.spec), halton(NULL), sig(NULL)
-{
- if (qpit.halton)
- halton = new HaltonSequence(*(qpit.halton));
- if (qpit.sig)
- sig = new ParameterSignal(qpit.spec->dimen());
-}
-
-@
-@<|qmcpit| destructor@>=
-qmcpit::~qmcpit()
-{
- if (halton)
- delete halton;
- if (sig)
- delete sig;
-}
-
-@
-@<|qmcpit::operator==| code@>=
-bool qmcpit::operator==(const qmcpit& qpit) const
-{
- return (spec == qpit.spec) &&
- ((halton == NULL && qpit.halton == NULL) ||
- (halton != NULL && qpit.halton != NULL && halton->getNum() == qpit.halton->getNum()));
-}
-
-@
-@<|qmcpit::operator=| code@>=
-const qmcpit& qmcpit::operator=(const qmcpit& qpit)
-{
- spec = qpit.spec;
- if (halton)
- delete halton;
- if (qpit.halton)
- halton = new HaltonSequence(*(qpit.halton));
- else
- halton = NULL;
- return *this;
-}
-
-
-@
-@<|qmcpit::operator++| code@>=
-qmcpit& qmcpit::operator++()
-{
- // todo: raise if |halton == null || qmcq == NULL|
- halton->increase();
- return *this;
-}
-
-@
-@<|qmcpit::weight| code@>=
-double qmcpit::weight() const
-{
- return 1.0/spec->level();
-}
-
-@
-@<|qmcnpit| empty constructor code@>=
-qmcnpit::qmcnpit()
- : qmcpit(), pnt(NULL)@+ {}
-
-@
-@<|qmcnpit| regular constructor code@>=
-qmcnpit::qmcnpit(const QMCSpecification& s, int n)
- : qmcpit(s, n), pnt(new Vector(s.dimen()))
-{
-}
-
-@
-@<|qmcnpit| copy constructor code@>=
-qmcnpit::qmcnpit(const qmcnpit& qpit)
- : qmcpit(qpit), pnt(NULL)
-{
- if (qpit.pnt)
- pnt = new Vector(*(qpit.pnt));
-}
-
-@
-@<|qmcnpit| destructor@>=
-qmcnpit::~qmcnpit()
-{
- if (pnt)
- delete pnt;
-}
-
-@
-@<|qmcnpit::operator=| code@>=
-const qmcnpit& qmcnpit::operator=(const qmcnpit& qpit)
-{
- qmcpit::operator=(qpit);
- if (pnt)
- delete pnt;
- if (qpit.pnt)
- pnt = new Vector(*(qpit.pnt));
- else
- pnt = NULL;
- return *this;
-}
-
-@ Here we inccrease a point in Halton sequence ant then store images
-of the points in |NormalICDF| function.
-
-@<|qmcnpit::operator++| code@>=
-qmcnpit& qmcnpit::operator++()
-{
- qmcpit::operator++();
- for (int i = 0; i < halton->point().length(); i++)
- (*pnt)[i] = NormalICDF::get(halton->point()[i]);
- return *this;
-}
-
-@ Clear from code.
-@<|WarnockPerScheme::permute| code@>=
-int WarnockPerScheme::permute(int i, int base, int c) const
-{
- return (c+i) % base;
-}
-
-@ Clear from code.
-@<|ReversePerScheme::permute| code@>=
-int ReversePerScheme::permute(int i, int base, int c) const
-{
- return (base-c) % base;
-}
-
-@ End of {\tt quasi\_mcarlo.cpp} file.
\ No newline at end of file
diff --git a/dynare++/integ/cc/quasi_mcarlo.hh b/dynare++/integ/cc/quasi_mcarlo.hh
new file mode 100644
index 0000000000000000000000000000000000000000..69e948930c6f1e7625b0eef5ce6d82689d7e1a58
--- /dev/null
+++ b/dynare++/integ/cc/quasi_mcarlo.hh
@@ -0,0 +1,334 @@
+// Copyright 2005, Ondra Kamenik
+
+// 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|. */
+
+#ifndef QUASI_MCARLO_H
+#define QUASI_MCARLO_H
+
+#include "int_sequence.h"
+#include "quadrature.hh"
+
+#include "Vector.h"
+
+#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$. */
+
+class PermutationScheme
+{
+public:
+ PermutationScheme()
+ {
+ }
+ virtual ~PermutationScheme()
+ {
+ }
+ 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|.
+
+ 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
+{
+ int num;
+ int base;
+ int maxn;
+ int j;
+ IntSequence coeff;
+public:
+ RadicalInverse(int n, int b, int mxn);
+ RadicalInverse(const RadicalInverse &ri)
+ : num(ri.num), base(ri.base), maxn(ri.maxn), j(ri.j), coeff(ri.coeff)
+ {
+ }
+ const RadicalInverse &
+ operator=(const RadicalInverse &radi)
+ {
+ num = radi.num; base = radi.base; maxn = radi.maxn;
+ j = radi.j; coeff = radi.coeff;
+ return *this;
+ }
+ double eval(const PermutationScheme &p) const;
+ void increase();
+ 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. */
+
+class HaltonSequence
+{
+private:
+ static int primes[];
+ static int num_primes;
+protected:
+ int num;
+ int maxn;
+ vector<RadicalInverse> ri;
+ const PermutationScheme &per;
+ Vector pt;
+public:
+ HaltonSequence(int n, int mxn, int dim, const PermutationScheme &p);
+ HaltonSequence(const HaltonSequence &hs)
+ : num(hs.num), maxn(hs.maxn), ri(hs.ri), per(hs.per), pt(hs.pt)
+ {
+ }
+ const HaltonSequence &operator=(const HaltonSequence &hs);
+ void increase();
+ const Vector &
+ point() const
+ {
+ return pt;
+ }
+ const int
+ getNum() const
+ {
+ return num;
+ }
+ void print() const;
+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. */
+
+class QMCSpecification
+{
+protected:
+ int dim;
+ int lev;
+ const PermutationScheme &per_scheme;
+public:
+ QMCSpecification(int d, int l, const PermutationScheme &p)
+ : dim(d), lev(l), per_scheme(p)
+ {
+ }
+ virtual ~QMCSpecification()
+ {
+ }
+ int
+ dimen() const
+ {
+ return dim;
+ }
+ int
+ level() const
+ {
+ return lev;
+ }
+ const PermutationScheme &
+ getPerScheme() const
+ {
+ return per_scheme;
+ }
+};
+
+/* 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. For
+ technical reasons, there is also an empty constructor; for that
+ reason, every member is a pointer. */
+
+class qmcpit
+{
+protected:
+ const QMCSpecification *spec;
+ HaltonSequence *halton;
+ ParameterSignal *sig;
+public:
+ qmcpit();
+ qmcpit(const QMCSpecification &s, int n);
+ qmcpit(const qmcpit &qpit);
+ ~qmcpit();
+ bool operator==(const qmcpit &qpit) const;
+ bool
+ operator!=(const qmcpit &qpit) const
+ {
+ return !operator==(qpit);
+ }
+ const qmcpit &operator=(const qmcpit &qpit);
+ qmcpit &operator++();
+ const ParameterSignal &
+ signal() const
+ {
+ return *sig;
+ }
+ const Vector &
+ point() const
+ {
+ return halton->point();
+ }
+ double weight() const;
+ void
+ print() const
+ {
+ halton->print();
+ }
+};
+
+/* 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
+{
+public:
+ QMCarloCubeQuadrature(int d, int l, const PermutationScheme &p)
+ : QuadratureImpl<qmcpit>(d), QMCSpecification(d, l, p)
+ {
+ }
+ virtual ~QMCarloCubeQuadrature()
+ {
+ }
+ int
+ numEvals(int l) const
+ {
+ return l;
+ }
+protected:
+ qmcpit
+ begin(int ti, int tn, int lev) const
+ {
+ return qmcpit(*this, ti*level()/tn + 1);
+ }
+};
+
+/* This is an iterator for |QMCarloNormalQuadrature|. It is equivalent
+ to an iterator for quasi Monte Carlo cube quadrature but a point. The
+ point is obtained from a point of |QMCarloCubeQuadrature| by a
+ transformation $\langle
+ 0,1\rangle\rightarrow\langle-\infty,\infty\rangle$ aplied to all
+ dimensions. The transformation yields a normal distribution from a
+ uniform distribution on $\langle0,1\rangle$. It is in fact
+ |NormalICDF|. */
+
+class qmcnpit : public qmcpit
+{
+protected:
+ Vector *pnt;
+public:
+ qmcnpit();
+ qmcnpit(const QMCSpecification &spec, int n);
+ qmcnpit(const qmcnpit &qpit);
+ ~qmcnpit();
+ bool
+ operator==(const qmcnpit &qpit) const
+ {
+ return qmcpit::operator==(qpit);
+ }
+ bool
+ operator!=(const qmcnpit &qpit) const
+ {
+ return !operator==(qpit);
+ }
+ const qmcnpit &operator=(const qmcnpit &qpit);
+ qmcnpit &operator++();
+ const ParameterSignal &
+ signal() const
+ {
+ return *sig;
+ }
+ const Vector &
+ point() const
+ {
+ return *pnt;
+ }
+ void
+ print() const
+ {
+ halton->print(); pnt->print();
+ }
+};
+
+/* This is an easy declaration of quasi Monte Carlo quadrature for a
+ cube. Everything important has been done in its iterator |qmcnpit|, so
+ we only inherit from general |Quadrature| and reimplement |begin| and
+ |numEvals|. */
+
+class QMCarloNormalQuadrature : public QuadratureImpl<qmcnpit>, public QMCSpecification
+{
+public:
+ QMCarloNormalQuadrature(int d, int l, const PermutationScheme &p)
+ : QuadratureImpl<qmcnpit>(d), QMCSpecification(d, l, p)
+ {
+ }
+ virtual ~QMCarloNormalQuadrature()
+ {
+ }
+ int
+ numEvals(int l) const
+ {
+ return l;
+ }
+protected:
+ qmcnpit
+ begin(int ti, int tn, int lev) const
+ {
+ return qmcnpit(*this, ti*level()/tn + 1);
+ }
+};
+
+/* Declares Warnock permutation scheme. */
+class WarnockPerScheme : public PermutationScheme
+{
+public:
+ int permute(int i, int base, int c) const;
+};
+
+/* Declares reverse permutation scheme. */
+class ReversePerScheme : public PermutationScheme
+{
+public:
+ int permute(int i, int base, int c) const;
+};
+
+/* Declares no permutation (identity) scheme. */
+class IdentityPerScheme : public PermutationScheme
+{
+public:
+ int
+ permute(int i, int base, int c) const
+ {
+ return c;
+ }
+};
+
+#endif
diff --git a/dynare++/integ/cc/quasi_mcarlo.hweb b/dynare++/integ/cc/quasi_mcarlo.hweb
deleted file mode 100644
index d720d6bed5c42aefbb7786d8cd8e9caad30b990a..0000000000000000000000000000000000000000
--- a/dynare++/integ/cc/quasi_mcarlo.hweb
+++ /dev/null
@@ -1,286 +0,0 @@
-@q $Id: quasi_mcarlo.hweb 431 2005-08-16 15:41:01Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@*2 Quasi Monte Carlo quadrature. This is {\tt quasi\_mcarlo.h} file.
-
-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|.
-
-@s PermutationScheme int
-@s RadicalInverse int
-@s HaltonSequence int
-@s QMCSpecification int
-@s qmcpit int
-@s QMCarloCubeQuadrature int
-@s qmcnpit int
-@s QMCarloNormalQuadrature int
-@s WarnockPerScheme int
-@s ReversePerScheme int
-@s IdentityPerScheme int
-
-@c
-#ifndef QUASI_MCARLO_H
-#define QUASI_MCARLO_H
-
-#include "int_sequence.h"
-#include "quadrature.h"
-
-#include "Vector.h"
-
-#include <vector>
-
-@<|PermutationScheme| class declaration@>;
-@<|RadicalInverse| class declaration@>;
-@<|HaltonSequence| class declaration@>;
-@<|QMCSpecification| class declaration@>;
-@<|qmcpit| class declaration@>;
-@<|QMCarloCubeQuadrature| class declaration@>;
-@<|qmcnpit| class declaration@>;
-@<|QMCarloNormalQuadrature| class declaration@>;
-@<|WarnockPerScheme| class declaration@>;
-@<|ReversePerScheme| class declaration@>;
-@<|IdentityPerScheme| class declaration@>;
-
-#endif
-
-@ 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$.
-
-@<|PermutationScheme| class declaration@>=
-class PermutationScheme {
-public:@;
- PermutationScheme()@+ {}
- virtual ~PermutationScheme()@+ {}
- 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|.
-
-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.
-
-@<|RadicalInverse| class declaration@>=
-class RadicalInverse {
- int num;
- int base;
- int maxn;
- int j;
- IntSequence coeff;
-public:@;
- RadicalInverse(int n, int b, int mxn);
- RadicalInverse(const RadicalInverse& ri)
- : num(ri.num), base(ri.base), maxn(ri.maxn), j(ri.j), coeff(ri.coeff)@+ {}
- const RadicalInverse& operator=(const RadicalInverse& radi)
- {
- num = radi.num; base = radi.base; maxn = radi.maxn;
- j = radi.j; coeff = radi.coeff;
- return *this;
- }
- double eval(const PermutationScheme& p) const;
- void increase();
- 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.
-
-@<|HaltonSequence| class declaration@>=
-class HaltonSequence {
-private:@;
- static int primes[];
- static int num_primes;
-protected:@;
- int num;
- int maxn;
- vector<RadicalInverse> ri;
- const PermutationScheme& per;
- Vector pt;
-public:@;
- HaltonSequence(int n, int mxn, int dim, const PermutationScheme& p);
- HaltonSequence(const HaltonSequence& hs)
- : num(hs.num), maxn(hs.maxn), ri(hs.ri), per(hs.per), pt(hs.pt)@+ {}
- const HaltonSequence& operator=(const HaltonSequence& hs);
- void increase();
- const Vector& point() const
- {@+ return pt;@+}
- const int getNum() const
- {@+ return num;@+}
- void print() const;
-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.
-
-@<|QMCSpecification| class declaration@>=
-class QMCSpecification {
-protected:@;
- int dim;
- int lev;
- const PermutationScheme& per_scheme;
-public:@;
- QMCSpecification(int d, int l, const PermutationScheme& p)
- : dim(d), lev(l), per_scheme(p)@+ {}
- virtual ~QMCSpecification() {}
- int dimen() const
- {@+ return dim;@+}
- int level() const
- {@+ return lev;@+}
- const PermutationScheme& getPerScheme() const
- {@+ return per_scheme;@+}
-};
-
-
-@ 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. For
-technical reasons, there is also an empty constructor; for that
-reason, every member is a pointer.
-
-@<|qmcpit| class declaration@>=
-class qmcpit {
-protected:@;
- const QMCSpecification* spec;
- HaltonSequence* halton;
- ParameterSignal* sig;
-public:@;
- qmcpit();
- qmcpit(const QMCSpecification& s, int n);
- qmcpit(const qmcpit& qpit);
- ~qmcpit();
- bool operator==(const qmcpit& qpit) const;
- bool operator!=(const qmcpit& qpit) const
- {@+ return ! operator==(qpit);@+}
- const qmcpit& operator=(const qmcpit& qpit);
- qmcpit& operator++();
- const ParameterSignal& signal() const
- {@+ return *sig;@+}
- const Vector& point() const
- {@+ return halton->point();@+}
- double weight() const;
- void print() const
- {@+ halton->print();@+}
-};
-
-@ 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|.
-
-@<|QMCarloCubeQuadrature| class declaration@>=
-class QMCarloCubeQuadrature : public QuadratureImpl<qmcpit>, public QMCSpecification {
-public:@;
- QMCarloCubeQuadrature(int d, int l, const PermutationScheme& p)
- : QuadratureImpl<qmcpit>(d), QMCSpecification(d, l, p)@+ {}
- virtual ~QMCarloCubeQuadrature()@+ {}
- int numEvals(int l) const
- {@+ return l;@+}
-protected:@;
- qmcpit begin(int ti, int tn, int lev) const
- {@+ return qmcpit(*this, ti*level()/tn + 1);@+}
-};
-
-@ This is an iterator for |QMCarloNormalQuadrature|. It is equivalent
-to an iterator for quasi Monte Carlo cube quadrature but a point. The
-point is obtained from a point of |QMCarloCubeQuadrature| by a
-transformation $\langle
-0,1\rangle\rightarrow\langle-\infty,\infty\rangle$ aplied to all
-dimensions. The transformation yields a normal distribution from a
-uniform distribution on $\langle0,1\rangle$. It is in fact
-|NormalICDF|.
-
-@<|qmcnpit| class declaration@>=
-class qmcnpit : public qmcpit {
-protected:@;
- Vector* pnt;
-public:@;
- qmcnpit();
- qmcnpit(const QMCSpecification& spec, int n);
- qmcnpit(const qmcnpit& qpit);
- ~qmcnpit();
- bool operator==(const qmcnpit& qpit) const
- {@+ return qmcpit::operator==(qpit);@+}
- bool operator!=(const qmcnpit& qpit) const
- {@+ return ! operator==(qpit);@+}
- const qmcnpit& operator=(const qmcnpit& qpit);
- qmcnpit& operator++();
- const ParameterSignal& signal() const
- {@+ return *sig;@+}
- const Vector& point() const
- {@+ return *pnt;@+}
- void print() const
- {@+ halton->print();pnt->print();@+}
-};
-
-@ This is an easy declaration of quasi Monte Carlo quadrature for a
-cube. Everything important has been done in its iterator |qmcnpit|, so
-we only inherit from general |Quadrature| and reimplement |begin| and
-|numEvals|.
-
-@<|QMCarloNormalQuadrature| class declaration@>=
-class QMCarloNormalQuadrature : public QuadratureImpl<qmcnpit>, public QMCSpecification {
-public:@;
- QMCarloNormalQuadrature(int d, int l, const PermutationScheme& p)
- : QuadratureImpl<qmcnpit>(d), QMCSpecification(d, l, p)@+ {}
- virtual ~QMCarloNormalQuadrature()@+ {}
- int numEvals(int l) const
- {@+ return l;@+}
-protected:@;
- qmcnpit begin(int ti, int tn, int lev) const
- {@+ return qmcnpit(*this, ti*level()/tn + 1);@+}
-};
-
-@ Declares Warnock permutation scheme.
-@<|WarnockPerScheme| class declaration@>=
-class WarnockPerScheme : public PermutationScheme {
-public:@;
- int permute(int i, int base, int c) const;
-};
-
-@ Declares reverse permutation scheme.
-@<|ReversePerScheme| class declaration@>=
-class ReversePerScheme : public PermutationScheme {
-public:@;
- int permute(int i, int base, int c) const;
-};
-
-@ Declares no permutation (identity) scheme.
-@<|IdentityPerScheme| class declaration@>=
-class IdentityPerScheme : public PermutationScheme {
-public:@;
- int permute(int i, int base, int c) const
- {@+ return c;@+}
-};
-
-@ End of {\tt quasi\_mcarlo.h} file
diff --git a/dynare++/integ/cc/smolyak.cc b/dynare++/integ/cc/smolyak.cc
new file mode 100644
index 0000000000000000000000000000000000000000..21fc850b1e98a4f88686c48e70e7e966947e5f63
--- /dev/null
+++ b/dynare++/integ/cc/smolyak.cc
@@ -0,0 +1,269 @@
+// Copyright 2005, Ondra Kamenik
+
+#include "smolyak.hh"
+#include "symmetry.h"
+
+smolpit::smolpit()
+ : smolq(NULL), isummand(0), jseq(NULL), sig(NULL), p(NULL)
+{
+}
+
+/* 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), jseq(new IntSequence(q.dimen(), 0)),
+ sig(new ParameterSignal(q.dimen())), p(new Vector(q.dimen()))
+{
+ if (isummand < q.numSummands())
+ {
+ setPointAndWeight();
+ }
+}
+
+smolpit::smolpit(const smolpit &spit)
+ : smolq(spit.smolq), isummand(spit.isummand), w(spit.w)
+{
+ if (spit.jseq)
+ jseq = new IntSequence(*(spit.jseq));
+ else
+ jseq = NULL;
+ if (spit.sig)
+ sig = new ParameterSignal(*(spit.sig));
+ else
+ sig = NULL;
+ if (spit.p)
+ p = new Vector(*(spit.p));
+ else
+ p = NULL;
+}
+
+smolpit::~smolpit()
+{
+ if (jseq)
+ delete jseq;
+ if (sig)
+ delete sig;
+ if (p)
+ delete p;
+}
+
+bool
+smolpit::operator==(const smolpit &spit) const
+{
+ bool ret = true;
+ ret = ret & smolq == spit.smolq;
+ ret = ret & isummand == spit.isummand;
+ ret = ret & ((jseq == NULL && spit.jseq == NULL)
+ || (jseq != NULL && spit.jseq != NULL && *jseq == *(spit.jseq)));
+ return ret;
+}
+
+const smolpit &
+smolpit::operator=(const smolpit &spit)
+{
+ smolq = spit.smolq;
+ isummand = spit.isummand;
+ w = spit.w;
+
+ if (jseq)
+ delete jseq;
+ if (sig)
+ delete sig;
+ if (p)
+ delete p;
+
+ if (spit.jseq)
+ jseq = new IntSequence(*(spit.jseq));
+ else
+ jseq = NULL;
+ if (spit.sig)
+ sig = new ParameterSignal(*(spit.sig));
+ else
+ sig = NULL;
+ if (spit.p)
+ p = new Vector(*(spit.p));
+ else
+ p = NULL;
+
+ return *this;
+}
+
+/* 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]++;
+ while (i >= 0 && (*jseq)[i] == levpts[i])
+ {
+ (*jseq)[i] = 0;
+ i--;
+ if (i >= 0)
+ (*jseq)[i]++;
+ }
+ sig->signalAfter(std::max(i, 0));
+
+ if (i < 0)
+ isummand++;
+
+ if (isummand < smolq->numSummands())
+ setPointAndWeight();
+
+ return *this;
+}
+
+/* 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()|
+ int l = smolq->level;
+ int d = smolq->dimen();
+ int sumk = (smolq->levels[isummand]).sum();
+ int m1exp = l + d - sumk - 1;
+ w = (2*(m1exp/2) == m1exp) ? 1.0 : -1.0;
+ w *= smolq->psc.noverk(d-1, sumk-l);
+ for (int i = 0; i < d; i++)
+ {
+ int ki = (smolq->levels[isummand])[i];
+ (*p)[i] = (smolq->uquad).point(ki, (*jseq)[i]);
+ w *= (smolq->uquad).weight(ki, (*jseq)[i]);
+ }
+}
+
+/* Debug print. */
+void
+smolpit::print() const
+{
+ printf("isum=%-3d: [", isummand);
+ for (int i = 0; i < smolq->dimen(); i++)
+ printf("%2d ", (smolq->levels[isummand])[i]);
+ printf("] j=[");
+ for (int i = 0; i < smolq->dimen(); i++)
+ printf("%2d ", (*jseq)[i]);
+ printf("] %+4.3f*(", w);
+ for (int i = 0; i < smolq->dimen()-1; i++)
+ printf("%+4.3f ", (*p)[i]);
+ printf("%+4.3f)\n", (*p)[smolq->dimen()-1]);
+}
+
+/* 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$. */
+
+SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature &uq)
+ : QuadratureImpl<smolpit>(d), level(l), uquad(uq), psc(d-1, d-1)
+{
+ // todo: check |l>1|, |l>=d|
+ // todo: check |l>=uquad.miLevel()|, |l<=uquad.maxLevel()|
+ int cum = 0;
+ SymmetrySet ss(l-1, d+1);
+ for (symiterator si(ss); !si.isEnd(); ++si)
+ {
+ if ((*si)[d] <= d-1)
+ {
+ IntSequence lev((const IntSequence &)*si, 0, d);
+ lev.add(1);
+ levels.push_back(lev);
+ IntSequence levpts(d);
+ for (int i = 0; i < d; i++)
+ levpts[i] = uquad.numPoints(lev[i]);
+ levpoints.push_back(levpts);
+ cum += levpts.mult();
+ cumevals.push_back(cum);
+ }
+ }
+}
+
+/* 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
+{
+ if (l != level)
+ return calcNumEvaluations(l);
+ else
+ 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|. */
+
+smolpit
+SmolyakQuadrature::begin(int ti, int tn, int l) const
+{
+ // todo: raise is |level!=l|
+ if (ti == tn)
+ return smolpit(*this, numSummands());
+
+ int totevals = cumevals[numSummands()-1];
+ int evals = (totevals*ti)/tn;
+ unsigned int isum = 0;
+ while (isum+1 < numSummands() && cumevals[isum+1] < evals)
+ isum++;
+ 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. */
+
+int
+SmolyakQuadrature::calcNumEvaluations(int lev) const
+{
+ int cum = 0;
+ SymmetrySet ss(lev-1, dim+1);
+ for (symiterator si(ss); !si.isEnd(); ++si)
+ {
+ if ((*si)[dim] <= dim-1)
+ {
+ IntSequence lev((const IntSequence &)*si, 0, dim);
+ lev.add(1);
+ IntSequence levpts(dim);
+ for (int i = 0; i < dim; i++)
+ levpts[i] = uquad.numPoints(lev[i]);
+ cum += levpts.mult();
+ }
+ }
+ return cum;
+}
+
+/* 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
+{
+ int last_evals;
+ evals = 1;
+ lev = 1;
+ do
+ {
+ lev++;
+ last_evals = evals;
+ evals = calcNumEvaluations(lev);
+ }
+ while (lev < uquad.numLevels() && evals <= max_evals);
+ lev--;
+ evals = last_evals;
+}
diff --git a/dynare++/integ/cc/smolyak.cweb b/dynare++/integ/cc/smolyak.cweb
deleted file mode 100644
index 8e0a57374cb931e0eac5056fe2d850d5c43eb783..0000000000000000000000000000000000000000
--- a/dynare++/integ/cc/smolyak.cweb
+++ /dev/null
@@ -1,294 +0,0 @@
-@q $Id: smolyak.cweb 1208 2007-03-19 21:33:12Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@ This is {\tt smolyak.cpp} file.
-
-@c
-#include "smolyak.h"
-#include "symmetry.h"
-
-@<|smolpit| empty constructor@>;
-@<|smolpit| regular constructor@>;
-@<|smolpit| copy constructor@>;
-@<|smolpit| destructor@>;
-@<|smolpit::operator==| code@>;
-@<|smolpit::operator=| code@>;
-@<|smolpit::operator++| code@>;
-@<|smolpit::setPointAndWeight| code@>;
-@<|smolpit::print| code@>;
-@<|SmolyakQuadrature| constructor@>;
-@<|SmolyakQuadrature::numEvals| code@>;
-@<|SmolyakQuadrature::begin| code@>;
-@<|SmolyakQuadrature::calcNumEvaluations| code@>;
-@<|SmolyakQuadrature::designLevelForEvals| code@>;
-
-@
-@<|smolpit| empty constructor@>=
-smolpit::smolpit()
- : smolq(NULL), isummand(0), jseq(NULL), sig(NULL), p(NULL)
-{
-}
-
-@ 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| regular constructor@>=
-smolpit::smolpit(const SmolyakQuadrature& q, unsigned int isum)
- : smolq(&q), isummand(isum), jseq(new IntSequence(q.dimen(), 0)),
- sig(new ParameterSignal(q.dimen())), p(new Vector(q.dimen()))
-{
- if (isummand < q.numSummands()) {
- setPointAndWeight();
- }
-}
-
-@
-@<|smolpit| copy constructor@>=
-smolpit::smolpit(const smolpit& spit)
- : smolq(spit.smolq), isummand(spit.isummand), w(spit.w)
-{
- if (spit.jseq)
- jseq = new IntSequence(*(spit.jseq));
- else
- jseq = NULL;
- if (spit.sig)
- sig = new ParameterSignal(*(spit.sig));
- else
- sig = NULL;
- if (spit.p)
- p = new Vector(*(spit.p));
- else
- p = NULL;
-}
-
-@
-@<|smolpit| destructor@>=
-smolpit::~smolpit()
-{
- if (jseq)
- delete jseq;
- if (sig)
- delete sig;
- if (p)
- delete p;
-}
-
-@
-@<|smolpit::operator==| code@>=
-bool smolpit::operator==(const smolpit& spit) const
-{
- bool ret = true;
- ret = ret & smolq == spit.smolq;
- ret = ret & isummand == spit.isummand;
- ret = ret & ((jseq==NULL && spit.jseq==NULL) ||
- (jseq!=NULL && spit.jseq!=NULL && *jseq == *(spit.jseq)));
- return ret;
-}
-
-@
-@<|smolpit::operator=| code@>=
-const smolpit& smolpit::operator=(const smolpit& spit)
-{
- smolq = spit.smolq;
- isummand = spit.isummand;
- w = spit.w;
-
- if (jseq)
- delete jseq;
- if (sig)
- delete sig;
- if (p)
- delete p;
-
- if (spit.jseq)
- jseq = new IntSequence(*(spit.jseq));
- else
- jseq = NULL;
- if (spit.sig)
- sig = new ParameterSignal(*(spit.sig));
- else
- sig = NULL;
- if (spit.p)
- p = new Vector(*(spit.p));
- else
- p = NULL;
-
- return *this;
-}
-
-@ 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::operator++| code@>=
-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]++;
- while (i >= 0 && (*jseq)[i] == levpts[i]) {
- (*jseq)[i] = 0;
- i--;
- if (i >= 0)
- (*jseq)[i]++;
- }
- sig->signalAfter(std::max(i,0));
-
- if (i < 0)
- isummand++;
-
- if (isummand < smolq->numSummands())
- setPointAndWeight();
-
- return *this;
-}
-
-
-@ Here we set the point coordinates according to |jseq| and
-|isummand|. Also the weight is set here.
-
-@<|smolpit::setPointAndWeight| code@>=
-void smolpit::setPointAndWeight()
-{
- // todo: raise if |smolq==NULL| or |jseq==NULL| or |sig==NULL| or
- // |p==NULL| or |isummand>=smolq->numSummands()|
- int l = smolq->level;
- int d = smolq->dimen();
- int sumk = (smolq->levels[isummand]).sum();
- int m1exp = l + d - sumk - 1;
- w = (2*(m1exp/2)==m1exp)? 1.0 : -1.0;
- w *= smolq->psc.noverk(d-1, sumk-l);
- for (int i = 0; i < d; i++) {
- int ki = (smolq->levels[isummand])[i];
- (*p)[i] = (smolq->uquad).point(ki, (*jseq)[i]);
- w *= (smolq->uquad).weight(ki, (*jseq)[i]);
- }
-}
-
-@ Debug print.
-@<|smolpit::print| code@>=
-void smolpit::print() const
-{
- printf("isum=%-3d: [", isummand);
- for (int i = 0; i < smolq->dimen(); i++)
- printf("%2d ", (smolq->levels[isummand])[i]);
- printf("] j=[");
- for (int i = 0; i < smolq->dimen(); i++)
- printf("%2d ", (*jseq)[i]);
- printf("] %+4.3f*(",w);
- for (int i = 0; i < smolq->dimen()-1; i++)
- printf("%+4.3f ", (*p)[i]);
- printf("%+4.3f)\n",(*p)[smolq->dimen()-1]);
-}
-
-@ 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$.
-
-@<|SmolyakQuadrature| constructor@>=
-SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature& uq)
- : QuadratureImpl<smolpit>(d), level(l), uquad(uq), psc(d-1,d-1)
-{
- // todo: check |l>1|, |l>=d|
- // todo: check |l>=uquad.miLevel()|, |l<=uquad.maxLevel()|
- int cum = 0;
- SymmetrySet ss(l-1, d+1);
- for (symiterator si(ss); !si.isEnd(); ++si) {
- if ((*si)[d] <= d-1) {
- IntSequence lev((const IntSequence&)*si, 0, d);
- lev.add(1);
- levels.push_back(lev);
- IntSequence levpts(d);
- for (int i = 0; i < d; i++)
- levpts[i] = uquad.numPoints(lev[i]);
- levpoints.push_back(levpts);
- cum += levpts.mult();
- cumevals.push_back(cum);
- }
- }
-}
-
-@ 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.
-
-@<|SmolyakQuadrature::numEvals| code@>=
-int SmolyakQuadrature::numEvals(int l) const
-{
- if (l != level)
- return calcNumEvaluations(l);
- else
- 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|.
-
-@<|SmolyakQuadrature::begin| code@>=
-smolpit SmolyakQuadrature::begin(int ti, int tn, int l) const
-{
- // todo: raise is |level!=l|
- if (ti == tn)
- return smolpit(*this, numSummands());
-
- int totevals = cumevals[numSummands()-1];
- int evals = (totevals*ti)/tn;
- unsigned int isum = 0;
- while (isum+1 < numSummands() && cumevals[isum+1] < evals)
- isum++;
- 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.
-
-@<|SmolyakQuadrature::calcNumEvaluations| code@>=
-int SmolyakQuadrature::calcNumEvaluations(int lev) const
-{
- int cum = 0;
- SymmetrySet ss(lev-1, dim+1);
- for (symiterator si(ss); !si.isEnd(); ++si) {
- if ((*si)[dim] <= dim-1) {
- IntSequence lev((const IntSequence&)*si, 0, dim);
- lev.add(1);
- IntSequence levpts(dim);
- for (int i = 0; i < dim; i++)
- levpts[i] = uquad.numPoints(lev[i]);
- cum += levpts.mult();
- }
- }
- return cum;
-}
-
-@ This returns a maximum level such that the number of evaluations is
-less than the given number.
-
-@<|SmolyakQuadrature::designLevelForEvals| code@>=
-void SmolyakQuadrature::designLevelForEvals(int max_evals, int& lev, int& evals) const
-{
- int last_evals;
- evals = 1;
- lev = 1;
- do {
- lev++;
- last_evals = evals;
- evals = calcNumEvaluations(lev);
- } while (lev < uquad.numLevels() && evals <= max_evals);
- lev--;
- evals = last_evals;
-}
-
-
-@ End of {\tt smolyak.cpp} file
diff --git a/dynare++/integ/cc/smolyak.hh b/dynare++/integ/cc/smolyak.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ed930814e3f633dd1e381cf6fc917a728aac64c0
--- /dev/null
+++ b/dynare++/integ/cc/smolyak.hh
@@ -0,0 +1,127 @@
+// Copyright 2005, Ondra Kamenik
+
+// Smolyak quadrature.
+
+/* This file defines Smolyak (sparse grid) multidimensional quadrature
+ for non-nested underlying one dimensional quadrature. Let $Q^1_l$ denote
+ the one dimensional quadrature of $l$ level. Let $n_l$ denote a
+ number of points in the $l$ level. Than the Smolyak quadrature can be
+ defined as
+ $$Q^df=\sum_{l\leq\vert k\vert\leq l+d-1}(-1)^{l+d-\vert k\vert-1}\left(\matrix{d-1\cr
+ \vert k\vert-l}\right)(Q_{k_1}^1\otimes\ldots\otimes Q_{k_d}^1)f,$$
+ where $d$ is the dimension, $k$ is $d$-dimensional sequence of
+ integers, and $\vert k\vert$ denotes a sum of the sequence.
+
+ Here we define |smolpit| as Smolyak iterator and |SmolyakQuadrature|. */
+
+#ifndef SMOLYAK_H
+#define SMOLYAK_H
+
+#include "int_sequence.h"
+#include "tl_static.h"
+#include "vector_function.hh"
+#include "quadrature.hh"
+
+/* Here we define the Smolyak point iterator. The Smolyak formula can
+ be broken to a sum of product quadratures with various combinations of
+ levels. The iterator follows this pattern. It maintains an index to a
+ summand and then a point coordinates within the summand (product
+ quadrature). The array of summands to which the |isummand| points is
+ maintained by the |SmolyakQuadrature| class to which the object knows
+ the pointer |smolq|.
+
+ We provide a constructor which points to the beginning of the given
+ summand. This constructor is used in |SmolyakQuadrature::begin| method
+ which approximately divideds all the iterators to subsets of equal
+ size. */
+
+class SmolyakQuadrature;
+
+class smolpit
+{
+protected:
+ const SmolyakQuadrature *smolq;
+ unsigned int isummand;
+ IntSequence *jseq;
+ ParameterSignal *sig;
+ Vector *p;
+ double w;
+public:
+ smolpit();
+ smolpit(const SmolyakQuadrature &q, unsigned int isum);
+ smolpit(const smolpit &spit);
+ ~smolpit();
+ bool operator==(const smolpit &spit) const;
+ bool
+ operator!=(const smolpit &spit) const
+ {
+ return !operator==(spit);
+ }
+ const smolpit &operator=(const smolpit &spit);
+ smolpit &operator++();
+ const ParameterSignal &
+ signal() const
+ {
+ return *sig;
+ }
+ const Vector &
+ point() const
+ {
+ return *p;
+ }
+ double
+ weight() const
+ {
+ return w;
+ }
+ void print() const;
+protected:
+ void setPointAndWeight();
+};
+
+/* Here we define the class |SmolyakQuadrature|. It maintains an array
+ of summands of the Smolyak quadrature formula:
+ $$\sum_{l\leq\vert k\vert\leq l+d-1}(-1)^{l+d-\vert
+ k\vert-1}\left(\matrix{d-1\cr
+ \vert k\vert-l}\right)(Q_{k_1}^1\otimes\ldots\otimes Q_{k_d}^1)f$$
+ Each summand is fully specified by sequence $k$. The summands are here
+ represented (besides $k$) also by sequence of number of points in each
+ level selected by $k$, and also by a cummulative number of
+ evaluations. The latter two are added only for conveniency.
+
+ The summands in the code are given by |levels|, which is a vector of
+ $k$ sequences, further by |levpoints| which is a vector of sequences
+ of nuber of points in each level, and by |cumevals| which is the
+ cumulative number of points, this is $\sum_k\prod_{i=1}^dn_{k_i}$,
+ where the sum is done through all $k$ before the current.
+
+ The |levels| and |levpoints| vectors are used by |smolpit|. */
+
+class SmolyakQuadrature : public QuadratureImpl<smolpit>
+{
+ friend class smolpit;
+ int level;
+ const OneDQuadrature &uquad;
+ vector<IntSequence> levels;
+ vector<IntSequence> levpoints;
+ vector<int> cumevals;
+ PascalTriangle psc;
+public:
+ SmolyakQuadrature(int d, int l, const OneDQuadrature &uq);
+ virtual ~SmolyakQuadrature()
+ {
+ }
+ virtual int numEvals(int level) const;
+ void designLevelForEvals(int max_eval, int &lev, int &evals) const;
+protected:
+ smolpit begin(int ti, int tn, int level) const;
+ unsigned int
+ numSummands() const
+ {
+ return levels.size();
+ }
+private:
+ int calcNumEvaluations(int level) const;
+};
+
+#endif
diff --git a/dynare++/integ/cc/smolyak.hweb b/dynare++/integ/cc/smolyak.hweb
deleted file mode 100644
index 4ce9c64fa01409e5ee46746534b225065e501d4f..0000000000000000000000000000000000000000
--- a/dynare++/integ/cc/smolyak.hweb
+++ /dev/null
@@ -1,123 +0,0 @@
-@q $Id: smolyak.hweb 431 2005-08-16 15:41:01Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@*2 Smolyak quadrature. This is {\tt smolyak.h} file
-
-This file defines Smolyak (sparse grid) multidimensional quadrature
-for non-nested underlying one dimensional quadrature. Let $Q^1_l$ denote
-the one dimensional quadrature of $l$ level. Let $n_l$ denote a
-number of points in the $l$ level. Than the Smolyak quadrature can be
-defined as
-$$Q^df=\sum_{l\leq\vert k\vert\leq l+d-1}(-1)^{l+d-\vert k\vert-1}\left(\matrix{d-1\cr
-\vert k\vert-l}\right)(Q_{k_1}^1\otimes\ldots\otimes Q_{k_d}^1)f,$$
-where $d$ is the dimension, $k$ is $d$-dimensional sequence of
-integers, and $\vert k\vert$ denotes a sum of the sequence.
-
-Here we define |smolpit| as Smolyak iterator and |SmolyakQuadrature|.
-
-@s smolpit int
-@s SmolyakQuadrature int
-@s PascalTriangle int
-@s SymmetrySet int
-@s symiterator int
-
-@c
-#ifndef SMOLYAK_H
-#define SMOLYAK_H
-
-#include "int_sequence.h"
-#include "tl_static.h"
-#include "vector_function.h"
-#include "quadrature.h"
-
-@<|smolpit| class declaration@>;
-@<|SmolyakQuadrature| class declaration@>;
-
-#endif
-
-@ Here we define the Smolyak point iterator. The Smolyak formula can
-be broken to a sum of product quadratures with various combinations of
-levels. The iterator follows this pattern. It maintains an index to a
-summand and then a point coordinates within the summand (product
-quadrature). The array of summands to which the |isummand| points is
-maintained by the |SmolyakQuadrature| class to which the object knows
-the pointer |smolq|.
-
-We provide a constructor which points to the beginning of the given
-summand. This constructor is used in |SmolyakQuadrature::begin| method
-which approximately divideds all the iterators to subsets of equal
-size.
-
-@<|smolpit| class declaration@>=
-class SmolyakQuadrature;
-
-class smolpit {
-protected:@;
- const SmolyakQuadrature* smolq;
- unsigned int isummand;
- IntSequence* jseq;
- ParameterSignal* sig;
- Vector* p;
- double w;
-public:@;
- smolpit();
- smolpit(const SmolyakQuadrature& q, unsigned int isum);
- smolpit(const smolpit& spit);
- ~smolpit();
- bool operator==(const smolpit& spit) const;
- bool operator!=(const smolpit& spit) const
- {@+ return ! operator==(spit);@+}
- const smolpit& operator=(const smolpit& spit);
- smolpit& operator++();
- const ParameterSignal& signal() const
- {@+ return *sig;@+}
- const Vector& point() const
- {@+ return *p;@+}
- double weight() const
- {@+ return w;@+}
- void print() const;
-protected:@;
- void setPointAndWeight();
-};
-
-@ Here we define the class |SmolyakQuadrature|. It maintains an array
-of summands of the Smolyak quadrature formula:
-$$\sum_{l\leq\vert k\vert\leq l+d-1}(-1)^{l+d-\vert
-k\vert-1}\left(\matrix{d-1\cr
-\vert k\vert-l}\right)(Q_{k_1}^1\otimes\ldots\otimes Q_{k_d}^1)f$$
-Each summand is fully specified by sequence $k$. The summands are here
-represented (besides $k$) also by sequence of number of points in each
-level selected by $k$, and also by a cummulative number of
-evaluations. The latter two are added only for conveniency.
-
-The summands in the code are given by |levels|, which is a vector of
-$k$ sequences, further by |levpoints| which is a vector of sequences
-of nuber of points in each level, and by |cumevals| which is the
-cumulative number of points, this is $\sum_k\prod_{i=1}^dn_{k_i}$,
-where the sum is done through all $k$ before the current.
-
-The |levels| and |levpoints| vectors are used by |smolpit|.
-
-@<|SmolyakQuadrature| class declaration@>=
-class SmolyakQuadrature : public QuadratureImpl<smolpit> {
- friend class smolpit;
- int level;
- const OneDQuadrature& uquad;
- vector<IntSequence> levels;
- vector<IntSequence> levpoints;
- vector<int> cumevals;
- PascalTriangle psc;
-public:@;
- SmolyakQuadrature(int d, int l, const OneDQuadrature& uq);
- virtual ~SmolyakQuadrature()@+ {}
- virtual int numEvals(int level) const;
- void designLevelForEvals(int max_eval, int& lev, int& evals) const;
-protected:@;
- smolpit begin(int ti, int tn, int level) const;
- unsigned int numSummands() const
- {@+ return levels.size();@+}
-private:@;
- int calcNumEvaluations(int level) const;
-};
-
-@ End of {\tt smolyak.h} file
diff --git a/dynare++/integ/cc/vector_function.cc b/dynare++/integ/cc/vector_function.cc
new file mode 100644
index 0000000000000000000000000000000000000000..aced635285a7a30bdb589e28ab374a1c0ddced43
--- /dev/null
+++ b/dynare++/integ/cc/vector_function.cc
@@ -0,0 +1,157 @@
+// Copyright 2005, Ondra Kamenik
+
+#include "vector_function.hh"
+
+#include <dynlapack.h>
+
+#include <cmath>
+
+#include <cstring>
+#include <algorithm>
+
+#ifdef __MINGW32__
+# define __CROSS_COMPILATION__
+#endif
+
+#ifdef __MINGW64__
+# define __CROSS_COMPILATION__
+#endif
+
+#ifdef __CROSS_COMPILATION__
+# define M_PI 3.14159265358979323846
+#endif
+
+/* Just an easy constructor of sequence of booleans defaulting to
+ change everywhere. */
+
+ParameterSignal::ParameterSignal(int n)
+ : data(new bool[n]), num(n)
+{
+ for (int i = 0; i < num; i++)
+ data[i] = true;
+}
+
+ParameterSignal::ParameterSignal(const ParameterSignal &sig)
+ : data(new bool[sig.num]), num(sig.num)
+{
+ memcpy(data, sig.data, num);
+}
+
+/* This sets |false| (no change) before a given parameter, and |true|
+ (change) after the given parameter (including). */
+
+void
+ParameterSignal::signalAfter(int l)
+{
+ for (int i = 0; i < std::min(l, num); i++)
+ data[i] = false;
+ for (int i = l; i < num; i++)
+ data[i] = true;
+}
+
+/* This constructs a function set hardcopying also the first. */
+VectorFunctionSet::VectorFunctionSet(const VectorFunction &f, int n)
+ : funcs(n), first_shallow(false)
+{
+ for (int i = 0; i < n; i++)
+ funcs[i] = f.clone();
+}
+
+/* This constructs a function set with shallow copy in the first and
+ hard copies in others. */
+VectorFunctionSet::VectorFunctionSet(VectorFunction &f, int n)
+ : funcs(n), first_shallow(true)
+{
+ if (n > 0)
+ funcs[0] = &f;
+ for (int i = 1; i < n; i++)
+ funcs[i] = f.clone();
+}
+
+/* This deletes the functions. The first is deleted only if it was not
+ a shallow copy. */
+
+VectorFunctionSet::~VectorFunctionSet()
+{
+ unsigned int start = first_shallow ? 1 : 0;
+ for (unsigned int i = start; i < funcs.size(); i++)
+ delete funcs[i];
+}
+
+/* Here we construct the object from the given function $f$ and given
+ variance-covariance matrix $\Sigma=$|vcov|. The matrix $A$ is
+ calculated as lower triangular and yields $\Sigma=AA^T$. */
+
+GaussConverterFunction::GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov)
+ : VectorFunction(f), func(&f), delete_flag(false), A(vcov.numRows(), vcov.numRows()),
+ multiplier(calcMultiplier())
+{
+ // todo: raise if |A.numRows() != indim()|
+ calcCholeskyFactor(vcov);
+}
+
+/* Here we construct the object in the same way, however we mark the
+ function as to be deleted. */
+
+GaussConverterFunction::GaussConverterFunction(VectorFunction *f, const GeneralMatrix &vcov)
+ : VectorFunction(*f), func(f), delete_flag(true), A(vcov.numRows(), vcov.numRows()),
+ multiplier(calcMultiplier())
+{
+ // todo: raise if |A.numRows() != indim()|
+ calcCholeskyFactor(vcov);
+}
+
+GaussConverterFunction::GaussConverterFunction(const GaussConverterFunction &f)
+ : VectorFunction(f), func(f.func->clone()), delete_flag(true), A(f.A),
+ multiplier(f.multiplier)
+{
+}
+
+/* Here we evaluate the function
+ $g(y)={1\over\sqrt{\pi^n}}f\left(\sqrt{2}Ay\right)$. Since the matrix $A$ is lower
+ triangular, the change signal for the function $f$ will look like
+ $(0,\ldots,0,1,\ldots,1)$ where the first $1$ is in the same position
+ as the first change in the given signal |sig| of the input
+ $y=$|point|. */
+
+void
+GaussConverterFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
+{
+ ParameterSignal s(sig);
+ int i = 0;
+ while (i < indim() && !sig[i])
+ i++;
+ s.signalAfter(i);
+
+ Vector x(indim());
+ x.zeros();
+ A.multaVec(x, point);
+ x.mult(sqrt(2.0));
+
+ func->eval(x, s, out);
+
+ out.mult(multiplier);
+}
+
+/* This returns $1\over\sqrt{\pi^n}$. */
+
+double
+GaussConverterFunction::calcMultiplier() const
+{
+ return sqrt(pow(M_PI, -1*indim()));
+}
+
+void
+GaussConverterFunction::calcCholeskyFactor(const GeneralMatrix &vcov)
+{
+ A = vcov;
+
+ lapack_int rows = A.numRows();
+ for (int i = 0; i < rows; i++)
+ for (int j = i+1; j < rows; j++)
+ A.get(i, j) = 0.0;
+
+ lapack_int info;
+ dpotrf("L", &rows, A.base(), &rows, &info);
+ // todo: raise if |info!=1|
+}
diff --git a/dynare++/integ/cc/vector_function.cweb b/dynare++/integ/cc/vector_function.cweb
deleted file mode 100644
index e2197f54145d518d77be16703df4f53e1fd8036c..0000000000000000000000000000000000000000
--- a/dynare++/integ/cc/vector_function.cweb
+++ /dev/null
@@ -1,190 +0,0 @@
-@q $Id: vector_function.cweb 431 2005-08-16 15:41:01Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@ This is {\tt vector\_function.cpp} file
-
-@c
-
-#include "vector_function.h"
-
-#include <dynlapack.h>
-
-#include <cmath>
-
-#include <cstring>
-#include <algorithm>
-
-#ifdef __MINGW32__
-#define __CROSS_COMPILATION__
-#endif
-
-#ifdef __MINGW64__
-#define __CROSS_COMPILATION__
-#endif
-
-#ifdef __CROSS_COMPILATION__
-#define M_PI 3.14159265358979323846
-#endif
-
-@<|ParameterSignal| constructor code@>;
-@<|ParameterSignal| copy constructor code@>;
-@<|ParameterSignal::signalAfter| code@>;
-@<|VectorFunctionSet| constructor 1 code@>;
-@<|VectorFunctionSet| constructor 2 code@>;
-@<|VectorFunctionSet| destructor code@>;
-@<|GaussConverterFunction| constructor code 1@>;
-@<|GaussConverterFunction| constructor code 2@>;
-@<|GaussConverterFunction| copy constructor code@>;
-@<|GaussConverterFunction::eval| code@>;
-@<|GaussConverterFunction::multiplier| code@>;
-@<|GaussConverterFunction::calcCholeskyFactor| code@>;
-
-@ Just an easy constructor of sequence of booleans defaulting to
-change everywhere.
-
-@<|ParameterSignal| constructor code@>=
-ParameterSignal::ParameterSignal(int n)
- : data(new bool[n]), num(n)
-{
- for (int i = 0; i < num; i++)
- data[i] = true;
-}
-
-@
-@<|ParameterSignal| copy constructor code@>=
-ParameterSignal::ParameterSignal(const ParameterSignal& sig)
- : data(new bool[sig.num]), num(sig.num)
-{
- memcpy(data, sig.data, num);
-}
-
-@ This sets |false| (no change) before a given parameter, and |true|
-(change) after the given parameter (including).
-
-@<|ParameterSignal::signalAfter| code@>=
-void ParameterSignal::signalAfter(int l)
-{
- for (int i = 0; i < std::min(l,num); i++)
- data[i] = false;
- for (int i = l; i < num; i++)
- data[i] = true;
-}
-
-@ This constructs a function set hardcopying also the first.
-@<|VectorFunctionSet| constructor 1 code@>=
-VectorFunctionSet::VectorFunctionSet(const VectorFunction& f, int n)
- : funcs(n), first_shallow(false)
-{
- for (int i = 0; i < n; i++)
- funcs[i] = f.clone();
-}
-
-@ This constructs a function set with shallow copy in the first and
-hard copies in others.
-
-@<|VectorFunctionSet| constructor 2 code@>=
-VectorFunctionSet::VectorFunctionSet(VectorFunction& f, int n)
- : funcs(n), first_shallow(true)
-{
- if (n > 0)
- funcs[0] = &f;
- for (int i = 1; i < n; i++)
- funcs[i] = f.clone();
-}
-
-@ This deletes the functions. The first is deleted only if it was not
-a shallow copy.
-
-@<|VectorFunctionSet| destructor code@>=
-VectorFunctionSet::~VectorFunctionSet()
-{
- unsigned int start = first_shallow ? 1 : 0;
- for (unsigned int i = start; i < funcs.size(); i++)
- delete funcs[i];
-}
-
-@ Here we construct the object from the given function $f$ and given
-variance-covariance matrix $\Sigma=$|vcov|. The matrix $A$ is
-calculated as lower triangular and yields $\Sigma=AA^T$.
-
-@<|GaussConverterFunction| constructor code 1@>=
-GaussConverterFunction::GaussConverterFunction(VectorFunction& f, const GeneralMatrix& vcov)
- : VectorFunction(f), func(&f), delete_flag(false), A(vcov.numRows(), vcov.numRows()),
- multiplier(calcMultiplier())
-{
- // todo: raise if |A.numRows() != indim()|
- calcCholeskyFactor(vcov);
-}
-
-@ Here we construct the object in the same way, however we mark the
-function as to be deleted.
-
-@<|GaussConverterFunction| constructor code 2@>=
-GaussConverterFunction::GaussConverterFunction(VectorFunction* f, const GeneralMatrix& vcov)
- : VectorFunction(*f), func(f), delete_flag(true), A(vcov.numRows(), vcov.numRows()),
- multiplier(calcMultiplier())
-{
- // todo: raise if |A.numRows() != indim()|
- calcCholeskyFactor(vcov);
-}
-
-
-@
-@<|GaussConverterFunction| copy constructor code@>=
-GaussConverterFunction::GaussConverterFunction(const GaussConverterFunction& f)
- : VectorFunction(f), func(f.func->clone()), delete_flag(true), A(f.A),
- multiplier(f.multiplier)
-{
-}
-
-@ Here we evaluate the function
-$g(y)={1\over\sqrt{\pi^n}}f\left(\sqrt{2}Ay\right)$. Since the matrix $A$ is lower
-triangular, the change signal for the function $f$ will look like
-$(0,\ldots,0,1,\ldots,1)$ where the first $1$ is in the same position
-as the first change in the given signal |sig| of the input
-$y=$|point|.
-
-@<|GaussConverterFunction::eval| code@>=
-void GaussConverterFunction::eval(const Vector& point, const ParameterSignal& sig, Vector& out)
-{
- ParameterSignal s(sig);
- int i = 0;
- while (i < indim() && !sig[i])
- i++;
- s.signalAfter(i);
-
- Vector x(indim());
- x.zeros();
- A.multaVec(x, point);
- x.mult(sqrt(2.0));
-
- func->eval(x, s, out);
-
- out.mult(multiplier);
-}
-
-@ This returns $1\over\sqrt{\pi^n}$.
-@<|GaussConverterFunction::multiplier| code@>=
-double GaussConverterFunction::calcMultiplier() const
-{
- return sqrt(pow(M_PI, -1*indim()));
-}
-
-@
-@<|GaussConverterFunction::calcCholeskyFactor| code@>=
-void GaussConverterFunction::calcCholeskyFactor(const GeneralMatrix& vcov)
-{
- A = vcov;
-
- lapack_int rows = A.numRows();
- for (int i = 0; i < rows; i++)
- for (int j = i+1; j < rows; j++)
- A.get(i,j) = 0.0;
-
- lapack_int info;
- dpotrf("L", &rows, A.base(), &rows, &info);
- // todo: raise if |info!=1|
-}
-
-
-@ End of {\tt vector\_function.cpp} file
diff --git a/dynare++/integ/cc/vector_function.hh b/dynare++/integ/cc/vector_function.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7ec856b82f8443eff1bc7ab448c92e2b25b19ea2
--- /dev/null
+++ b/dynare++/integ/cc/vector_function.hh
@@ -0,0 +1,174 @@
+// Copyright 2005, Ondra Kamenik
+
+// Vector function.
+
+/* This file defines interface for functions taking a vector as an input
+ and returning a vector (with a different size) as an output. We are
+ also introducing a parameter signalling; it is a boolean vector which
+ tracks parameters which were changed from the previous call. The
+ |VectorFunction| implementation can exploit this information and
+ evaluate the function more efficiently. The information can be
+ completely ignored.
+
+ From the signalling reason, and from other reasons, the function
+ evaluation is not |const|. */
+
+#ifndef VECTOR_FUNCTION_H
+#define VECTOR_FUNCTION_H
+
+#include "Vector.h"
+#include "GeneralMatrix.h"
+
+#include <vector>
+
+/* This is a simple class representing a vector of booleans. The items
+ night be retrieved or changed, or can be set |true| after some
+ point. This is useful when we multiply the vector with lower
+ triangular matrix.
+
+ |true| means that a parameter was changed. */
+
+class ParameterSignal
+{
+protected:
+ bool *data;
+ int num;
+public:
+ ParameterSignal(int n);
+ ParameterSignal(const ParameterSignal &sig);
+ ~ParameterSignal()
+ {
+ delete [] data;
+ }
+ void signalAfter(int l);
+ const bool &
+ operator[](int i) const
+ {
+ return data[i];
+ }
+ bool &
+ operator[](int i)
+ {
+ return data[i];
+ }
+};
+
+/* This is the abstract class for vector function. At this level of
+ abstraction we only need to know size of input vector and a size of
+ output vector.
+
+ The important thing here is a clone method, we will need to make hard
+ copies of vector functions since the evaluations are not |const|. The
+ hardcopies apply for parallelization. */
+
+class VectorFunction
+{
+protected:
+ int in_dim;
+ int out_dim;
+public:
+ VectorFunction(int idim, int odim)
+ : in_dim(idim), out_dim(odim)
+ {
+ }
+ VectorFunction(const VectorFunction &func)
+ : in_dim(func.in_dim), out_dim(func.out_dim)
+ {
+ }
+ virtual ~VectorFunction()
+ {
+ }
+ virtual VectorFunction *clone() const = 0;
+ virtual void eval(const Vector &point, const ParameterSignal &sig, Vector &out) = 0;
+ int
+ indim() const
+ {
+ return in_dim;
+ }
+ int
+ outdim() const
+ {
+ return out_dim;
+ }
+};
+
+/* This makes |n| copies of |VectorFunction|. The first constructor
+ make exactly |n| new copies, the second constructor copies only the
+ pointer to the first and others are hard (real) copies.
+
+ The class is useful for making a given number of copies at once, and
+ this set can be reused many times if we need mupliple copis of the
+ function (for example for paralelizing the code). */
+
+class VectorFunctionSet
+{
+protected:
+ std::vector<VectorFunction *> funcs;
+ bool first_shallow;
+public:
+ VectorFunctionSet(const VectorFunction &f, int n);
+ VectorFunctionSet(VectorFunction &f, int n);
+ ~VectorFunctionSet();
+ VectorFunction &
+ getFunc(int i)
+ {
+ return *(funcs[i]);
+ }
+ int
+ getNum() const
+ {
+ return funcs.size();
+ }
+};
+
+/* This class wraps another |VectorFunction| to allow integration of a
+ function through normally distributed inputs. Namely, if one wants to
+ integrate
+ $${1\over\sqrt{(2\pi)^n\vert\Sigma\vert}}\int f(x)e^{-{1\over2}x^T\Sigma^{-1}x}{\rm d}x$$
+ then if we write $\Sigma=AA^T$ and $x=\sqrt{2}Ay$, we get integral
+ $${1\over\sqrt{(2\pi)^n\vert\Sigma\vert}}
+ \int f\left(\sqrt{2}Ay\right)e^{-y^Ty}\sqrt{2^n}\vert A\vert{\rm d}y=
+ {1\over\sqrt{\pi^n}}\int f\left(\sqrt{2}Ay\right)e^{-y^Ty}{\rm d}y,$$
+ which means that a given function $f$ we have to wrap to yield a function
+ $$g(y)={1\over\sqrt{\pi^n}}f\left(\sqrt{2}Ay\right).$$
+ This is exactly what this class is doing. This transformation is
+ useful since the Gauss--Hermite points and weights are defined for
+ weighting function $e^{-y^2}$, so this transformation allows using
+ Gauss--Hermite quadratures seemlessly in a context of integration through
+ normally distributed inputs.
+
+ The class maintains a pointer to the function $f$. When the object is
+ constructed by the first constructor, the $f$ is not copied. If the
+ object of this class is copied, then $f$ is copied and we need to
+ remember to destroy it in the desctructor; hence |delete_flag|. The
+ second constructor takes a pointer to the function and differs from
+ the first only by setting |delete_flag| to |true|. */
+
+class GaussConverterFunction : public VectorFunction
+{
+protected:
+ VectorFunction *func;
+ bool delete_flag;
+ GeneralMatrix A;
+ double multiplier;
+public:
+ GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov);
+ GaussConverterFunction(VectorFunction *f, const GeneralMatrix &vcov);
+ GaussConverterFunction(const GaussConverterFunction &f);
+ virtual ~GaussConverterFunction()
+ {
+ if (delete_flag)
+ delete func;
+ }
+ virtual VectorFunction *
+ clone() const
+ {
+ return new GaussConverterFunction(*this);
+ }
+ virtual void eval(const Vector &point, const ParameterSignal &sig, Vector &out);
+private:
+ double calcMultiplier() const;
+ void calcCholeskyFactor(const GeneralMatrix &vcov);
+};
+
+#endif
diff --git a/dynare++/integ/cc/vector_function.hweb b/dynare++/integ/cc/vector_function.hweb
deleted file mode 100644
index 840b5892b093fa95983b0406d06c104005fe0e29..0000000000000000000000000000000000000000
--- a/dynare++/integ/cc/vector_function.hweb
+++ /dev/null
@@ -1,156 +0,0 @@
-@q $Id: vector_function.hweb 431 2005-08-16 15:41:01Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@*2 Vector function. This is {\tt vector\_function.h} file
-
-This file defines interface for functions taking a vector as an input
-and returning a vector (with a different size) as an output. We are
-also introducing a parameter signalling; it is a boolean vector which
-tracks parameters which were changed from the previous call. The
-|VectorFunction| implementation can exploit this information and
-evaluate the function more efficiently. The information can be
-completely ignored.
-
-From the signalling reason, and from other reasons, the function
-evaluation is not |const|.
-
-@s ParameterSignal int
-@s VectorFunction int
-@s VectorFunctionSet int
-@s GaussConverterFunction int
-
-@c
-#ifndef VECTOR_FUNCTION_H
-#define VECTOR_FUNCTION_H
-
-#include "Vector.h"
-#include "GeneralMatrix.h"
-
-#include <vector>
-
-@<|ParameterSignal| class declaration@>;
-@<|VectorFunction| class declaration@>;
-@<|VectorFunctionSet| class declaration@>;
-@<|GaussConverterFunction| class declaration@>;
-
-#endif
-
-@ This is a simple class representing a vector of booleans. The items
-night be retrieved or changed, or can be set |true| after some
-point. This is useful when we multiply the vector with lower
-triangular matrix.
-
-|true| means that a parameter was changed.
-
-@<|ParameterSignal| class declaration@>=
-class ParameterSignal {
-protected:@;
- bool* data;
- int num;
-public:@;
- ParameterSignal(int n);
- ParameterSignal(const ParameterSignal& sig);
- ~ParameterSignal()
- {@+ delete [] data;@+}
- void signalAfter(int l);
- const bool& operator[](int i) const
- {@+ return data[i];@+}
- bool& operator[](int i)
- {@+ return data[i];@+}
-};
-
-@ This is the abstract class for vector function. At this level of
-abstraction we only need to know size of input vector and a size of
-output vector.
-
-The important thing here is a clone method, we will need to make hard
-copies of vector functions since the evaluations are not |const|. The
-hardcopies apply for parallelization.
-
-@<|VectorFunction| class declaration@>=
-class VectorFunction {
-protected:@;
- int in_dim;
- int out_dim;
-public:@;
- VectorFunction(int idim, int odim)
- : in_dim(idim), out_dim(odim)@+ {}
- VectorFunction(const VectorFunction& func)
- : in_dim(func.in_dim), out_dim(func.out_dim)@+ {}
- virtual ~VectorFunction()@+ {}
- virtual VectorFunction* clone() const =0;
- virtual void eval(const Vector& point, const ParameterSignal& sig, Vector& out) =0;
- int indim() const
- {@+ return in_dim;@+}
- int outdim() const
- {@+ return out_dim;@+}
-};
-
-@ This makes |n| copies of |VectorFunction|. The first constructor
-make exactly |n| new copies, the second constructor copies only the
-pointer to the first and others are hard (real) copies.
-
-The class is useful for making a given number of copies at once, and
-this set can be reused many times if we need mupliple copis of the
-function (for example for paralelizing the code).
-
-@<|VectorFunctionSet| class declaration@>=
-class VectorFunctionSet {
-protected:@;
- std::vector<VectorFunction*> funcs;
- bool first_shallow;
-public:@;
- VectorFunctionSet(const VectorFunction& f, int n);
- VectorFunctionSet(VectorFunction& f, int n);
- ~VectorFunctionSet();
- VectorFunction& getFunc(int i)
- {@+ return *(funcs[i]);@+}
- int getNum() const
- {@+ return funcs.size(); @+}
-};
-
-@ This class wraps another |VectorFunction| to allow integration of a
-function through normally distributed inputs. Namely, if one wants to
-integrate
-$${1\over\sqrt{(2\pi)^n\vert\Sigma\vert}}\int f(x)e^{-{1\over2}x^T\Sigma^{-1}x}{\rm d}x$$
-then if we write $\Sigma=AA^T$ and $x=\sqrt{2}Ay$, we get integral
-$${1\over\sqrt{(2\pi)^n\vert\Sigma\vert}}
-\int f\left(\sqrt{2}Ay\right)e^{-y^Ty}\sqrt{2^n}\vert A\vert{\rm d}y=
-{1\over\sqrt{\pi^n}}\int f\left(\sqrt{2}Ay\right)e^{-y^Ty}{\rm d}y,$$
-which means that a given function $f$ we have to wrap to yield a function
-$$g(y)={1\over\sqrt{\pi^n}}f\left(\sqrt{2}Ay\right).$$
-This is exactly what this class is doing. This transformation is
-useful since the Gauss--Hermite points and weights are defined for
-weighting function $e^{-y^2}$, so this transformation allows using
-Gauss--Hermite quadratures seemlessly in a context of integration through
-normally distributed inputs.
-
-The class maintains a pointer to the function $f$. When the object is
-constructed by the first constructor, the $f$ is not copied. If the
-object of this class is copied, then $f$ is copied and we need to
-remember to destroy it in the desctructor; hence |delete_flag|. The
-second constructor takes a pointer to the function and differs from
-the first only by setting |delete_flag| to |true|.
-
-@<|GaussConverterFunction| class declaration@>=
-class GaussConverterFunction : public VectorFunction {
-protected:@;
- VectorFunction* func;
- bool delete_flag;
- GeneralMatrix A;
- double multiplier;
-public:@;
- GaussConverterFunction(VectorFunction& f, const GeneralMatrix& vcov);
- GaussConverterFunction(VectorFunction* f, const GeneralMatrix& vcov);
- GaussConverterFunction(const GaussConverterFunction& f);
- virtual ~GaussConverterFunction()
- {@+ if (delete_flag) delete func; @+}
- virtual VectorFunction* clone() const
- {@+ return new GaussConverterFunction(*this);@+}
- virtual void eval(const Vector& point, const ParameterSignal& sig, Vector& out);
-private:@;
- double calcMultiplier() const;
- void calcCholeskyFactor(const GeneralMatrix& vcov);
-};
-
-@ End of {\tt vector\_function.h} file
diff --git a/dynare++/integ/src/Makefile.am b/dynare++/integ/src/Makefile.am
index d964d2eff4436079e6d6c49eb16306206d71a3ce..8c55df19e3005630f28399b676356993d437bb97 100644
--- a/dynare++/integ/src/Makefile.am
+++ b/dynare++/integ/src/Makefile.am
@@ -1,6 +1,6 @@
noinst_PROGRAMS = quadrature-points
-quadrature_points_SOURCES = quadrature-points.cpp
+quadrature_points_SOURCES = quadrature-points.cc
quadrature_points_CPPFLAGS = -I../.. -I../../sylv/cc -I../../integ/cc -I../../tl/cc
quadrature_points_CXXFLAGS = $(AM_CXXFLAGS) $(PTHREAD_CFLAGS)
quadrature_points_LDADD = ../cc/libinteg.a ../../tl/cc/libtl.a ../../parser/cc/libparser.a ../../sylv/cc/libsylv.a ../../utils/cc/libutils.a $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(PTHREAD_LIBS)
diff --git a/dynare++/integ/src/quadrature-points.cc b/dynare++/integ/src/quadrature-points.cc
new file mode 100644
index 0000000000000000000000000000000000000000..7d8358aa604726fdbdd33740653b86b4bf2b20e7
--- /dev/null
+++ b/dynare++/integ/src/quadrature-points.cc
@@ -0,0 +1,212 @@
+// Copyright (C) 2008-2011, Ondra Kamenik
+
+#include "parser/cc/matrix_parser.h"
+#include "utils/cc/memory_file.h"
+#include "utils/cc/exception.h"
+#include "sylv/cc/GeneralMatrix.h"
+#include "sylv/cc/Vector.h"
+#include "sylv/cc/SymSchurDecomp.h"
+#include "sylv/cc/SylvException.h"
+#include "integ/cc/quadrature.hh"
+#include "integ/cc/smolyak.hh"
+#include "integ/cc/product.hh"
+
+#include <getopt.h>
+#include <cstdio>
+
+#include <cmath>
+
+struct QuadParams
+{
+ const char *outname;
+ const char *vcovname;
+ int max_level;
+ double discard_weight;
+ QuadParams(int argc, char **argv);
+ void check_consistency() const;
+private:
+ enum {opt_max_level, opt_discard_weight, opt_vcov};
+};
+
+QuadParams::QuadParams(int argc, char **argv)
+ : outname(NULL), vcovname(NULL), max_level(3), discard_weight(0.0)
+{
+ if (argc == 1)
+ {
+ // print the help and exit
+ exit(1);
+ }
+
+ outname = argv[argc-1];
+ argc--;
+
+ struct option const opts [] = {
+ {"max-level", required_argument, NULL, opt_max_level},
+ {"discard-weight", required_argument, NULL, opt_discard_weight},
+ {"vcov", required_argument, NULL, opt_vcov},
+ {NULL, 0, NULL, 0}
+ };
+
+ int ret;
+ int index;
+ while (-1 != (ret = getopt_long(argc, argv, "", opts, &index)))
+ {
+ switch (ret)
+ {
+ case opt_max_level:
+ if (1 != sscanf(optarg, "%d", &max_level))
+ fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
+ break;
+ case opt_discard_weight:
+ if (1 != sscanf(optarg, "%lf", &discard_weight))
+ fprintf(stderr, "Couldn't parse float %s, ignored\n", optarg);
+ break;
+ case opt_vcov:
+ vcovname = optarg;
+ break;
+ }
+ }
+
+ check_consistency();
+}
+
+void
+QuadParams::check_consistency() const
+{
+ if (outname == NULL)
+ {
+ fprintf(stderr, "Error: output name not set\n");
+ exit(1);
+ }
+
+ if (vcovname == NULL)
+ {
+ fprintf(stderr, "Error: vcov file name not set\n");
+ exit(1);
+ }
+}
+
+/** Utility class for ordering pointers to vectors according their
+ * ordering. */
+struct OrderVec
+{
+ bool
+ operator()(const Vector *a, const Vector *b) const
+ {
+ return *a < *b;
+ }
+};
+
+int
+main(int argc, char **argv)
+{
+ QuadParams params(argc, argv);
+
+ // open output file for writing
+ FILE *fout;
+ if (NULL == (fout = fopen(params.outname, "w")))
+ {
+ fprintf(stderr, "Could not open %s for writing\n", params.outname);
+ exit(1);
+ }
+
+ try
+ {
+
+ // open memory file for vcov
+ ogu::MemoryFile vcov_mf(params.vcovname);
+
+ // parse the vcov matrix
+ ogp::MatrixParser mp;
+ mp.parse(vcov_mf.length(), vcov_mf.base());
+ if (mp.nrows() != mp.ncols())
+ throw ogu::Exception(__FILE__, __LINE__,
+ "VCOV matrix not square");
+ // and put to the GeneralMatrix
+ GeneralMatrix vcov(mp.nrows(), mp.ncols());
+ vcov.zeros();
+ for (ogp::MPIterator it = mp.begin(); it != mp.end(); ++it)
+ vcov.get(it.row(), it.col()) = *it;
+
+ // calculate the factor A of vcov, so that A*A^T=VCOV
+ GeneralMatrix A(vcov.numRows(), vcov.numRows());
+ SymSchurDecomp ssd(vcov);
+ ssd.getFactor(A);
+
+ // construct Gauss-Hermite quadrature
+ GaussHermite ghq;
+ // construct Smolyak quadrature
+ int level = params.max_level;
+ SmolyakQuadrature sq(vcov.numRows(), level, ghq);
+
+ printf("Dimension: %d\n", vcov.numRows());
+ printf("Maximum level: %d\n", level);
+ printf("Total number of nodes: %d\n", sq.numEvals(level));
+
+ // put the points to the vector
+ std::vector<Vector *> points;
+ for (smolpit qit = sq.start(level); qit != sq.end(level); ++qit)
+ points.push_back(new Vector((const Vector &) qit.point()));
+ // sort and uniq
+ OrderVec ordvec;
+ std::sort(points.begin(), points.end(), ordvec);
+ std::vector<Vector *>::iterator new_end = std::unique(points.begin(), points.end());
+ for (std::vector<Vector *>::iterator it = new_end; it != points.end(); ++it)
+ delete *it;
+ points.erase(new_end, points.end());
+
+ printf("Duplicit nodes removed: %lu\n", (unsigned long) (sq.numEvals(level)-points.size()));
+
+ // calculate weights and mass
+ double mass = 0.0;
+ std::vector<double> weights;
+ for (int i = 0; i < (int) points.size(); i++)
+ {
+ weights.push_back(std::exp(-points[i]->dot(*(points[i]))));
+ mass += weights.back();
+ }
+
+ // calculate discarded mass
+ double discard_mass = 0.0;
+ for (int i = 0; i < (int) weights.size(); i++)
+ if (weights[i]/mass < params.discard_weight)
+ discard_mass += weights[i];
+
+ printf("Total mass discarded: %f\n", discard_mass/mass);
+
+ // dump the results
+ int npoints = 0;
+ double upscale_weight = 1/(mass-discard_mass);
+ Vector x(vcov.numRows());
+ for (int i = 0; i < (int) weights.size(); i++)
+ if (weights[i]/mass >= params.discard_weight)
+ {
+ // print the upscaled weight
+ fprintf(fout, "%20.16g", upscale_weight*weights[i]);
+ // multiply point with the factor A and sqrt(2)
+ A.multVec(0.0, x, std::sqrt(2.), *(points[i]));
+ // print the coordinates
+ for (int j = 0; j < x.length(); j++)
+ fprintf(fout, " %20.16g", x[j]);
+ fprintf(fout, "\n");
+ npoints++;
+ }
+
+ printf("Final number of points: %d\n", npoints);
+
+ fclose(fout);
+
+ }
+ catch (const SylvException &e)
+ {
+ e.printMessage();
+ return 1;
+ }
+ catch (const ogu::Exception &e)
+ {
+ e.print();
+ return 1;
+ }
+
+ return 0;
+}
diff --git a/dynare++/integ/src/quadrature-points.cpp b/dynare++/integ/src/quadrature-points.cpp
deleted file mode 100644
index d1dd477f3cab452ddf1d7b8c8a600ebcdca0952d..0000000000000000000000000000000000000000
--- a/dynare++/integ/src/quadrature-points.cpp
+++ /dev/null
@@ -1,192 +0,0 @@
-// Copyright (C) 2008-2011, Ondra Kamenik
-
-#include "parser/cc/matrix_parser.h"
-#include "utils/cc/memory_file.h"
-#include "utils/cc/exception.h"
-#include "sylv/cc/GeneralMatrix.h"
-#include "sylv/cc/Vector.h"
-#include "sylv/cc/SymSchurDecomp.h"
-#include "sylv/cc/SylvException.h"
-#include "integ/cc/quadrature.h"
-#include "integ/cc/smolyak.h"
-#include "integ/cc/product.h"
-
-#include <getopt.h>
-#include <cstdio>
-
-#include <cmath>
-
-struct QuadParams {
- const char* outname;
- const char* vcovname;
- int max_level;
- double discard_weight;
- QuadParams(int argc, char** argv);
- void check_consistency() const;
-private:
- enum {opt_max_level, opt_discard_weight, opt_vcov};
-};
-
-QuadParams::QuadParams(int argc, char** argv)
- : outname(NULL), vcovname(NULL), max_level(3), discard_weight(0.0)
-{
- if (argc == 1) {
- // print the help and exit
- exit(1);
- }
-
- outname = argv[argc-1];
- argc--;
-
- struct option const opts [] = {
- {"max-level", required_argument, NULL, opt_max_level},
- {"discard-weight", required_argument, NULL, opt_discard_weight},
- {"vcov", required_argument, NULL, opt_vcov},
- {NULL, 0, NULL, 0}
- };
-
- int ret;
- int index;
- while (-1 != (ret = getopt_long(argc, argv, "", opts, &index))) {
- switch (ret) {
- case opt_max_level:
- if (1 != sscanf(optarg, "%d", &max_level))
- fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
- break;
- case opt_discard_weight:
- if (1 != sscanf(optarg, "%lf", &discard_weight))
- fprintf(stderr, "Couldn't parse float %s, ignored\n", optarg);
- break;
- case opt_vcov:
- vcovname = optarg;
- break;
- }
- }
-
- check_consistency();
-}
-
-void QuadParams::check_consistency() const
-{
- if (outname == NULL) {
- fprintf(stderr, "Error: output name not set\n");
- exit(1);
- }
-
- if (vcovname == NULL) {
- fprintf(stderr, "Error: vcov file name not set\n");
- exit(1);
- }
-}
-
-/** Utility class for ordering pointers to vectors according their
- * ordering. */
-struct OrderVec {
- bool operator()(const Vector* a, const Vector* b) const
- {return *a < *b;}
-};
-
-int main(int argc, char** argv)
-{
- QuadParams params(argc, argv);
-
- // open output file for writing
- FILE* fout;
- if (NULL == (fout=fopen(params.outname, "w"))) {
- fprintf(stderr, "Could not open %s for writing\n", params.outname);
- exit(1);
- }
-
- try {
-
- // open memory file for vcov
- ogu::MemoryFile vcov_mf(params.vcovname);
-
- // parse the vcov matrix
- ogp::MatrixParser mp;
- mp.parse(vcov_mf.length(), vcov_mf.base());
- if (mp.nrows() != mp.ncols())
- throw ogu::Exception(__FILE__,__LINE__,
- "VCOV matrix not square");
- // and put to the GeneralMatrix
- GeneralMatrix vcov(mp.nrows(), mp.ncols());
- vcov.zeros();
- for (ogp::MPIterator it = mp.begin(); it != mp.end(); ++it)
- vcov.get(it.row(), it.col()) = *it;
-
- // calculate the factor A of vcov, so that A*A^T=VCOV
- GeneralMatrix A(vcov.numRows(), vcov.numRows());
- SymSchurDecomp ssd(vcov);
- ssd.getFactor(A);
-
- // construct Gauss-Hermite quadrature
- GaussHermite ghq;
- // construct Smolyak quadrature
- int level = params.max_level;
- SmolyakQuadrature sq(vcov.numRows(), level, ghq);
-
- printf("Dimension: %d\n", vcov.numRows());
- printf("Maximum level: %d\n", level);
- printf("Total number of nodes: %d\n", sq.numEvals(level));
-
- // put the points to the vector
- std::vector<Vector*> points;
- for (smolpit qit = sq.start(level); qit != sq.end(level); ++qit)
- points.push_back(new Vector((const Vector&)qit.point()));
- // sort and uniq
- OrderVec ordvec;
- std::sort(points.begin(), points.end(), ordvec);
- std::vector<Vector*>::iterator new_end = std::unique(points.begin(), points.end());
- for (std::vector<Vector*>::iterator it = new_end; it != points.end(); ++it)
- delete *it;
- points.erase(new_end, points.end());
-
- printf("Duplicit nodes removed: %lu\n", (unsigned long) (sq.numEvals(level)-points.size()));
-
- // calculate weights and mass
- double mass = 0.0;
- std::vector<double> weights;
- for (int i = 0; i < (int)points.size(); i++) {
- weights.push_back(std::exp(-points[i]->dot(*(points[i]))));
- mass += weights.back();
- }
-
- // calculate discarded mass
- double discard_mass = 0.0;
- for (int i = 0; i < (int)weights.size(); i++)
- if (weights[i]/mass < params.discard_weight)
- discard_mass += weights[i];
-
- printf("Total mass discarded: %f\n", discard_mass/mass);
-
- // dump the results
- int npoints = 0;
- double upscale_weight = 1/(mass-discard_mass);
- Vector x(vcov.numRows());
- for (int i = 0; i < (int)weights.size(); i++)
- if (weights[i]/mass >= params.discard_weight) {
- // print the upscaled weight
- fprintf(fout, "%20.16g", upscale_weight*weights[i]);
- // multiply point with the factor A and sqrt(2)
- A.multVec(0.0, x, std::sqrt(2.), *(points[i]));
- // print the coordinates
- for (int j = 0; j < x.length(); j++)
- fprintf(fout, " %20.16g", x[j]);
- fprintf(fout, "\n");
- npoints++;
- }
-
- printf("Final number of points: %d\n", npoints);
-
- fclose(fout);
-
- } catch (const SylvException& e) {
- e.printMessage();
- return 1;
- } catch (const ogu::Exception& e) {
- e.print();
- return 1;
- }
-
- return 0;
-}
diff --git a/dynare++/integ/testing/Makefile.am b/dynare++/integ/testing/Makefile.am
index c30ef377a9f657e9325fe574cde8512bf3ceef42..3838bcb888b7fac82bf080bf515337b2a56f5bec 100644
--- a/dynare++/integ/testing/Makefile.am
+++ b/dynare++/integ/testing/Makefile.am
@@ -1,6 +1,6 @@
check_PROGRAMS = tests
-tests_SOURCES = tests.cpp
+tests_SOURCES = tests.cc
tests_CPPFLAGS = -I../cc -I../../tl/cc -I../../sylv/cc -I$(top_srcdir)/mex/sources
tests_CXXFLAGS = $(AM_CXXFLAGS) $(PTHREAD_CFLAGS)
tests_LDFLAGS = $(AM_LDFLAGS) $(LDFLAGS_MATIO)
diff --git a/dynare++/integ/testing/tests.cc b/dynare++/integ/testing/tests.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d10e06393cc9ed6802621b9209edee7f34458318
--- /dev/null
+++ b/dynare++/integ/testing/tests.cc
@@ -0,0 +1,637 @@
+/* $Id: tests.cpp 431 2005-08-16 15:41:01Z kamenik $ */
+/* Copyright 2005, Ondra Kamenik */
+
+#include "GeneralMatrix.h"
+#include <dynlapack.h>
+#include "SylvException.h"
+
+#include "rfs_tensor.h"
+#include "normal_moments.h"
+
+#include "vector_function.h"
+#include "quadrature.h"
+#include "smolyak.h"
+#include "product.h"
+#include "quasi_mcarlo.h"
+
+#include <cstdio>
+#include <cstring>
+#include <sys/time.h>
+#include <cmath>
+
+const int num_threads = 2; // does nothing if DEBUG defined
+
+// evaluates unfolded (Dx)^k power, where x is a vector, D is a
+// Cholesky factor (lower triangular)
+class MomentFunction : public VectorFunction
+{
+ GeneralMatrix D;
+ int k;
+public:
+ MomentFunction(const GeneralMatrix &inD, int kk)
+ : VectorFunction(inD.numRows(), UFSTensor::calcMaxOffset(inD.numRows(), kk)),
+ D(inD), k(kk)
+ {
+ }
+ MomentFunction(const MomentFunction &func)
+ : VectorFunction(func), D(func.D), k(func.k)
+ {
+ }
+ VectorFunction *
+ clone() const
+ {
+ return new MomentFunction(*this);
+ }
+ void eval(const Vector &point, const ParameterSignal &sig, Vector &out);
+};
+
+void
+MomentFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
+{
+ if (point.length() != indim() || out.length() != outdim())
+ {
+ printf("Wrong length of vectors in MomentFunction::eval\n");
+ exit(1);
+ }
+ Vector y(point);
+ y.zeros();
+ D.multaVec(y, point);
+ URSingleTensor ypow(y, k);
+ out.zeros();
+ out.add(1.0, ypow.getData());
+}
+
+class TensorPower : public VectorFunction
+{
+ int k;
+public:
+ TensorPower(int nvar, int kk)
+ : VectorFunction(nvar, UFSTensor::calcMaxOffset(nvar, kk)), k(kk)
+ {
+ }
+ TensorPower(const TensorPower &func)
+ : VectorFunction(func), k(func.k)
+ {
+ }
+ VectorFunction *
+ clone() const
+ {
+ return new TensorPower(*this);
+ }
+ void eval(const Vector &point, const ParameterSignal &sig, Vector &out);
+};
+
+void
+TensorPower::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
+{
+ if (point.length() != indim() || out.length() != outdim())
+ {
+ printf("Wrong length of vectors in TensorPower::eval\n");
+ exit(1);
+ }
+ URSingleTensor ypow(point, k);
+ out.zeros();
+ out.add(1.0, ypow.getData());
+}
+
+// evaluates (1+1/d)^d*(x_1*...*x_d)^(1/d), its integral over <0,1>^d
+// is 1.0, and its variation grows exponetially
+// d = dim
+class Function1 : public VectorFunction
+{
+ int dim;
+public:
+ Function1(int d)
+ : VectorFunction(d, 1), dim(d)
+ {
+ }
+ Function1(const Function1 &f)
+ : VectorFunction(f.indim(), f.outdim()), dim(f.dim)
+ {
+ }
+ VectorFunction *
+ clone() const
+ {
+ return new Function1(*this);
+ }
+ virtual void eval(const Vector &point, const ParameterSignal &sig, Vector &out);
+};
+
+void
+Function1::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
+{
+ if (point.length() != dim || out.length() != 1)
+ {
+ printf("Wrong length of vectors in Function1::eval\n");
+ exit(1);
+ }
+ double r = 1;
+ for (int i = 0; i < dim; i++)
+ r *= point[i];
+ r = pow(r, 1.0/dim);
+ r *= pow(1.0 + 1.0/dim, (double) dim);
+ out[0] = r;
+}
+
+// evaluates Function1 but with transformation x_i=0.5(y_i+1)
+// this makes the new function integrate over <-1,1>^d to 1.0
+class Function1Trans : public Function1
+{
+public:
+ Function1Trans(int d)
+ : Function1(d)
+ {
+ }
+ Function1Trans(const Function1Trans &func)
+ : Function1(func)
+ {
+ }
+ VectorFunction *
+ clone() const
+ {
+ return new Function1Trans(*this);
+ }
+ virtual void eval(const Vector &point, const ParameterSignal &sig, Vector &out);
+};
+
+void
+Function1Trans::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
+{
+ Vector p(point.length());
+ for (int i = 0; i < p.length(); i++)
+ p[i] = 0.5*(point[i]+1);
+ Function1::eval(p, sig, out);
+ out.mult(pow(0.5, indim()));
+}
+
+// WallTimer class. Constructor saves the wall time, destructor
+// cancels the current time from the saved, and prints the message
+// with time information
+class WallTimer
+{
+ char mes[100];
+ struct timeval start;
+ bool new_line;
+public:
+ WallTimer(const char *m, bool nl = true)
+ {
+ strcpy(mes, m); new_line = nl; gettimeofday(&start, NULL);
+ }
+ ~WallTimer()
+ {
+ struct timeval end;
+ gettimeofday(&end, NULL);
+ printf("%s%8.4g", mes,
+ end.tv_sec-start.tv_sec + (end.tv_usec-start.tv_usec)*1.0e-6);
+ if (new_line)
+ printf("\n");
+ }
+};
+
+/****************************************************/
+/* declaration of TestRunnable class */
+/****************************************************/
+class TestRunnable
+{
+ char name[100];
+public:
+ int dim; // dimension of the solved problem
+ int nvar; // number of variable of the solved problem
+ TestRunnable(const char *n, int d, int nv)
+ : dim(d), nvar(nv)
+ {
+ strncpy(name, n, 100);
+ }
+ bool test() const;
+ virtual bool run() const = 0;
+ const char *
+ getName() const
+ {
+ return name;
+ }
+protected:
+ static bool smolyak_normal_moments(const GeneralMatrix &m, int imom, int level);
+ static bool product_normal_moments(const GeneralMatrix &m, int imom, int level);
+ static bool qmc_normal_moments(const GeneralMatrix &m, int imom, int level);
+ static bool smolyak_product_cube(const VectorFunction &func, const Vector &res,
+ double tol, int level);
+ static bool qmc_cube(const VectorFunction &func, double res, double tol, int level);
+};
+
+bool
+TestRunnable::test() const
+{
+ printf("Running test <%s>\n", name);
+ bool passed;
+ {
+ WallTimer tim("Wall clock time ", false);
+ passed = run();
+ }
+ if (passed)
+ {
+ printf("............................ passed\n\n");
+ return passed;
+ }
+ else
+ {
+ printf("............................ FAILED\n\n");
+ return passed;
+ }
+}
+
+/****************************************************/
+/* definition of TestRunnable static methods */
+/****************************************************/
+bool
+TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level)
+{
+ // first make m*m' and then Cholesky factor
+ GeneralMatrix mtr(m, "transpose");
+ GeneralMatrix msq(m, mtr);
+
+ // make vector function
+ int dim = m.numRows();
+ TensorPower tp(dim, imom);
+ GaussConverterFunction func(tp, msq);
+
+ // smolyak quadrature
+ Vector smol_out(UFSTensor::calcMaxOffset(dim, imom));
+ {
+ WallTimer tim("\tSmolyak quadrature time: ");
+ GaussHermite gs;
+ SmolyakQuadrature quad(dim, level, gs);
+ quad.integrate(func, level, num_threads, smol_out);
+ printf("\tNumber of Smolyak evaluations: %d\n", quad.numEvals(level));
+ }
+
+ // check against theoretical moments
+ UNormalMoments moments(imom, msq);
+ smol_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
+ printf("\tError: %16.12g\n", smol_out.getMax());
+ return smol_out.getMax() < 1.e-7;
+}
+
+bool
+TestRunnable::product_normal_moments(const GeneralMatrix &m, int imom, int level)
+{
+ // first make m*m' and then Cholesky factor
+ GeneralMatrix mtr(m, "transpose");
+ GeneralMatrix msq(m, mtr);
+
+ // make vector function
+ int dim = m.numRows();
+ TensorPower tp(dim, imom);
+ GaussConverterFunction func(tp, msq);
+
+ // product quadrature
+ Vector prod_out(UFSTensor::calcMaxOffset(dim, imom));
+ {
+ WallTimer tim("\tProduct quadrature time: ");
+ GaussHermite gs;
+ ProductQuadrature quad(dim, gs);
+ quad.integrate(func, level, num_threads, prod_out);
+ printf("\tNumber of product evaluations: %d\n", quad.numEvals(level));
+ }
+
+ // check against theoretical moments
+ UNormalMoments moments(imom, msq);
+ prod_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
+ printf("\tError: %16.12g\n", prod_out.getMax());
+ return prod_out.getMax() < 1.e-7;
+}
+
+bool
+TestRunnable::qmc_normal_moments(const GeneralMatrix &m, int imom, int level)
+{
+ // first make m*m' and then Cholesky factor
+ GeneralMatrix mtr(m, "transpose");
+ GeneralMatrix msq(m, mtr);
+ GeneralMatrix mchol(msq);
+ int rows = mchol.numRows();
+ for (int i = 0; i < rows; i++)
+ for (int j = i+1; j < rows; j++)
+ mchol.get(i, j) = 0.0;
+ int info;
+ dpotrf("L", &rows, mchol.base(), &rows, &info);
+
+ // make vector function
+ MomentFunction func(mchol, imom);
+
+ // permutation schemes
+ WarnockPerScheme wps;
+ ReversePerScheme rps;
+ IdentityPerScheme ips;
+ PermutationScheme *scheme[] = {&wps, &rps, &ips};
+ const char *labs[] = {"Warnock", "Reverse", "Identity"};
+
+ // theoretical result
+ int dim = mchol.numRows();
+ UNormalMoments moments(imom, msq);
+ Vector res((const Vector &)((moments.get(Symmetry(imom)))->getData()));
+
+ // quasi monte carlo normal quadrature
+ double max_error = 0.0;
+ Vector qmc_out(UFSTensor::calcMaxOffset(dim, imom));
+ for (int i = 0; i < 3; i++)
+ {
+ {
+ char mes[100];
+ sprintf(mes, "\tQMC normal quadrature time %8s: ", labs[i]);
+ WallTimer tim(mes);
+ QMCarloNormalQuadrature quad(dim, level, *(scheme[i]));
+ quad.integrate(func, level, num_threads, qmc_out);
+ }
+ qmc_out.add(-1.0, res);
+ printf("\tError %8s: %16.12g\n", labs[i], qmc_out.getMax());
+ if (qmc_out.getMax() > max_error)
+ {
+ max_error = qmc_out.getMax();
+ }
+ }
+
+ return max_error < 1.e-7;
+}
+
+bool
+TestRunnable::smolyak_product_cube(const VectorFunction &func, const Vector &res,
+ double tol, int level)
+{
+ if (res.length() != func.outdim())
+ {
+ fprintf(stderr, "Incompatible dimensions of check value and function.\n");
+ exit(1);
+ }
+
+ GaussLegendre glq;
+ Vector out(func.outdim());
+ double smol_error;
+ double prod_error;
+ {
+ WallTimer tim("\tSmolyak quadrature time: ");
+ SmolyakQuadrature quad(func.indim(), level, glq);
+ quad.integrate(func, level, num_threads, out);
+ out.add(-1.0, res);
+ smol_error = out.getMax();
+ printf("\tNumber of Smolyak evaluations: %d\n", quad.numEvals(level));
+ printf("\tError: %16.12g\n", smol_error);
+ }
+ {
+ WallTimer tim("\tProduct quadrature time: ");
+ ProductQuadrature quad(func.indim(), glq);
+ quad.integrate(func, level, num_threads, out);
+ out.add(-1.0, res);
+ prod_error = out.getMax();
+ printf("\tNumber of product evaluations: %d\n", quad.numEvals(level));
+ printf("\tError: %16.12g\n", prod_error);
+ }
+
+ return smol_error < tol && prod_error < tol;
+}
+
+bool
+TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int level)
+{
+ Vector r(1);
+ double error1;
+ {
+ WallTimer tim("\tQuasi-Monte Carlo (Warnock scrambling) time: ");
+ WarnockPerScheme wps;
+ QMCarloCubeQuadrature qmc(func.indim(), level, wps);
+ // qmc.savePoints("warnock.txt", level);
+ qmc.integrate(func, level, num_threads, r);
+ error1 = std::max(res - r[0], r[0] - res);
+ printf("\tQuasi-Monte Carlo (Warnock scrambling) error: %16.12g\n",
+ error1);
+ }
+ double error2;
+ {
+ WallTimer tim("\tQuasi-Monte Carlo (reverse scrambling) time: ");
+ ReversePerScheme rps;
+ QMCarloCubeQuadrature qmc(func.indim(), level, rps);
+ // qmc.savePoints("reverse.txt", level);
+ qmc.integrate(func, level, num_threads, r);
+ error2 = std::max(res - r[0], r[0] - res);
+ printf("\tQuasi-Monte Carlo (reverse scrambling) error: %16.12g\n",
+ error2);
+ }
+ double error3;
+ {
+ WallTimer tim("\tQuasi-Monte Carlo (no scrambling) time: ");
+ IdentityPerScheme ips;
+ QMCarloCubeQuadrature qmc(func.indim(), level, ips);
+ // qmc.savePoints("identity.txt", level);
+ qmc.integrate(func, level, num_threads, r);
+ error3 = std::max(res - r[0], r[0] - res);
+ printf("\tQuasi-Monte Carlo (no scrambling) error: %16.12g\n",
+ error3);
+ }
+
+ return error1 < tol && error2 < tol && error3 < tol;
+}
+
+/****************************************************/
+/* definition of TestRunnable subclasses */
+/****************************************************/
+class SmolyakNormalMom1 : public TestRunnable
+{
+public:
+ SmolyakNormalMom1()
+ : TestRunnable("Smolyak normal moments (dim=2, level=4, order=4)", 4, 2)
+ {
+ }
+
+ bool
+ run() const
+ {
+ GeneralMatrix m(2, 2);
+ m.zeros(); m.get(0, 0) = 1; m.get(1, 1) = 1;
+ return smolyak_normal_moments(m, 4, 4);
+ }
+};
+
+class SmolyakNormalMom2 : public TestRunnable
+{
+public:
+ SmolyakNormalMom2()
+ : TestRunnable("Smolyak normal moments (dim=3, level=8, order=8)", 8, 3)
+ {
+ }
+
+ bool
+ run() const
+ {
+ GeneralMatrix m(3, 3);
+ m.zeros();
+ m.get(0, 0) = 1; m.get(0, 2) = 0.5; m.get(1, 1) = 1;
+ m.get(1, 0) = 0.5; m.get(2, 2) = 2; m.get(2, 1) = 4;
+ return smolyak_normal_moments(m, 8, 8);
+ }
+};
+
+class ProductNormalMom1 : public TestRunnable
+{
+public:
+ ProductNormalMom1()
+ : TestRunnable("Product normal moments (dim=2, level=4, order=4)", 4, 2)
+ {
+ }
+
+ bool
+ run() const
+ {
+ GeneralMatrix m(2, 2);
+ m.zeros(); m.get(0, 0) = 1; m.get(1, 1) = 1;
+ return product_normal_moments(m, 4, 4);
+ }
+};
+
+class ProductNormalMom2 : public TestRunnable
+{
+public:
+ ProductNormalMom2()
+ : TestRunnable("Product normal moments (dim=3, level=8, order=8)", 8, 3)
+ {
+ }
+
+ bool
+ run() const
+ {
+ GeneralMatrix m(3, 3);
+ m.zeros();
+ m.get(0, 0) = 1; m.get(0, 2) = 0.5; m.get(1, 1) = 1;
+ m.get(1, 0) = 0.5; m.get(2, 2) = 2; m.get(2, 1) = 4;
+ return product_normal_moments(m, 8, 8);
+ }
+};
+
+class QMCNormalMom1 : public TestRunnable
+{
+public:
+ QMCNormalMom1()
+ : TestRunnable("QMC normal moments (dim=2, level=1000, order=4)", 4, 2)
+ {
+ }
+
+ bool
+ run() const
+ {
+ GeneralMatrix m(2, 2);
+ m.zeros(); m.get(0, 0) = 1; m.get(1, 1) = 1;
+ return qmc_normal_moments(m, 4, 1000);
+ }
+};
+
+class QMCNormalMom2 : public TestRunnable
+{
+public:
+ QMCNormalMom2()
+ : TestRunnable("QMC normal moments (dim=3, level=10000, order=8)", 8, 3)
+ {
+ }
+
+ bool
+ run() const
+ {
+ GeneralMatrix m(3, 3);
+ m.zeros();
+ m.get(0, 0) = 1; m.get(0, 2) = 0.5; m.get(1, 1) = 1;
+ m.get(1, 0) = 0.5; m.get(2, 2) = 2; m.get(2, 1) = 4;
+ return qmc_normal_moments(m, 8, 10000);
+ }
+};
+
+// note that here we pass 1,1 to tls since smolyak has its own PascalTriangle
+class F1GaussLegendre : public TestRunnable
+{
+public:
+ F1GaussLegendre()
+ : TestRunnable("Function1 Gauss-Legendre (dim=6, level=13", 1, 1)
+ {
+ }
+
+ bool
+ run() const
+ {
+ Function1Trans f1(6);
+ Vector res(1); res[0] = 1.0;
+ return smolyak_product_cube(f1, res, 1e-2, 13);
+ }
+};
+
+class F1QuasiMCarlo : public TestRunnable
+{
+public:
+ F1QuasiMCarlo()
+ : TestRunnable("Function1 Quasi-Monte Carlo (dim=6, level=1000000)", 1, 1)
+ {
+ }
+
+ bool
+ run() const
+ {
+ Function1 f1(6);
+ return qmc_cube(f1, 1.0, 1.e-4, 1000000);
+ }
+};
+
+int
+main()
+{
+ TestRunnable *all_tests[50];
+ // fill in vector of all tests
+ int num_tests = 0;
+ all_tests[num_tests++] = new SmolyakNormalMom1();
+ all_tests[num_tests++] = new SmolyakNormalMom2();
+ all_tests[num_tests++] = new ProductNormalMom1();
+ all_tests[num_tests++] = new ProductNormalMom2();
+ all_tests[num_tests++] = new QMCNormalMom1();
+ all_tests[num_tests++] = new QMCNormalMom2();
+ /*
+ all_tests[num_tests++] = new F1GaussLegendre();
+ all_tests[num_tests++] = new F1QuasiMCarlo();
+ */
+ // find maximum dimension and maximum nvar
+ int dmax = 0;
+ int nvmax = 0;
+ for (int i = 0; i < num_tests; i++)
+ {
+ if (dmax < all_tests[i]->dim)
+ dmax = all_tests[i]->dim;
+ if (nvmax < all_tests[i]->nvar)
+ nvmax = all_tests[i]->nvar;
+ }
+ tls.init(dmax, nvmax); // initialize library
+ THREAD_GROUP::max_parallel_threads = num_threads;
+
+ // launch the tests
+ int success = 0;
+ for (int i = 0; i < num_tests; i++)
+ {
+ try
+ {
+ if (all_tests[i]->test())
+ success++;
+ }
+ catch (const TLException &e)
+ {
+ printf("Caugth TL exception in <%s>:\n", all_tests[i]->getName());
+ e.print();
+ }
+ catch (SylvException &e)
+ {
+ printf("Caught Sylv exception in <%s>:\n", all_tests[i]->getName());
+ e.printMessage();
+ }
+ }
+
+ printf("There were %d tests that failed out of %d tests run.\n",
+ num_tests - success, num_tests);
+
+ // destroy
+ for (int i = 0; i < num_tests; i++)
+ {
+ delete all_tests[i];
+ }
+
+ return 0;
+}
diff --git a/dynare++/integ/testing/tests.cpp b/dynare++/integ/testing/tests.cpp
deleted file mode 100644
index 7862a08373cfc714eadbe58d91c0b9ffe41daef4..0000000000000000000000000000000000000000
--- a/dynare++/integ/testing/tests.cpp
+++ /dev/null
@@ -1,544 +0,0 @@
-/* $Id: tests.cpp 431 2005-08-16 15:41:01Z kamenik $ */
-/* Copyright 2005, Ondra Kamenik */
-
-#include "GeneralMatrix.h"
-#include <dynlapack.h>
-#include "SylvException.h"
-
-#include "rfs_tensor.h"
-#include "normal_moments.h"
-
-#include "vector_function.h"
-#include "quadrature.h"
-#include "smolyak.h"
-#include "product.h"
-#include "quasi_mcarlo.h"
-
-#include <cstdio>
-#include <cstring>
-#include <sys/time.h>
-#include <cmath>
-
-const int num_threads = 2; // does nothing if DEBUG defined
-
-// evaluates unfolded (Dx)^k power, where x is a vector, D is a
-// Cholesky factor (lower triangular)
-class MomentFunction : public VectorFunction {
- GeneralMatrix D;
- int k;
-public:
- MomentFunction(const GeneralMatrix& inD, int kk)
- : VectorFunction(inD.numRows(), UFSTensor::calcMaxOffset(inD.numRows(), kk)),
- D(inD), k(kk) {}
- MomentFunction(const MomentFunction& func)
- : VectorFunction(func), D(func.D), k(func.k) {}
- VectorFunction* clone() const
- {return new MomentFunction(*this);}
- void eval(const Vector& point, const ParameterSignal& sig, Vector& out);
-};
-
-void MomentFunction::eval(const Vector& point, const ParameterSignal& sig, Vector& out)
-{
- if (point.length() != indim() || out.length() != outdim()) {
- printf("Wrong length of vectors in MomentFunction::eval\n");
- exit(1);
- }
- Vector y(point);
- y.zeros();
- D.multaVec(y, point);
- URSingleTensor ypow(y, k);
- out.zeros();
- out.add(1.0, ypow.getData());
-}
-
-class TensorPower : public VectorFunction {
- int k;
-public:
- TensorPower(int nvar, int kk)
- : VectorFunction(nvar, UFSTensor::calcMaxOffset(nvar, kk)), k(kk) {}
- TensorPower(const TensorPower& func)
- : VectorFunction(func), k(func.k) {}
- VectorFunction* clone() const
- {return new TensorPower(*this);}
- void eval(const Vector& point, const ParameterSignal& sig, Vector& out);
-};
-
-void TensorPower::eval(const Vector& point, const ParameterSignal& sig, Vector& out)
-{
- if (point.length() != indim() || out.length() != outdim()) {
- printf("Wrong length of vectors in TensorPower::eval\n");
- exit(1);
- }
- URSingleTensor ypow(point, k);
- out.zeros();
- out.add(1.0, ypow.getData());
-}
-
-
-// evaluates (1+1/d)^d*(x_1*...*x_d)^(1/d), its integral over <0,1>^d
-// is 1.0, and its variation grows exponetially
-// d = dim
-class Function1 : public VectorFunction {
- int dim;
-public:
- Function1(int d)
- : VectorFunction(d, 1), dim(d) {}
- Function1(const Function1& f)
- : VectorFunction(f.indim(), f.outdim()), dim(f.dim) {}
- VectorFunction* clone() const
- {return new Function1(*this);}
- virtual void eval(const Vector& point, const ParameterSignal& sig, Vector& out);
-};
-
-void Function1::eval(const Vector& point, const ParameterSignal& sig, Vector& out)
-{
- if (point.length() != dim || out.length() != 1) {
- printf("Wrong length of vectors in Function1::eval\n");
- exit(1);
- }
- double r = 1;
- for (int i = 0; i < dim; i++)
- r *= point[i];
- r = pow(r, 1.0/dim);
- r *= pow(1.0 + 1.0/dim, (double)dim);
- out[0] = r;
-}
-
-// evaluates Function1 but with transformation x_i=0.5(y_i+1)
-// this makes the new function integrate over <-1,1>^d to 1.0
-class Function1Trans : public Function1 {
-public:
- Function1Trans(int d)
- : Function1(d) {}
- Function1Trans(const Function1Trans& func)
- : Function1(func) {}
- VectorFunction* clone() const
- {return new Function1Trans(*this);}
- virtual void eval(const Vector& point, const ParameterSignal& sig, Vector& out);
-};
-
-void Function1Trans::eval(const Vector& point, const ParameterSignal& sig, Vector& out)
-{
- Vector p(point.length());
- for (int i = 0; i < p.length(); i++)
- p[i] = 0.5*(point[i]+1);
- Function1::eval(p, sig, out);
- out.mult(pow(0.5,indim()));
-}
-
-
-// WallTimer class. Constructor saves the wall time, destructor
-// cancels the current time from the saved, and prints the message
-// with time information
-class WallTimer {
- char mes[100];
- struct timeval start;
- bool new_line;
-public:
- WallTimer(const char* m, bool nl = true)
- {strcpy(mes, m);new_line = nl; gettimeofday(&start, NULL);}
- ~WallTimer()
- {
- struct timeval end;
- gettimeofday(&end, NULL);
- printf("%s%8.4g", mes,
- end.tv_sec-start.tv_sec + (end.tv_usec-start.tv_usec)*1.0e-6);
- if (new_line)
- printf("\n");
- }
-};
-
-/****************************************************/
-/* declaration of TestRunnable class */
-/****************************************************/
-class TestRunnable {
- char name[100];
-public:
- int dim; // dimension of the solved problem
- int nvar; // number of variable of the solved problem
- TestRunnable(const char* n, int d, int nv)
- : dim(d), nvar(nv)
- {strncpy(name, n, 100);}
- bool test() const;
- virtual bool run() const =0;
- const char* getName() const
- {return name;}
-protected:
- static bool smolyak_normal_moments(const GeneralMatrix& m, int imom, int level);
- static bool product_normal_moments(const GeneralMatrix& m, int imom, int level);
- static bool qmc_normal_moments(const GeneralMatrix& m, int imom, int level);
- static bool smolyak_product_cube(const VectorFunction& func, const Vector& res,
- double tol, int level);
- static bool qmc_cube(const VectorFunction& func, double res, double tol, int level);
-};
-
-bool TestRunnable::test() const
-{
- printf("Running test <%s>\n",name);
- bool passed;
- {
- WallTimer tim("Wall clock time ", false);
- passed = run();
- }
- if (passed) {
- printf("............................ passed\n\n");
- return passed;
- } else {
- printf("............................ FAILED\n\n");
- return passed;
- }
-}
-
-
-/****************************************************/
-/* definition of TestRunnable static methods */
-/****************************************************/
-bool TestRunnable::smolyak_normal_moments(const GeneralMatrix& m, int imom, int level)
-{
- // first make m*m' and then Cholesky factor
- GeneralMatrix mtr(m, "transpose");
- GeneralMatrix msq(m, mtr);
-
- // make vector function
- int dim = m.numRows();
- TensorPower tp(dim, imom);
- GaussConverterFunction func(tp, msq);
-
- // smolyak quadrature
- Vector smol_out(UFSTensor::calcMaxOffset(dim, imom));
- {
- WallTimer tim("\tSmolyak quadrature time: ");
- GaussHermite gs;
- SmolyakQuadrature quad(dim, level, gs);
- quad.integrate(func, level, num_threads, smol_out);
- printf("\tNumber of Smolyak evaluations: %d\n", quad.numEvals(level));
- }
-
- // check against theoretical moments
- UNormalMoments moments(imom, msq);
- smol_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
- printf("\tError: %16.12g\n", smol_out.getMax());
- return smol_out.getMax() < 1.e-7;
-}
-
-bool TestRunnable::product_normal_moments(const GeneralMatrix& m, int imom, int level)
-{
- // first make m*m' and then Cholesky factor
- GeneralMatrix mtr(m, "transpose");
- GeneralMatrix msq(m, mtr);
-
- // make vector function
- int dim = m.numRows();
- TensorPower tp(dim, imom);
- GaussConverterFunction func(tp, msq);
-
- // product quadrature
- Vector prod_out(UFSTensor::calcMaxOffset(dim, imom));
- {
- WallTimer tim("\tProduct quadrature time: ");
- GaussHermite gs;
- ProductQuadrature quad(dim, gs);
- quad.integrate(func, level, num_threads, prod_out);
- printf("\tNumber of product evaluations: %d\n", quad.numEvals(level));
- }
-
- // check against theoretical moments
- UNormalMoments moments(imom, msq);
- prod_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
- printf("\tError: %16.12g\n", prod_out.getMax());
- return prod_out.getMax() < 1.e-7;
-}
-
-bool TestRunnable::qmc_normal_moments(const GeneralMatrix& m, int imom, int level)
-{
- // first make m*m' and then Cholesky factor
- GeneralMatrix mtr(m, "transpose");
- GeneralMatrix msq(m, mtr);
- GeneralMatrix mchol(msq);
- int rows = mchol.numRows();
- for (int i = 0; i < rows; i++)
- for (int j = i+1; j < rows; j++)
- mchol.get(i,j) = 0.0;
- int info;
- dpotrf("L", &rows, mchol.base(), &rows, &info);
-
- // make vector function
- MomentFunction func(mchol, imom);
-
- // permutation schemes
- WarnockPerScheme wps;
- ReversePerScheme rps;
- IdentityPerScheme ips;
- PermutationScheme* scheme[] = {&wps, &rps, &ips};
- const char* labs[] = {"Warnock", "Reverse", "Identity"};
-
- // theoretical result
- int dim = mchol.numRows();
- UNormalMoments moments(imom, msq);
- Vector res((const Vector&)((moments.get(Symmetry(imom)))->getData()));
-
- // quasi monte carlo normal quadrature
- double max_error = 0.0;
- Vector qmc_out(UFSTensor::calcMaxOffset(dim, imom));
- for (int i = 0; i < 3; i++) {
- {
- char mes[100];
- sprintf(mes, "\tQMC normal quadrature time %8s: ", labs[i]);
- WallTimer tim(mes);
- QMCarloNormalQuadrature quad(dim, level, *(scheme[i]));
- quad.integrate(func, level, num_threads, qmc_out);
- }
- qmc_out.add(-1.0, res);
- printf("\tError %8s: %16.12g\n", labs[i], qmc_out.getMax());
- if (qmc_out.getMax() > max_error) {
- max_error = qmc_out.getMax();
- }
- }
-
- return max_error < 1.e-7;
-}
-
-
-bool TestRunnable::smolyak_product_cube(const VectorFunction& func, const Vector& res,
- double tol, int level)
-{
- if (res.length() != func.outdim()) {
- fprintf(stderr, "Incompatible dimensions of check value and function.\n");
- exit(1);
- }
-
- GaussLegendre glq;
- Vector out(func.outdim());
- double smol_error;
- double prod_error;
- {
- WallTimer tim("\tSmolyak quadrature time: ");
- SmolyakQuadrature quad(func.indim(), level, glq);
- quad.integrate(func, level, num_threads, out);
- out.add(-1.0, res);
- smol_error = out.getMax();
- printf("\tNumber of Smolyak evaluations: %d\n", quad.numEvals(level));
- printf("\tError: %16.12g\n", smol_error);
- }
- {
- WallTimer tim("\tProduct quadrature time: ");
- ProductQuadrature quad(func.indim(), glq);
- quad.integrate(func, level, num_threads, out);
- out.add(-1.0, res);
- prod_error = out.getMax();
- printf("\tNumber of product evaluations: %d\n", quad.numEvals(level));
- printf("\tError: %16.12g\n", prod_error);
- }
-
- return smol_error < tol && prod_error < tol;
-}
-
-bool TestRunnable::qmc_cube(const VectorFunction& func, double res, double tol, int level)
-{
- Vector r(1);
- double error1;
- {
- WallTimer tim("\tQuasi-Monte Carlo (Warnock scrambling) time: ");
- WarnockPerScheme wps;
- QMCarloCubeQuadrature qmc(func.indim(), level, wps);
-// qmc.savePoints("warnock.txt", level);
- qmc.integrate(func, level, num_threads, r);
- error1 = std::max(res - r[0], r[0] - res);
- printf("\tQuasi-Monte Carlo (Warnock scrambling) error: %16.12g\n",
- error1);
- }
- double error2;
- {
- WallTimer tim("\tQuasi-Monte Carlo (reverse scrambling) time: ");
- ReversePerScheme rps;
- QMCarloCubeQuadrature qmc(func.indim(), level, rps);
-// qmc.savePoints("reverse.txt", level);
- qmc.integrate(func, level, num_threads, r);
- error2 = std::max(res - r[0], r[0] - res);
- printf("\tQuasi-Monte Carlo (reverse scrambling) error: %16.12g\n",
- error2);
- }
- double error3;
- {
- WallTimer tim("\tQuasi-Monte Carlo (no scrambling) time: ");
- IdentityPerScheme ips;
- QMCarloCubeQuadrature qmc(func.indim(), level, ips);
-// qmc.savePoints("identity.txt", level);
- qmc.integrate(func, level, num_threads, r);
- error3 = std::max(res - r[0], r[0] - res);
- printf("\tQuasi-Monte Carlo (no scrambling) error: %16.12g\n",
- error3);
- }
-
-
- return error1 < tol && error2 < tol && error3 < tol;
-}
-
-/****************************************************/
-/* definition of TestRunnable subclasses */
-/****************************************************/
-class SmolyakNormalMom1 : public TestRunnable {
-public:
- SmolyakNormalMom1()
- : TestRunnable("Smolyak normal moments (dim=2, level=4, order=4)", 4, 2) {}
-
- bool run() const
- {
- GeneralMatrix m(2,2);
- m.zeros(); m.get(0,0)=1; m.get(1,1)=1;
- return smolyak_normal_moments(m, 4, 4);
- }
-};
-
-class SmolyakNormalMom2 : public TestRunnable {
-public:
- SmolyakNormalMom2()
- : TestRunnable("Smolyak normal moments (dim=3, level=8, order=8)", 8, 3) {}
-
- bool run() const
- {
- GeneralMatrix m(3,3);
- m.zeros();
- m.get(0,0)=1; m.get(0,2)=0.5; m.get(1,1)=1;
- m.get(1,0)=0.5;m.get(2,2)=2;m.get(2,1)=4;
- return smolyak_normal_moments(m, 8, 8);
- }
-};
-
-class ProductNormalMom1 : public TestRunnable {
-public:
- ProductNormalMom1()
- : TestRunnable("Product normal moments (dim=2, level=4, order=4)", 4, 2) {}
-
- bool run() const
- {
- GeneralMatrix m(2,2);
- m.zeros(); m.get(0,0)=1; m.get(1,1)=1;
- return product_normal_moments(m, 4, 4);
- }
-};
-
-class ProductNormalMom2 : public TestRunnable {
-public:
- ProductNormalMom2()
- : TestRunnable("Product normal moments (dim=3, level=8, order=8)", 8, 3) {}
-
- bool run() const
- {
- GeneralMatrix m(3,3);
- m.zeros();
- m.get(0,0)=1; m.get(0,2)=0.5; m.get(1,1)=1;
- m.get(1,0)=0.5;m.get(2,2)=2;m.get(2,1)=4;
- return product_normal_moments(m, 8, 8);
- }
-};
-
-class QMCNormalMom1 : public TestRunnable {
-public:
- QMCNormalMom1()
- : TestRunnable("QMC normal moments (dim=2, level=1000, order=4)", 4, 2) {}
-
- bool run() const
- {
- GeneralMatrix m(2,2);
- m.zeros(); m.get(0,0)=1; m.get(1,1)=1;
- return qmc_normal_moments(m, 4, 1000);
- }
-};
-
-class QMCNormalMom2 : public TestRunnable {
-public:
- QMCNormalMom2()
- : TestRunnable("QMC normal moments (dim=3, level=10000, order=8)", 8, 3) {}
-
- bool run() const
- {
- GeneralMatrix m(3,3);
- m.zeros();
- m.get(0,0)=1; m.get(0,2)=0.5; m.get(1,1)=1;
- m.get(1,0)=0.5;m.get(2,2)=2;m.get(2,1)=4;
- return qmc_normal_moments(m, 8, 10000);
- }
-};
-
-
-
-// note that here we pass 1,1 to tls since smolyak has its own PascalTriangle
-class F1GaussLegendre : public TestRunnable {
-public:
- F1GaussLegendre()
- : TestRunnable("Function1 Gauss-Legendre (dim=6, level=13", 1, 1) {}
-
- bool run() const
- {
- Function1Trans f1(6);
- Vector res(1); res[0] = 1.0;
- return smolyak_product_cube(f1, res, 1e-2, 13);
- }
-};
-
-
-class F1QuasiMCarlo : public TestRunnable {
-public:
- F1QuasiMCarlo()
- : TestRunnable("Function1 Quasi-Monte Carlo (dim=6, level=1000000)", 1, 1) {}
-
- bool run() const
- {
- Function1 f1(6);
- return qmc_cube(f1, 1.0, 1.e-4, 1000000);
- }
-};
-
-int main()
-{
- TestRunnable* all_tests[50];
- // fill in vector of all tests
- int num_tests = 0;
- all_tests[num_tests++] = new SmolyakNormalMom1();
- all_tests[num_tests++] = new SmolyakNormalMom2();
- all_tests[num_tests++] = new ProductNormalMom1();
- all_tests[num_tests++] = new ProductNormalMom2();
- all_tests[num_tests++] = new QMCNormalMom1();
- all_tests[num_tests++] = new QMCNormalMom2();
-/*
- all_tests[num_tests++] = new F1GaussLegendre();
- all_tests[num_tests++] = new F1QuasiMCarlo();
-*/
- // find maximum dimension and maximum nvar
- int dmax=0;
- int nvmax = 0;
- for (int i = 0; i < num_tests; i++) {
- if (dmax < all_tests[i]->dim)
- dmax = all_tests[i]->dim;
- if (nvmax < all_tests[i]->nvar)
- nvmax = all_tests[i]->nvar;
- }
- tls.init(dmax, nvmax); // initialize library
- THREAD_GROUP::max_parallel_threads = num_threads;
-
- // launch the tests
- int success = 0;
- for (int i = 0; i < num_tests; i++) {
- try {
- if (all_tests[i]->test())
- success++;
- } catch (const TLException& e) {
- printf("Caugth TL exception in <%s>:\n", all_tests[i]->getName());
- e.print();
- } catch (SylvException& e) {
- printf("Caught Sylv exception in <%s>:\n", all_tests[i]->getName());
- e.printMessage();
- }
- }
-
- printf("There were %d tests that failed out of %d tests run.\n",
- num_tests - success, num_tests);
-
- // destroy
- for (int i = 0; i < num_tests; i++) {
- delete all_tests[i];
- }
-
- return 0;
-}
diff --git a/dynare++/kord/global_check.cc b/dynare++/kord/global_check.cc
index b2549000972e7a0dba1cd7c79f7ed0ebb981fd67..d0a5556ef194d01198db5240bdfbda018ee7ad70 100644
--- a/dynare++/kord/global_check.cc
+++ b/dynare++/kord/global_check.cc
@@ -4,9 +4,9 @@
#include "global_check.hh"
-#include "smolyak.h"
-#include "product.h"
-#include "quasi_mcarlo.h"
+#include "smolyak.hh"
+#include "product.hh"
+#include "quasi_mcarlo.hh"
#ifdef __MINGW32__
# define __CROSS_COMPILATION__
diff --git a/dynare++/kord/global_check.hh b/dynare++/kord/global_check.hh
index f11b9a49a0f8af512ddfd5881663b13c34892d18..ed1b1a1a6c625d4685744f7d511c7750e317e1e9 100644
--- a/dynare++/kord/global_check.hh
+++ b/dynare++/kord/global_check.hh
@@ -37,8 +37,8 @@
#include <matio.h>
-#include "vector_function.h"
-#include "quadrature.h"
+#include "vector_function.hh"
+#include "quadrature.hh"
#include "dynamic_model.hh"
#include "journal.hh"
diff --git a/mex/build/libdynare++.am b/mex/build/libdynare++.am
index 5f76584fc8f29e8eb5b59dbee456db222564f2db..f65b31eadfa792267212fda5881ed283e5f4689e 100644
--- a/mex/build/libdynare++.am
+++ b/mex/build/libdynare++.am
@@ -76,11 +76,16 @@ TL_SRCS = \
$(TOPDIR)/tl/cc/tl_static.cpp
INTEG_SRCS = \
- $(TOPDIR)/integ/cc/product.cpp \
- $(TOPDIR)/integ/cc/quadrature.cpp \
- $(TOPDIR)/integ/cc/quasi_mcarlo.cpp \
- $(TOPDIR)/integ/cc/smolyak.cpp \
- $(TOPDIR)/integ/cc/vector_function.cpp
+ $(TOPDIR)/integ/cc/quadrature.cc \
+ $(TOPDIR)/integ/cc/quadrature.hh \
+ $(TOPDIR)/integ/cc/quasi_mcarlo.cc \
+ $(TOPDIR)/integ/cc/quasi_mcarlo.hh \
+ $(TOPDIR)/integ/cc/product.cc \
+ $(TOPDIR)/integ/cc/product.hh \
+ $(TOPDIR)/integ/cc/smolyak.cc \
+ $(TOPDIR)/integ/cc/smolyak.hh \
+ $(TOPDIR)/integ/cc/vector_function.cc \
+ $(TOPDIR)/integ/cc/vector_function.hh
nodist_libdynare___a_SOURCES = \
$(KORD_SRCS) \
diff --git a/windows/dynare.nsi b/windows/dynare.nsi
index aee97f89fcc443315baa4b7641c25bbfc91e6721..cda1a0cd5b2a88ebb78d87c9ad1e725529873373 100644
--- a/windows/dynare.nsi
+++ b/windows/dynare.nsi
@@ -146,7 +146,7 @@ Section "Documentation and examples (Dynare and Dynare++)"
File ..\doc\dynare.html\*.html ..\doc\dynare.html\*.png
SetOutPath $INSTDIR\doc\dynare++
- File ..\dynare++\doc\dynare++-tutorial.pdf ..\dynare++\doc\dynare++-ramsey.pdf ..\dynare++\sylv\sylvester.pdf ..\dynare++\tl\cc\tl.pdf ..\dynare++\integ\cc\integ.pdf ..\dynare++\kord\kord.pdf
+ File ..\dynare++\doc\dynare++-tutorial.pdf ..\dynare++\doc\dynare++-ramsey.pdf ..\dynare++\sylv\sylvester.pdf ..\dynare++\tl\cc\tl.pdf
CreateShortcut "${SMLOC}\Documentation.lnk" "$INSTDIR\doc"