diff --git a/dynare++/integ/cc/product.cc b/dynare++/integ/cc/product.cc index 891af46dc..d1f368c1a 100644 --- a/dynare++/integ/cc/product.cc +++ b/dynare++/integ/cc/product.cc @@ -3,111 +3,47 @@ #include "product.hh" #include "symmetry.hh" -prodpit::prodpit() - -{ -} +#include +#include /* 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())) + : prodq(q), level(l), npoints(q.uquad.numPoints(l)), + jseq{q.dimen(), 0}, + end_flag(false), + sig{q.dimen()}, + p{q.dimen()} { if (j0 < npoints) { - (*jseq)[0] = j0; + 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 = 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; + end_flag = true; } 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 == 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; + return &prodq == &ppit.prodq && end_flag == ppit.end_flag && jseq == ppit.jseq; } 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) + int i = prodq.dimen()-1; + jseq[i]++; + while (i >= 0 && jseq[i] == npoints) { - (*jseq)[i] = 0; + jseq[i] = 0; i--; if (i >= 0) - (*jseq)[i]++; + jseq[i]++; } - sig->signalAfter(std::max(i, 0)); + sig.signalAfter(std::max(i, 0)); if (i == -1) end_flag = true; @@ -126,10 +62,10 @@ 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++) + for (int i = 0; i < prodq.dimen(); i++) { - (*p)[i] = (prodq->uquad).point(level, (*jseq)[i]); - w *= (prodq->uquad).weight(level, (*jseq)[i]); + p[i] = (prodq.uquad).point(level, jseq[i]); + w *= (prodq.uquad).weight(level, jseq[i]); } } @@ -138,13 +74,16 @@ prodpit::setPointAndWeight() 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]); + auto ff = std::cout.flags(); + std::cout << "j=["; + for (int i = 0; i < prodq.dimen(); i++) + std::cout << std::setw(2) << jseq[i]; + std::cout << std::showpos << std::fixed << std::setprecision(3) + << "] " << std::setw(4) << w << "*("; + 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) diff --git a/dynare++/integ/cc/product.hh b/dynare++/integ/cc/product.hh index a381c8973..061e1d962 100644 --- a/dynare++/integ/cc/product.hh +++ b/dynare++/integ/cc/product.hh @@ -38,36 +38,36 @@ class ProductQuadrature; class prodpit { protected: - const ProductQuadrature *prodq{nullptr}; + const ProductQuadrature &prodq; int level{0}; int npoints{0}; - IntSequence *jseq{nullptr}; + IntSequence jseq; bool end_flag{true}; - ParameterSignal *sig{nullptr}; - Vector *p{nullptr}; + ParameterSignal sig; + Vector p; double w; public: - prodpit(); + prodpit() = default; prodpit(const ProductQuadrature &q, int j0, int l); - prodpit(const prodpit &ppit); - ~prodpit(); + prodpit(const prodpit &ppit) = default; + ~prodpit() = default; bool operator==(const prodpit &ppit) const; bool operator!=(const prodpit &ppit) const { return !operator==(ppit); } - const prodpit &operator=(const prodpit &spit); + prodpit &operator=(const prodpit &spit) = delete; prodpit &operator++(); const ParameterSignal & signal() const { - return *sig; + return sig; } const Vector & point() const { - return *p; + return p; } double weight() const @@ -93,8 +93,7 @@ class ProductQuadrature : public QuadratureImpl const OneDQuadrature &uquad; public: ProductQuadrature(int d, const OneDQuadrature &uq); - ~ProductQuadrature() - override = default; + ~ProductQuadrature() override = default; int numEvals(int l) const override { diff --git a/dynare++/integ/cc/quadrature.hh b/dynare++/integ/cc/quadrature.hh index 0049caf54..aae21a274 100644 --- a/dynare++/integ/cc/quadrature.hh +++ b/dynare++/integ/cc/quadrature.hh @@ -30,7 +30,10 @@ #ifndef QUADRATURE_H #define QUADRATURE_H -#include +#include +#include +#include + #include "vector_function.hh" #include "int_sequence.hh" #include "sthread.hh" @@ -43,8 +46,7 @@ class OneDQuadrature { public: - virtual ~OneDQuadrature() - = default; + virtual ~OneDQuadrature() = default; virtual int numLevels() const = 0; virtual int numPoints(int level) const = 0; virtual double point(int level, int i) const = 0; @@ -72,8 +74,7 @@ public: Quadrature(int d) : dim(d) { } - virtual ~Quadrature() - = default; + virtual ~Quadrature() = default; int dimen() const { @@ -183,25 +184,26 @@ public: /* Just for debugging. */ void - savePoints(const char *fname, int level) const + savePoints(const string &fname, int level) const { - FILE *fd; - if (nullptr == (fd = fopen(fname, "w"))) + ofstream fd{fname, std::ios::out | std::ios::trunc}; + if (fd.fail()) { // todo: raise - fprintf(stderr, "Cannot open file %s for writing.\n", fname); - exit(1); + std::cerr << "Cannot open file " << fname << " for writing." << std::endl; + exit(EXIT_FAILURE); } _Tpit beg = begin(0, 1, level); _Tpit end = begin(1, 1, level); + fd << std::setprecision(12); 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++) - fprintf(fd, "\t%16.12g", run.point()[i]); - fprintf(fd, "\n"); + fd << '\t' << std::setw(16) << run.point()[i]; + fd << std::endl; } - fclose(fd); + fd.close(); } _Tpit @@ -244,8 +246,7 @@ public: { calcOffsets(); } - ~OneDPrecalcQuadrature() - override = default; + ~OneDPrecalcQuadrature() override = default; int numLevels() const override { diff --git a/dynare++/integ/cc/quasi_mcarlo.cc b/dynare++/integ/cc/quasi_mcarlo.cc index 4d0167a93..5eacd7f90 100644 --- a/dynare++/integ/cc/quasi_mcarlo.cc +++ b/dynare++/integ/cc/quasi_mcarlo.cc @@ -3,6 +3,8 @@ #include "quasi_mcarlo.hh" #include +#include +#include /* Here in the constructor, we have to calculate a maximum length of |coeff| array for a given |base| and given maximum |maxn|. After @@ -10,7 +12,7 @@ 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) + coeff(static_cast(floor(log(static_cast(maxn))/log(static_cast(b)))+2), 0) { int nr = num; j = -1; @@ -75,15 +77,14 @@ RadicalInverse::increase() void RadicalInverse::print() const { - printf("n=%d b=%d c=", num, base); + std::cout << "n=" << num << " b=" << base << " c="; 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[] = { +std::array 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, @@ -116,18 +117,6 @@ HaltonSequence::HaltonSequence(int n, int mxn, int dim, const PermutationScheme 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 calls |eval|. */ @@ -153,113 +142,48 @@ HaltonSequence::eval() void HaltonSequence::print() const { + auto ff = std::cout.flags(); for (const auto & i : ri) i.print(); - printf("point=[ "); + std::cout << "point=[ " + << std::fixed << std::setprecision(6); for (unsigned int i = 0; i < ri.size(); i++) - printf("%7.6f ", pt[i]); - printf("]\n"); -} - -qmcpit::qmcpit() - -{ + std::cout << std::setw(7) << pt[i] << ' '; + std::cout << ']' << std::endl; + std::cout.flags(ff); } qmcpit::qmcpit(const QMCSpecification &s, int n) - : spec(&s), halton(new HaltonSequence(n, s.level(), s.dimen(), s.getPerScheme())), - sig(new ParameterSignal(s.dimen())) + : spec(s), halton{n, s.level(), s.dimen(), s.getPerScheme()}, + 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 qmcpit::operator==(const qmcpit &qpit) const { - return (spec == qpit.spec) - && ((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; + return &spec == &qpit.spec && halton.getNum() == qpit.halton.getNum(); } qmcpit & qmcpit::operator++() { // todo: raise if |halton == null || qmcq == NULL| - halton->increase(); + halton.increase(); return *this; } double qmcpit::weight() const { - return 1.0/spec->level(); -} - -qmcnpit::qmcnpit() - : qmcpit() -{ + return 1.0/spec.level(); } 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 of the points in |NormalICDF| function. */ @@ -267,8 +191,8 @@ qmcnpit & qmcnpit::operator++() { qmcpit::operator++(); - for (int i = 0; i < halton->point().length(); i++) - (*pnt)[i] = NormalICDF::get(halton->point()[i]); + for (int i = 0; i < halton.point().length(); i++) + pnt[i] = NormalICDF::get(halton.point()[i]); return *this; } diff --git a/dynare++/integ/cc/quasi_mcarlo.hh b/dynare++/integ/cc/quasi_mcarlo.hh index 6e2477905..3118c8089 100644 --- a/dynare++/integ/cc/quasi_mcarlo.hh +++ b/dynare++/integ/cc/quasi_mcarlo.hh @@ -39,11 +39,9 @@ class PermutationScheme { public: - PermutationScheme() - = default; - virtual ~PermutationScheme() - = default; - virtual int permute(int i, int base, int c) const = 0; + PermutationScheme() = default; + virtual ~PermutationScheme() = default; + virtual int permute(int i, int base, int c) const = 0; }; /* This class represents an integer number |num| as @@ -65,12 +63,8 @@ class RadicalInverse IntSequence coeff; public: RadicalInverse(int n, int b, int mxn); - RadicalInverse(const RadicalInverse &ri) - - = default; - RadicalInverse & - operator=(const RadicalInverse &radi) - = default; + RadicalInverse(const RadicalInverse &ri) = default; + RadicalInverse &operator=(const RadicalInverse &radi) = default; double eval(const PermutationScheme &p) const; void increase(); void print() const; @@ -85,8 +79,7 @@ public: class HaltonSequence { private: - static int primes[]; - static int num_primes; + static std::array primes; protected: int num; int maxn; @@ -95,10 +88,8 @@ protected: Vector pt; public: HaltonSequence(int n, int mxn, int dim, const PermutationScheme &p); - HaltonSequence(const HaltonSequence &hs) - - = default; - const HaltonSequence &operator=(const HaltonSequence &hs); + HaltonSequence(const HaltonSequence &hs) = default; + HaltonSequence &operator=(const HaltonSequence &hs) = delete; void increase(); const Vector & point() const @@ -131,8 +122,7 @@ public: : dim(d), lev(l), per_scheme(p) { } - virtual ~QMCSpecification() - = default; + virtual ~QMCSpecification() = default; int dimen() const { @@ -153,44 +143,41 @@ public: /* This is an iterator for quasi Monte Carlo over a cube |QMCarloCubeQuadrature|. The iterator maintains |HaltonSequence| of the same dimension as given by the specification. An iterator can be - constructed from a given number |n|, or by a copy constructor. For - technical reasons, there is also an empty constructor; for that - reason, every member is a pointer. */ + constructed from a given number |n|, or by a copy constructor. */ class qmcpit { protected: - const QMCSpecification *spec{nullptr}; - HaltonSequence *halton{nullptr}; - ParameterSignal *sig{nullptr}; + const QMCSpecification &spec; + HaltonSequence halton; + ParameterSignal sig; public: - qmcpit(); qmcpit(const QMCSpecification &s, int n); - qmcpit(const qmcpit &qpit); - ~qmcpit(); + qmcpit(const qmcpit &qpit) = default; + virtual ~qmcpit() = default; bool operator==(const qmcpit &qpit) const; bool operator!=(const qmcpit &qpit) const { return !operator==(qpit); } - const qmcpit &operator=(const qmcpit &qpit); + qmcpit &operator=(const qmcpit &qpit) = delete; qmcpit &operator++(); const ParameterSignal & signal() const { - return *sig; + return sig; } const Vector & point() const { - return halton->point(); + return halton.point(); } double weight() const; void print() const { - halton->print(); + halton.print(); } }; @@ -206,8 +193,7 @@ public: : QuadratureImpl(d), QMCSpecification(d, l, p) { } - ~QMCarloCubeQuadrature() - override = default; + ~QMCarloCubeQuadrature() override = default; int numEvals(int l) const override { @@ -233,12 +219,11 @@ protected: class qmcnpit : public qmcpit { protected: - Vector *pnt{nullptr}; + Vector pnt; public: - qmcnpit(); qmcnpit(const QMCSpecification &spec, int n); - qmcnpit(const qmcnpit &qpit); - ~qmcnpit(); + qmcnpit(const qmcnpit &qpit) = default; + ~qmcnpit() override = default; bool operator==(const qmcnpit &qpit) const { @@ -249,22 +234,23 @@ public: { return !operator==(qpit); } - const qmcnpit &operator=(const qmcnpit &qpit); + qmcnpit &operator=(const qmcnpit &qpit) = delete; qmcnpit &operator++(); const ParameterSignal & signal() const { - return *sig; + return sig; } const Vector & point() const { - return *pnt; + return pnt; } void print() const { - halton->print(); pnt->print(); + halton.print(); + pnt.print(); } }; @@ -280,8 +266,7 @@ public: : QuadratureImpl(d), QMCSpecification(d, l, p) { } - ~QMCarloNormalQuadrature() - override = default; + ~QMCarloNormalQuadrature() override = default; int numEvals(int l) const override { diff --git a/dynare++/integ/cc/smolyak.cc b/dynare++/integ/cc/smolyak.cc index b764c38e5..c155e64e7 100644 --- a/dynare++/integ/cc/smolyak.cc +++ b/dynare++/integ/cc/smolyak.cc @@ -3,91 +3,27 @@ #include "smolyak.hh" #include "symmetry.hh" -smolpit::smolpit() - -{ -} +#include +#include /* This constructs a beginning of |isum| summand in |smolq|. We must be careful here, since |isum| can be past-the-end, so no reference to vectors in |smolq| by |isum| must be done in this case. */ 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())) + : smolq(q), isummand(isum), + jseq{q.dimen(), 0}, + sig{q.dimen()}, + p{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 = 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; + setPointAndWeight(); } bool smolpit::operator==(const smolpit &spit) const { - bool ret = true; - 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; + return &smolq == &spit.smolq && isummand == spit.isummand && jseq == spit.jseq; } /* We first try to increase index within the current summand. If we are @@ -98,22 +34,22 @@ 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]) + const IntSequence &levpts = smolq.levpoints[isummand]; + int i = smolq.dimen()-1; + jseq[i]++; + while (i >= 0 && jseq[i] == levpts[i]) { - (*jseq)[i] = 0; + jseq[i] = 0; i--; if (i >= 0) - (*jseq)[i]++; + jseq[i]++; } - sig->signalAfter(std::max(i, 0)); + sig.signalAfter(std::max(i, 0)); if (i < 0) isummand++; - if (isummand < smolq->numSummands()) + if (isummand < smolq.numSummands()) setPointAndWeight(); return *this; @@ -126,18 +62,18 @@ 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(); + // |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); + 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]); + int ki = (smolq.levels[isummand])[i]; + p[i] = (smolq.uquad).point(ki, jseq[i]); + w *= (smolq.uquad).weight(ki, jseq[i]); } } @@ -145,16 +81,19 @@ smolpit::setPointAndWeight() 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]); + auto ff = std::cout.flags(); + std::cout << "isum=" << std::left << std::setw(3) << isummand << std::right << ": ["; + for (int i = 0; i < smolq.dimen(); i++) + std::cout << std::setw(2) << (smolq.levels[isummand])[i] << ' '; + std::cout << "] j=["; + for (int i = 0; i < smolq.dimen(); i++) + std::cout << std::setw(2) << jseq[i] << ' '; + std::cout << std::showpos << std::fixed << std::setprecision(3) + << "] " << std::setw(4) << w << "*("; + 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 diff --git a/dynare++/integ/cc/smolyak.hh b/dynare++/integ/cc/smolyak.hh index dbfbfd0cf..83d0103e6 100644 --- a/dynare++/integ/cc/smolyak.hh +++ b/dynare++/integ/cc/smolyak.hh @@ -40,34 +40,33 @@ class SmolyakQuadrature; class smolpit { protected: - const SmolyakQuadrature *smolq{nullptr}; + const SmolyakQuadrature &smolq; unsigned int isummand{0}; - IntSequence *jseq{nullptr}; - ParameterSignal *sig{nullptr}; - Vector *p{nullptr}; + IntSequence jseq; + ParameterSignal sig; + Vector p; double w; public: - smolpit(); smolpit(const SmolyakQuadrature &q, unsigned int isum); - smolpit(const smolpit &spit); - ~smolpit(); + smolpit(const smolpit &spit) = default; + ~smolpit() = default; bool operator==(const smolpit &spit) const; bool operator!=(const smolpit &spit) const { return !operator==(spit); } - const smolpit &operator=(const smolpit &spit); + smolpit &operator=(const smolpit &spit) = delete; smolpit &operator++(); const ParameterSignal & signal() const { - return *sig; + return sig; } const Vector & point() const { - return *p; + return p; } double weight() const @@ -108,8 +107,7 @@ class SmolyakQuadrature : public QuadratureImpl PascalTriangle psc; public: SmolyakQuadrature(int d, int l, const OneDQuadrature &uq); - ~SmolyakQuadrature() - override = default; + ~SmolyakQuadrature() override = default; int numEvals(int level) const override; void designLevelForEvals(int max_eval, int &lev, int &evals) const; protected: diff --git a/dynare++/integ/cc/vector_function.cc b/dynare++/integ/cc/vector_function.cc index cc30f8123..92afd6ee7 100644 --- a/dynare++/integ/cc/vector_function.cc +++ b/dynare++/integ/cc/vector_function.cc @@ -5,24 +5,14 @@ #include #include - -#include #include /* Just an easy constructor of sequence of booleans defaulting to change everywhere. */ 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| @@ -31,47 +21,44 @@ ParameterSignal::ParameterSignal(const ParameterSignal &sig) void 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; - for (int i = l; i < num; i++) + for (size_t i = l; i < data.size(); 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) + : funcs(n) { 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 hard copies in others. */ VectorFunctionSet::VectorFunctionSet(VectorFunction &f, int n) - : funcs(n), first_shallow(true) + : funcs(n) { if (n > 0) funcs[0] = &f; 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 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()), + : VectorFunction(f), func(&f), A(vcov.numRows(), vcov.numRows()), multiplier(calcMultiplier()) { // 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 function as to be deleted. */ - -GaussConverterFunction::GaussConverterFunction(VectorFunction *f, const GeneralMatrix &vcov) - : VectorFunction(*f), func(f), delete_flag(true), A(vcov.numRows(), vcov.numRows()), +GaussConverterFunction::GaussConverterFunction(std::unique_ptr f, const GeneralMatrix &vcov) + : VectorFunction(*f), func_storage{move(f)}, func{func_storage.get()}, A(vcov.numRows(), vcov.numRows()), multiplier(calcMultiplier()) { // todo: raise if |A.numRows() != indim()| @@ -90,7 +76,7 @@ GaussConverterFunction::GaussConverterFunction(VectorFunction *f, const GeneralM } 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) { } diff --git a/dynare++/integ/cc/vector_function.hh b/dynare++/integ/cc/vector_function.hh index 8f9dab812..cbdf8cbe1 100644 --- a/dynare++/integ/cc/vector_function.hh +++ b/dynare++/integ/cc/vector_function.hh @@ -20,6 +20,7 @@ #include "GeneralMatrix.hh" #include +#include /* This is a simple class representing a vector of booleans. The items night be retrieved or changed, or can be set |true| after some @@ -31,22 +32,18 @@ class ParameterSignal { protected: - bool *data; - int num; + std::vector data; public: ParameterSignal(int n); - ParameterSignal(const ParameterSignal &sig); - ~ParameterSignal() - { - delete [] data; - } + ParameterSignal(const ParameterSignal &sig) = default; + ~ParameterSignal() = default; void signalAfter(int l); - const bool & + bool operator[](int i) const { return data[i]; } - bool & + std::vector::reference operator[](int i) { return data[i]; @@ -71,12 +68,9 @@ public: : in_dim(idim), out_dim(odim) { } - VectorFunction(const VectorFunction &func) - - = default; - virtual ~VectorFunction() - = default; - virtual VectorFunction *clone() const = 0; + VectorFunction(const VectorFunction &func) = default; + virtual ~VectorFunction() = default; + virtual std::unique_ptr clone() const = 0; virtual void eval(const Vector &point, const ParameterSignal &sig, Vector &out) = 0; int indim() const @@ -100,13 +94,15 @@ public: class VectorFunctionSet { +private: + // Stores the hard copies made by the class + std::vector> func_copies; protected: std::vector funcs; - bool first_shallow; public: VectorFunctionSet(const VectorFunction &f, int n); VectorFunctionSet(VectorFunction &f, int n); - ~VectorFunctionSet(); + ~VectorFunctionSet() = default; VectorFunction & getFunc(int i) { @@ -136,32 +132,28 @@ public: 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|. */ + constructed by the first constructor, the $f$ is assumed to be owned by the + caller. If the object of this class is copied, then $f$ is copied and hence + stored in a std::unique_ptr. The second constructor takes a smart pointer to + the function and in that case the class takes ownership of $f$. */ class GaussConverterFunction : public VectorFunction { +private: + std::unique_ptr func_storage; protected: VectorFunction *func; - bool delete_flag; GeneralMatrix A; double multiplier; public: GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov); - GaussConverterFunction(VectorFunction *f, const GeneralMatrix &vcov); + GaussConverterFunction(std::unique_ptr f, const GeneralMatrix &vcov); GaussConverterFunction(const GaussConverterFunction &f); - ~GaussConverterFunction() override - { - if (delete_flag) - delete func; - } - VectorFunction * + ~GaussConverterFunction() override = default; + std::unique_ptr clone() const override { - return new GaussConverterFunction(*this); + return std::make_unique(*this); } void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override; private: diff --git a/dynare++/integ/src/quadrature-points.cc b/dynare++/integ/src/quadrature-points.cc index d51434846..10cc32257 100644 --- a/dynare++/integ/src/quadrature-points.cc +++ b/dynare++/integ/src/quadrature-points.cc @@ -11,19 +11,21 @@ #include "integ/cc/product.hh" #include -#include #include - +#include +#include #include #include +#include +#include struct QuadParams { - const char *outname; - const char *vcovname; - int max_level; - double discard_weight; + string outname; + string vcovname; + int max_level{3}; + double discard_weight{0.0}; QuadParams(int argc, char **argv); void check_consistency() const; private: @@ -31,12 +33,12 @@ private: }; QuadParams::QuadParams(int argc, char **argv) - : outname(nullptr), vcovname(nullptr), max_level(3), discard_weight(0.0) { if (argc == 1) { // 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]; @@ -56,12 +58,24 @@ QuadParams::QuadParams(int argc, char **argv) switch (ret) { case opt_max_level: - if (1 != sscanf(optarg, "%d", &max_level)) - fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg); + try + { + max_level = std::stoi(string{optarg}); + } + catch (const std::invalid_argument &e) + { + std::cerr << "Couldn't parse integer " << optarg << ", ignored" << std::endl; + } break; case opt_discard_weight: - if (1 != sscanf(optarg, "%lf", &discard_weight)) - fprintf(stderr, "Couldn't parse float %s, ignored\n", optarg); + try + { + discard_weight = std::stod(string{optarg}); + } + catch (const std::invalid_argument &e) + { + std::cerr << "Couldn't parse float " << optarg << ", ignored" << std::endl; + } break; case opt_vcov: vcovname = optarg; @@ -75,41 +89,30 @@ QuadParams::QuadParams(int argc, char **argv) void QuadParams::check_consistency() const { - if (outname == nullptr) + if (outname.empty()) { - fprintf(stderr, "Error: output name not set\n"); - exit(1); + std::cerr << "Error: output name not set" << std::endl; + std::exit(EXIT_FAILURE); } - if (vcovname == nullptr) + if (vcovname.empty()) { - fprintf(stderr, "Error: vcov file name not set\n"); - exit(1); + std::cerr << "Error: vcov file name not set" << std::endl; + 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 main(int argc, char **argv) { QuadParams params(argc, argv); // open output file for writing - FILE *fout; - if (nullptr == (fout = fopen(params.outname, "w"))) + std::ofstream fout{params.outname, std::ios::out | std::ios::trunc}; + if (fout.fail()) { - fprintf(stderr, "Could not open %s for writing\n", params.outname); - exit(1); + std::cerr << "Could not open " << params.outname << " for writing" << std::endl; + std::exit(EXIT_FAILURE); } try @@ -142,23 +145,21 @@ main(int argc, char **argv) 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)); + std::cout << "Dimension: " << vcov.numRows() << std::endl + << "Maximum level: " << level << std::endl + << "Total number of nodes: " << sq.numEvals(level) << std::endl; // put the points to the vector - std::vector points; + std::vector> points; 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((const Vector &) qit.point())); // sort and uniq - OrderVec ordvec; - std::sort(points.begin(), points.end(), ordvec); + std::sort(points.begin(), points.end(), [](auto &a, auto &b) { return a.get() < b.get(); }); 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()); - 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 double mass = 0.0; @@ -175,41 +176,41 @@ main(int argc, char **argv) if (weight/mass < params.discard_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 int npoints = 0; double upscale_weight = 1/(mass-discard_mass); Vector x(vcov.numRows()); + fout << std::setprecision(16); 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]); + fout << std::setw(20) << 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"); + fout << ' ' << std::setw(20) << x[j]; + fout << std::endl; npoints++; } - printf("Final number of points: %d\n", npoints); - - fclose(fout); + std::cout << "Final number of points: " << npoints << std::endl; + fout.close(); } catch (const SylvException &e) { e.printMessage(); - return 1; + return EXIT_FAILURE; } catch (const ogu::Exception &e) { e.print(); - return 1; + return EXIT_FAILURE; } - return 0; + return EXIT_SUCCESS; } diff --git a/dynare++/integ/testing/tests.cc b/dynare++/integ/testing/tests.cc index fb88c5146..f3d375338 100644 --- a/dynare++/integ/testing/tests.cc +++ b/dynare++/integ/testing/tests.cc @@ -14,10 +14,14 @@ #include "product.hh" #include "quasi_mcarlo.hh" -#include -#include -#include +#include +#include #include +#include +#include +#include +#include +#include const int num_threads = 2; // does nothing if DEBUG defined @@ -33,13 +37,12 @@ public: D(inD), k(kk) { } - MomentFunction(const MomentFunction &func) - - = default; - VectorFunction * + MomentFunction(const MomentFunction &func) = default; + + std::unique_ptr clone() const override { - return new MomentFunction(*this); + return std::make_unique(*this); } 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()) { - printf("Wrong length of vectors in MomentFunction::eval\n"); - exit(1); + std::cerr << "Wrong length of vectors in MomentFunction::eval" << std::endl; + std::exit(EXIT_FAILURE); } Vector y(point); y.zeros(); @@ -68,13 +71,12 @@ public: : VectorFunction(nvar, UFSTensor::calcMaxOffset(nvar, kk)), k(kk) { } - TensorPower(const TensorPower &func) - - = default; - VectorFunction * + TensorPower(const TensorPower &func) = default; + + std::unique_ptr clone() const override { - return new TensorPower(*this); + return std::make_unique(*this); } 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()) { - printf("Wrong length of vectors in TensorPower::eval\n"); - exit(1); + std::cerr << "Wrong length of vectors in TensorPower::eval" << std::endl; + std::exit(EXIT_FAILURE); } URSingleTensor ypow(point, k); out.zeros(); @@ -107,10 +109,10 @@ public: : VectorFunction(f.indim(), f.outdim()), dim(f.dim) { } - VectorFunction * + std::unique_ptr clone() const override { - return new Function1(*this); + return std::make_unique(*this); } 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) { - printf("Wrong length of vectors in Function1::eval\n"); - exit(1); + std::cerr << "Wrong length of vectors in Function1::eval" << std::endl; + std::exit(EXIT_FAILURE); } double r = 1; for (int i = 0; i < dim; i++) @@ -140,13 +142,12 @@ public: : Function1(d) { } - Function1Trans(const Function1Trans &func) - - = default; - VectorFunction * + Function1Trans(const Function1Trans &func) = default; + + std::unique_ptr clone() const override { - return new Function1Trans(*this); + return std::make_unique(*this); } 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 class WallTimer { - char mes[100]; - struct timeval start; + string mes; + std::chrono::time_point start; bool new_line; 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() { - struct timeval end; - gettimeofday(&end, nullptr); - printf("%s%8.4g", mes, - end.tv_sec-start.tv_sec + (end.tv_usec-start.tv_usec)*1.0e-6); + auto end = std::chrono::high_resolution_clock::now(); + std::chrono::duration duration = end - start; + std::cout << mes << std::setw(8) << std::setprecision(4) << duration.count(); if (new_line) - printf("\n"); + std::cout << std::endl; } }; @@ -190,22 +190,17 @@ public: /****************************************************/ class TestRunnable { - char name[100]; public: + const string name; 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) + TestRunnable(string name_arg, int d, int nv) + : name{move(name_arg)}, dim(d), nvar(nv) { - strncpy(name, n, 100); } + virtual ~TestRunnable() = default; 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); @@ -218,7 +213,7 @@ protected: bool TestRunnable::test() const { - printf("Running test <%s>\n", name); + cout << "Running test <" << name << ">" << endl; bool passed; { WallTimer tim("Wall clock time ", false); @@ -226,12 +221,12 @@ TestRunnable::test() const } if (passed) { - printf("............................ passed\n\n"); + std::cout << "............................ passed" << std::endl << std::endl; return passed; } else { - printf("............................ FAILED\n\n"); + std::cout << "............................ FAILED" << std::endl << std::endl; return passed; } } @@ -258,13 +253,13 @@ TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level 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)); + std::cout << "\tNumber of Smolyak evaluations: " << quad.numEvals(level) << std::endl; } // 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()); + std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << smol_out.getMax() << std::endl; return smol_out.getMax() < 1.e-7; } @@ -287,13 +282,13 @@ TestRunnable::product_normal_moments(const GeneralMatrix &m, int imom, int level GaussHermite gs; ProductQuadrature quad(dim, gs); 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 UNormalMoments moments(imom, msq); 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; } @@ -318,8 +313,8 @@ TestRunnable::qmc_normal_moments(const GeneralMatrix &m, int imom, int level) WarnockPerScheme wps; ReversePerScheme rps; IdentityPerScheme ips; - PermutationScheme *scheme[] = {&wps, &rps, &ips}; - const char *labs[] = {"Warnock", "Reverse", "Identity"}; + array scheme = {&wps, &rps, &ips}; + array labs = {"Warnock", "Reverse", "Identity"}; // theoretical result int dim = mchol.numRows(); @@ -329,21 +324,19 @@ TestRunnable::qmc_normal_moments(const GeneralMatrix &m, int imom, int level) // quasi monte carlo normal quadrature double max_error = 0.0; 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]; - sprintf(mes, "\tQMC normal quadrature time %8s: ", labs[i]); - WallTimer tim(mes); + ostringstream mes; + mes << "\tQMC normal quadrature time " << std::setw(8) << labs[i] << ": "; + WallTimer tim(mes.str()); 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()); + 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) - { - max_error = qmc_out.getMax(); - } + max_error = qmc_out.getMax(); } 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()) { - fprintf(stderr, "Incompatible dimensions of check value and function.\n"); - exit(1); + std::cerr << "Incompatible dimensions of check value and function." << std::endl; + std::exit(EXIT_FAILURE); } GaussLegendre glq; @@ -369,8 +362,8 @@ TestRunnable::smolyak_product_cube(const VectorFunction &func, const Vector &res 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); + std::cout << "\tNumber of Smolyak evaluations: " << quad.numEvals(level) << std::endl; + std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << smol_error << std::endl; } { 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); 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); + std::cout << "\tNumber of product evaluations: " << quad.numEvals(level) << std::endl; + std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << prod_error << std::endl; } 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.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); + std::cout << "\tQuasi-Monte Carlo (Warnock scrambling) error: " << std::setw(16) << std::setprecision(12) << error1 << std::endl; } double error2; { @@ -408,8 +400,7 @@ TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int l // 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); + std::cout << "\tQuasi-Monte Carlo (reverse scrambling) error: " << std::setw(16) << std::setprecision(12) << error2 << std::endl; } double error3; { @@ -419,8 +410,7 @@ TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int l // 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); + std::cout << "\tQuasi-Monte Carlo (no scrambling) error: " << std::setw(16) << std::setprecision(12) << error3 << std::endl; } return error1 < tol && error2 < tol && error3 < tol; @@ -574,61 +564,54 @@ public: int main() { - TestRunnable *all_tests[50]; + std::vector> all_tests; // 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.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); + all_tests.push_back(std::make_unique()); /* - all_tests[num_tests++] = new F1GaussLegendre(); - all_tests[num_tests++] = new F1QuasiMCarlo(); + all_tests.push_back(make_unique()); + all_tests.push_back(make_unique()); */ // find maximum dimension and maximum nvar int dmax = 0; int nvmax = 0; - for (int i = 0; i < num_tests; i++) + for (const auto &test : all_tests) { - if (dmax < all_tests[i]->dim) - dmax = all_tests[i]->dim; - if (nvmax < all_tests[i]->nvar) - nvmax = all_tests[i]->nvar; + if (dmax < test->dim) + dmax = test->dim; + if (nvmax < test->nvar) + nvmax = test->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++) + for (const auto &test : all_tests) { try { - if (all_tests[i]->test()) + if (test->test()) success++; } 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(); } 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(); } } - printf("There were %d tests that failed out of %d tests run.\n", - num_tests - success, num_tests); + std::cout << "There were " << all_tests.size() - success << " tests that failed out of " + << all_tests.size() << " tests run." << std::endl; - // destroy - for (int i = 0; i < num_tests; i++) - { - delete all_tests[i]; - } - - return 0; + return EXIT_SUCCESS; } diff --git a/dynare++/kord/global_check.hh b/dynare++/kord/global_check.hh index 4c88d6a9d..e49d23923 100644 --- a/dynare++/kord/global_check.hh +++ b/dynare++/kord/global_check.hh @@ -37,6 +37,8 @@ #include +#include + #include "vector_function.hh" #include "quadrature.hh" @@ -76,10 +78,10 @@ public: ResidFunction(const ResidFunction &rf); ~ResidFunction() override; - VectorFunction * + std::unique_ptr clone() const override { - return new ResidFunction(*this); + return std::make_unique(*this); } void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override; void setYU(const Vector &ys, const Vector &xx); @@ -91,18 +93,15 @@ class GResidFunction : public GaussConverterFunction { public: GResidFunction(const Approximation &app) - : GaussConverterFunction(new ResidFunction(app), app.getModel().getVcov()) + : GaussConverterFunction(std::make_unique(app), app.getModel().getVcov()) { } - GResidFunction(const GResidFunction &rf) - - = default; - ~GResidFunction() - override = default; - VectorFunction * + GResidFunction(const GResidFunction &rf) = default; + ~GResidFunction() override = default; + std::unique_ptr clone() const override { - return new GResidFunction(*this); + return std::make_unique(*this); } void setYU(const Vector &ys, const Vector &xx)