diff --git a/dynare++/kord/approximation.cc b/dynare++/kord/approximation.cc
index 12c0dd0c9eece64ba7c1573984e8ee3960a1e6d4..e04807abeadbb28756749cc83f226b476485acab 100644
--- a/dynare++/kord/approximation.cc
+++ b/dynare++/kord/approximation.cc
@@ -81,13 +81,13 @@ Approximation::approxAtSteady()
model.getVcov(), journal);
korder.switchToFolded();
for (int k = 2; k <= model.order(); k++)
- korder.performStep<KOrder::fold>(k);
+ korder.performStep<Storage::fold>(k);
saveRuleDerivs(korder.getFoldDers());
}
else
{
- FirstOrderDerivs<KOrder::fold> fo_ders(fo);
+ FirstOrderDerivs<Storage::fold> fo_ders(fo);
saveRuleDerivs(fo_ders);
}
check(0.0);
@@ -148,7 +148,7 @@ Approximation::walkStochSteady()
/* We form the |DRFixPoint| object from the last rule with
$\sigma=dsigma$. Then we save the steady state to |ss|. The new steady
is also put to |model.getSteady()|. */
- DRFixPoint<KOrder::fold> fp(*rule_ders, ypart, model.getSteady(), dsigma);
+ DRFixPoint<Storage::fold> fp(*rule_ders, ypart, model.getSteady(), dsigma);
bool converged = fp.calcFixPoint(DecisionRule::emethod::horner, model.getSteady());
JournalRecord rec(journal);
rec << "Fix point calcs: iter=" << fp.getNumIter() << ", newton_iter="
@@ -171,7 +171,7 @@ Approximation::walkStochSteady()
Vector dy(model.getSteady());
dy.add(-1.0, last_steady);
- StochForwardDerivs<KOrder::fold> hh(ypart, model.nexog(), *rule_ders_ss, mom, dy,
+ StochForwardDerivs<Storage::fold> hh(ypart, model.nexog(), *rule_ders_ss, mom, dy,
dsigma, sigma_so_far);
JournalRecord rec1(journal);
rec1 << "Calculation of g** expectations done" << endrec;
@@ -183,7 +183,7 @@ Approximation::walkStochSteady()
KOrderStoch korder_stoch(ypart, model.nexog(), model.getModelDerivatives(),
hh, journal);
for (int d = 1; d <= model.order(); d++)
- korder_stoch.performStep<KOrder::fold>(d);
+ korder_stoch.performStep<Storage::fold>(d);
saveRuleDerivs(korder_stoch.getFoldDers());
@@ -198,7 +198,7 @@ Approximation::walkStochSteady()
if (steps == 0 && dr_centralize)
{
// centralize decision rule for zero steps
- DRFixPoint<KOrder::fold> fp(*rule_ders, ypart, model.getSteady(), 1.0);
+ DRFixPoint<Storage::fold> fp(*rule_ders, ypart, model.getSteady(), 1.0);
bool converged = fp.calcFixPoint(DecisionRule::emethod::horner, model.getSteady());
JournalRecord rec(journal);
rec << "Fix point calcs: iter=" << fp.getNumIter() << ", newton_iter="
diff --git a/dynare++/kord/decision_rule.cc b/dynare++/kord/decision_rule.cc
index 10de34f239d88b010b35aefb59f8cb8b003dfc36..deee95946dea01cdcce3d11e26765f93476a65ed 100644
--- a/dynare++/kord/decision_rule.cc
+++ b/dynare++/kord/decision_rule.cc
@@ -15,20 +15,20 @@
// |FoldDecisionRule| conversion from |UnfoldDecisionRule|
FoldDecisionRule::FoldDecisionRule(const UnfoldDecisionRule &udr)
- : DecisionRuleImpl<KOrder::fold>(ctraits<KOrder::fold>::Tpol(udr.nrows(), udr.nvars()),
+ : DecisionRuleImpl<Storage::fold>(ctraits<Storage::fold>::Tpol(udr.nrows(), udr.nvars()),
udr.ypart, udr.nu, udr.ysteady)
{
for (const auto &it : udr)
- insert(std::make_unique<ctraits<KOrder::fold>::Ttensym>(*(it.second)));
+ insert(std::make_unique<ctraits<Storage::fold>::Ttensym>(*(it.second)));
}
// |UnfoldDecisionRule| conversion from |FoldDecisionRule|
UnfoldDecisionRule::UnfoldDecisionRule(const FoldDecisionRule &fdr)
- : DecisionRuleImpl<KOrder::unfold>(ctraits<KOrder::unfold>::Tpol(fdr.nrows(), fdr.nvars()),
+ : DecisionRuleImpl<Storage::unfold>(ctraits<Storage::unfold>::Tpol(fdr.nrows(), fdr.nvars()),
fdr.ypart, fdr.nu, fdr.ysteady)
{
for (const auto &it : fdr)
- insert(std::make_unique<ctraits<KOrder::unfold>::Ttensym>(*(it.second)));
+ insert(std::make_unique<ctraits<Storage::unfold>::Ttensym>(*(it.second)));
}
/* This runs simulations with an output to journal file. Note that we
diff --git a/dynare++/kord/decision_rule.hh b/dynare++/kord/decision_rule.hh
index ed0165c7f827608b18d8acbc0fe51058fb8a7f8a..ff201fd062dd5a641f83a30fc6c68b4b20c25a77 100644
--- a/dynare++/kord/decision_rule.hh
+++ b/dynare++/kord/decision_rule.hh
@@ -90,25 +90,28 @@ public:
\left[\matrix{y^*_{t-1}-\tilde y^*\cr u_t}\right]^{\alpha_m}.$$
We provide a method and a constructor to transform a rule to the centralized form.
- The class is templated, the template argument is either |KOrder::fold|
- or |KOrder::unfold|. So, there are two implementations of |DecisionRule|
+ The class is templated, the template argument is either |Storage::fold|
+ or |Storage::unfold|. So, there are two implementations of |DecisionRule|
interface. */
-template <int t>
+template <Storage t>
class DecisionRuleImpl : public ctraits<t>::Tpol, public DecisionRule
{
protected:
- using _Tparent = typename ctraits<t>::Tpol;
+ using _Tpol = typename ctraits<t>::Tpol;
+ using _Tg = typename ctraits<t>::Tg;
+ using _Ttensor = typename ctraits<t>::Ttensor;
+ using _Ttensym = typename ctraits<t>::Ttensym;
const Vector ysteady;
const PartitionY ypart;
const int nu;
public:
- DecisionRuleImpl(const _Tparent &pol, const PartitionY &yp, int nuu,
+ DecisionRuleImpl(const _Tpol &pol, const PartitionY &yp, int nuu,
const ConstVector &ys)
: ctraits<t>::Tpol(pol), ysteady(ys), ypart(yp), nu(nuu)
{
}
- DecisionRuleImpl(_Tparent &pol, const PartitionY &yp, int nuu,
+ DecisionRuleImpl(_Tpol &pol, const PartitionY &yp, int nuu,
const ConstVector &ys)
: ctraits<t>::Tpol(0, yp.ny(), pol), ysteady(ys), ypart(yp),
nu(nuu)
@@ -190,7 +193,7 @@ protected:
So we go through $i+j=d=0\ldots q$ and in each loop we form the fully
symmetric tensor $[g_{(yu)^l}]$ and insert it to the container. */
-template <int t>
+template <Storage t>
void
DecisionRuleImpl<t>::fillTensors(const _Tg &g, double sigma)
{
@@ -244,7 +247,7 @@ DecisionRuleImpl<t>::fillTensors(const _Tg &g, double sigma)
$$dstate=\left[\matrix{\tilde y^*-\bar y^*\cr 0}\right],$$
where $\bar y$ is the steady state of the original rule |dr|. */
-template <int t>
+template <Storage t>
void
DecisionRuleImpl<t>::centralize(const DecisionRuleImpl &dr)
{
@@ -278,7 +281,7 @@ DecisionRuleImpl<t>::centralize(const DecisionRuleImpl &dr)
|ysteady| is canceled from |ystart|, we simulate, and at the end
|ysteady| is added to all columns of the result. */
-template <int t>
+template <Storage t>
TwoDMatrix
DecisionRuleImpl<t>::simulate(emethod em, int np, const ConstVector &ystart,
ShockRealization &sr) const
@@ -347,7 +350,7 @@ DecisionRuleImpl<t>::simulate(emethod em, int np, const ConstVector &ystart,
that the steady state (fix point) is cancelled and added once. Hence
we have two special methods. */
-template <int t>
+template <Storage t>
void
DecisionRuleImpl<t>::evaluate(emethod em, Vector &out, const ConstVector &ys,
const ConstVector &u) const
@@ -370,7 +373,7 @@ DecisionRuleImpl<t>::evaluate(emethod em, Vector &out, const ConstVector &ys,
/* This is easy. We just return the newly created copy using the
centralized constructor. */
-template <int t>
+template <Storage t>
std::unique_ptr<DecisionRule>
DecisionRuleImpl<t>::centralizedClone(const Vector &fixpoint) const
{
@@ -380,19 +383,19 @@ DecisionRuleImpl<t>::centralizedClone(const Vector &fixpoint) const
/* Here we only encapsulate two implementations to one, deciding
according to the parameter. */
-template <int t>
+template <Storage t>
void
DecisionRuleImpl<t>::eval(emethod em, Vector &out, const ConstVector &v) const
{
if (em == emethod::horner)
- _Tparent::evalHorner(out, v);
+ _Tpol::evalHorner(out, v);
else
- _Tparent::evalTrad(out, v);
+ _Tpol::evalTrad(out, v);
}
/* Write the decision rule and steady state to the MAT file. */
-template <int t>
+template <Storage t>
void
DecisionRuleImpl<t>::writeMat(mat_t *fd, const std::string &prefix) const
{
@@ -402,63 +405,63 @@ DecisionRuleImpl<t>::writeMat(mat_t *fd, const std::string &prefix) const
ConstTwoDMatrix(dum).writeMat(fd, prefix + "_ss");
}
-/* This is exactly the same as |DecisionRuleImpl<KOrder::fold>|. The
+/* This is exactly the same as |DecisionRuleImpl<Storage::fold>|. The
only difference is that we have a conversion from
|UnfoldDecisionRule|, which is exactly
- |DecisionRuleImpl<KOrder::unfold>|. */
+ |DecisionRuleImpl<Storage::unfold>|. */
class UnfoldDecisionRule;
-class FoldDecisionRule : public DecisionRuleImpl<KOrder::fold>
+class FoldDecisionRule : public DecisionRuleImpl<Storage::fold>
{
friend class UnfoldDecisionRule;
public:
- FoldDecisionRule(const ctraits<KOrder::fold>::Tpol &pol, const PartitionY &yp, int nuu,
+ FoldDecisionRule(const ctraits<Storage::fold>::Tpol &pol, const PartitionY &yp, int nuu,
const ConstVector &ys)
- : DecisionRuleImpl<KOrder::fold>(pol, yp, nuu, ys)
+ : DecisionRuleImpl<Storage::fold>(pol, yp, nuu, ys)
{
}
- FoldDecisionRule(ctraits<KOrder::fold>::Tpol &pol, const PartitionY &yp, int nuu,
+ FoldDecisionRule(ctraits<Storage::fold>::Tpol &pol, const PartitionY &yp, int nuu,
const ConstVector &ys)
- : DecisionRuleImpl<KOrder::fold>(pol, yp, nuu, ys)
+ : DecisionRuleImpl<Storage::fold>(pol, yp, nuu, ys)
{
}
- FoldDecisionRule(const ctraits<KOrder::fold>::Tg &g, const PartitionY &yp, int nuu,
+ FoldDecisionRule(const ctraits<Storage::fold>::Tg &g, const PartitionY &yp, int nuu,
const ConstVector &ys, double sigma)
- : DecisionRuleImpl<KOrder::fold>(g, yp, nuu, ys, sigma)
+ : DecisionRuleImpl<Storage::fold>(g, yp, nuu, ys, sigma)
{
}
- FoldDecisionRule(const DecisionRuleImpl<KOrder::fold> &dr, const ConstVector &fixpoint)
- : DecisionRuleImpl<KOrder::fold>(dr, fixpoint)
+ FoldDecisionRule(const DecisionRuleImpl<Storage::fold> &dr, const ConstVector &fixpoint)
+ : DecisionRuleImpl<Storage::fold>(dr, fixpoint)
{
}
FoldDecisionRule(const UnfoldDecisionRule &udr);
};
-/* This is exactly the same as |DecisionRuleImpl<KOrder::unfold>|, but
+/* This is exactly the same as |DecisionRuleImpl<Storage::unfold>|, but
with a conversion from |FoldDecisionRule|, which is exactly
- |DecisionRuleImpl<KOrder::fold>|. */
+ |DecisionRuleImpl<Storage::fold>|. */
-class UnfoldDecisionRule : public DecisionRuleImpl<KOrder::unfold>
+class UnfoldDecisionRule : public DecisionRuleImpl<Storage::unfold>
{
friend class FoldDecisionRule;
public:
- UnfoldDecisionRule(const ctraits<KOrder::unfold>::Tpol &pol, const PartitionY &yp, int nuu,
+ UnfoldDecisionRule(const ctraits<Storage::unfold>::Tpol &pol, const PartitionY &yp, int nuu,
const ConstVector &ys)
- : DecisionRuleImpl<KOrder::unfold>(pol, yp, nuu, ys)
+ : DecisionRuleImpl<Storage::unfold>(pol, yp, nuu, ys)
{
}
- UnfoldDecisionRule(ctraits<KOrder::unfold>::Tpol &pol, const PartitionY &yp, int nuu,
+ UnfoldDecisionRule(ctraits<Storage::unfold>::Tpol &pol, const PartitionY &yp, int nuu,
const ConstVector &ys)
- : DecisionRuleImpl<KOrder::unfold>(pol, yp, nuu, ys)
+ : DecisionRuleImpl<Storage::unfold>(pol, yp, nuu, ys)
{
}
- UnfoldDecisionRule(const ctraits<KOrder::unfold>::Tg &g, const PartitionY &yp, int nuu,
+ UnfoldDecisionRule(const ctraits<Storage::unfold>::Tg &g, const PartitionY &yp, int nuu,
const ConstVector &ys, double sigma)
- : DecisionRuleImpl<KOrder::unfold>(g, yp, nuu, ys, sigma)
+ : DecisionRuleImpl<Storage::unfold>(g, yp, nuu, ys, sigma)
{
}
- UnfoldDecisionRule(const DecisionRuleImpl<KOrder::unfold> &dr, const ConstVector &fixpoint)
- : DecisionRuleImpl<KOrder::unfold>(dr, fixpoint)
+ UnfoldDecisionRule(const DecisionRuleImpl<Storage::unfold> &dr, const ConstVector &fixpoint)
+ : DecisionRuleImpl<Storage::unfold>(dr, fixpoint)
{
}
UnfoldDecisionRule(const FoldDecisionRule &udr);
@@ -477,18 +480,21 @@ public:
Jacobian of the solved system is given by derivatives stored in
|bigfder|. */
-template <int t>
+template <Storage t>
class DRFixPoint : public ctraits<t>::Tpol
{
- using _Tparent = typename ctraits<t>::Tpol;
+ using _Tpol = typename ctraits<t>::Tpol;
+ using _Tg = typename ctraits<t>::Tg;
+ using _Ttensor = typename ctraits<t>::Ttensor;
+ using _Ttensym = typename ctraits<t>::Ttensym;
constexpr static int max_iter = 10000;
constexpr static int max_newton_iter = 50;
constexpr static int newton_pause = 100;
constexpr static double tol = 1e-10;
const Vector ysteady;
const PartitionY ypart;
- std::unique_ptr<_Tparent> bigf;
- std::unique_ptr<_Tparent> bigfder;
+ std::unique_ptr<_Tpol> bigf;
+ std::unique_ptr<_Tpol> bigfder;
public:
using emethod = typename DecisionRule::emethod;
DRFixPoint(const _Tg &g, const PartitionY &yp,
@@ -526,19 +532,19 @@ private:
(|Symmetry(1)|). Then the derivative of the $F$ polynomial is
calculated. */
-template <int t>
+template <Storage t>
DRFixPoint<t>::DRFixPoint(const _Tg &g, const PartitionY &yp,
const Vector &ys, double sigma)
: ctraits<t>::Tpol(yp.ny(), yp.nys()),
ysteady(ys), ypart(yp)
{
fillTensors(g, sigma);
- _Tparent yspol(ypart.nstat, ypart.nys(), *this);
- bigf = std::make_unique<_Tparent>(const_cast<const _Tparent &>(yspol));
+ _Tpol yspol(ypart.nstat, ypart.nys(), *this);
+ bigf = std::make_unique<_Tpol>(const_cast<const _Tpol &>(yspol));
_Ttensym &frst = bigf->get(Symmetry{1});
for (int i = 0; i < ypart.nys(); i++)
frst.get(i, i) = frst.get(i, i) - 1;
- bigfder = std::make_unique<_Tparent>(*bigf, 0);
+ bigfder = std::make_unique<_Tpol>(*bigf, 0);
}
/* Here we fill the tensors for the |DRFixPoint| class. We ignore the
@@ -547,7 +553,7 @@ DRFixPoint<t>::DRFixPoint(const _Tg &g, const PartitionY &yp,
dimension and |d|, and add ${\sigma^k\over d!k!}g_{y^d\sigma^k}$ to
the tensor $g_{(y)^d}$. */
-template <int t>
+template <Storage t>
void
DRFixPoint<t>::fillTensors(const _Tg &g, double sigma)
{
@@ -584,7 +590,7 @@ DRFixPoint<t>::fillTensors(const _Tg &g, double sigma)
|urelax| is less that |urelax_threshold|, we stop searching and stop
the Newton. */
-template <int t>
+template <Storage t>
bool
DRFixPoint<t>::solveNewton(Vector &y)
{
@@ -661,7 +667,7 @@ DRFixPoint<t>::solveNewton(Vector &y)
The |out| vector is not touched if the algorithm has not convered. */
-template <int t>
+template <Storage t>
bool
DRFixPoint<t>::calcFixPoint(emethod em, Vector &out)
{
@@ -695,7 +701,7 @@ DRFixPoint<t>::calcFixPoint(emethod em, Vector &out)
if (converged)
{
- _Tparent::evalHorner(out, ystar);
+ _Tpol::evalHorner(out, ystar);
out.add(1.0, ysteady);
}
diff --git a/dynare++/kord/dynamic_model.cc b/dynare++/kord/dynamic_model.cc
index 6766a706015ed8014b4a229239c76d338bdf2396..118e7f1251a2ec3c4114b0bbf8eba9ac6a4a564e 100644
--- a/dynare++/kord/dynamic_model.cc
+++ b/dynare++/kord/dynamic_model.cc
@@ -2,31 +2,31 @@
#include "dynamic_model.hh"
-#include <cstring>
+#include <iostream>
+#include <algorithm>
void
NameList::print() const
{
for (int i = 0; i < getNum(); i++)
- printf("%s\n", getName(i));
+ std::cout << getName(i) << '\n';
}
void
-NameList::writeMat(mat_t *fd, const char *vname) const
+NameList::writeMat(mat_t *fd, const std::string &vname) const
{
int maxlen = 0;
for (int i = 0; i < getNum(); i++)
- if (maxlen < (int) strlen(getName(i)))
- maxlen = (int) strlen(getName(i));
+ maxlen = std::max(maxlen, static_cast<int>(getName(i).size()));
if (maxlen == 0)
return;
- auto *m = new char[getNum()*maxlen];
+ std::vector<char> m(getNum()*maxlen);
for (int i = 0; i < getNum(); i++)
for (int j = 0; j < maxlen; j++)
- if (j < (int) strlen(getName(i)))
+ if (j < static_cast<int>(getName(i).size()))
m[j*getNum()+i] = getName(i)[j];
else
m[j*getNum()+i] = ' ';
@@ -35,23 +35,20 @@ NameList::writeMat(mat_t *fd, const char *vname) const
dims[0] = getNum();
dims[1] = maxlen;
- matvar_t *v = Mat_VarCreate(vname, MAT_C_CHAR, MAT_T_UINT8, 2, dims, m, 0);
+ matvar_t *v = Mat_VarCreate(vname.c_str(), MAT_C_CHAR, MAT_T_UINT8, 2, dims, m.data(), 0);
Mat_VarWrite(fd, v, MAT_COMPRESSION_NONE);
Mat_VarFree(v);
- delete[] m;
}
void
-NameList::writeMatIndices(mat_t *fd, const char *prefix) const
+NameList::writeMatIndices(mat_t *fd, const std::string &prefix) const
{
- char tmp[100];
TwoDMatrix aux(1, 1);
for (int i = 0; i < getNum(); i++)
{
- sprintf(tmp, "%s_i_%s", prefix, getName(i));
aux.get(0, 0) = i+1;
- aux.writeMat(fd, tmp);
+ aux.writeMat(fd, prefix + "_i_" + getName(i));
}
}
diff --git a/dynare++/kord/dynamic_model.hh b/dynare++/kord/dynamic_model.hh
index 58021f01d05997d2605a49384971d6f51d36df7c..df8609b99df831f2d32e735893b79a1bb57dc2ab 100644
--- a/dynare++/kord/dynamic_model.hh
+++ b/dynare++/kord/dynamic_model.hh
@@ -2,7 +2,7 @@
// Dynamic model abstraction
-/* This file only defines a generic interface to an SDGE model. The model
+/* This file only defines a generic interface to a DSGE model. The model
takes the form:
$$E_t\left[f(g^{**}(g^*(y,u_t),u_{t+1}),g(y,u),y,u_t)\right]=0$$
The interface is defined via pure virtual class |DynamicModel|. */
@@ -15,22 +15,24 @@
#include "Vector.hh"
+#include <memory>
+#include <string>
+
/* The class is a virtual pure class which provides an access to names
of the variables. */
class NameList
{
public:
- virtual ~NameList()
- = default;
+ virtual ~NameList() = default;
virtual int getNum() const = 0;
- virtual const char *getName(int i) const = 0;
+ virtual const std::string &getName(int i) const = 0;
void print() const;
- void writeMat(mat_t *fd, const char *vname) const;
- void writeMatIndices(mat_t *fd, const char *prefix) const;
+ void writeMat(mat_t *fd, const std::string &vname) const;
+ void writeMatIndices(mat_t *fd, const std::string &prefix) const;
};
-/* This is the interface to an information on a generic SDGE
+/* This is the interface to an information on a generic DSGE
model. It is sufficient for calculations of policy rule Taylor
approximations at some (not necessarily deterministic) steady state.
@@ -82,9 +84,8 @@ public:
class DynamicModel
{
public:
- virtual DynamicModel *clone() const = 0;
- virtual ~DynamicModel()
- = default;
+ virtual std::unique_ptr<DynamicModel> clone() const = 0;
+ virtual ~DynamicModel() = default;
virtual int nstat() const = 0;
virtual int nboth() const = 0;
@@ -98,13 +99,13 @@ public:
return nstat()+nboth()+npred()+nforw();
}
- virtual const NameList&getAllEndoNames() const = 0;
- virtual const NameList&getStateNames() const = 0;
- virtual const NameList&getExogNames() const = 0;
- virtual const TwoDMatrix&getVcov() const = 0;
- virtual const TensorContainer<FSSparseTensor>&getModelDerivatives() const = 0;
- virtual const Vector&getSteady() const = 0;
- virtual Vector&getSteady() = 0;
+ virtual const NameList &getAllEndoNames() const = 0;
+ virtual const NameList &getStateNames() const = 0;
+ virtual const NameList &getExogNames() const = 0;
+ virtual const TwoDMatrix &getVcov() const = 0;
+ virtual const TensorContainer<FSSparseTensor> &getModelDerivatives() const = 0;
+ virtual const Vector &getSteady() const = 0;
+ virtual Vector &getSteady() = 0;
virtual void solveDeterministicSteady() = 0;
virtual void evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) = 0;
diff --git a/dynare++/kord/faa_di_bruno.cc b/dynare++/kord/faa_di_bruno.cc
index d303d2b945ac79c6017defdd85dd3f43775c86f4..5cde7551bf502375df389162636166ddd338d7e4 100644
--- a/dynare++/kord/faa_di_bruno.cc
+++ b/dynare++/kord/faa_di_bruno.cc
@@ -5,8 +5,6 @@
#include <cmath>
-double FaaDiBruno::magic_mult = 1.5;
-
// |FaaDiBruno::calculate| folded sparse code
/* We take an opportunity to refine the stack container to avoid
allocation of more memory than available. */
@@ -121,19 +119,19 @@ FaaDiBruno::estimRefinment(const TensorDimens &tdims, int nr, int l,
int &avmem_mb, int &tmpmem_mb)
{
int nthreads = sthread::detach_thread_group::max_parallel_threads;
- long int per_size1 = tdims.calcUnfoldMaxOffset();
- auto per_size2 = (long int) pow((double) tdims.getNVS().getMax(), l);
+ long per_size1 = tdims.calcUnfoldMaxOffset();
+ long per_size2 = static_cast<long>(std::pow(tdims.getNVS().getMax(), l));
double lambda = 0.0;
- long int per_size = sizeof(double)*nr
- *(long int) (lambda*per_size1+(1-lambda)*per_size2);
- long int mem = SystemResources::availableMemory();
+ long per_size = sizeof(double)*nr
+ *static_cast<long>(lambda*per_size1+(1-lambda)*per_size2);
+ long mem = SystemResources::availableMemory();
int max = 0;
- double num_cols = ((double) (mem-magic_mult*nthreads*per_size))
+ double num_cols = static_cast<double>(mem-magic_mult*nthreads*per_size)
/nthreads/sizeof(double)/nr;
if (num_cols > 0)
{
- double maxd = pow(num_cols, ((double) 1)/l);
- max = (int) floor(maxd);
+ double maxd = std::pow(num_cols, 1.0/l);
+ max = static_cast<int>(std::floor(maxd));
}
if (max == 0)
{
@@ -145,6 +143,6 @@ FaaDiBruno::estimRefinment(const TensorDimens &tdims, int nr, int l,
rec << endrec;
}
avmem_mb = mem/1024/1024;
- tmpmem_mb = (nthreads*per_size)/1024/1024;
+ tmpmem_mb = nthreads*per_size/1024/1024;
return max;
}
diff --git a/dynare++/kord/faa_di_bruno.hh b/dynare++/kord/faa_di_bruno.hh
index 07913e2b8baf2ed43ff2c32d044d7ded89698bca..aabcff0276036d8db0f1a2bb071b282ca5bcf815 100644
--- a/dynare++/kord/faa_di_bruno.hh
+++ b/dynare++/kord/faa_di_bruno.hh
@@ -38,7 +38,7 @@ public:
UGSTensor &out);
protected:
int estimRefinment(const TensorDimens &tdims, int nr, int l, int &avmem_mb, int &tmpmem_mb);
- static double magic_mult;
+ constexpr static double magic_mult = 1.5;
};
#endif
diff --git a/dynare++/kord/first_order.cc b/dynare++/kord/first_order.cc
index 21de0808bd151f5d46f0238969c42359e798cd02..d5672342937c85c42257d58965daf029baf01f71 100644
--- a/dynare++/kord/first_order.cc
+++ b/dynare++/kord/first_order.cc
@@ -5,16 +5,17 @@
#include <dynlapack.h>
-double qz_criterium = 1.000001;
+double FirstOrder::qz_criterium_global;
+std::mutex FirstOrder::mut;
/* This is a function which selects the eigenvalues pair used by
|dgges|. See documentation to DGGES for details. Here we want
to select (return true) the pairs for which $\alpha<\beta$. */
lapack_int
-order_eigs(const double *alphar, const double *alphai, const double *beta)
+FirstOrder::order_eigs(const double *alphar, const double *alphai, const double *beta)
{
- return (*alphar **alphar + *alphai **alphai < *beta **beta * qz_criterium * qz_criterium);
+ return (*alphar **alphar + *alphai **alphai < *beta **beta * qz_criterium_global * qz_criterium_global);
}
/* Here we solve the linear approximation. The result are the matrices
@@ -36,8 +37,6 @@ FirstOrder::solve(const TwoDMatrix &fd)
JournalRecordPair pa(journal);
pa << "Recovering first order derivatives " << endrec;
- ::qz_criterium = FirstOrder::qz_criterium;
-
// **********************
// solve derivatives |gy|
// **********************
@@ -164,21 +163,22 @@ FirstOrder::solve(const TwoDMatrix &fd)
lapack_int ldvsl = vsl.getLD(), ldvsr = vsr.getLD();
lapack_int lwork = 100*n+16;
Vector work(lwork);
- auto *bwork = new lapack_int[n];
+ std::vector<lapack_int> bwork(n);
lapack_int info;
lapack_int sdim2 = sdim;
- dgges("N", "V", "S", order_eigs, &n, matE.getData().base(), &lda,
+ {
+ std::lock_guard<std::mutex> lk{mut};
+ qz_criterium_global = qz_criterium;
+ dgges("N", "V", "S", order_eigs, &n, matE.getData().base(), &lda,
matD.getData().base(), &ldb, &sdim2, alphar.base(), alphai.base(),
beta.base(), vsl.getData().base(), &ldvsl, vsr.getData().base(), &ldvsr,
- work.base(), &lwork, bwork, &info);
+ work.base(), &lwork, bwork.data(), &info);
+ }
if (info)
- {
- throw KordException(__FILE__, __LINE__,
- "DGGES returns an error in FirstOrder::solve");
- }
+ throw KordException(__FILE__, __LINE__,
+ "DGGES returns an error in FirstOrder::solve");
sdim = sdim2;
bk_cond = (sdim == ypart.nys());
- delete[] bwork;
// make submatrices of right space
/* Here we setup submatrices of the matrix $Z$. */
@@ -214,7 +214,7 @@ FirstOrder::solve(const TwoDMatrix &fd)
gy.place(fder, ypart.nstat+ypart.nys(), 0);
// check difference for derivatives of both
- GeneralMatrix bder((const GeneralMatrix &)sfder, ypart.nstat, 0, ypart.nboth, ypart.nys());
+ GeneralMatrix bder(const_cast<const GeneralMatrix &>(sfder), ypart.nstat, 0, ypart.nboth, ypart.nys());
GeneralMatrix bder2(preder, ypart.npred, 0, ypart.nboth, ypart.nys());
bder.add(-1, bder2);
b_error = bder.getData().getMax();
@@ -248,10 +248,8 @@ FirstOrder::solve(const TwoDMatrix &fd)
journalEigs();
if (!gy.isFinite() || !gu.isFinite())
- {
- throw KordException(__FILE__, __LINE__,
- "NaN or Inf asserted in first order derivatives in FirstOrder::solve");
- }
+ throw KordException(__FILE__, __LINE__,
+ "NaN or Inf asserted in first order derivatives in FirstOrder::solve");
}
void
@@ -269,24 +267,22 @@ FirstOrder::journalEigs()
<< " " << "npred=" << ypart.nys() << endrec;
}
if (!bk_cond)
- {
- for (int i = 0; i < alphar.length(); i++)
- {
- if (i == sdim || i == ypart.nys())
- {
- JournalRecord jr(journal);
- jr << "---------------------------------------------------- ";
- if (i == sdim)
- jr << "sdim";
- else
- jr << "npred";
- jr << endrec;
- }
- JournalRecord jr(journal);
- double mod = sqrt(alphar[i]*alphar[i]+alphai[i]*alphai[i]);
- mod = mod/round(100000*std::abs(beta[i]))*100000;
- jr << i << "\t(" << alphar[i] << "," << alphai[i] << ") / " << beta[i]
- << " \t" << mod << endrec;
- }
- }
+ for (int i = 0; i < alphar.length(); i++)
+ {
+ if (i == sdim || i == ypart.nys())
+ {
+ JournalRecord jr(journal);
+ jr << "---------------------------------------------------- ";
+ if (i == sdim)
+ jr << "sdim";
+ else
+ jr << "npred";
+ jr << endrec;
+ }
+ JournalRecord jr(journal);
+ double mod = std::sqrt(alphar[i]*alphar[i]+alphai[i]*alphai[i]);
+ mod = mod/std::round(100000*std::abs(beta[i]))*100000;
+ jr << i << "\t(" << alphar[i] << "," << alphai[i] << ") / " << beta[i]
+ << " \t" << mod << endrec;
+ }
}
diff --git a/dynare++/kord/first_order.hh b/dynare++/kord/first_order.hh
index 79045a9378570928f0ac985c24c6a21ac9091482..03326927677f28d41dadbb4d9b6f43833045fc89 100644
--- a/dynare++/kord/first_order.hh
+++ b/dynare++/kord/first_order.hh
@@ -7,11 +7,13 @@
#include "korder.hh"
-template<int>
+#include <mutex>
+
+template<Storage>
class FirstOrderDerivs;
class FirstOrder
{
- template <int>
+ template <Storage>
friend class FirstOrderDerivs;
PartitionY ypart;
int nu;
@@ -25,6 +27,17 @@ class FirstOrder
Vector beta;
double qz_criterium;
Journal &journal;
+
+ // Passed to LAPACK's DGGES
+ static lapack_int order_eigs(const double *alphar, const double *alphai, const double *beta);
+
+ // The value of qz_criterium_global used by the order_eigs function
+ /* NB: we have no choice but to use a global variable, since LAPACK won't
+ take a closure */
+ static double qz_criterium_global;
+
+ // Protects the static qz_criterium_global
+ static std::mutex mut;
public:
FirstOrder(int num_stat, int num_pred, int num_both, int num_forw,
int num_u, const FSSparseTensor &f, Journal &jr, double qz_crit)
@@ -63,7 +76,7 @@ protected:
/* This class only converts the derivatives $g_{y^*}$ and $g_u$ to a
folded or unfolded container. */
-template <int t>
+template <Storage t>
class FirstOrderDerivs : public ctraits<t>::Tg
{
public:
@@ -71,11 +84,11 @@ public:
: ctraits<t>::Tg(4)
{
IntSequence nvs{fo.ypart.nys(), fo.nu, fo.nu, 1};
- auto ten = std::make_unique<_Ttensor>(fo.ypart.ny(), TensorDimens(Symmetry{1, 0, 0, 0}, nvs));
+ auto ten = std::make_unique<typename ctraits<t>::Ttensor>(fo.ypart.ny(), TensorDimens(Symmetry{1, 0, 0, 0}, nvs));
ten->zeros();
ten->add(1.0, fo.gy);
this->insert(std::move(ten));
- ten = std::make_unique<_Ttensor>(fo.ypart.ny(), TensorDimens(Symmetry{0, 1, 0, 0}, nvs));
+ ten = std::make_unique<typename ctraits<t>::Ttensor>(fo.ypart.ny(), TensorDimens(Symmetry{0, 1, 0, 0}, nvs));
ten->zeros();
ten->add(1.0, fo.gu);
this->insert(std::move(ten));
diff --git a/dynare++/kord/global_check.cc b/dynare++/kord/global_check.cc
index 43027d5e5e05eea0ea6dd7c4bb627ed42417fe75..80f66eeef7ddb454e065abaf6133cca4ab8dadd6 100644
--- a/dynare++/kord/global_check.cc
+++ b/dynare++/kord/global_check.cc
@@ -9,6 +9,7 @@
#include "product.hh"
#include "quasi_mcarlo.hh"
+#include <utility>
#include <cmath>
/* Here we just set a reference to the approximation, and create a new
@@ -16,37 +17,21 @@
ResidFunction::ResidFunction(const Approximation &app)
: VectorFunction(app.getModel().nexog(), app.getModel().numeq()), approx(app),
- model(app.getModel().clone()),
- yplus(nullptr), ystar(nullptr), u(nullptr), hss(nullptr)
+ model(app.getModel().clone())
{
}
ResidFunction::ResidFunction(const ResidFunction &rf)
- : VectorFunction(rf), approx(rf.approx), model(rf.model->clone()),
- yplus(nullptr), ystar(nullptr), u(nullptr), hss(nullptr)
+ : VectorFunction(rf), approx(rf.approx), model(rf.model->clone())
{
if (rf.yplus)
- yplus = new Vector(*(rf.yplus));
+ yplus = std::make_unique<Vector>(*(rf.yplus));
if (rf.ystar)
- ystar = new Vector(*(rf.ystar));
+ ystar = std::make_unique<Vector>(*(rf.ystar));
if (rf.u)
- u = new Vector(*(rf.u));
+ u = std::make_unique<Vector>(*(rf.u));
if (rf.hss)
- hss = new FTensorPolynomial(*(rf.hss));
-}
-
-ResidFunction::~ResidFunction()
-{
- delete model;
- // delete |y| and |u| dependent data
- if (yplus)
- delete yplus;
- if (ystar)
- delete ystar;
- if (u)
- delete u;
- if (hss)
- delete hss;
+ hss = std::make_unique<FTensorPolynomial>(*(rf.hss));
}
/* This sets $y^*$ and $u$. We have to create |ystar|, |u|, |yplus| and
@@ -55,58 +40,41 @@ ResidFunction::~ResidFunction()
void
ResidFunction::setYU(const ConstVector &ys, const ConstVector &xx)
{
- // delete |y| and |u| dependent data
- /* NB: code shared with the destructor */
- if (yplus)
- delete yplus;
- if (ystar)
- delete ystar;
- if (u)
- delete u;
- if (hss)
- delete hss;
-
- ystar = new Vector(ys);
- u = new Vector(xx);
- yplus = new Vector(model->numeq());
+ ystar = std::make_unique<Vector>(ys);
+ u = std::make_unique<Vector>(xx);
+ yplus = std::make_unique<Vector>(model->numeq());
approx.getFoldDecisionRule().evaluate(DecisionRule::emethod::horner,
*yplus, *ystar, *u);
// make a tensor polynomial of in-place subtensors from decision rule
- /* Here we use a dirty tricky of converting |const| to non-|const| to
- obtain a polynomial of subtensor corresponding to non-predetermined
- variables. However, this new non-|const| polynomial will be used for a
+ /* Note that the non-|const| polynomial will be used for a
construction of |hss| and will be used in |const| context. So this
- dirty thing is safe.
+ const_cast is safe.
Note, that there is always a folded decision rule in |Approximation|. */
- union {const FoldDecisionRule *c; FoldDecisionRule *n;} dr;
- dr.c = &(approx.getFoldDecisionRule());
+ const FoldDecisionRule &dr = approx.getFoldDecisionRule();
FTensorPolynomial dr_ss(model->nstat()+model->npred(), model->nboth()+model->nforw(),
- *(dr.n));
+ const_cast<FoldDecisionRule &>(dr));
// make |ytmp_star| be a difference of |yplus| from steady
Vector ytmp_star(ConstVector(*yplus, model->nstat(), model->npred()+model->nboth()));
- ConstVector ysteady_star(dr.c->getSteady(), model->nstat(),
+ ConstVector ysteady_star(dr.getSteady(), model->nstat(),
model->npred()+model->nboth());
ytmp_star.add(-1.0, ysteady_star);
// make |hss| and add steady to it
/* Here is the |const| context of |dr_ss|. */
- hss = new FTensorPolynomial(dr_ss, ytmp_star);
- ConstVector ysteady_ss(dr.c->getSteady(), model->nstat()+model->npred(),
+ hss = std::make_unique<FTensorPolynomial>(dr_ss, ytmp_star);
+ ConstVector ysteady_ss(dr.getSteady(), model->nstat()+model->npred(),
model->nboth()+model->nforw());
if (hss->check(Symmetry{0}))
- {
- hss->get(Symmetry{0}).getData().add(1.0, ysteady_ss);
- }
+ hss->get(Symmetry{0}).getData().add(1.0, ysteady_ss);
else
{
auto ten = std::make_unique<FFSTensor>(hss->nrows(), hss->nvars(), 0);
ten->getData() = ysteady_ss;
hss->insert(std::move(ten));
}
-
}
/* Here we evaluate the residual $F(y^*,u,u')$. We have to evaluate |hss|
@@ -133,7 +101,7 @@ GlobalChecker::check(const Quadrature &quad, int level,
const ConstVector &ys, const ConstVector &x, Vector &out)
{
for (int ifunc = 0; ifunc < vfs.getNum(); ifunc++)
- ((GResidFunction &) (vfs.getFunc(ifunc))).setYU(ys, x);
+ dynamic_cast<GResidFunction &>(vfs.getFunc(ifunc)).setYU(ys, x);
quad.integrate(vfs, level, out);
}
@@ -168,13 +136,13 @@ GlobalChecker::check(int max_evals, const ConstTwoDMatrix &y,
bool take_smolyak = (smol_evals < prod_evals) && (smol_level >= prod_level-1);
- Quadrature *quad;
+ std::unique_ptr<Quadrature> quad;
int lev;
// create the quadrature and report the decision
if (take_smolyak)
{
- quad = new SmolyakQuadrature(model.nexog(), smol_level, gh);
+ quad = std::make_unique<SmolyakQuadrature>(model.nexog(), smol_level, gh);
lev = smol_level;
JournalRecord rec(journal);
rec << "Selected Smolyak (level,evals)=(" << smol_level << ","
@@ -183,7 +151,7 @@ GlobalChecker::check(int max_evals, const ConstTwoDMatrix &y,
}
else
{
- quad = new ProductQuadrature(model.nexog(), gh);
+ quad = std::make_unique<ProductQuadrature>(model.nexog(), gh);
lev = prod_level;
JournalRecord rec(journal);
rec << "Selected product (level,evals)=(" << prod_level << ","
@@ -201,8 +169,6 @@ GlobalChecker::check(int max_evals, const ConstTwoDMatrix &y,
Vector outj{out.getCol(j)};
check(*quad, lev, yj, xj, outj);
}
-
- delete quad;
}
/* This method checks an error of the approximation by evaluating
@@ -211,7 +177,7 @@ GlobalChecker::check(int max_evals, const ConstTwoDMatrix &y,
$-mult\cdot\sigma$ to $mult\cdot\sigma$ in |m| steps. */
void
-GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const char *prefix,
+GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const std::string &prefix,
int m, double mult, int max_evals)
{
JournalRecordPair pa(journal);
@@ -219,15 +185,15 @@ GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const char *prefix,
<< mult << " std errors, granularity " << m << endrec;
// setup |y_mat| of steady states for checking
- TwoDMatrix y_mat(model.numeq(), 2*m *model.nexog()+1);
+ TwoDMatrix y_mat(model.numeq(), 2*m*model.nexog()+1);
for (int j = 0; j < 2*m*model.nexog()+1; j++)
{
Vector yj{y_mat.getCol(j)};
- yj = (const Vector &) model.getSteady();
+ yj = model.getSteady();
}
// setup |exo_mat| for checking
- TwoDMatrix exo_mat(model.nexog(), 2*m *model.nexog()+1);
+ TwoDMatrix exo_mat(model.nexog(), 2*m*model.nexog()+1);
exo_mat.zeros();
for (int ishock = 0; ishock < model.nexog(); ishock++)
{
@@ -235,12 +201,11 @@ GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const char *prefix,
for (int j = 0; j < 2*m; j++)
{
int jmult = (j < m) ? j-m : j-m+1;
- exo_mat.get(ishock, 1+2*m*ishock+j)
- = mult*jmult*max_sigma/m;
+ exo_mat.get(ishock, 1+2*m*ishock+j) = mult*jmult*max_sigma/m;
}
}
- TwoDMatrix errors(model.numeq(), 2*m *model.nexog()+1);
+ TwoDMatrix errors(model.numeq(), 2*m*model.nexog()+1);
check(max_evals, y_mat, exo_mat, errors);
// report errors along shock and save them
@@ -248,12 +213,9 @@ GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const char *prefix,
JournalRecord rec(journal);
rec << "Shock value error" << endrec;
ConstVector err0{errors.getCol(0)};
- char shock[9];
- char erbuf[17];
for (int ishock = 0; ishock < model.nexog(); ishock++)
{
TwoDMatrix err_out(model.numeq(), 2*m+1);
- sprintf(shock, "%-8s", model.getExogNames().getName(ishock));
for (int j = 0; j < 2*m+1; j++)
{
int jj;
@@ -273,13 +235,12 @@ GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const char *prefix,
error = err0;
}
JournalRecord rec1(journal);
- sprintf(erbuf, "%12.7g ", error.getMax());
- rec1 << shock << " " << exo_mat.get(ishock, jj)
- << "\t" << erbuf << endrec;
+ std::string shockname{model.getExogNames().getName(ishock)};
+ shockname.resize(8, ' ');
+ rec1 << shockname << ' ' << exo_mat.get(ishock, jj)
+ << "\t" << error.getMax() << endrec;
}
- char tmp[100];
- sprintf(tmp, "%s_shock_%s_errors", prefix, model.getExogNames().getName(ishock));
- err_out.writeMat(fd, tmp);
+ err_out.writeMat(fd, prefix + "_shock_" + model.getExogNames().getName(ishock) + "_errors");
}
}
@@ -298,7 +259,7 @@ GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const char *prefix,
the results. */
void
-GlobalChecker::checkOnEllipseAndSave(mat_t *fd, const char *prefix,
+GlobalChecker::checkOnEllipseAndSave(mat_t *fd, const std::string &prefix,
int m, double mult, int max_evals)
{
JournalRecordPair pa(journal);
@@ -310,7 +271,7 @@ GlobalChecker::checkOnEllipseAndSave(mat_t *fd, const char *prefix,
a submatrix of covariances of all endogenous variables. The submatrix
corresponds to state variables (predetermined plus both). */
TwoDMatrix ycov{approx.calcYCov()};
- TwoDMatrix ycovpred((const TwoDMatrix &) ycov, model.nstat(), model.nstat(),
+ TwoDMatrix ycovpred(const_cast<const TwoDMatrix &>(ycov), model.nstat(), model.nstat(),
model.npred()+model.nboth(), model.npred()+model.nboth());
SymSchurDecomp ssd(ycovpred);
ssd.correctDefinitness(1.e-05);
@@ -381,18 +342,15 @@ GlobalChecker::checkOnEllipseAndSave(mat_t *fd, const char *prefix,
TwoDMatrix out(model.numeq(), ymat.ncols());
check(max_evals, ymat, umat, out);
- char tmp[100];
- sprintf(tmp, "%s_ellipse_points", prefix);
- ymat.writeMat(fd, tmp);
- sprintf(tmp, "%s_ellipse_errors", prefix);
- out.writeMat(fd, tmp);
+ ymat.writeMat(fd, prefix + "_ellipse_points");
+ out.writeMat(fd, prefix + "_ellipse_errors");
}
/* Here we check the errors along a simulation. We simulate, then set
|x| to zeros, check and save results. */
void
-GlobalChecker::checkAlongSimulationAndSave(mat_t *fd, const char *prefix,
+GlobalChecker::checkAlongSimulationAndSave(mat_t *fd, const std::string &prefix,
int m, int max_evals)
{
JournalRecordPair pa(journal);
@@ -406,9 +364,6 @@ GlobalChecker::checkAlongSimulationAndSave(mat_t *fd, const char *prefix,
TwoDMatrix out(model.numeq(), m);
check(max_evals, y, x, out);
- char tmp[100];
- sprintf(tmp, "%s_simul_points", prefix);
- y.writeMat(fd, tmp);
- sprintf(tmp, "%s_simul_errors", prefix);
- out.writeMat(fd, tmp);
+ y.writeMat(fd, prefix + "_simul_points");
+ out.writeMat(fd, prefix + "_simul_errors");
}
diff --git a/dynare++/kord/global_check.hh b/dynare++/kord/global_check.hh
index e35be844e65f63d9dc46d0c3e355eba6917b1ba6..c25536fa61c476b70c7539596576659ded5f05d1 100644
--- a/dynare++/kord/global_check.hh
+++ b/dynare++/kord/global_check.hh
@@ -68,16 +68,13 @@ class ResidFunction : public VectorFunction
{
protected:
const Approximation ≈
- DynamicModel *model;
- Vector *yplus;
- Vector *ystar;
- Vector *u;
- FTensorPolynomial *hss;
+ std::unique_ptr<DynamicModel> model;
+ std::unique_ptr<Vector> yplus, ystar, u;
+ std::unique_ptr<FTensorPolynomial> hss;
public:
ResidFunction(const Approximation &app);
ResidFunction(const ResidFunction &rf);
- ~ResidFunction() override;
std::unique_ptr<VectorFunction>
clone() const override
{
@@ -96,8 +93,6 @@ public:
: GaussConverterFunction(std::make_unique<ResidFunction>(app), app.getModel().getVcov())
{
}
- GResidFunction(const GResidFunction &rf) = default;
- ~GResidFunction() override = default;
std::unique_ptr<VectorFunction>
clone() const override
{
@@ -106,7 +101,7 @@ public:
void
setYU(const ConstVector &ys, const ConstVector &xx)
{
- ((ResidFunction *) func)->setYU(ys, xx);
+ dynamic_cast<ResidFunction *>(func)->setYU(ys, xx);
}
};
@@ -137,13 +132,13 @@ public:
}
void check(int max_evals, const ConstTwoDMatrix &y,
const ConstTwoDMatrix &x, TwoDMatrix &out);
- void checkAlongShocksAndSave(mat_t *fd, const char *prefix,
+ void checkAlongShocksAndSave(mat_t *fd, const std::string &prefix,
int m, double mult, int max_evals);
- void checkOnEllipseAndSave(mat_t *fd, const char *prefix,
+ void checkOnEllipseAndSave(mat_t *fd, const std::string &prefix,
int m, double mult, int max_evals);
- void checkAlongSimulationAndSave(mat_t *fd, const char *prefix,
+ void checkAlongSimulationAndSave(mat_t *fd, const std::string &prefix,
int m, int max_evals);
- void checkUnconditionalAndSave(mat_t *fd, const char *prefix,
+ void checkUnconditionalAndSave(mat_t *fd, const std::string &prefix,
int m, int max_evals);
protected:
void check(const Quadrature &quad, int level,
diff --git a/dynare++/kord/journal.cc b/dynare++/kord/journal.cc
index 5f42e3384bcc5a1fa20260e10265b033092606e4..c38128bf668e719708bd428cd16bf0b2af9b385c 100644
--- a/dynare++/kord/journal.cc
+++ b/dynare++/kord/journal.cc
@@ -6,6 +6,7 @@
#include <iomanip>
#include <cmath>
#include <ctime>
+#include <limits>
#ifndef __MINGW32__
# include <sys/time.h> // For getrusage()
@@ -89,8 +90,8 @@ SystemResources::SystemResources()
idrss = rus.ru_idrss;
majflt = rus.ru_majflt;
#else
- utime = -1.0;
- stime = -1.0;
+ utime = std::numeric_limits<double>::quiet_NaN();
+ stime = std::numeric_limits<double>::quiet_NaN();
idrss = -1;
majflt = -1;
#endif
@@ -98,7 +99,7 @@ SystemResources::SystemResources()
#ifndef __MINGW32__
getloadavg(&load_avg, 1);
#else
- load_avg = -1.0;
+ load_avg = std::numeric_limits<double>::quiet_NaN();
#endif
pg_avail = availablePhysicalPages();
diff --git a/dynare++/kord/korder.cc b/dynare++/kord/korder.cc
index fec3b0ac8efae3a014ec5731339cecb4f34733af..162025e9a3916afcd1c477ce86ea9d23a846718e 100644
--- a/dynare++/kord/korder.cc
+++ b/dynare++/kord/korder.cc
@@ -3,12 +3,6 @@
#include "kord_exception.hh"
#include "korder.hh"
-PLUMatrix::PLUMatrix(const PLUMatrix &plu)
- : TwoDMatrix(plu), inv(plu.inv), ipiv(new lapack_int[nrows()])
-{
- memcpy(ipiv, plu.ipiv, nrows()*sizeof(lapack_int));
-}
-
/* Here we set |ipiv| and |inv| members of the |PLUMatrix| depending on
its content. It is assumed that subclasses will call this method at
the end of their constructors. */
@@ -18,8 +12,8 @@ PLUMatrix::calcPLU()
{
lapack_int info;
lapack_int rows = nrows(), lda = ld;
- inv = (const Vector &) getData();
- dgetrf(&rows, &rows, inv.base(), &lda, ipiv, &info);
+ inv = const_cast<const Vector &>(getData());
+ dgetrf(&rows, &rows, inv.base(), &lda, ipiv.data(), &info);
}
/* Here we just call the LAPACK machinery to multiply by the inverse. */
@@ -34,9 +28,8 @@ PLUMatrix::multInv(TwoDMatrix &m) const
lapack_int mcols = m.ncols();
lapack_int mrows = m.nrows();
lapack_int ldb = m.getLD();
- double *mbase = m.getData().base();
- dgetrs("N", &mrows, &mcols, inv.base(), &lda, ipiv,
- mbase, &ldb, &info);
+ dgetrs("N", &mrows, &mcols, inv.base(), &lda, ipiv.data(),
+ m.getData().base(), &ldb, &info);
KORD_RAISE_IF(info != 0,
"Info!=0 in PLUMatrix::multInv");
}
@@ -101,116 +94,116 @@ MatrixS::MatrixS(const FSSparseTensor &f, const IntSequence &ss,
interesting here. */
template<>
-ctraits<KOrder::unfold>::Tg& KOrder::g<KOrder::unfold>()
+ctraits<Storage::unfold>::Tg& KOrder::g<Storage::unfold>()
{
return _ug;
}
template<>
-const ctraits<KOrder::unfold>::Tg& KOrder::g<KOrder::unfold>() const
+const ctraits<Storage::unfold>::Tg& KOrder::g<Storage::unfold>() const
{ return _ug;}
template<>
-ctraits<KOrder::fold>::Tg& KOrder::g<KOrder::fold>()
+ctraits<Storage::fold>::Tg& KOrder::g<Storage::fold>()
{
return _fg;
}
template<>
-const ctraits<KOrder::fold>::Tg& KOrder::g<KOrder::fold>() const
+const ctraits<Storage::fold>::Tg& KOrder::g<Storage::fold>() const
{ return _fg;}
template<>
-ctraits<KOrder::unfold>::Tgs& KOrder::gs<KOrder::unfold>()
+ctraits<Storage::unfold>::Tgs& KOrder::gs<Storage::unfold>()
{
return _ugs;
}
template<>
-const ctraits<KOrder::unfold>::Tgs& KOrder::gs<KOrder::unfold>() const
+const ctraits<Storage::unfold>::Tgs& KOrder::gs<Storage::unfold>() const
{ return _ugs;}
template<>
-ctraits<KOrder::fold>::Tgs& KOrder::gs<KOrder::fold>()
+ctraits<Storage::fold>::Tgs& KOrder::gs<Storage::fold>()
{
return _fgs;
}
template<>
-const ctraits<KOrder::fold>::Tgs& KOrder::gs<KOrder::fold>() const
+const ctraits<Storage::fold>::Tgs& KOrder::gs<Storage::fold>() const
{ return _fgs;}
template<>
-ctraits<KOrder::unfold>::Tgss& KOrder::gss<KOrder::unfold>()
+ctraits<Storage::unfold>::Tgss& KOrder::gss<Storage::unfold>()
{
return _ugss;
}
template<>
-const ctraits<KOrder::unfold>::Tgss& KOrder::gss<KOrder::unfold>() const
+const ctraits<Storage::unfold>::Tgss& KOrder::gss<Storage::unfold>() const
{ return _ugss;}
template<>
-ctraits<KOrder::fold>::Tgss& KOrder::gss<KOrder::fold>()
+ctraits<Storage::fold>::Tgss& KOrder::gss<Storage::fold>()
{
return _fgss;
}
template<>
-const ctraits<KOrder::fold>::Tgss& KOrder::gss<KOrder::fold>() const
+const ctraits<Storage::fold>::Tgss& KOrder::gss<Storage::fold>() const
{ return _fgss;}
template<>
-ctraits<KOrder::unfold>::TG& KOrder::G<KOrder::unfold>()
+ctraits<Storage::unfold>::TG& KOrder::G<Storage::unfold>()
{
return _uG;
}
template<>
-const ctraits<KOrder::unfold>::TG& KOrder::G<KOrder::unfold>() const
+const ctraits<Storage::unfold>::TG& KOrder::G<Storage::unfold>() const
{ return _uG;}
template<>
-ctraits<KOrder::fold>::TG& KOrder::G<KOrder::fold>()
+ctraits<Storage::fold>::TG& KOrder::G<Storage::fold>()
{
return _fG;
}
template<>
-const ctraits<KOrder::fold>::TG& KOrder::G<KOrder::fold>() const
+const ctraits<Storage::fold>::TG& KOrder::G<Storage::fold>() const
{ return _fG;}
template<>
-ctraits<KOrder::unfold>::TZstack& KOrder::Zstack<KOrder::unfold>()
+ctraits<Storage::unfold>::TZstack& KOrder::Zstack<Storage::unfold>()
{
return _uZstack;
}
template<>
-const ctraits<KOrder::unfold>::TZstack& KOrder::Zstack<KOrder::unfold>() const
+const ctraits<Storage::unfold>::TZstack& KOrder::Zstack<Storage::unfold>() const
{ return _uZstack;}
template<>
-ctraits<KOrder::fold>::TZstack& KOrder::Zstack<KOrder::fold>()
+ctraits<Storage::fold>::TZstack& KOrder::Zstack<Storage::fold>()
{
return _fZstack;
}
template<>
-const ctraits<KOrder::fold>::TZstack& KOrder::Zstack<KOrder::fold>() const
+const ctraits<Storage::fold>::TZstack& KOrder::Zstack<Storage::fold>() const
{ return _fZstack;}
template<>
-ctraits<KOrder::unfold>::TGstack& KOrder::Gstack<KOrder::unfold>()
+ctraits<Storage::unfold>::TGstack& KOrder::Gstack<Storage::unfold>()
{
return _uGstack;
}
template<>
-const ctraits<KOrder::unfold>::TGstack& KOrder::Gstack<KOrder::unfold>() const
+const ctraits<Storage::unfold>::TGstack& KOrder::Gstack<Storage::unfold>() const
{ return _uGstack;}
template<>
-ctraits<KOrder::fold>::TGstack& KOrder::Gstack<KOrder::fold>()
+ctraits<Storage::fold>::TGstack& KOrder::Gstack<Storage::fold>()
{
return _fGstack;
}
template<>
-const ctraits<KOrder::fold>::TGstack& KOrder::Gstack<KOrder::fold>() const
+const ctraits<Storage::fold>::TGstack& KOrder::Gstack<Storage::fold>() const
{ return _fGstack;}
template<>
-ctraits<KOrder::unfold>::Tm& KOrder::m<KOrder::unfold>()
+ctraits<Storage::unfold>::Tm& KOrder::m<Storage::unfold>()
{
return _um;
}
template<>
-const ctraits<KOrder::unfold>::Tm& KOrder::m<KOrder::unfold>() const
+const ctraits<Storage::unfold>::Tm& KOrder::m<Storage::unfold>() const
{ return _um;}
template<>
-ctraits<KOrder::fold>::Tm& KOrder::m<KOrder::fold>()
+ctraits<Storage::fold>::Tm& KOrder::m<Storage::fold>()
{
return _fm;
}
template<>
-const ctraits<KOrder::fold>::Tm& KOrder::m<KOrder::fold>() const
+const ctraits<Storage::fold>::Tm& KOrder::m<Storage::fold>() const
{ return _fm;}
/* Here is the constructor of the |KOrder| class. We pass what we have
@@ -264,19 +257,19 @@ KOrder::KOrder(int num_stat, int num_pred, int num_both, int num_forw,
the container. We insert a new physical copies. */
auto tgy = std::make_unique<UGSTensor>(ny, TensorDimens(Symmetry{1, 0, 0, 0}, nvs));
tgy->getData() = gy.getData();
- insertDerivative<unfold>(std::move(tgy));
+ insertDerivative<Storage::unfold>(std::move(tgy));
auto tgu = std::make_unique<UGSTensor>(ny, TensorDimens(Symmetry{0, 1, 0, 0}, nvs));
tgu->getData() = gu.getData();
- insertDerivative<unfold>(std::move(tgu));
+ insertDerivative<Storage::unfold>(std::move(tgu));
// put $G_y$, $G_u$ and $G_{u'}$ to the container
/* Also note that since $g_\sigma$ is zero, so $G_\sigma$. */
- auto tGy = faaDiBrunoG<unfold>(Symmetry{1, 0, 0, 0});
- G<unfold>().insert(std::move(tGy));
- auto tGu = faaDiBrunoG<unfold>(Symmetry{0, 1, 0, 0});
- G<unfold>().insert(std::move(tGu));
- auto tGup = faaDiBrunoG<unfold>(Symmetry{0, 0, 1, 0});
- G<unfold>().insert(std::move(tGup));
+ auto tGy = faaDiBrunoG<Storage::unfold>(Symmetry{1, 0, 0, 0});
+ G<Storage::unfold>().insert(std::move(tGy));
+ auto tGu = faaDiBrunoG<Storage::unfold>(Symmetry{0, 1, 0, 0});
+ G<Storage::unfold>().insert(std::move(tGu));
+ auto tGup = faaDiBrunoG<Storage::unfold>(Symmetry{0, 0, 1, 0});
+ G<Storage::unfold>().insert(std::move(tGup));
}
// |KOrder::sylvesterSolve| unfolded specialization
@@ -295,7 +288,7 @@ KOrder::KOrder(int num_stat, int num_pred, int num_both, int num_forw,
template<>
void
-KOrder::sylvesterSolve<KOrder::unfold>(ctraits<unfold>::Ttensor &der) const
+KOrder::sylvesterSolve<Storage::unfold>(ctraits<Storage::unfold>::Ttensor &der) const
{
JournalRecordPair pa(journal);
pa << "Sylvester equation for dimension = " << der.getSym()[0] << endrec;
@@ -303,7 +296,7 @@ KOrder::sylvesterSolve<KOrder::unfold>(ctraits<unfold>::Ttensor &der) const
{
KORD_RAISE_IF(!der.isFinite(),
"RHS of Sylverster is not finite");
- TwoDMatrix gs_y(gs<unfold>().get(Symmetry{1, 0, 0, 0}));
+ TwoDMatrix gs_y(gs<Storage::unfold>().get(Symmetry{1, 0, 0, 0}));
GeneralSylvester sylv(der.getSym()[0], ny, ypart.nys(),
ypart.nstat+ypart.npred,
matA.getData(), matB.getData(),
@@ -323,11 +316,11 @@ KOrder::sylvesterSolve<KOrder::unfold>(ctraits<unfold>::Ttensor &der) const
template<>
void
-KOrder::sylvesterSolve<KOrder::fold>(ctraits<fold>::Ttensor &der) const
+KOrder::sylvesterSolve<Storage::fold>(ctraits<Storage::fold>::Ttensor &der) const
{
- ctraits<unfold>::Ttensor tmp(der);
- sylvesterSolve<unfold>(tmp);
- ctraits<fold>::Ttensor ftmp(tmp);
+ ctraits<Storage::unfold>::Ttensor tmp(der);
+ sylvesterSolve<Storage::unfold>(tmp);
+ ctraits<Storage::fold>::Ttensor ftmp(tmp);
der.getData() = (const Vector &) (ftmp.getData());
}
@@ -337,28 +330,28 @@ KOrder::switchToFolded()
JournalRecordPair pa(journal);
pa << "Switching from unfolded to folded" << endrec;
- int maxdim = g<unfold>().getMaxDim();
+ int maxdim = g<Storage::unfold>().getMaxDim();
for (int dim = 1; dim <= maxdim; dim++)
for (auto &si : SymmetrySet(dim, 4))
{
- if (si[2] == 0 && g<unfold>().check(si))
+ if (si[2] == 0 && g<Storage::unfold>().check(si))
{
- auto ft = std::make_unique<FGSTensor>(g<unfold>().get(si));
- insertDerivative<fold>(std::move(ft));
+ auto ft = std::make_unique<FGSTensor>(g<Storage::unfold>().get(si));
+ insertDerivative<Storage::fold>(std::move(ft));
if (dim > 1)
{
- gss<unfold>().remove(si);
- gs<unfold>().remove(si);
- g<unfold>().remove(si);
+ gss<Storage::unfold>().remove(si);
+ gs<Storage::unfold>().remove(si);
+ g<Storage::unfold>().remove(si);
}
}
- if (G<unfold>().check(si))
+ if (G<Storage::unfold>().check(si))
{
- auto ft = std::make_unique<FGSTensor>(G<unfold>().get(si));
- G<fold>().insert(std::move(ft));
+ auto ft = std::make_unique<FGSTensor>(G<Storage::unfold>().get(si));
+ G<Storage::fold>().insert(std::move(ft));
if (dim > 1)
{
- G<fold>().remove(si);
+ G<Storage::fold>().remove(si);
}
}
}
diff --git a/dynare++/kord/korder.hh b/dynare++/kord/korder.hh
index 11cb6ebe9706eb496a4639d89fd929b781a82394..da1e122148cc9962bc53a1c203c4016e04fe31e4 100644
--- a/dynare++/kord/korder.hh
+++ b/dynare++/kord/korder.hh
@@ -44,50 +44,34 @@
#include <cmath>
#include <type_traits>
-#define TYPENAME typename
-
-#define _Ttensor TYPENAME ctraits<t>::Ttensor
-#define _Ttensym TYPENAME ctraits<t>::Ttensym
-#define _Tg TYPENAME ctraits<t>::Tg
-#define _Tgs TYPENAME ctraits<t>::Tgs
-#define _Tgss TYPENAME ctraits<t>::Tgss
-#define _TG TYPENAME ctraits<t>::TG
-#define _TZstack TYPENAME ctraits<t>::TZstack
-#define _TGstack TYPENAME ctraits<t>::TGstack
-#define _TZXstack TYPENAME ctraits<t>::TZXstack
-#define _TGXstack TYPENAME ctraits<t>::TGXstack
-#define __Tm TYPENAME ctraits<t>::Tm
-#define _Tpol TYPENAME ctraits<t>::Tpol
-
-/* Here we use a classical IF template, and in |ctraits| we define a
- number of types. We have a type for tensor |Ttensor|, and types for
- each pair of folded/unfolded containers used as a member in |KOrder|.
-
- Note that we have enumeration |fold| and |unfold|. These must have the
- same value as the same enumeration in |KOrder|. */
+// The enum class passed as template parameter for many data structures
+enum class Storage { fold, unfold };
+
+/* In |ctraits| we define a number of types. We have a type for tensor
+ |Ttensor|, and types for each pair of folded/unfolded containers used as a
+ member in |KOrder|. */
class FoldedZXContainer;
class UnfoldedZXContainer;
class FoldedGXContainer;
class UnfoldedGXContainer;
-template <int type>
+template <Storage type>
class ctraits
{
public:
- enum { fold, unfold };
- using Ttensor = std::conditional_t<type == fold, FGSTensor, UGSTensor>;
- using Ttensym = std::conditional_t<type == fold, FFSTensor, UFSTensor>;
- using Tg = std::conditional_t<type == fold, FGSContainer, UGSContainer>;
- using Tgs = std::conditional_t<type == fold, FGSContainer, UGSContainer>;
- using Tgss = std::conditional_t<type == fold, FGSContainer, UGSContainer>;
- using TG = std::conditional_t<type == fold, FGSContainer, UGSContainer>;
- using TZstack = std::conditional_t<type == fold, FoldedZContainer, UnfoldedZContainer>;
- using TGstack = std::conditional_t<type == fold, FoldedGContainer, UnfoldedGContainer>;
- using Tm = std::conditional_t<type == fold, FNormalMoments, UNormalMoments>;
- using Tpol = std::conditional_t<type == fold, FTensorPolynomial, UTensorPolynomial>;
- using TZXstack = std::conditional_t<type == fold, FoldedZXContainer, UnfoldedZXContainer>;
- using TGXstack = std::conditional_t<type == fold, FoldedGXContainer, UnfoldedGXContainer>;
+ using Ttensor = std::conditional_t<type == Storage::fold, FGSTensor, UGSTensor>;
+ using Ttensym = std::conditional_t<type == Storage::fold, FFSTensor, UFSTensor>;
+ using Tg = std::conditional_t<type == Storage::fold, FGSContainer, UGSContainer>;
+ using Tgs = std::conditional_t<type == Storage::fold, FGSContainer, UGSContainer>;
+ using Tgss = std::conditional_t<type == Storage::fold, FGSContainer, UGSContainer>;
+ using TG = std::conditional_t<type == Storage::fold, FGSContainer, UGSContainer>;
+ using TZstack = std::conditional_t<type == Storage::fold, FoldedZContainer, UnfoldedZContainer>;
+ using TGstack = std::conditional_t<type == Storage::fold, FoldedGContainer, UnfoldedGContainer>;
+ using Tm = std::conditional_t<type == Storage::fold, FNormalMoments, UNormalMoments>;
+ using Tpol = std::conditional_t<type == Storage::fold, FTensorPolynomial, UTensorPolynomial>;
+ using TZXstack = std::conditional_t<type == Storage::fold, FoldedZXContainer, UnfoldedZXContainer>;
+ using TGXstack = std::conditional_t<type == Storage::fold, FoldedGXContainer, UnfoldedGXContainer>;
};
/* The |PartitionY| class defines the partitioning of state variables
@@ -148,18 +132,13 @@ public:
PLUMatrix(int n)
: TwoDMatrix(n, n),
inv(nrows()*ncols()),
- ipiv(new lapack_int[nrows()])
- {
- }
- PLUMatrix(const PLUMatrix &plu);
- ~PLUMatrix() override
+ ipiv(nrows())
{
- delete [] ipiv;
}
void multInv(TwoDMatrix &m) const;
private:
Vector inv;
- lapack_int *ipiv;
+ std::vector<lapack_int> ipiv;
protected:
void calcPLU();
};
@@ -258,9 +237,7 @@ public:
Most of the code is templated, and template types are calculated in
|ctraits|. So all templated methods get a template argument |T|, which
- can be either |fold|, or |unfold|. To shorten a reference to a type
- calculated by |ctraits| for a particular |t|, we define the following
- macros.
+ can be either |Storage::fold|, or |Storage::unfold|.
*/
class KOrder
@@ -296,34 +273,40 @@ protected:
/* These are the declarations of the template functions accessing the
containers. */
- template<int t>
- _Tg&g();
- template<int t>
- const _Tg&g() const;
- template<int t>
- _Tgs&gs();
- template<int t>
- const _Tgs&gs() const;
- template<int t>
- _Tgss&gss();
- template<int t>
- const _Tgss&gss() const;
- template<int t>
- _TG&G();
- template<int t>
- const _TG&G() const;
- template<int t>
- _TZstack&Zstack();
- template<int t>
- const _TZstack&Zstack() const;
- template<int t>
- _TGstack&Gstack();
- template<int t>
- const _TGstack&Gstack() const;
- template<int t>
- __Tm&m();
- template<int t>
- const __Tm&m() const;
+ template<Storage t>
+ typename ctraits<t>::Tg &g();
+ template<Storage t>
+ const typename ctraits<t>::Tg &g() const;
+
+ template<Storage t>
+ typename ctraits<t>::Tgs &gs();
+ template<Storage t>
+ const typename ctraits<t>::Tgs &gs() const;
+
+ template<Storage t>
+ typename ctraits<t>::Tgss &gss();
+ template<Storage t>
+ const typename ctraits<t>::Tgss &gss() const;
+
+ template<Storage t>
+ typename ctraits<t>::TG &G();
+ template<Storage t>
+ const typename ctraits<t>::TG &G() const;
+
+ template<Storage t>
+ typename ctraits<t>::TZstack &Zstack();
+ template<Storage t>
+ const typename ctraits<t>::TZstack &Zstack() const;
+
+ template<Storage t>
+ typename ctraits<t>::TGstack &Gstack();
+ template<Storage t>
+ const typename ctraits<t>::TGstack &Gstack() const;
+
+ template<Storage t>
+ typename ctraits<t>::Tm &m();
+ template<Storage t>
+ const typename ctraits<t>::Tm &m() const;
Journal &journal;
public:
@@ -331,13 +314,12 @@ public:
const TensorContainer<FSSparseTensor> &fcont,
const TwoDMatrix &gy, const TwoDMatrix &gu, const TwoDMatrix &v,
Journal &jr);
- enum { fold, unfold };
- template <int t>
+ template <Storage t>
void performStep(int order);
- template <int t>
+ template <Storage t>
double check(int dim) const;
- template <int t>
- Vector *calcStochShift(int order, double sigma) const;
+ template <Storage t>
+ Vector calcStochShift(int order, double sigma) const;
void switchToFolded();
const PartitionY &
getPartY() const
@@ -357,59 +339,59 @@ public:
static bool
is_even(int i)
{
- return (i/2)*2 == i;
+ return i % 2 == 0;
}
protected:
- template <int t>
- void insertDerivative(std::unique_ptr<_Ttensor> der);
- template<int t>
- void sylvesterSolve(_Ttensor &der) const;
+ template <Storage t>
+ void insertDerivative(std::unique_ptr<typename ctraits<t>::Ttensor> der);
+ template<Storage t>
+ void sylvesterSolve(typename ctraits<t>::Ttensor &der) const;
- template <int t>
- std::unique_ptr<_Ttensor> faaDiBrunoZ(const Symmetry &sym) const;
- template <int t>
- std::unique_ptr<_Ttensor> faaDiBrunoG(const Symmetry &sym) const;
+ template <Storage t>
+ std::unique_ptr<typename ctraits<t>::Ttensor> faaDiBrunoZ(const Symmetry &sym) const;
+ template <Storage t>
+ std::unique_ptr<typename ctraits<t>::Ttensor> faaDiBrunoG(const Symmetry &sym) const;
- template<int t>
+ template<Storage t>
void recover_y(int i);
- template <int t>
+ template <Storage t>
void recover_yu(int i, int j);
- template <int t>
+ template <Storage t>
void recover_ys(int i, int j);
- template <int t>
+ template <Storage t>
void recover_yus(int i, int j, int k);
- template <int t>
+ template <Storage t>
void recover_s(int i);
- template<int t>
+ template<Storage t>
void fillG(int i, int j, int k);
- template <int t>
- _Ttensor *calcD_ijk(int i, int j, int k) const;
- template <int t>
- _Ttensor *calcD_ik(int i, int k) const;
- template <int t>
- _Ttensor *calcD_k(int k) const;
-
- template <int t>
- _Ttensor *calcE_ijk(int i, int j, int k) const;
- template <int t>
- _Ttensor *calcE_ik(int i, int k) const;
- template <int t>
- _Ttensor *calcE_k(int k) const;
+ template <Storage t>
+ typename ctraits<t>::Ttensor calcD_ijk(int i, int j, int k) const;
+ template <Storage t>
+ typename ctraits<t>::Ttensor calcD_ik(int i, int k) const;
+ template <Storage t>
+ typename ctraits<t>::Ttensor calcD_k(int k) const;
+
+ template <Storage t>
+ typename ctraits<t>::Ttensor calcE_ijk(int i, int j, int k) const;
+ template <Storage t>
+ typename ctraits<t>::Ttensor calcE_ik(int i, int k) const;
+ template <Storage t>
+ typename ctraits<t>::Ttensor calcE_k(int k) const;
};
/* Here we insert the result to the container. Along the insertion, we
also create subtensors and insert as well. */
-template <int t>
+template <Storage t>
void
-KOrder::insertDerivative(std::unique_ptr<_Ttensor> der)
+KOrder::insertDerivative(std::unique_ptr<typename ctraits<t>::Ttensor> der)
{
auto der_ptr = der.get();
g<t>().insert(std::move(der));
- gs<t>().insert(std::make_unique<_Ttensor>(ypart.nstat, ypart.nys(), *der_ptr));
- gss<t>().insert(std::make_unique<_Ttensor>(ypart.nstat+ypart.npred,
- ypart.nyss(), *der_ptr));
+ gs<t>().insert(std::make_unique<typename ctraits<t>::Ttensor>(ypart.nstat, ypart.nys(), *der_ptr));
+ gss<t>().insert(std::make_unique<typename ctraits<t>::Ttensor>(ypart.nstat+ypart.npred,
+ ypart.nyss(), *der_ptr));
}
/* Here we implement Faa Di Bruno formula
@@ -419,13 +401,13 @@ KOrder::insertDerivative(std::unique_ptr<_Ttensor> der)
where $s$ is a given outer symmetry and $k$ is the dimension of the
symmetry. */
-template <int t>
-std::unique_ptr<_Ttensor>
+template <Storage t>
+std::unique_ptr<typename ctraits<t>::Ttensor>
KOrder::faaDiBrunoZ(const Symmetry &sym) const
{
JournalRecordPair pa(journal);
pa << "Faa Di Bruno Z container for " << sym << endrec;
- auto res = std::make_unique<_Ttensor>(ny, TensorDimens(sym, nvs));
+ auto res = std::make_unique<typename ctraits<t>::Ttensor>(ny, TensorDimens(sym, nvs));
FaaDiBruno bruno(journal);
bruno.calculate(Zstack<t>(), f, *res);
return res;
@@ -434,14 +416,14 @@ KOrder::faaDiBrunoZ(const Symmetry &sym) const
/* The same as |@<|KOrder::faaDiBrunoZ| templated code@>|, but for
$g^{**}$ and $G$ stack. */
-template <int t>
-std::unique_ptr<_Ttensor>
+template <Storage t>
+std::unique_ptr<typename ctraits<t>::Ttensor>
KOrder::faaDiBrunoG(const Symmetry &sym) const
{
JournalRecordPair pa(journal);
pa << "Faa Di Bruno G container for " << sym << endrec;
TensorDimens tdims(sym, nvs);
- auto res = std::make_unique<_Ttensor>(ypart.nyss(), tdims);
+ auto res = std::make_unique<typename ctraits<t>::Ttensor>(ypart.nyss(), tdims);
FaaDiBruno bruno(journal);
bruno.calculate(Gstack<t>(), gss<t>(), *res);
return res;
@@ -459,7 +441,7 @@ KOrder::faaDiBrunoG(const Symmetry &sym) const
{\bf Provides:} $g_{y^i}$, and $G_{y^i}$ */
-template<int t>
+template<Storage t>
void
KOrder::recover_y(int i)
{
@@ -478,9 +460,9 @@ KOrder::recover_y(int i)
insertDerivative<t>(std::move(g_yi));
- _Ttensor &gss_y = gss<t>().get(Symmetry{1, 0, 0, 0});
+ auto &gss_y = gss<t>().get(Symmetry{1, 0, 0, 0});
gs<t>().multAndAdd(gss_y, *G_yi_ptr);
- _Ttensor &gss_yi = gss<t>().get(sym);
+ auto &gss_yi = gss<t>().get(sym);
gs<t>().multAndAdd(gss_yi, *G_yi_ptr);
}
@@ -495,7 +477,7 @@ KOrder::recover_y(int i)
{\bf Provides:} $g_{y^iu^j}$, and $G_{y^iu^j}$ */
-template <int t>
+template <Storage t>
void
KOrder::recover_yu(int i, int j)
{
@@ -535,7 +517,7 @@ KOrder::recover_yu(int i, int j)
$G_{y^iu'^m\sigma^{j-m}}$ for $m=1,\ldots,j$. The latter is calculated
by |fillG| before the actual calculation. */
-template <int t>
+template <Storage t>
void
KOrder::recover_ys(int i, int j)
{
@@ -554,16 +536,14 @@ KOrder::recover_ys(int i, int j)
auto g_yisj = faaDiBrunoZ<t>(sym);
{
- _Ttensor *D_ij = calcD_ik<t>(i, j);
- g_yisj->add(1.0, *D_ij);
- delete D_ij;
+ auto D_ij = calcD_ik<t>(i, j);
+ g_yisj->add(1.0, D_ij);
}
if (j >= 3)
{
- _Ttensor *E_ij = calcE_ik<t>(i, j);
- g_yisj->add(1.0, *E_ij);
- delete E_ij;
+ auto E_ij = calcE_ik<t>(i, j);
+ g_yisj->add(1.0, E_ij);
}
g_yisj->mult(-1.0);
@@ -596,7 +576,7 @@ KOrder::recover_ys(int i, int j)
{\bf Provides:} $g_{y^iu^j\sigma^k}$, $G_{y^iu^j\sigma^k}$, and
$G_{y^iu^ju'^m\sigma^{k-m}}$ for $m=1,\ldots, k$ */
-template <int t>
+template <Storage t>
void
KOrder::recover_yus(int i, int j, int k)
{
@@ -615,16 +595,14 @@ KOrder::recover_yus(int i, int j, int k)
auto g_yiujsk = faaDiBrunoZ<t>(sym);
{
- _Ttensor *D_ijk = calcD_ijk<t>(i, j, k);
- g_yiujsk->add(1.0, *D_ijk);
- delete D_ijk;
+ auto D_ijk = calcD_ijk<t>(i, j, k);
+ g_yiujsk->add(1.0, D_ijk);
}
if (k >= 3)
{
- _Ttensor *E_ijk = calcE_ijk<t>(i, j, k);
- g_yiujsk->add(1.0, *E_ijk);
- delete E_ijk;
+ auto E_ijk = calcE_ijk<t>(i, j, k);
+ g_yiujsk->add(1.0, E_ijk);
}
g_yiujsk->mult(-1.0);
@@ -664,7 +642,7 @@ KOrder::recover_yus(int i, int j, int k)
{\bf Provides:} $g_{\sigma^i}$, $G_{\sigma^i}$, and
$G_{u'^m\sigma^{i-m}}$ for $m=1,\ldots,i$ */
-template <int t>
+template <Storage t>
void
KOrder::recover_s(int i)
{
@@ -683,16 +661,14 @@ KOrder::recover_s(int i)
auto g_si = faaDiBrunoZ<t>(sym);
{
- _Ttensor *D_i = calcD_k<t>(i);
- g_si->add(1.0, *D_i);
- delete D_i;
+ auto D_i = calcD_k<t>(i);
+ g_si->add(1.0, D_i);
}
if (i >= 3)
{
- _Ttensor *E_i = calcE_k<t>(i);
- g_si->add(1.0, *E_i);
- delete E_i;
+ auto E_i = calcE_k<t>(i);
+ g_si->add(1.0, E_i);
}
g_si->mult(-1.0);
@@ -709,18 +685,16 @@ KOrder::recover_s(int i)
$m=1,\ldots, k$. The derivatives are inserted only for $k-m$ being
even. */
-template<int t>
+template<Storage t>
void
KOrder::fillG(int i, int j, int k)
{
for (int m = 1; m <= k; m++)
- {
- if (is_even(k-m))
- {
- auto G_yiujupms = faaDiBrunoG<t>(Symmetry{i, j, m, k-m});
- G<t>().insert(std::move(G_yiujupms));
- }
- }
+ if (is_even(k-m))
+ {
+ auto G_yiujupms = faaDiBrunoG<t>(Symmetry{i, j, m, k-m});
+ G<t>().insert(std::move(G_yiujupms));
+ }
}
/* Here we calculate
@@ -730,16 +704,16 @@ KOrder::fillG(int i, int j, int k)
\left[\Sigma\right]^{\gamma_1\ldots\gamma_k}$$
So it is non zero only for even $k$. */
-template <int t>
-_Ttensor *
+template <Storage t>
+typename ctraits<t>::Ttensor
KOrder::calcD_ijk(int i, int j, int k) const
{
- _Ttensor *res = new _Ttensor(ny, TensorDimens(Symmetry{i, j, 0, 0}, nvs));
- res->zeros();
+ typename ctraits<t>::Ttensor res(ny, TensorDimens(Symmetry{i, j, 0, 0}, nvs));
+ res.zeros();
if (is_even(k))
{
auto tmp = faaDiBrunoZ<t>(Symmetry{i, j, k, 0});
- tmp->contractAndAdd(2, *res, m<t>().get(Symmetry{k}));
+ tmp->contractAndAdd(2, res, m<t>().get(Symmetry{k}));
}
return res;
}
@@ -751,44 +725,44 @@ KOrder::calcD_ijk(int i, int j, int k) const
\left[\Sigma\right]^{\gamma_1\ldots\gamma_m}$$
The sum can sum only for even $m$. */
-template <int t>
-_Ttensor *
+template <Storage t>
+typename ctraits<t>::Ttensor
KOrder::calcE_ijk(int i, int j, int k) const
{
- _Ttensor *res = new _Ttensor(ny, TensorDimens(Symmetry{i, j, 0, 0}, nvs));
- res->zeros();
+ typename ctraits<t>::Ttensor res(ny, TensorDimens(Symmetry{i, j, 0, 0}, nvs));
+ res.zeros();
for (int n = 2; n <= k-1; n += 2)
{
auto tmp = faaDiBrunoZ<t>(Symmetry{i, j, n, k-n});
- tmp->mult((double) (PascalTriangle::noverk(k, n)));
- tmp->contractAndAdd(2, *res, m<t>().get(Symmetry{n}));
+ tmp->mult(static_cast<double>(PascalTriangle::noverk(k, n)));
+ tmp->contractAndAdd(2, res, m<t>().get(Symmetry{n}));
}
return res;
}
-template <int t>
-_Ttensor *
+template <Storage t>
+typename ctraits<t>::Ttensor
KOrder::calcD_ik(int i, int k) const
{
return calcD_ijk<t>(i, 0, k);
}
-template <int t>
-_Ttensor *
+template <Storage t>
+typename ctraits<t>::Ttensor
KOrder::calcD_k(int k) const
{
return calcD_ijk<t>(0, 0, k);
}
-template <int t>
-_Ttensor *
+template <Storage t>
+typename ctraits<t>::Ttensor
KOrder::calcE_ik(int i, int k) const
{
return calcE_ijk<t>(i, 0, k);
}
-template <int t>
-_Ttensor *
+template <Storage t>
+typename ctraits<t>::Ttensor
KOrder::calcE_k(int k) const
{
return calcE_ijk<t>(0, 0, k);
@@ -816,7 +790,7 @@ KOrder::calcE_k(int k) const
through all the recovering methods, he should find out that also all
$G$ are provided. */
-template <int t>
+template <Storage t>
void
KOrder::performStep(int order)
{
@@ -828,23 +802,18 @@ KOrder::performStep(int order)
recover_y<t>(order);
for (int i = 0; i < order; i++)
- {
- recover_yu<t>(i, order-i);
- }
+ recover_yu<t>(i, order-i);
for (int j = 1; j < order; j++)
{
for (int i = j-1; i >= 1; i--)
- {
- recover_yus<t>(order-j, i, j-i);
- }
+ recover_yus<t>(order-j, i, j-i);
recover_ys<t>(order-j, j);
}
for (int i = order-1; i >= 1; i--)
- {
- recover_yus<t>(0, i, order-i);
- }
+ recover_yus<t>(0, i, order-i);
+
recover_s<t>(order);
}
@@ -852,7 +821,7 @@ KOrder::performStep(int order)
order. The method returns the largest residual size. Each check simply
evaluates the equation. */
-template <int t>
+template <Storage t>
double
KOrder::check(int dim) const
{
@@ -884,12 +853,10 @@ KOrder::check(int dim) const
{
Symmetry sym{i, j, 0, k};
auto r = faaDiBrunoZ<t>(sym);
- _Ttensor *D_ijk = calcD_ijk<t>(i, j, k);
- r->add(1.0, *D_ijk);
- delete D_ijk;
- _Ttensor *E_ijk = calcE_ijk<t>(i, j, k);
- r->add(1.0, *E_ijk);
- delete E_ijk;
+ auto D_ijk = calcD_ijk<t>(i, j, k);
+ r->add(1.0, D_ijk);
+ auto E_ijk = calcE_ijk<t>(i, j, k);
+ r->add(1.0, E_ijk);
double err = r->getData().getMax();
JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
}
@@ -897,12 +864,10 @@ KOrder::check(int dim) const
// check for $F_{\sigma^i}+D_i+E_i=0
auto r = faaDiBrunoZ<t>(Symmetry{0, 0, 0, dim});
- _Ttensor *D_k = calcD_k<t>(dim);
- r->add(1.0, *D_k);
- delete D_k;
- _Ttensor *E_k = calcE_k<t>(dim);
- r->add(1.0, *E_k);
- delete E_k;
+ auto D_k = calcD_k<t>(dim);
+ r->add(1.0, D_k);
+ auto E_k = calcE_k<t>(dim);
+ r->add(1.0, E_k);
double err = r->getData().getMax();
Symmetry sym{0, 0, 0, dim};
JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
@@ -912,19 +877,18 @@ KOrder::check(int dim) const
return maxerror;
}
-template <int t>
-Vector *
+template <Storage t>
+Vector
KOrder::calcStochShift(int order, double sigma) const
{
- auto *res = new Vector(ny);
- res->zeros();
+ Vector res(ny);
+ res.zeros();
int jfac = 1;
for (int j = 1; j <= order; j++, jfac *= j)
if (is_even(j))
{
- _Ttensor *ten = calcD_k<t>(j);
- res->add(std::pow(sigma, j)/jfac, ten->getData());
- delete ten;
+ auto ten = calcD_k<t>(j);
+ res.add(std::pow(sigma, j)/jfac, ten.getData());
}
return res;
}
diff --git a/dynare++/kord/korder_stoch.cc b/dynare++/kord/korder_stoch.cc
index 924f797a7cb811793f08d045750fc30566a7a0f4..1bdc404b7756879d71885fcd991e4c1e4cb29bfe 100644
--- a/dynare++/kord/korder_stoch.cc
+++ b/dynare++/kord/korder_stoch.cc
@@ -11,7 +11,7 @@ MatrixAA::MatrixAA(const FSSparseTensor &f, const IntSequence &ss,
{
zeros();
- IntSequence c(1); c[0] = 1;
+ IntSequence c{1};
FGSTensor f_y(f, ss, c, TensorDimens(ss, c));
add(1.0, f_y);
@@ -59,88 +59,112 @@ KOrderStoch::KOrderStoch(const PartitionY &yp, int nu,
// |KOrderStoch| convenience method specializations
template<>
-ctraits<KOrder::unfold>::Tg& KOrderStoch::g<KOrder::unfold>()
+ctraits<Storage::unfold>::Tg &KOrderStoch::g<Storage::unfold>()
{
return _ug;
}
template<>
-const ctraits<KOrder::unfold>::Tg& KOrderStoch::g<KOrder::unfold>() const
-{ return _ug;}
+const ctraits<Storage::unfold>::Tg &KOrderStoch::g<Storage::unfold>() const
+{
+ return _ug;
+}
template<>
-ctraits<KOrder::fold>::Tg& KOrderStoch::g<KOrder::fold>()
+ctraits<Storage::fold>::Tg &KOrderStoch::g<Storage::fold>()
{
return _fg;
}
template<>
-const ctraits<KOrder::fold>::Tg& KOrderStoch::g<KOrder::fold>() const
-{ return _fg;}
+const ctraits<Storage::fold>::Tg &KOrderStoch::g<Storage::fold>() const
+{
+ return _fg;
+}
template<>
-ctraits<KOrder::unfold>::Tgs& KOrderStoch::gs<KOrder::unfold>()
+ctraits<Storage::unfold>::Tgs &KOrderStoch::gs<Storage::unfold>()
{
return _ugs;
}
template<>
-const ctraits<KOrder::unfold>::Tgs& KOrderStoch::gs<KOrder::unfold>() const
-{ return _ugs;}
+const ctraits<Storage::unfold>::Tgs &KOrderStoch::gs<Storage::unfold>() const
+{
+ return _ugs;
+}
template<>
-ctraits<KOrder::fold>::Tgs& KOrderStoch::gs<KOrder::fold>()
+ctraits<Storage::fold>::Tgs &KOrderStoch::gs<Storage::fold>()
{
return _fgs;
}
template<>
-const ctraits<KOrder::fold>::Tgs& KOrderStoch::gs<KOrder::fold>() const
-{ return _fgs;}
+const ctraits<Storage::fold>::Tgs &KOrderStoch::gs<Storage::fold>() const
+{
+ return _fgs;
+}
template<>
-const ctraits<KOrder::unfold>::Tgss& KOrderStoch::h<KOrder::unfold>() const
-{ return *_uh;}
+const ctraits<Storage::unfold>::Tgss &KOrderStoch::h<Storage::unfold>() const
+{
+ return *_uh;
+}
template<>
-const ctraits<KOrder::fold>::Tgss& KOrderStoch::h<KOrder::fold>() const
-{ return *_fh;}
+const ctraits<Storage::fold>::Tgss &KOrderStoch::h<Storage::fold>() const
+{
+ return *_fh;
+}
template<>
-ctraits<KOrder::unfold>::TG& KOrderStoch::G<KOrder::unfold>()
+ctraits<Storage::unfold>::TG &KOrderStoch::G<Storage::unfold>()
{
return _uG;
}
template<>
-const ctraits<KOrder::unfold>::TG& KOrderStoch::G<KOrder::unfold>() const
-{ return _uG;}
+const ctraits<Storage::unfold>::TG &KOrderStoch::G<Storage::unfold>() const
+{
+ return _uG;
+}
template<>
-ctraits<KOrder::fold>::TG& KOrderStoch::G<KOrder::fold>()
+ctraits<Storage::fold>::TG &KOrderStoch::G<Storage::fold>()
{
return _fG;
}
template<>
-const ctraits<KOrder::fold>::TG& KOrderStoch::G<KOrder::fold>() const
-{ return _fG;}
+const ctraits<Storage::fold>::TG& KOrderStoch::G<Storage::fold>() const
+{
+ return _fG;
+}
template<>
-ctraits<KOrder::unfold>::TZXstack& KOrderStoch::Zstack<KOrder::unfold>()
+ctraits<Storage::unfold>::TZXstack &KOrderStoch::Zstack<Storage::unfold>()
{
return _uZstack;
}
template<>
-const ctraits<KOrder::unfold>::TZXstack& KOrderStoch::Zstack<KOrder::unfold>() const
-{ return _uZstack;}
+const ctraits<Storage::unfold>::TZXstack &KOrderStoch::Zstack<Storage::unfold>() const
+{
+ return _uZstack;
+}
template<>
-ctraits<KOrder::fold>::TZXstack& KOrderStoch::Zstack<KOrder::fold>()
+ctraits<Storage::fold>::TZXstack &KOrderStoch::Zstack<Storage::fold>()
{
return _fZstack;
}
template<>
-const ctraits<KOrder::fold>::TZXstack& KOrderStoch::Zstack<KOrder::fold>() const
-{ return _fZstack;}
+const ctraits<Storage::fold>::TZXstack &KOrderStoch::Zstack<Storage::fold>() const
+{
+ return _fZstack;
+}
template<>
-ctraits<KOrder::unfold>::TGXstack& KOrderStoch::Gstack<KOrder::unfold>()
+ctraits<Storage::unfold>::TGXstack &KOrderStoch::Gstack<Storage::unfold>()
{
return _uGstack;
}
template<>
-const ctraits<KOrder::unfold>::TGXstack& KOrderStoch::Gstack<KOrder::unfold>() const
-{ return _uGstack;}
+const ctraits<Storage::unfold>::TGXstack &KOrderStoch::Gstack<Storage::unfold>() const
+{
+ return _uGstack;
+}
template<>
-ctraits<KOrder::fold>::TGXstack& KOrderStoch::Gstack<KOrder::fold>()
+ctraits<Storage::fold>::TGXstack &KOrderStoch::Gstack<Storage::fold>()
{
return _fGstack;
}
template<>
-const ctraits<KOrder::fold>::TGXstack& KOrderStoch::Gstack<KOrder::fold>() const
-{ return _fGstack;}
+const ctraits<Storage::fold>::TGXstack &KOrderStoch::Gstack<Storage::fold>() const
+{
+ return _fGstack;
+}
diff --git a/dynare++/kord/korder_stoch.hh b/dynare++/kord/korder_stoch.hh
index 38f0cc040937bba9fc60419b7cfaeb2388b34389..2adf37cd0a2f8b57486c116779ab39936c74c42c 100644
--- a/dynare++/kord/korder_stoch.hh
+++ b/dynare++/kord/korder_stoch.hh
@@ -32,12 +32,12 @@
/* This class is a container, which has a specialized constructor
integrating the policy rule at given $\sigma$. */
-template <int t>
+template <Storage t>
class IntegDerivs : public ctraits<t>::Tgss
{
public:
- IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const __Tm &mom,
- double at_sigma);
+ IntegDerivs(int r, const IntSequence &nvs, const typename ctraits<t>::Tgss &g,
+ const typename ctraits<t>::Tm &mom, double at_sigma);
};
/* This constructor integrates a rule (namely its $g^{**}$ part) with
@@ -79,9 +79,9 @@ public:
\Sigma^{m+n}{\tilde\sigma}^m$$
and this is exactly what the code does. */
-template <int t>
-IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const __Tm &mom,
- double at_sigma)
+template <Storage t>
+IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const typename ctraits<t>::Tgss &g,
+ const typename ctraits<t>::Tm &mom, double at_sigma)
: ctraits<t>::Tgss(4)
{
int maxd = g.getMaxDim();
@@ -91,7 +91,7 @@ IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const
{
int p = d-i;
Symmetry sym{i, 0, 0, p};
- auto ten = std::make_unique<_Ttensor>(r, TensorDimens(sym, nvs));
+ auto ten = std::make_unique<typename ctraits<t>::Ttensor>(r, TensorDimens(sym, nvs));
// calculate derivative $h_{y^i\sigma^p}$
/* This code calculates
@@ -113,7 +113,7 @@ IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const
ten->add(mult, g.get(sym_mn));
if (m+n > 0 && KOrder::is_even(m+n) && g.check(sym_mn))
{
- _Ttensor gtmp(g.get(sym_mn));
+ typename ctraits<t>::Ttensor gtmp(g.get(sym_mn));
gtmp.mult(mult);
gtmp.contractAndAdd(1, *ten, mom.get(Symmetry{m+n}));
}
@@ -133,12 +133,12 @@ IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const
y^*,\bar\sigma)$. The derivatives are extrapolated based on
derivatives at $(\tilde y^*,\tilde\sigma)$. */
-template <int t>
+template <Storage t>
class StochForwardDerivs : public ctraits<t>::Tgss
{
public:
StochForwardDerivs(const PartitionY &ypart, int nu,
- const _Tgss &g, const __Tm &m,
+ const typename ctraits<t>::Tgss &g, const typename ctraits<t>::Tm &m,
const Vector &ydelta, double sdelta,
double at_sigma);
};
@@ -160,9 +160,10 @@ public:
\li Recover general symmetry tensors from the (full symmetric) polynomial
\endorderedlist */
-template <int t>
+template <Storage t>
StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
- const _Tgss &g, const __Tm &m,
+ const typename ctraits<t>::Tgss &g,
+ const typename ctraits<t>::Tm &m,
const Vector &ydelta, double sdelta,
double at_sigma)
: ctraits<t>::Tgss(4)
@@ -181,10 +182,10 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
// make |g_int_sym| be full symmetric polynomial from |g_int|
/* Here we just form a polynomial whose unique variable corresponds to
$\left[\matrix{y^*\cr\sigma}\right]$ stack. */
- _Tpol g_int_sym(r, ypart.nys()+1);
+ typename ctraits<t>::Tpol g_int_sym(r, ypart.nys()+1);
for (int d = 1; d <= maxd; d++)
{
- auto ten = std::make_unique<_Ttensym>(r, ypart.nys()+1, d);
+ auto ten = std::make_unique<typename ctraits<t>::Ttensym>(r, ypart.nys()+1, d);
ten->zeros();
for (int i = 0; i <= d; i++)
{
@@ -207,7 +208,7 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
ConstVector dy_in(ydelta, ypart.nstat, ypart.nys());
dy = dy_in;
delta[ypart.nys()] = sdelta;
- _Tpol g_int_cent(r, ypart.nys()+1);
+ typename ctraits<t>::Tpol g_int_cent(r, ypart.nys()+1);
for (int d = 1; d <= maxd; d++)
{
g_int_sym.derivative(d-1);
@@ -224,19 +225,16 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
true_nvs[1] = nu;
true_nvs[2] = nu;
for (int d = 1; d <= maxd; d++)
- {
- if (g_int_cent.check(Symmetry{d}))
+ if (g_int_cent.check(Symmetry{d}))
+ for (int i = 0; i <= d; i++)
{
- for (int i = 0; i <= d; i++)
- {
- Symmetry sym{i, 0, 0, d-i};
- IntSequence coor(pp.unfold(sym));
- auto ten = std::make_unique<_Ttensor>(g_int_cent.get(Symmetry{d}), ss, coor,
- TensorDimens(sym, true_nvs));
- this->insert(std::move(ten));
- }
+ Symmetry sym{i, 0, 0, d-i};
+ IntSequence coor(pp.unfold(sym));
+ auto ten = std::make_unique<typename ctraits<t>::Ttensor>(g_int_cent.get(Symmetry{d}),
+ ss, coor,
+ TensorDimens(sym, true_nvs));
+ this->insert(std::move(ten));
}
- }
}
/* This container corresponds to $h(g^*(y,u,\sigma),\sigma)$. Note that
@@ -425,7 +423,7 @@ public:
const FGSContainer &hh, Journal &jr);
KOrderStoch(const PartitionY &ypart, int nu, const TensorContainer<FSSparseTensor> &fcont,
const UGSContainer &hh, Journal &jr);
- template <int t>
+ template <Storage t>
void performStep(int order);
const FGSContainer &
getFoldDers() const
@@ -438,46 +436,51 @@ public:
return _ug;
}
protected:
- template <int t>
- std::unique_ptr<_Ttensor> faaDiBrunoZ(const Symmetry &sym) const;
- template <int t>
- std::unique_ptr<_Ttensor> faaDiBrunoG(const Symmetry &sym) const;
+ template <Storage t>
+ std::unique_ptr<typename ctraits<t>::Ttensor> faaDiBrunoZ(const Symmetry &sym) const;
+ template <Storage t>
+ std::unique_ptr<typename ctraits<t>::Ttensor> faaDiBrunoG(const Symmetry &sym) const;
// convenience access methods
- template<int t>
- _Tg&g();
- template<int t>
- const _Tg&g() const;
- template<int t>
- _Tgs&gs();
- template<int t>
- const _Tgs&gs() const;
- template<int t>
- const _Tgss&h() const;
- template<int t>
- _TG&G();
- template<int t>
- const _TG&G() const;
- template<int t>
- _TZXstack&Zstack();
- template<int t>
- const _TZXstack&Zstack() const;
- template<int t>
- _TGXstack&Gstack();
- template<int t>
- const _TGXstack&Gstack() const;
+ template<Storage t>
+ typename ctraits<t>::Tg &g();
+ template<Storage t>
+ const typename ctraits<t>::Tg &g() const;
+
+ template<Storage t>
+ typename ctraits<t>::Tgs &gs();
+ template<Storage t>
+ const typename ctraits<t>::Tgs &gs() const;
+
+ template<Storage t>
+ const typename ctraits<t>::Tgss &h() const;
+
+ template<Storage t>
+ typename ctraits<t>::TG &G();
+ template<Storage t>
+ const typename ctraits<t>::TG &G() const;
+
+ template<Storage t>
+ typename ctraits<t>::TZXstack &Zstack();
+ template<Storage t>
+ const typename ctraits<t>::TZXstack &Zstack() const;
+
+ template<Storage t>
+ typename ctraits<t>::TGXstack &Gstack();
+ template<Storage t>
+ const typename ctraits<t>::TGXstack &Gstack() const;
};
/* This calculates a derivative of $f(G(y,u,\sigma),g(y,u,\sigma),y,u)$
of a given symmetry. */
-template <int t>
-std::unique_ptr<_Ttensor>
+template <Storage t>
+std::unique_ptr<typename ctraits<t>::Ttensor>
KOrderStoch::faaDiBrunoZ(const Symmetry &sym) const
{
JournalRecordPair pa(journal);
pa << "Faa Di Bruno ZX container for " << sym << endrec;
- auto res = std::make_unique<_Ttensor>(ypart.ny(), TensorDimens(sym, nvs));
+ auto res = std::make_unique<typename ctraits<t>::Ttensor>(ypart.ny(), TensorDimens(sym, nvs));
FaaDiBruno bruno(journal);
bruno.calculate(Zstack<t>(), f, *res);
return res;
@@ -486,14 +489,14 @@ KOrderStoch::faaDiBrunoZ(const Symmetry &sym) const
/* This calculates a derivative of
$G(y,u,\sigma)=h(g^*(y,u,\sigma),\sigma)$ of a given symmetry. */
-template <int t>
-std::unique_ptr<_Ttensor>
+template <Storage t>
+std::unique_ptr<typename ctraits<t>::Ttensor>
KOrderStoch::faaDiBrunoG(const Symmetry &sym) const
{
JournalRecordPair pa(journal);
pa << "Faa Di Bruno GX container for " << sym << endrec;
TensorDimens tdims(sym, nvs);
- auto res = std::make_unique<_Ttensor>(ypart.nyss(), tdims);
+ auto res = std::make_unique<typename ctraits<t>::Ttensor>(ypart.nyss(), tdims);
FaaDiBruno bruno(journal);
bruno.calculate(Gstack<t>(), h<t>(), *res);
return res;
@@ -510,7 +513,7 @@ KOrderStoch::faaDiBrunoG(const Symmetry &sym) const
$$g_s=-matA^{-1}\cdot RHS.$$ Finally we have to update $G_s$ by
calling |Gstack<t>().multAndAdd(1, h<t>(), *G_sym)|. */
-template <int t>
+template <Storage t>
void
KOrderStoch::performStep(int order)
{
@@ -518,24 +521,23 @@ KOrderStoch::performStep(int order)
KORD_RAISE_IF(order-1 != maxd && (order != 1 || maxd != -1),
"Wrong order for KOrderStoch::performStep");
for (auto &si : SymmetrySet(order, 4))
- {
- if (si[2] == 0)
- {
- JournalRecordPair pa(journal);
- pa << "Recovering symmetry " << si << endrec;
-
- auto G_sym = faaDiBrunoG<t>(si);
- auto G_sym_ptr = G_sym.get();
- G<t>().insert(std::move(G_sym));
-
- auto g_sym = faaDiBrunoZ<t>(si);
- auto g_sym_ptr = g_sym.get();
- g_sym->mult(-1.0);
- matA.multInv(*g_sym);
- g<t>().insert(std::move(g_sym));
- gs<t>().insert(std::make_unique<_Ttensor>(ypart.nstat, ypart.nys(), *g_sym_ptr));
-
- Gstack<t>().multAndAdd(1, h<t>(), *G_sym_ptr);
- }
- }
+ if (si[2] == 0)
+ {
+ JournalRecordPair pa(journal);
+ pa << "Recovering symmetry " << si << endrec;
+
+ auto G_sym = faaDiBrunoG<t>(si);
+ auto G_sym_ptr = G_sym.get();
+ G<t>().insert(std::move(G_sym));
+
+ auto g_sym = faaDiBrunoZ<t>(si);
+ auto g_sym_ptr = g_sym.get();
+ g_sym->mult(-1.0);
+ matA.multInv(*g_sym);
+ g<t>().insert(std::move(g_sym));
+ gs<t>().insert(std::make_unique<typename ctraits<t>::Ttensor>(ypart.nstat, ypart.nys(),
+ *g_sym_ptr));
+
+ Gstack<t>().multAndAdd(1, h<t>(), *G_sym_ptr);
+ }
}
diff --git a/dynare++/kord/tests.cc b/dynare++/kord/tests.cc
index 8a626733aa051e9db71b0661200fb204ad84943f..26db77f44d352b1ad68c4dbb50a46ad2954f8913 100644
--- a/dynare++/kord/tests.cc
+++ b/dynare++/kord/tests.cc
@@ -298,12 +298,12 @@ TestRunnable::korder_unfold_fold(int maxdim, int unfold_dim,
for (int d = 2; d <= unfold_dim; d++)
{
clock_t pertime = clock();
- kord.performStep<KOrder::unfold>(d);
+ kord.performStep<Storage::unfold>(d);
pertime = clock()-pertime;
std::cout << "\ttime for unfolded step dim=" << d << ": " << std::setprecision(4)
<< static_cast<double>(pertime)/CLOCKS_PER_SEC << std::endl;
clock_t checktime = clock();
- double err = kord.check<KOrder::unfold>(d);
+ double err = kord.check<Storage::unfold>(d);
checktime = clock()-checktime;
std::cout << "\ttime for step check dim=" << d << ": " << std::setprecision(4)
<< static_cast<double>(checktime)/CLOCKS_PER_SEC << '\n'
@@ -324,12 +324,12 @@ TestRunnable::korder_unfold_fold(int maxdim, int unfold_dim,
for (int d = unfold_dim+1; d <= maxdim; d++)
{
clock_t pertime = clock();
- kord.performStep<KOrder::fold>(d);
+ kord.performStep<Storage::fold>(d);
pertime = clock()-pertime;
std::cout << "\ttime for folded step dim=" << d << ": " << std::setprecision(4)
<< static_cast<double>(pertime)/CLOCKS_PER_SEC << std::endl;
clock_t checktime = clock();
- double err = kord.check<KOrder::fold>(d);
+ double err = kord.check<Storage::fold>(d);
checktime = clock()-checktime;
std::cout << "\ttime for step check dim=" << d << ": " << std::setprecision(4)
<< static_cast<double>(checktime)/CLOCKS_PER_SEC << '\n'
diff --git a/dynare++/src/dynare3.cc b/dynare++/src/dynare3.cc
index 973e36f0e7a0a1194675364c6032d246ccc19ecd..36f3f65a5b3661e44ef19d07bf88ce5c01fb413f 100644
--- a/dynare++/src/dynare3.cc
+++ b/dynare++/src/dynare3.cc
@@ -20,13 +20,13 @@
/* DynareNameList class */
/**************************************************************************************/
std::vector<int>
-DynareNameList::selectIndices(const std::vector<const char *> &ns) const
+DynareNameList::selectIndices(const std::vector<std::string> &ns) const
{
std::vector<int> res;
- for (auto n : ns)
+ for (const auto &n : ns)
{
int j = 0;
- while (j < getNum() && strcmp(getName(j), n) != 0)
+ while (j < getNum() && n != getName(j))
j++;
if (j == getNum())
throw DynareException(__FILE__, __LINE__,
diff --git a/dynare++/src/dynare3.hh b/dynare++/src/dynare3.hh
index 6a7fa48463f751b870d6add4b049297154807d2a..f270010af682eb1c7be0c81144246b7b35de740b 100644
--- a/dynare++/src/dynare3.hh
+++ b/dynare++/src/dynare3.hh
@@ -20,15 +20,15 @@ class Dynare;
class DynareNameList : public NameList
{
- std::vector<const char *> names;
+ std::vector<std::string> names;
public:
DynareNameList(const Dynare &dynare);
int
getNum() const override
{
- return (int) names.size();
+ return static_cast<int>(names.size());
}
- const char *
+ const std::string &
getName(int i) const override
{
return names[i];
@@ -36,20 +36,20 @@ public:
/** This for each string of the input vector calculates its index
* in the names. And returns the resulting vector of indices. If
* the name cannot be found, then an exception is raised. */
- std::vector<int> selectIndices(const std::vector<const char *> &ns) const;
+ std::vector<int> selectIndices(const std::vector<std::string> &ns) const;
};
class DynareExogNameList : public NameList
{
- std::vector<const char *> names;
+ std::vector<std::string> names;
public:
DynareExogNameList(const Dynare &dynare);
int
getNum() const override
{
- return (int) names.size();
+ return static_cast<int>(names.size());
}
- const char *
+ const std::string &
getName(int i) const override
{
return names[i];
@@ -58,16 +58,16 @@ public:
class DynareStateNameList : public NameList
{
- std::vector<const char *> names;
+ std::vector<std::string> names;
public:
DynareStateNameList(const Dynare &dynare, const DynareNameList &dnl,
const DynareExogNameList &denl);
int
getNum() const override
{
- return (int) names.size();
+ return static_cast<int>(names.size());
}
- const char *
+ const std::string &
getName(int i) const override
{
return names[i];
@@ -105,10 +105,10 @@ public:
double sstol, Journal &jr);
/** Makes a deep copy of the object. */
Dynare(const Dynare &dyn);
- DynamicModel *
+ std::unique_ptr<DynamicModel>
clone() const override
{
- return new Dynare(*this);
+ return std::make_unique<Dynare>(*this);
}
~Dynare() override;
diff --git a/dynare++/src/dynare_params.hh b/dynare++/src/dynare_params.hh
index 3e5ebe6789bf3946992f4acd9ecf1fbcf1596e46..e8aef5c438700ecb939b9aab857cf011aef7171c 100644
--- a/dynare++/src/dynare_params.hh
+++ b/dynare++/src/dynare_params.hh
@@ -40,7 +40,7 @@ struct DynareParams
/** Flag for doing IRFs even if the irf_list is empty. */
bool do_irfs_all;
/** List of shocks for which IRF will be calculated. */
- std::vector<const char *> irf_list;
+ std::vector<std::string> irf_list;
bool do_centralize;
double qz_criterium;
bool help;
diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc
index 700f61fc6b1ba4fc62e49a9863146b0cc95e69e1..00e0c2d78f813cff2d370a0b55c55d6a583a42ed 100644
--- a/mex/sources/k_order_perturbation/k_ord_dynare.cc
+++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc
@@ -340,7 +340,7 @@ KordpDynare::ReorderDynareJacobianIndices() noexcept(false)
/* DynareNameList class */
/**************************************************************************************/
-DynareNameList::DynareNameList(const KordpDynare &dynare, const std::vector<std::string> &names_arg) : names(names_arg)
+DynareNameList::DynareNameList(const KordpDynare &dynare, std::vector<std::string> names_arg) : names(std::move(names_arg))
{
}
@@ -348,7 +348,7 @@ DynareStateNameList::DynareStateNameList(const KordpDynare &dynare, const Dynare
const DynareNameList &denl)
{
for (int i = 0; i < dynare.nys(); i++)
- names.push_back(std::string{dnl.getName(i+dynare.nstat())});
+ names.emplace_back(dnl.getName(i+dynare.nstat()));
for (int i = 0; i < dynare.nexog(); i++)
- names.push_back(std::string{denl.getName(i)});
+ names.emplace_back(denl.getName(i));
}
diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.hh b/mex/sources/k_order_perturbation/k_ord_dynare.hh
index 49cc1d96ec58ca3e15b214782281e3a6ee372ef2..f8325083c73ad3c5a3e606a5a305fe45487882fd 100644
--- a/mex/sources/k_order_perturbation/k_ord_dynare.hh
+++ b/mex/sources/k_order_perturbation/k_ord_dynare.hh
@@ -44,16 +44,16 @@ class DynareNameList : public NameList
{
std::vector<std::string> names;
public:
- DynareNameList(const KordpDynare &dynare, const std::vector<std::string> &names_arg);
+ DynareNameList(const KordpDynare &dynare, std::vector<std::string> names_arg);
int
getNum() const
{
- return (int) names.size();
+ return static_cast<int>(names.size());
}
- const char *
+ const std::string &
getName(int i) const
{
- return names[i].c_str();
+ return names[i];
}
};
@@ -66,12 +66,12 @@ public:
int
getNum() const
{
- return (int) names.size();
+ return static_cast<int>(names.size());
}
- const char *
+ const std::string &
getName(int i) const
{
- return names[i].c_str();
+ return names[i];
}
};
/*********************************************/
@@ -235,7 +235,7 @@ public:
const ConstVector &yyp, const Vector &xx) noexcept(false);
void calcDerivativesAtSteady();
std::unique_ptr<DynamicModelAC> dynamicModelFile;
- DynamicModel *
+ std::unique_ptr<DynamicModel>
clone() const
{
std::cerr << "KordpDynare::clone() not implemented" << std::endl;