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

Dynare++: various modernizations of numerical integration library

time-shift
Sébastien Villemot 2019-01-14 16:09:49 +01:00
parent 4a72266d05
commit b8791e9f13
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
12 changed files with 339 additions and 593 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

119
dynare++/integ/cc/product.cc
View File

@ -3,111 +3,47 @@
#include "product.hh" #include "product.hh"
#include "symmetry.hh" #include "symmetry.hh"
prodpit::prodpit() #include <iostream>
#include <iomanip>
{
}
/* This constructs a product iterator corresponding to index $(j0,0\ldots,0)$. */ /* This constructs a product iterator corresponding to index $(j0,0\ldots,0)$. */
prodpit::prodpit(const ProductQuadrature &q, int j0, int l) prodpit::prodpit(const ProductQuadrature &q, int j0, int l)
: prodq(&q), level(l), npoints(q.uquad.numPoints(l)), jseq(new IntSequence(q.dimen(), 0)), : prodq(q), level(l), npoints(q.uquad.numPoints(l)),
end_flag(false), sig(new ParameterSignal(q.dimen())), p(new Vector(q.dimen())) jseq{q.dimen(), 0},
end_flag(false),
sig{q.dimen()},
p{q.dimen()}
{ {
if (j0 < npoints) if (j0 < npoints)
{ {
(*jseq)[0] = j0; jseq[0] = j0;
setPointAndWeight(); setPointAndWeight();
} }
else else
{ end_flag = true;
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 = nullptr;
if (ppit.sig)
sig = new ParameterSignal(*(ppit.sig));
else
sig = nullptr;
if (ppit.p)
p = new Vector(*(ppit.p));
else
p = nullptr;
}
prodpit::~prodpit()
{
if (jseq)
delete jseq;
if (sig)
delete sig;
if (p)
delete p;
} }
bool bool
prodpit::operator==(const prodpit &ppit) const prodpit::operator==(const prodpit &ppit) const
{ {
bool ret = true; return &prodq == &ppit.prodq && end_flag == ppit.end_flag && jseq == ppit.jseq;
ret = ret & prodq == ppit.prodq;
ret = ret & end_flag == ppit.end_flag;
ret = ret & ((jseq == nullptr && ppit.jseq == nullptr)
|| (jseq != nullptr && ppit.jseq != nullptr && *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 = nullptr;
if (ppit.sig)
sig = new ParameterSignal(*(ppit.sig));
else
sig = nullptr;
if (ppit.p)
p = new Vector(*(ppit.p));
else
p = nullptr;
return *this;
} }
prodpit & prodpit &
prodpit::operator++() prodpit::operator++()
{ {
// todo: throw if |prodq==NULL| or |jseq==NULL| or |sig==NULL| or |end_flag==true| // todo: throw if |prodq==NULL| or |jseq==NULL| or |sig==NULL| or |end_flag==true|
int i = prodq->dimen()-1; int i = prodq.dimen()-1;
(*jseq)[i]++; jseq[i]++;
while (i >= 0 && (*jseq)[i] == npoints) while (i >= 0 && jseq[i] == npoints)
{ {
(*jseq)[i] = 0; jseq[i] = 0;
i--; i--;
if (i >= 0) if (i >= 0)
(*jseq)[i]++; jseq[i]++;
} }
sig->signalAfter(std::max(i, 0)); sig.signalAfter(std::max(i, 0));
if (i == -1) if (i == -1)
end_flag = true; end_flag = true;
@ -126,10 +62,10 @@ prodpit::setPointAndWeight()
// todo: raise if |prodq==NULL| or |jseq==NULL| or |sig==NULL| or // todo: raise if |prodq==NULL| or |jseq==NULL| or |sig==NULL| or
// |p==NULL| or |end_flag==true| // |p==NULL| or |end_flag==true|
w = 1.0; w = 1.0;
for (int i = 0; i < prodq->dimen(); i++) for (int i = 0; i < prodq.dimen(); i++)
{ {
(*p)[i] = (prodq->uquad).point(level, (*jseq)[i]); p[i] = (prodq.uquad).point(level, jseq[i]);
w *= (prodq->uquad).weight(level, (*jseq)[i]); w *= (prodq.uquad).weight(level, jseq[i]);
} }
} }
@ -138,13 +74,16 @@ prodpit::setPointAndWeight()
void void
prodpit::print() const prodpit::print() const
{ {
printf("j=["); auto ff = std::cout.flags();
for (int i = 0; i < prodq->dimen(); i++) std::cout << "j=[";
printf("%2d ", (*jseq)[i]); for (int i = 0; i < prodq.dimen(); i++)
printf("] %+4.3f*(", w); std::cout << std::setw(2) << jseq[i];
for (int i = 0; i < prodq->dimen()-1; i++) std::cout << std::showpos << std::fixed << std::setprecision(3)
printf("%+4.3f ", (*p)[i]); << "] " << std::setw(4) << w << "*(";
printf("%+4.3f)\n", (*p)[prodq->dimen()-1]); for (int i = 0; i < prodq.dimen()-1; i++)
std::cout << std::setw(4) << p[i] << ' ';
std::cout << std::setw(4) << p[prodq.dimen()-1] << ')' << std::endl;
std::cout.flags(ff);
} }
ProductQuadrature::ProductQuadrature(int d, const OneDQuadrature &uq) ProductQuadrature::ProductQuadrature(int d, const OneDQuadrature &uq)

23
dynare++/integ/cc/product.hh
View File

@ -38,36 +38,36 @@ class ProductQuadrature;
class prodpit class prodpit
{ {
protected: protected:
const ProductQuadrature *prodq{nullptr}; const ProductQuadrature &prodq;
int level{0}; int level{0};
int npoints{0}; int npoints{0};
IntSequence *jseq{nullptr}; IntSequence jseq;
bool end_flag{true}; bool end_flag{true};
ParameterSignal *sig{nullptr}; ParameterSignal sig;
Vector *p{nullptr}; Vector p;
double w; double w;
public: public:
prodpit(); prodpit() = default;
prodpit(const ProductQuadrature &q, int j0, int l); prodpit(const ProductQuadrature &q, int j0, int l);
prodpit(const prodpit &ppit); prodpit(const prodpit &ppit) = default;
~prodpit(); ~prodpit() = default;
bool operator==(const prodpit &ppit) const; bool operator==(const prodpit &ppit) const;
bool bool
operator!=(const prodpit &ppit) const operator!=(const prodpit &ppit) const
{ {
return !operator==(ppit); return !operator==(ppit);
} }
const prodpit &operator=(const prodpit &spit); prodpit &operator=(const prodpit &spit) = delete;
prodpit &operator++(); prodpit &operator++();
const ParameterSignal & const ParameterSignal &
signal() const signal() const
{ {
return *sig; return sig;
} }
const Vector & const Vector &
point() const point() const
{ {
return *p; return p;
} }
double double
weight() const weight() const
@ -93,8 +93,7 @@ class ProductQuadrature : public QuadratureImpl<prodpit>
const OneDQuadrature &uquad; const OneDQuadrature &uquad;
public: public:
ProductQuadrature(int d, const OneDQuadrature &uq); ProductQuadrature(int d, const OneDQuadrature &uq);
~ProductQuadrature() ~ProductQuadrature() override = default;
override = default;
int int
numEvals(int l) const override numEvals(int l) const override
{ {

33
dynare++/integ/cc/quadrature.hh
View File

@ -30,7 +30,10 @@
#ifndef QUADRATURE_H #ifndef QUADRATURE_H
#define QUADRATURE_H #define QUADRATURE_H
#include <cstdlib> #include <iostream>
#include <fstream>
#include <iomanip>
#include "vector_function.hh" #include "vector_function.hh"
#include "int_sequence.hh" #include "int_sequence.hh"
#include "sthread.hh" #include "sthread.hh"
@ -43,8 +46,7 @@
class OneDQuadrature class OneDQuadrature
{ {
public: public:
virtual ~OneDQuadrature() virtual ~OneDQuadrature() = default;
= default;
virtual int numLevels() const = 0; virtual int numLevels() const = 0;
virtual int numPoints(int level) const = 0; virtual int numPoints(int level) const = 0;
virtual double point(int level, int i) const = 0; virtual double point(int level, int i) const = 0;
@ -72,8 +74,7 @@ public:
Quadrature(int d) : dim(d) Quadrature(int d) : dim(d)
{ {
} }
virtual ~Quadrature() virtual ~Quadrature() = default;
= default;
int int
dimen() const dimen() const
{ {
@ -183,25 +184,26 @@ public:
/* Just for debugging. */ /* Just for debugging. */
void void
savePoints(const char *fname, int level) const savePoints(const string &fname, int level) const
{ {
FILE *fd; ofstream fd{fname, std::ios::out | std::ios::trunc};
if (nullptr == (fd = fopen(fname, "w"))) if (fd.fail())
{ {
// todo: raise // todo: raise
fprintf(stderr, "Cannot open file %s for writing.\n", fname); std::cerr << "Cannot open file " << fname << " for writing." << std::endl;
exit(1); exit(EXIT_FAILURE);
} }
_Tpit beg = begin(0, 1, level); _Tpit beg = begin(0, 1, level);
_Tpit end = begin(1, 1, level); _Tpit end = begin(1, 1, level);
fd << std::setprecision(12);
for (_Tpit run = beg; run != end; ++run) for (_Tpit run = beg; run != end; ++run)
{ {
fprintf(fd, "%16.12g", run.weight()); fd << std::setw(16) << run.weight();
for (int i = 0; i < dimen(); i++) for (int i = 0; i < dimen(); i++)
fprintf(fd, "\t%16.12g", run.point()[i]); fd << '\t' << std::setw(16) << run.point()[i];
fprintf(fd, "\n"); fd << std::endl;
} }
fclose(fd); fd.close();
} }
_Tpit _Tpit
@ -244,8 +246,7 @@ public:
{ {
calcOffsets(); calcOffsets();
} }
~OneDPrecalcQuadrature() ~OneDPrecalcQuadrature() override = default;
override = default;
int int
numLevels() const override numLevels() const override
{ {

114
dynare++/integ/cc/quasi_mcarlo.cc
View File

@ -3,6 +3,8 @@
#include "quasi_mcarlo.hh" #include "quasi_mcarlo.hh"
#include <cmath> #include <cmath>
#include <iostream>
#include <iomanip>
/* Here in the constructor, we have to calculate a maximum length of /* Here in the constructor, we have to calculate a maximum length of
|coeff| array for a given |base| and given maximum |maxn|. After |coeff| array for a given |base| and given maximum |maxn|. After
@ -10,7 +12,7 @@
RadicalInverse::RadicalInverse(int n, int b, int mxn) RadicalInverse::RadicalInverse(int n, int b, int mxn)
: num(n), base(b), maxn(mxn), : num(n), base(b), maxn(mxn),
coeff((int) (floor(log((double)maxn)/log((double)b))+2), 0) coeff(static_cast<int>(floor(log(static_cast<double>(maxn))/log(static_cast<double>(b)))+2), 0)
{ {
int nr = num; int nr = num;
j = -1; j = -1;
@ -75,15 +77,14 @@ RadicalInverse::increase()
void void
RadicalInverse::print() const RadicalInverse::print() const
{ {
printf("n=%d b=%d c=", num, base); std::cout << "n=" << num << " b=" << base << " c=";
coeff.print(); coeff.print();
} }
/* Here we have the first 170 primes. This means that we are not able /* Here we have the first 170 primes. This means that we are not able
to integrate dimensions greater than 170. */ to integrate dimensions greater than 170. */
int HaltonSequence::num_primes = 170; std::array<int, 170> HaltonSequence::primes = {
int HaltonSequence::primes[] = {
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
@ -116,18 +117,6 @@ HaltonSequence::HaltonSequence(int n, int mxn, int dim, const PermutationScheme
eval(); eval();
} }
const HaltonSequence &
HaltonSequence::operator=(const HaltonSequence &hs)
{
num = hs.num;
maxn = hs.maxn;
ri.clear();
for (const auto & i : hs.ri)
ri.emplace_back(i);
pt = hs.pt;
return *this;
}
/* This calls |RadicalInverse::increase| for all radical inverses and /* This calls |RadicalInverse::increase| for all radical inverses and
calls |eval|. */ calls |eval|. */
@ -153,113 +142,48 @@ HaltonSequence::eval()
void void
HaltonSequence::print() const HaltonSequence::print() const
{ {
auto ff = std::cout.flags();
for (const auto & i : ri) for (const auto & i : ri)
i.print(); i.print();
printf("point=[ "); std::cout << "point=[ "
<< std::fixed << std::setprecision(6);
for (unsigned int i = 0; i < ri.size(); i++) for (unsigned int i = 0; i < ri.size(); i++)
printf("%7.6f ", pt[i]); std::cout << std::setw(7) << pt[i] << ' ';
printf("]\n"); std::cout << ']' << std::endl;
} std::cout.flags(ff);
qmcpit::qmcpit()
{
} }
qmcpit::qmcpit(const QMCSpecification &s, int n) qmcpit::qmcpit(const QMCSpecification &s, int n)
: spec(&s), halton(new HaltonSequence(n, s.level(), s.dimen(), s.getPerScheme())), : spec(s), halton{n, s.level(), s.dimen(), s.getPerScheme()},
sig(new ParameterSignal(s.dimen())) sig{s.dimen()}
{ {
} }
qmcpit::qmcpit(const qmcpit &qpit)
: spec(qpit.spec), halton(nullptr), sig(nullptr)
{
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 bool
qmcpit::operator==(const qmcpit &qpit) const qmcpit::operator==(const qmcpit &qpit) const
{ {
return (spec == qpit.spec) return &spec == &qpit.spec && halton.getNum() == qpit.halton.getNum();
&& ((halton == nullptr && qpit.halton == nullptr)
|| (halton != nullptr && qpit.halton != nullptr && 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 = nullptr;
return *this;
} }
qmcpit & qmcpit &
qmcpit::operator++() qmcpit::operator++()
{ {
// todo: raise if |halton == null || qmcq == NULL| // todo: raise if |halton == null || qmcq == NULL|
halton->increase(); halton.increase();
return *this; return *this;
} }
double double
qmcpit::weight() const qmcpit::weight() const
{ {
return 1.0/spec->level(); return 1.0/spec.level();
}
qmcnpit::qmcnpit()
: qmcpit()
{
} }
qmcnpit::qmcnpit(const QMCSpecification &s, int n) qmcnpit::qmcnpit(const QMCSpecification &s, int n)
: qmcpit(s, n), pnt(new Vector(s.dimen())) : qmcpit(s, n), pnt{s.dimen()}
{ {
} }
qmcnpit::qmcnpit(const qmcnpit &qpit)
: qmcpit(qpit), pnt(nullptr)
{
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 = nullptr;
return *this;
}
/* Here we inccrease a point in Halton sequence ant then store images /* Here we inccrease a point in Halton sequence ant then store images
of the points in |NormalICDF| function. */ of the points in |NormalICDF| function. */
@ -267,8 +191,8 @@ qmcnpit &
qmcnpit::operator++() qmcnpit::operator++()
{ {
qmcpit::operator++(); qmcpit::operator++();
for (int i = 0; i < halton->point().length(); i++) for (int i = 0; i < halton.point().length(); i++)
(*pnt)[i] = NormalICDF::get(halton->point()[i]); pnt[i] = NormalICDF::get(halton.point()[i]);
return *this; return *this;
} }

73
dynare++/integ/cc/quasi_mcarlo.hh
View File

@ -39,11 +39,9 @@
class PermutationScheme class PermutationScheme
{ {
public: public:
PermutationScheme() PermutationScheme() = default;
= default; virtual ~PermutationScheme() = default;
virtual ~PermutationScheme() virtual int permute(int i, int base, int c) const = 0;
= default;
virtual int permute(int i, int base, int c) const = 0;
}; };
/* This class represents an integer number |num| as /* This class represents an integer number |num| as
@ -65,12 +63,8 @@ class RadicalInverse
IntSequence coeff; IntSequence coeff;
public: public:
RadicalInverse(int n, int b, int mxn); RadicalInverse(int n, int b, int mxn);
RadicalInverse(const RadicalInverse &ri) RadicalInverse(const RadicalInverse &ri) = default;
RadicalInverse &operator=(const RadicalInverse &radi) = default;
= default;
RadicalInverse &
operator=(const RadicalInverse &radi)
= default;
double eval(const PermutationScheme &p) const; double eval(const PermutationScheme &p) const;
void increase(); void increase();
void print() const; void print() const;
@ -85,8 +79,7 @@ public:
class HaltonSequence class HaltonSequence
{ {
private: private:
static int primes[]; static std::array<int, 170> primes;
static int num_primes;
protected: protected:
int num; int num;
int maxn; int maxn;
@ -95,10 +88,8 @@ protected:
Vector pt; Vector pt;
public: public:
HaltonSequence(int n, int mxn, int dim, const PermutationScheme &p); HaltonSequence(int n, int mxn, int dim, const PermutationScheme &p);
HaltonSequence(const HaltonSequence &hs) HaltonSequence(const HaltonSequence &hs) = default;
HaltonSequence &operator=(const HaltonSequence &hs) = delete;
= default;
const HaltonSequence &operator=(const HaltonSequence &hs);
void increase(); void increase();
const Vector & const Vector &
point() const point() const
@ -131,8 +122,7 @@ public:
: dim(d), lev(l), per_scheme(p) : dim(d), lev(l), per_scheme(p)
{ {
} }
virtual ~QMCSpecification() virtual ~QMCSpecification() = default;
= default;
int int
dimen() const dimen() const
{ {
@ -153,44 +143,41 @@ public:
/* This is an iterator for quasi Monte Carlo over a cube /* This is an iterator for quasi Monte Carlo over a cube
|QMCarloCubeQuadrature|. The iterator maintains |HaltonSequence| of |QMCarloCubeQuadrature|. The iterator maintains |HaltonSequence| of
the same dimension as given by the specification. An iterator can be the same dimension as given by the specification. An iterator can be
constructed from a given number |n|, or by a copy constructor. For constructed from a given number |n|, or by a copy constructor. */
technical reasons, there is also an empty constructor; for that
reason, every member is a pointer. */
class qmcpit class qmcpit
{ {
protected: protected:
const QMCSpecification *spec{nullptr}; const QMCSpecification &spec;
HaltonSequence *halton{nullptr}; HaltonSequence halton;
ParameterSignal *sig{nullptr}; ParameterSignal sig;
public: public:
qmcpit();
qmcpit(const QMCSpecification &s, int n); qmcpit(const QMCSpecification &s, int n);
qmcpit(const qmcpit &qpit); qmcpit(const qmcpit &qpit) = default;
~qmcpit(); virtual ~qmcpit() = default;
bool operator==(const qmcpit &qpit) const; bool operator==(const qmcpit &qpit) const;
bool bool
operator!=(const qmcpit &qpit) const operator!=(const qmcpit &qpit) const
{ {
return !operator==(qpit); return !operator==(qpit);
} }
const qmcpit &operator=(const qmcpit &qpit); qmcpit &operator=(const qmcpit &qpit) = delete;
qmcpit &operator++(); qmcpit &operator++();
const ParameterSignal & const ParameterSignal &
signal() const signal() const
{ {
return *sig; return sig;
} }
const Vector & const Vector &
point() const point() const
{ {
return halton->point(); return halton.point();
} }
double weight() const; double weight() const;
void void
print() const print() const
{ {
halton->print(); halton.print();
} }
}; };
@ -206,8 +193,7 @@ public:
: QuadratureImpl<qmcpit>(d), QMCSpecification(d, l, p) : QuadratureImpl<qmcpit>(d), QMCSpecification(d, l, p)
{ {
} }
~QMCarloCubeQuadrature() ~QMCarloCubeQuadrature() override = default;
override = default;
int int
numEvals(int l) const override numEvals(int l) const override
{ {
@ -233,12 +219,11 @@ protected:
class qmcnpit : public qmcpit class qmcnpit : public qmcpit
{ {
protected: protected:
Vector *pnt{nullptr}; Vector pnt;
public: public:
qmcnpit();
qmcnpit(const QMCSpecification &spec, int n); qmcnpit(const QMCSpecification &spec, int n);
qmcnpit(const qmcnpit &qpit); qmcnpit(const qmcnpit &qpit) = default;
~qmcnpit(); ~qmcnpit() override = default;
bool bool
operator==(const qmcnpit &qpit) const operator==(const qmcnpit &qpit) const
{ {
@ -249,22 +234,23 @@ public:
{ {
return !operator==(qpit); return !operator==(qpit);
} }
const qmcnpit &operator=(const qmcnpit &qpit); qmcnpit &operator=(const qmcnpit &qpit) = delete;
qmcnpit &operator++(); qmcnpit &operator++();
const ParameterSignal & const ParameterSignal &
signal() const signal() const
{ {
return *sig; return sig;
} }
const Vector & const Vector &
point() const point() const
{ {
return *pnt; return pnt;
} }
void void
print() const print() const
{ {
halton->print(); pnt->print(); halton.print();
pnt.print();
} }
}; };
@ -280,8 +266,7 @@ public:
: QuadratureImpl<qmcnpit>(d), QMCSpecification(d, l, p) : QuadratureImpl<qmcnpit>(d), QMCSpecification(d, l, p)
{ {
} }
~QMCarloNormalQuadrature() ~QMCarloNormalQuadrature() override = default;
override = default;
int int
numEvals(int l) const override numEvals(int l) const override
{ {

135
dynare++/integ/cc/smolyak.cc
View File

@ -3,91 +3,27 @@
#include "smolyak.hh" #include "smolyak.hh"
#include "symmetry.hh" #include "symmetry.hh"
smolpit::smolpit() #include <iostream>
#include <iomanip>
{
}
/* This constructs a beginning of |isum| summand in |smolq|. We must be /* 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 careful here, since |isum| can be past-the-end, so no reference to
vectors in |smolq| by |isum| must be done in this case. */ vectors in |smolq| by |isum| must be done in this case. */
smolpit::smolpit(const SmolyakQuadrature &q, unsigned int isum) smolpit::smolpit(const SmolyakQuadrature &q, unsigned int isum)
: smolq(&q), isummand(isum), jseq(new IntSequence(q.dimen(), 0)), : smolq(q), isummand(isum),
sig(new ParameterSignal(q.dimen())), p(new Vector(q.dimen())) jseq{q.dimen(), 0},
sig{q.dimen()},
p{q.dimen()}
{ {
if (isummand < q.numSummands()) if (isummand < q.numSummands())
{ setPointAndWeight();
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 = nullptr;
if (spit.sig)
sig = new ParameterSignal(*(spit.sig));
else
sig = nullptr;
if (spit.p)
p = new Vector(*(spit.p));
else
p = nullptr;
}
smolpit::~smolpit()
{
if (jseq)
delete jseq;
if (sig)
delete sig;
if (p)
delete p;
} }
bool bool
smolpit::operator==(const smolpit &spit) const smolpit::operator==(const smolpit &spit) const
{ {
bool ret = true; return &smolq == &spit.smolq && isummand == spit.isummand && jseq == spit.jseq;
ret = ret & smolq == spit.smolq;
ret = ret & isummand == spit.isummand;
ret = ret & ((jseq == nullptr && spit.jseq == nullptr)
|| (jseq != nullptr && spit.jseq != nullptr && *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 = nullptr;
if (spit.sig)
sig = new ParameterSignal(*(spit.sig));
else
sig = nullptr;
if (spit.p)
p = new Vector(*(spit.p));
else
p = nullptr;
return *this;
} }
/* We first try to increase index within the current summand. If we are /* We first try to increase index within the current summand. If we are
@ -98,22 +34,22 @@ smolpit &
smolpit::operator++() smolpit::operator++()
{ {
// todo: throw if |smolq==NULL| or |jseq==NULL| or |sig==NULL| // todo: throw if |smolq==NULL| or |jseq==NULL| or |sig==NULL|
const IntSequence &levpts = smolq->levpoints[isummand]; const IntSequence &levpts = smolq.levpoints[isummand];
int i = smolq->dimen()-1; int i = smolq.dimen()-1;
(*jseq)[i]++; jseq[i]++;
while (i >= 0 && (*jseq)[i] == levpts[i]) while (i >= 0 && jseq[i] == levpts[i])
{ {
(*jseq)[i] = 0; jseq[i] = 0;
i--; i--;
if (i >= 0) if (i >= 0)
(*jseq)[i]++; jseq[i]++;
} }
sig->signalAfter(std::max(i, 0)); sig.signalAfter(std::max(i, 0));
if (i < 0) if (i < 0)
isummand++; isummand++;
if (isummand < smolq->numSummands()) if (isummand < smolq.numSummands())
setPointAndWeight(); setPointAndWeight();
return *this; return *this;
@ -126,18 +62,18 @@ void
smolpit::setPointAndWeight() smolpit::setPointAndWeight()
{ {
// todo: raise if |smolq==NULL| or |jseq==NULL| or |sig==NULL| or // todo: raise if |smolq==NULL| or |jseq==NULL| or |sig==NULL| or
// |p==NULL| or |isummand>=smolq->numSummands()| // |p==NULL| or |isummand>=smolq.numSummands()|
int l = smolq->level; int l = smolq.level;
int d = smolq->dimen(); int d = smolq.dimen();
int sumk = (smolq->levels[isummand]).sum(); int sumk = (smolq.levels[isummand]).sum();
int m1exp = l + d - sumk - 1; int m1exp = l + d - sumk - 1;
w = (2*(m1exp/2) == m1exp) ? 1.0 : -1.0; w = (2*(m1exp/2) == m1exp) ? 1.0 : -1.0;
w *= smolq->psc.noverk(d-1, sumk-l); w *= smolq.psc.noverk(d-1, sumk-l);
for (int i = 0; i < d; i++) for (int i = 0; i < d; i++)
{ {
int ki = (smolq->levels[isummand])[i]; int ki = (smolq.levels[isummand])[i];
(*p)[i] = (smolq->uquad).point(ki, (*jseq)[i]); p[i] = (smolq.uquad).point(ki, jseq[i]);
w *= (smolq->uquad).weight(ki, (*jseq)[i]); w *= (smolq.uquad).weight(ki, jseq[i]);
} }
} }
@ -145,16 +81,19 @@ smolpit::setPointAndWeight()
void void
smolpit::print() const smolpit::print() const
{ {
printf("isum=%-3d: [", isummand); auto ff = std::cout.flags();
for (int i = 0; i < smolq->dimen(); i++) std::cout << "isum=" << std::left << std::setw(3) << isummand << std::right << ": [";
printf("%2d ", (smolq->levels[isummand])[i]); for (int i = 0; i < smolq.dimen(); i++)
printf("] j=["); std::cout << std::setw(2) << (smolq.levels[isummand])[i] << ' ';
for (int i = 0; i < smolq->dimen(); i++) std::cout << "] j=[";
printf("%2d ", (*jseq)[i]); for (int i = 0; i < smolq.dimen(); i++)
printf("] %+4.3f*(", w); std::cout << std::setw(2) << jseq[i] << ' ';
for (int i = 0; i < smolq->dimen()-1; i++) std::cout << std::showpos << std::fixed << std::setprecision(3)
printf("%+4.3f ", (*p)[i]); << "] " << std::setw(4) << w << "*(";
printf("%+4.3f)\n", (*p)[smolq->dimen()-1]); for (int i = 0; i < smolq.dimen()-1; i++)
std::cout << std::setw(4) << p[i] << ' ';
std::cout << std::setw(4) << p[smolq.dimen()-1] << ')' << std::endl;
std::cout.flags(ff);
} }
/* Here is the constructor of |SmolyakQuadrature|. We have to setup /* Here is the constructor of |SmolyakQuadrature|. We have to setup

22
dynare++/integ/cc/smolyak.hh
View File

@ -40,34 +40,33 @@ class SmolyakQuadrature;
class smolpit class smolpit
{ {
protected: protected:
const SmolyakQuadrature *smolq{nullptr}; const SmolyakQuadrature &smolq;
unsigned int isummand{0}; unsigned int isummand{0};
IntSequence *jseq{nullptr}; IntSequence jseq;
ParameterSignal *sig{nullptr}; ParameterSignal sig;
Vector *p{nullptr}; Vector p;
double w; double w;
public: public:
smolpit();
smolpit(const SmolyakQuadrature &q, unsigned int isum); smolpit(const SmolyakQuadrature &q, unsigned int isum);
smolpit(const smolpit &spit); smolpit(const smolpit &spit) = default;
~smolpit(); ~smolpit() = default;
bool operator==(const smolpit &spit) const; bool operator==(const smolpit &spit) const;
bool bool
operator!=(const smolpit &spit) const operator!=(const smolpit &spit) const
{ {
return !operator==(spit); return !operator==(spit);
} }
const smolpit &operator=(const smolpit &spit); smolpit &operator=(const smolpit &spit) = delete;
smolpit &operator++(); smolpit &operator++();
const ParameterSignal & const ParameterSignal &
signal() const signal() const
{ {
return *sig; return sig;
} }
const Vector & const Vector &
point() const point() const
{ {
return *p; return p;
} }
double double
weight() const weight() const
@ -108,8 +107,7 @@ class SmolyakQuadrature : public QuadratureImpl<smolpit>
PascalTriangle psc; PascalTriangle psc;
public: public:
SmolyakQuadrature(int d, int l, const OneDQuadrature &uq); SmolyakQuadrature(int d, int l, const OneDQuadrature &uq);
~SmolyakQuadrature() ~SmolyakQuadrature() override = default;
override = default;
int numEvals(int level) const override; int numEvals(int level) const override;
void designLevelForEvals(int max_eval, int &lev, int &evals) const; void designLevelForEvals(int max_eval, int &lev, int &evals) const;
protected: protected:

48
dynare++/integ/cc/vector_function.cc
View File

@ -5,24 +5,14 @@
#include <dynlapack.h> #include <dynlapack.h>
#include <cmath> #include <cmath>
#include <cstring>
#include <algorithm> #include <algorithm>
/* Just an easy constructor of sequence of booleans defaulting to /* Just an easy constructor of sequence of booleans defaulting to
change everywhere. */ change everywhere. */
ParameterSignal::ParameterSignal(int n) ParameterSignal::ParameterSignal(int n)
: data(new bool[n]), num(n) : data(n, true)
{ {
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| /* This sets |false| (no change) before a given parameter, and |true|
@ -31,47 +21,44 @@ ParameterSignal::ParameterSignal(const ParameterSignal &sig)
void void
ParameterSignal::signalAfter(int l) ParameterSignal::signalAfter(int l)
{ {
for (int i = 0; i < std::min(l, num); i++) for (size_t i = 0; i < std::min((size_t) l, data.size()); i++)
data[i] = false; data[i] = false;
for (int i = l; i < num; i++) for (size_t i = l; i < data.size(); i++)
data[i] = true; data[i] = true;
} }
/* This constructs a function set hardcopying also the first. */ /* This constructs a function set hardcopying also the first. */
VectorFunctionSet::VectorFunctionSet(const VectorFunction &f, int n) VectorFunctionSet::VectorFunctionSet(const VectorFunction &f, int n)
: funcs(n), first_shallow(false) : funcs(n)
{ {
for (int i = 0; i < n; i++) for (int i = 0; i < n; i++)
funcs[i] = f.clone(); {
func_copies.push_back(f.clone());
funcs[i] = func_copies.back().get();
}
} }
/* This constructs a function set with shallow copy in the first and /* This constructs a function set with shallow copy in the first and
hard copies in others. */ hard copies in others. */
VectorFunctionSet::VectorFunctionSet(VectorFunction &f, int n) VectorFunctionSet::VectorFunctionSet(VectorFunction &f, int n)
: funcs(n), first_shallow(true) : funcs(n)
{ {
if (n > 0) if (n > 0)
funcs[0] = &f; funcs[0] = &f;
for (int i = 1; i < n; i++) for (int i = 1; i < n; i++)
funcs[i] = f.clone(); {
func_copies.push_back(f.clone());
funcs[i] = func_copies.back().get();
}
} }
/* 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 /* Here we construct the object from the given function $f$ and given
variance-covariance matrix $\Sigma=$|vcov|. The matrix $A$ is variance-covariance matrix $\Sigma=$|vcov|. The matrix $A$ is
calculated as lower triangular and yields $\Sigma=AA^T$. */ calculated as lower triangular and yields $\Sigma=AA^T$. */
GaussConverterFunction::GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov) GaussConverterFunction::GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov)
: VectorFunction(f), func(&f), delete_flag(false), A(vcov.numRows(), vcov.numRows()), : VectorFunction(f), func(&f), A(vcov.numRows(), vcov.numRows()),
multiplier(calcMultiplier()) multiplier(calcMultiplier())
{ {
// todo: raise if |A.numRows() != indim()| // todo: raise if |A.numRows() != indim()|
@ -80,9 +67,8 @@ GaussConverterFunction::GaussConverterFunction(VectorFunction &f, const GeneralM
/* Here we construct the object in the same way, however we mark the /* Here we construct the object in the same way, however we mark the
function as to be deleted. */ function as to be deleted. */
GaussConverterFunction::GaussConverterFunction(std::unique_ptr<VectorFunction> f, const GeneralMatrix &vcov)
GaussConverterFunction::GaussConverterFunction(VectorFunction *f, const GeneralMatrix &vcov) : VectorFunction(*f), func_storage{move(f)}, func{func_storage.get()}, A(vcov.numRows(), vcov.numRows()),
: VectorFunction(*f), func(f), delete_flag(true), A(vcov.numRows(), vcov.numRows()),
multiplier(calcMultiplier()) multiplier(calcMultiplier())
{ {
// todo: raise if |A.numRows() != indim()| // todo: raise if |A.numRows() != indim()|
@ -90,7 +76,7 @@ GaussConverterFunction::GaussConverterFunction(VectorFunction *f, const GeneralM
} }
GaussConverterFunction::GaussConverterFunction(const GaussConverterFunction &f) GaussConverterFunction::GaussConverterFunction(const GaussConverterFunction &f)
: VectorFunction(f), func(f.func->clone()), delete_flag(true), A(f.A), : VectorFunction(f), func_storage{f.func->clone()}, func{func_storage.get()}, A(f.A),
multiplier(f.multiplier) multiplier(f.multiplier)
{ {
} }

54
dynare++/integ/cc/vector_function.hh
View File

@ -20,6 +20,7 @@
#include "GeneralMatrix.hh" #include "GeneralMatrix.hh"
#include <vector> #include <vector>
#include <memory>
/* This is a simple class representing a vector of booleans. The items /* This is a simple class representing a vector of booleans. The items
night be retrieved or changed, or can be set |true| after some night be retrieved or changed, or can be set |true| after some
@ -31,22 +32,18 @@
class ParameterSignal class ParameterSignal
{ {
protected: protected:
bool *data; std::vector<bool> data;
int num;
public: public:
ParameterSignal(int n); ParameterSignal(int n);
ParameterSignal(const ParameterSignal &sig); ParameterSignal(const ParameterSignal &sig) = default;
~ParameterSignal() ~ParameterSignal() = default;
{
delete [] data;
}
void signalAfter(int l); void signalAfter(int l);
const bool & bool
operator[](int i) const operator[](int i) const
{ {
return data[i]; return data[i];
} }
bool & std::vector<bool>::reference
operator[](int i) operator[](int i)
{ {
return data[i]; return data[i];
@ -71,12 +68,9 @@ public:
: in_dim(idim), out_dim(odim) : in_dim(idim), out_dim(odim)
{ {
} }
VectorFunction(const VectorFunction &func) VectorFunction(const VectorFunction &func) = default;
virtual ~VectorFunction() = default;
= default; virtual std::unique_ptr<VectorFunction> clone() const = 0;
virtual ~VectorFunction()
= default;
virtual VectorFunction *clone() const = 0;
virtual void eval(const Vector &point, const ParameterSignal &sig, Vector &out) = 0; virtual void eval(const Vector &point, const ParameterSignal &sig, Vector &out) = 0;
int int
indim() const indim() const
@ -100,13 +94,15 @@ public:
class VectorFunctionSet class VectorFunctionSet
{ {
private:
// Stores the hard copies made by the class
std::vector<std::unique_ptr<VectorFunction>> func_copies;
protected: protected:
std::vector<VectorFunction *> funcs; std::vector<VectorFunction *> funcs;
bool first_shallow;
public: public:
VectorFunctionSet(const VectorFunction &f, int n); VectorFunctionSet(const VectorFunction &f, int n);
VectorFunctionSet(VectorFunction &f, int n); VectorFunctionSet(VectorFunction &f, int n);
~VectorFunctionSet(); ~VectorFunctionSet() = default;
VectorFunction & VectorFunction &
getFunc(int i) getFunc(int i)
{ {
@ -136,32 +132,28 @@ public:
normally distributed inputs. normally distributed inputs.
The class maintains a pointer to the function $f$. When the object is 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 constructed by the first constructor, the $f$ is assumed to be owned by the
object of this class is copied, then $f$ is copied and we need to caller. If the object of this class is copied, then $f$ is copied and hence
remember to destroy it in the desctructor; hence |delete_flag|. The stored in a std::unique_ptr. The second constructor takes a smart pointer to
second constructor takes a pointer to the function and differs from the function and in that case the class takes ownership of $f$. */
the first only by setting |delete_flag| to |true|. */
class GaussConverterFunction : public VectorFunction class GaussConverterFunction : public VectorFunction
{ {
private:
std::unique_ptr<VectorFunction> func_storage;
protected: protected:
VectorFunction *func; VectorFunction *func;
bool delete_flag;
GeneralMatrix A; GeneralMatrix A;
double multiplier; double multiplier;
public: public:
GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov); GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov);
GaussConverterFunction(VectorFunction *f, const GeneralMatrix &vcov); GaussConverterFunction(std::unique_ptr<VectorFunction> f, const GeneralMatrix &vcov);
GaussConverterFunction(const GaussConverterFunction &f); GaussConverterFunction(const GaussConverterFunction &f);
~GaussConverterFunction() override ~GaussConverterFunction() override = default;
{ std::unique_ptr<VectorFunction>
if (delete_flag)
delete func;
}
VectorFunction *
clone() const override clone() const override
{ {
return new GaussConverterFunction(*this); return std::make_unique<GaussConverterFunction>(*this);
} }
void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override; void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
private: private:

107
dynare++/integ/src/quadrature-points.cc
View File

@ -11,19 +11,21 @@
#include "integ/cc/product.hh" #include "integ/cc/product.hh"
#include <getopt.h> #include <getopt.h>
#include <cstdio>
#include <cmath> #include <cmath>
#include <cstdlib>
#include <iostream>
#include <fstream> #include <fstream>
#include <sstream> #include <sstream>
#include <memory>
#include <string>
struct QuadParams struct QuadParams
{ {
const char *outname; string outname;
const char *vcovname; string vcovname;
int max_level; int max_level{3};
double discard_weight; double discard_weight{0.0};
QuadParams(int argc, char **argv); QuadParams(int argc, char **argv);
void check_consistency() const; void check_consistency() const;
private: private:
@ -31,12 +33,12 @@ private:
}; };
QuadParams::QuadParams(int argc, char **argv) QuadParams::QuadParams(int argc, char **argv)
: outname(nullptr), vcovname(nullptr), max_level(3), discard_weight(0.0)
{ {
if (argc == 1) if (argc == 1)
{ {
// print the help and exit // print the help and exit
exit(1); std::cerr << "Usage: " << argv[0] << " [--max-level INTEGER] [--discard-weight FLOAT] [--vcov FILENAME] OUTPUT_FILENAME" << std::endl;
std::exit(EXIT_FAILURE);
} }
outname = argv[argc-1]; outname = argv[argc-1];
@ -56,12 +58,24 @@ QuadParams::QuadParams(int argc, char **argv)
switch (ret) switch (ret)
{ {
case opt_max_level: case opt_max_level:
if (1 != sscanf(optarg, "%d", &max_level)) try
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg); {
max_level = std::stoi(string{optarg});
}
catch (const std::invalid_argument &e)
{
std::cerr << "Couldn't parse integer " << optarg << ", ignored" << std::endl;
}
break; break;
case opt_discard_weight: case opt_discard_weight:
if (1 != sscanf(optarg, "%lf", &discard_weight)) try
fprintf(stderr, "Couldn't parse float %s, ignored\n", optarg); {
discard_weight = std::stod(string{optarg});
}
catch (const std::invalid_argument &e)
{
std::cerr << "Couldn't parse float " << optarg << ", ignored" << std::endl;
}
break; break;
case opt_vcov: case opt_vcov:
vcovname = optarg; vcovname = optarg;
@ -75,41 +89,30 @@ QuadParams::QuadParams(int argc, char **argv)
void void
QuadParams::check_consistency() const QuadParams::check_consistency() const
{ {
if (outname == nullptr) if (outname.empty())
{ {
fprintf(stderr, "Error: output name not set\n"); std::cerr << "Error: output name not set" << std::endl;
exit(1); std::exit(EXIT_FAILURE);
} }
if (vcovname == nullptr) if (vcovname.empty())
{ {
fprintf(stderr, "Error: vcov file name not set\n"); std::cerr << "Error: vcov file name not set" << std::endl;
exit(1); std::exit(EXIT_FAILURE);
} }
} }
/** 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 int
main(int argc, char **argv) main(int argc, char **argv)
{ {
QuadParams params(argc, argv); QuadParams params(argc, argv);
// open output file for writing // open output file for writing
FILE *fout; std::ofstream fout{params.outname, std::ios::out | std::ios::trunc};
if (nullptr == (fout = fopen(params.outname, "w"))) if (fout.fail())
{ {
fprintf(stderr, "Could not open %s for writing\n", params.outname); std::cerr << "Could not open " << params.outname << " for writing" << std::endl;
exit(1); std::exit(EXIT_FAILURE);
} }
try try
@ -142,23 +145,21 @@ main(int argc, char **argv)
int level = params.max_level; int level = params.max_level;
SmolyakQuadrature sq(vcov.numRows(), level, ghq); SmolyakQuadrature sq(vcov.numRows(), level, ghq);
printf("Dimension: %d\n", vcov.numRows()); std::cout << "Dimension: " << vcov.numRows() << std::endl
printf("Maximum level: %d\n", level); << "Maximum level: " << level << std::endl
printf("Total number of nodes: %d\n", sq.numEvals(level)); << "Total number of nodes: " << sq.numEvals(level) << std::endl;
// put the points to the vector // put the points to the vector
std::vector<Vector *> points; std::vector<std::unique_ptr<Vector>> points;
for (smolpit qit = sq.start(level); qit != sq.end(level); ++qit) for (smolpit qit = sq.start(level); qit != sq.end(level); ++qit)
points.push_back(new Vector((const Vector &) qit.point())); points.push_back(std::make_unique<Vector>((const Vector &) qit.point()));
// sort and uniq // sort and uniq
OrderVec ordvec; std::sort(points.begin(), points.end(), [](auto &a, auto &b) { return a.get() < b.get(); });
std::sort(points.begin(), points.end(), ordvec);
auto new_end = std::unique(points.begin(), points.end()); auto new_end = std::unique(points.begin(), points.end());
for (auto it = new_end; it != points.end(); ++it)
delete *it;
points.erase(new_end, points.end()); points.erase(new_end, points.end());
printf("Duplicit nodes removed: %lu\n", (unsigned long) (sq.numEvals(level)-points.size())); std::cout << "Duplicit nodes removed: " << (unsigned long) (sq.numEvals(level)-points.size())
<< std::endl;
// calculate weights and mass // calculate weights and mass
double mass = 0.0; double mass = 0.0;
@ -175,41 +176,41 @@ main(int argc, char **argv)
if (weight/mass < params.discard_weight) if (weight/mass < params.discard_weight)
discard_mass += weight; discard_mass += weight;
printf("Total mass discarded: %f\n", discard_mass/mass); std::cout << "Total mass discarded: " << std::fixed << discard_mass/mass << std::endl;
// dump the results // dump the results
int npoints = 0; int npoints = 0;
double upscale_weight = 1/(mass-discard_mass); double upscale_weight = 1/(mass-discard_mass);
Vector x(vcov.numRows()); Vector x(vcov.numRows());
fout << std::setprecision(16);
for (int i = 0; i < (int) weights.size(); i++) for (int i = 0; i < (int) weights.size(); i++)
if (weights[i]/mass >= params.discard_weight) if (weights[i]/mass >= params.discard_weight)
{ {
// print the upscaled weight // print the upscaled weight
fprintf(fout, "%20.16g", upscale_weight*weights[i]); fout << std::setw(20) << upscale_weight*weights[i];
// multiply point with the factor A and sqrt(2) // multiply point with the factor A and sqrt(2)
A.multVec(0.0, x, std::sqrt(2.), *(points[i])); A.multVec(0.0, x, std::sqrt(2.), *(points[i]));
// print the coordinates // print the coordinates
for (int j = 0; j < x.length(); j++) for (int j = 0; j < x.length(); j++)
fprintf(fout, " %20.16g", x[j]); fout << ' ' << std::setw(20) << x[j];
fprintf(fout, "\n"); fout << std::endl;
npoints++; npoints++;
} }
printf("Final number of points: %d\n", npoints); std::cout << "Final number of points: " << npoints << std::endl;
fclose(fout);
fout.close();
} }
catch (const SylvException &e) catch (const SylvException &e)
{ {
e.printMessage(); e.printMessage();
return 1; return EXIT_FAILURE;
} }
catch (const ogu::Exception &e) catch (const ogu::Exception &e)
{ {
e.print(); e.print();
return 1; return EXIT_FAILURE;
} }
return 0; return EXIT_SUCCESS;
} }

185
dynare++/integ/testing/tests.cc
View File

@ -14,10 +14,14 @@
#include "product.hh" #include "product.hh"
#include "quasi_mcarlo.hh" #include "quasi_mcarlo.hh"
#include <cstdio> #include <iomanip>
#include <cstring> #include <chrono>
#include <sys/time.h>
#include <cmath> #include <cmath>
#include <iostream>
#include <utility>
#include <array>
#include <memory>
#include <cstdlib>
const int num_threads = 2; // does nothing if DEBUG defined const int num_threads = 2; // does nothing if DEBUG defined
@ -33,13 +37,12 @@ public:
D(inD), k(kk) D(inD), k(kk)
{ {
} }
MomentFunction(const MomentFunction &func) MomentFunction(const MomentFunction &func) = default;
= default; std::unique_ptr<VectorFunction>
VectorFunction *
clone() const override clone() const override
{ {
return new MomentFunction(*this); return std::make_unique<MomentFunction>(*this);
} }
void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override; void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
}; };
@ -49,8 +52,8 @@ MomentFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &ou
{ {
if (point.length() != indim() || out.length() != outdim()) if (point.length() != indim() || out.length() != outdim())
{ {
printf("Wrong length of vectors in MomentFunction::eval\n"); std::cerr << "Wrong length of vectors in MomentFunction::eval" << std::endl;
exit(1); std::exit(EXIT_FAILURE);
} }
Vector y(point); Vector y(point);
y.zeros(); y.zeros();
@ -68,13 +71,12 @@ public:
: VectorFunction(nvar, UFSTensor::calcMaxOffset(nvar, kk)), k(kk) : VectorFunction(nvar, UFSTensor::calcMaxOffset(nvar, kk)), k(kk)
{ {
} }
TensorPower(const TensorPower &func) TensorPower(const TensorPower &func) = default;
= default; std::unique_ptr<VectorFunction>
VectorFunction *
clone() const override clone() const override
{ {
return new TensorPower(*this); return std::make_unique<TensorPower>(*this);
} }
void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override; void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
}; };
@ -84,8 +86,8 @@ TensorPower::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{ {
if (point.length() != indim() || out.length() != outdim()) if (point.length() != indim() || out.length() != outdim())
{ {
printf("Wrong length of vectors in TensorPower::eval\n"); std::cerr << "Wrong length of vectors in TensorPower::eval" << std::endl;
exit(1); std::exit(EXIT_FAILURE);
} }
URSingleTensor ypow(point, k); URSingleTensor ypow(point, k);
out.zeros(); out.zeros();
@ -107,10 +109,10 @@ public:
: VectorFunction(f.indim(), f.outdim()), dim(f.dim) : VectorFunction(f.indim(), f.outdim()), dim(f.dim)
{ {
} }
VectorFunction * std::unique_ptr<VectorFunction>
clone() const override clone() const override
{ {
return new Function1(*this); return std::make_unique<Function1>(*this);
} }
void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override; void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
}; };
@ -120,8 +122,8 @@ Function1::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{ {
if (point.length() != dim || out.length() != 1) if (point.length() != dim || out.length() != 1)
{ {
printf("Wrong length of vectors in Function1::eval\n"); std::cerr << "Wrong length of vectors in Function1::eval" << std::endl;
exit(1); std::exit(EXIT_FAILURE);
} }
double r = 1; double r = 1;
for (int i = 0; i < dim; i++) for (int i = 0; i < dim; i++)
@ -140,13 +142,12 @@ public:
: Function1(d) : Function1(d)
{ {
} }
Function1Trans(const Function1Trans &func) Function1Trans(const Function1Trans &func) = default;
= default; std::unique_ptr<VectorFunction>
VectorFunction *
clone() const override clone() const override
{ {
return new Function1Trans(*this); return std::make_unique<Function1Trans>(*this);
} }
void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override; void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
}; };
@ -166,22 +167,21 @@ Function1Trans::eval(const Vector &point, const ParameterSignal &sig, Vector &ou
// with time information // with time information
class WallTimer class WallTimer
{ {
char mes[100]; string mes;
struct timeval start; std::chrono::time_point<std::chrono::high_resolution_clock> start;
bool new_line; bool new_line;
public: public:
WallTimer(const char *m, bool nl = true) WallTimer(string m, bool nl = true)
: mes{m}, start{std::chrono::high_resolution_clock::now()}, new_line{nl}
{ {
strcpy(mes, m); new_line = nl; gettimeofday(&start, nullptr);
} }
~WallTimer() ~WallTimer()
{ {
struct timeval end; auto end = std::chrono::high_resolution_clock::now();
gettimeofday(&end, nullptr); std::chrono::duration<double> duration = end - start;
printf("%s%8.4g", mes, std::cout << mes << std::setw(8) << std::setprecision(4) << duration.count();
end.tv_sec-start.tv_sec + (end.tv_usec-start.tv_usec)*1.0e-6);
if (new_line) if (new_line)
printf("\n"); std::cout << std::endl;
} }
}; };
@ -190,22 +190,17 @@ public:
/****************************************************/ /****************************************************/
class TestRunnable class TestRunnable
{ {
char name[100];
public: public:
const string name;
int dim; // dimension of the solved problem int dim; // dimension of the solved problem
int nvar; // number of variable of the solved problem int nvar; // number of variable of the solved problem
TestRunnable(const char *n, int d, int nv) TestRunnable(string name_arg, int d, int nv)
: dim(d), nvar(nv) : name{move(name_arg)}, dim(d), nvar(nv)
{ {
strncpy(name, n, 100);
} }
virtual ~TestRunnable() = default;
bool test() const; bool test() const;
virtual bool run() const = 0; virtual bool run() const = 0;
const char *
getName() const
{
return name;
}
protected: protected:
static bool smolyak_normal_moments(const GeneralMatrix &m, int imom, int level); 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 product_normal_moments(const GeneralMatrix &m, int imom, int level);
@ -218,7 +213,7 @@ protected:
bool bool
TestRunnable::test() const TestRunnable::test() const
{ {
printf("Running test <%s>\n", name); cout << "Running test <" << name << ">" << endl;
bool passed; bool passed;
{ {
WallTimer tim("Wall clock time ", false); WallTimer tim("Wall clock time ", false);
@ -226,12 +221,12 @@ TestRunnable::test() const
} }
if (passed) if (passed)
{ {
printf("............................ passed\n\n"); std::cout << "............................ passed" << std::endl << std::endl;
return passed; return passed;
} }
else else
{ {
printf("............................ FAILED\n\n"); std::cout << "............................ FAILED" << std::endl << std::endl;
return passed; return passed;
} }
} }
@ -258,13 +253,13 @@ TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level
GaussHermite gs; GaussHermite gs;
SmolyakQuadrature quad(dim, level, gs); SmolyakQuadrature quad(dim, level, gs);
quad.integrate(func, level, num_threads, smol_out); quad.integrate(func, level, num_threads, smol_out);
printf("\tNumber of Smolyak evaluations: %d\n", quad.numEvals(level)); std::cout << "\tNumber of Smolyak evaluations: " << quad.numEvals(level) << std::endl;
} }
// check against theoretical moments // check against theoretical moments
UNormalMoments moments(imom, msq); UNormalMoments moments(imom, msq);
smol_out.add(-1.0, (moments.get(Symmetry(imom)))->getData()); smol_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
printf("\tError: %16.12g\n", smol_out.getMax()); std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << smol_out.getMax() << std::endl;
return smol_out.getMax() < 1.e-7; return smol_out.getMax() < 1.e-7;
} }
@ -287,13 +282,13 @@ TestRunnable::product_normal_moments(const GeneralMatrix &m, int imom, int level
GaussHermite gs; GaussHermite gs;
ProductQuadrature quad(dim, gs); ProductQuadrature quad(dim, gs);
quad.integrate(func, level, num_threads, prod_out); quad.integrate(func, level, num_threads, prod_out);
printf("\tNumber of product evaluations: %d\n", quad.numEvals(level)); std::cout << "\tNumber of product evaluations: " << quad.numEvals(level) << std::endl;
} }
// check against theoretical moments // check against theoretical moments
UNormalMoments moments(imom, msq); UNormalMoments moments(imom, msq);
prod_out.add(-1.0, (moments.get(Symmetry(imom)))->getData()); prod_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
printf("\tError: %16.12g\n", prod_out.getMax()); std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << prod_out.getMax() << endl;
return prod_out.getMax() < 1.e-7; return prod_out.getMax() < 1.e-7;
} }
@ -318,8 +313,8 @@ TestRunnable::qmc_normal_moments(const GeneralMatrix &m, int imom, int level)
WarnockPerScheme wps; WarnockPerScheme wps;
ReversePerScheme rps; ReversePerScheme rps;
IdentityPerScheme ips; IdentityPerScheme ips;
PermutationScheme *scheme[] = {&wps, &rps, &ips}; array<PermutationScheme *, 3> scheme = {&wps, &rps, &ips};
const char *labs[] = {"Warnock", "Reverse", "Identity"}; array<string, 3> labs = {"Warnock", "Reverse", "Identity"};
// theoretical result // theoretical result
int dim = mchol.numRows(); int dim = mchol.numRows();
@ -329,21 +324,19 @@ TestRunnable::qmc_normal_moments(const GeneralMatrix &m, int imom, int level)
// quasi monte carlo normal quadrature // quasi monte carlo normal quadrature
double max_error = 0.0; double max_error = 0.0;
Vector qmc_out(UFSTensor::calcMaxOffset(dim, imom)); Vector qmc_out(UFSTensor::calcMaxOffset(dim, imom));
for (int i = 0; i < 3; i++) for (int i = 0; i < (int) scheme.size(); i++)
{ {
{ {
char mes[100]; ostringstream mes;
sprintf(mes, "\tQMC normal quadrature time %8s: ", labs[i]); mes << "\tQMC normal quadrature time " << std::setw(8) << labs[i] << ": ";
WallTimer tim(mes); WallTimer tim(mes.str());
QMCarloNormalQuadrature quad(dim, level, *(scheme[i])); QMCarloNormalQuadrature quad(dim, level, *(scheme[i]));
quad.integrate(func, level, num_threads, qmc_out); quad.integrate(func, level, num_threads, qmc_out);
} }
qmc_out.add(-1.0, res); qmc_out.add(-1.0, res);
printf("\tError %8s: %16.12g\n", labs[i], qmc_out.getMax()); std::cout << "\tError " << std::setw(8) << labs[i] << ": " << std::setw(16) << std::setprecision(12) << qmc_out.getMax() << std::endl;
if (qmc_out.getMax() > max_error) if (qmc_out.getMax() > max_error)
{ max_error = qmc_out.getMax();
max_error = qmc_out.getMax();
}
} }
return max_error < 1.e-7; return max_error < 1.e-7;
@ -355,8 +348,8 @@ TestRunnable::smolyak_product_cube(const VectorFunction &func, const Vector &res
{ {
if (res.length() != func.outdim()) if (res.length() != func.outdim())
{ {
fprintf(stderr, "Incompatible dimensions of check value and function.\n"); std::cerr << "Incompatible dimensions of check value and function." << std::endl;
exit(1); std::exit(EXIT_FAILURE);
} }
GaussLegendre glq; GaussLegendre glq;
@ -369,8 +362,8 @@ TestRunnable::smolyak_product_cube(const VectorFunction &func, const Vector &res
quad.integrate(func, level, num_threads, out); quad.integrate(func, level, num_threads, out);
out.add(-1.0, res); out.add(-1.0, res);
smol_error = out.getMax(); smol_error = out.getMax();
printf("\tNumber of Smolyak evaluations: %d\n", quad.numEvals(level)); std::cout << "\tNumber of Smolyak evaluations: " << quad.numEvals(level) << std::endl;
printf("\tError: %16.12g\n", smol_error); std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << smol_error << std::endl;
} }
{ {
WallTimer tim("\tProduct quadrature time: "); WallTimer tim("\tProduct quadrature time: ");
@ -378,8 +371,8 @@ TestRunnable::smolyak_product_cube(const VectorFunction &func, const Vector &res
quad.integrate(func, level, num_threads, out); quad.integrate(func, level, num_threads, out);
out.add(-1.0, res); out.add(-1.0, res);
prod_error = out.getMax(); prod_error = out.getMax();
printf("\tNumber of product evaluations: %d\n", quad.numEvals(level)); std::cout << "\tNumber of product evaluations: " << quad.numEvals(level) << std::endl;
printf("\tError: %16.12g\n", prod_error); std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << prod_error << std::endl;
} }
return smol_error < tol && prod_error < tol; return smol_error < tol && prod_error < tol;
@ -397,8 +390,7 @@ TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int l
// qmc.savePoints("warnock.txt", level); // qmc.savePoints("warnock.txt", level);
qmc.integrate(func, level, num_threads, r); qmc.integrate(func, level, num_threads, r);
error1 = std::max(res - r[0], r[0] - res); error1 = std::max(res - r[0], r[0] - res);
printf("\tQuasi-Monte Carlo (Warnock scrambling) error: %16.12g\n", std::cout << "\tQuasi-Monte Carlo (Warnock scrambling) error: " << std::setw(16) << std::setprecision(12) << error1 << std::endl;
error1);
} }
double error2; double error2;
{ {
@ -408,8 +400,7 @@ TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int l
// qmc.savePoints("reverse.txt", level); // qmc.savePoints("reverse.txt", level);
qmc.integrate(func, level, num_threads, r); qmc.integrate(func, level, num_threads, r);
error2 = std::max(res - r[0], r[0] - res); error2 = std::max(res - r[0], r[0] - res);
printf("\tQuasi-Monte Carlo (reverse scrambling) error: %16.12g\n", std::cout << "\tQuasi-Monte Carlo (reverse scrambling) error: " << std::setw(16) << std::setprecision(12) << error2 << std::endl;
error2);
} }
double error3; double error3;
{ {
@ -419,8 +410,7 @@ TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int l
// qmc.savePoints("identity.txt", level); // qmc.savePoints("identity.txt", level);
qmc.integrate(func, level, num_threads, r); qmc.integrate(func, level, num_threads, r);
error3 = std::max(res - r[0], r[0] - res); error3 = std::max(res - r[0], r[0] - res);
printf("\tQuasi-Monte Carlo (no scrambling) error: %16.12g\n", std::cout << "\tQuasi-Monte Carlo (no scrambling) error: " << std::setw(16) << std::setprecision(12) << error3 << std::endl;
error3);
} }
return error1 < tol && error2 < tol && error3 < tol; return error1 < tol && error2 < tol && error3 < tol;
@ -574,61 +564,54 @@ public:
int int
main() main()
{ {
TestRunnable *all_tests[50]; std::vector<std::unique_ptr<TestRunnable>> all_tests;
// fill in vector of all tests // fill in vector of all tests
int num_tests = 0; all_tests.push_back(std::make_unique<SmolyakNormalMom1>());
all_tests[num_tests++] = new SmolyakNormalMom1(); all_tests.push_back(std::make_unique<SmolyakNormalMom2>());
all_tests[num_tests++] = new SmolyakNormalMom2(); all_tests.push_back(std::make_unique<ProductNormalMom1>());
all_tests[num_tests++] = new ProductNormalMom1(); all_tests.push_back(std::make_unique<ProductNormalMom2>());
all_tests[num_tests++] = new ProductNormalMom2(); all_tests.push_back(std::make_unique<QMCNormalMom1>());
all_tests[num_tests++] = new QMCNormalMom1(); all_tests.push_back(std::make_unique<QMCNormalMom2>());
all_tests[num_tests++] = new QMCNormalMom2();
/* /*
all_tests[num_tests++] = new F1GaussLegendre(); all_tests.push_back(make_unique<F1GaussLegendre>());
all_tests[num_tests++] = new F1QuasiMCarlo(); all_tests.push_back(make_unique<F1QuasiMCarlo>());
*/ */
// find maximum dimension and maximum nvar // find maximum dimension and maximum nvar
int dmax = 0; int dmax = 0;
int nvmax = 0; int nvmax = 0;
for (int i = 0; i < num_tests; i++) for (const auto &test : all_tests)
{ {
if (dmax < all_tests[i]->dim) if (dmax < test->dim)
dmax = all_tests[i]->dim; dmax = test->dim;
if (nvmax < all_tests[i]->nvar) if (nvmax < test->nvar)
nvmax = all_tests[i]->nvar; nvmax = test->nvar;
} }
tls.init(dmax, nvmax); // initialize library tls.init(dmax, nvmax); // initialize library
THREAD_GROUP::max_parallel_threads = num_threads; THREAD_GROUP::max_parallel_threads = num_threads;
// launch the tests // launch the tests
int success = 0; int success = 0;
for (int i = 0; i < num_tests; i++) for (const auto &test : all_tests)
{ {
try try
{ {
if (all_tests[i]->test()) if (test->test())
success++; success++;
} }
catch (const TLException &e) catch (const TLException &e)
{ {
printf("Caugth TL exception in <%s>:\n", all_tests[i]->getName()); std::cout << "Caught TL exception in <" << test->name << ">:" << std::endl;
e.print(); e.print();
} }
catch (SylvException &e) catch (SylvException &e)
{ {
printf("Caught Sylv exception in <%s>:\n", all_tests[i]->getName()); std::cout << "Caught Sylv exception in <" << test->name << ">:" << std::endl;
e.printMessage(); e.printMessage();
} }
} }
printf("There were %d tests that failed out of %d tests run.\n", std::cout << "There were " << all_tests.size() - success << " tests that failed out of "
num_tests - success, num_tests); << all_tests.size() << " tests run." << std::endl;
// destroy return EXIT_SUCCESS;
for (int i = 0; i < num_tests; i++)
{
delete all_tests[i];
}
return 0;
} }

19
dynare++/kord/global_check.hh
View File

@ -37,6 +37,8 @@
#include <matio.h> #include <matio.h>
#include <memory>
#include "vector_function.hh" #include "vector_function.hh"
#include "quadrature.hh" #include "quadrature.hh"
@ -76,10 +78,10 @@ public:
ResidFunction(const ResidFunction &rf); ResidFunction(const ResidFunction &rf);
~ResidFunction() override; ~ResidFunction() override;
VectorFunction * std::unique_ptr<VectorFunction>
clone() const override clone() const override
{ {
return new ResidFunction(*this); return std::make_unique<ResidFunction>(*this);
} }
void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override; void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
void setYU(const Vector &ys, const Vector &xx); void setYU(const Vector &ys, const Vector &xx);
@ -91,18 +93,15 @@ class GResidFunction : public GaussConverterFunction
{ {
public: public:
GResidFunction(const Approximation &app) GResidFunction(const Approximation &app)
: GaussConverterFunction(new ResidFunction(app), app.getModel().getVcov()) : GaussConverterFunction(std::make_unique<ResidFunction>(app), app.getModel().getVcov())
{ {
} }
GResidFunction(const GResidFunction &rf) GResidFunction(const GResidFunction &rf) = default;
~GResidFunction() override = default;
= default; std::unique_ptr<VectorFunction>
~GResidFunction()
override = default;
VectorFunction *
clone() const override clone() const override
{ {
return new GResidFunction(*this); return std::make_unique<GResidFunction>(*this);
} }
void void
setYU(const Vector &ys, const Vector &xx) setYU(const Vector &ys, const Vector &xx)