diff --git a/dynare++/extern/matlab/dynare_simul.cc b/dynare++/extern/matlab/dynare_simul.cc
index 6036090a2db84caf7d12ed9e4f398c4a075a92ef..56486d608fea5fe6ea144e956aaa152f6c3176e1 100644
--- a/dynare++/extern/matlab/dynare_simul.cc
+++ b/dynare++/extern/matlab/dynare_simul.cc
@@ -113,12 +113,10 @@ extern "C" {
ConstTwoDMatrix vcov_mat(nexog, nexog, ConstVector{vcov});
GenShockRealization sr(vcov_mat, shocks_mat, seed);
// simulate and copy the results
- TwoDMatrix *res_mat
- = dr.simulate(DecisionRule::horner, nper,
- ConstVector{ystart}, sr);
+ TwoDMatrix res_mat{dr.simulate(DecisionRule::emethod::horner, nper,
+ ConstVector{ystart}, sr)};
TwoDMatrix res_tmp_mat{ny, nper, Vector{res}};
- res_tmp_mat = (const TwoDMatrix &) (*res_mat);
- delete res_mat;
+ res_tmp_mat = const_cast<const TwoDMatrix &>(res_mat);
plhs[1] = res;
}
catch (const KordException &e)
diff --git a/dynare++/kord/approximation.cc b/dynare++/kord/approximation.cc
index b4eb7975700002895d9af23dc4c0a6ce0293556e..12c0dd0c9eece64ba7c1573984e8ee3960a1e6d4 100644
--- a/dynare++/kord/approximation.cc
+++ b/dynare++/kord/approximation.cc
@@ -10,8 +10,7 @@
ZAuxContainer::ZAuxContainer(const _Ctype *gss, int ngss, int ng, int ny, int nu)
: StackContainer<FGSTensor>(4, 1)
{
- stack_sizes[0] = ngss; stack_sizes[1] = ng;
- stack_sizes[2] = ny; stack_sizes[3] = nu;
+ stack_sizes = { ngss, ng, ny, nu };
conts[0] = gss;
calculateOffsets();
}
@@ -32,34 +31,21 @@ ZAuxContainer::getType(int i, const Symmetry &s) const
}
Approximation::Approximation(DynamicModel &m, Journal &j, int ns, bool dr_centr, double qz_crit)
- : model(m), journal(j), rule_ders(nullptr), rule_ders_ss(nullptr), fdr(nullptr), udr(nullptr),
+ : model(m), journal(j),
ypart(model.nstat(), model.npred(), model.nboth(), model.nforw()),
- mom(UNormalMoments(model.order(), model.getVcov())), nvs(4), steps(ns),
+ mom(UNormalMoments(model.order(), model.getVcov())),
+ nvs{ypart.nys(), model.nexog(), model.nexog(), 1 },
+ steps(ns),
dr_centralize(dr_centr), qz_criterium(qz_crit), ss(ypart.ny(), steps+1)
{
- nvs[0] = ypart.nys(); nvs[1] = model.nexog();
- nvs[2] = model.nexog(); nvs[3] = 1;
-
ss.nans();
}
-Approximation::~Approximation()
-{
- if (rule_ders_ss)
- delete rule_ders_ss;
- if (rule_ders)
- delete rule_ders;
- if (fdr)
- delete fdr;
- if (udr)
- delete udr;
-}
-
/* This just returns |fdr| with a check that it is created. */
const FoldDecisionRule &
Approximation::getFoldDecisionRule() const
{
- KORD_RAISE_IF(fdr == nullptr,
+ KORD_RAISE_IF(!fdr,
"Folded decision rule has not been created in Approximation::getFoldDecisionRule");
return *fdr;
}
@@ -68,7 +54,7 @@ Approximation::getFoldDecisionRule() const
const UnfoldDecisionRule &
Approximation::getUnfoldDecisionRule() const
{
- KORD_RAISE_IF(udr == nullptr,
+ KORD_RAISE_IF(!udr,
"Unfolded decision rule has not been created in Approximation::getUnfoldDecisionRule");
return *udr;
}
@@ -147,7 +133,7 @@ Approximation::walkStochSteady()
model.solveDeterministicSteady();
approxAtSteady();
Vector steady0{ss.getCol(0)};
- steady0 = (const Vector &) model.getSteady();
+ steady0 = model.getSteady();
double sigma_so_far = 0.0;
double dsigma = (steps == 0) ? 0.0 : 1.0/steps;
@@ -156,14 +142,14 @@ Approximation::walkStochSteady()
JournalRecordPair pa(journal);
pa << "Approximation about stochastic steady for sigma=" << sigma_so_far+dsigma << endrec;
- Vector last_steady((const Vector &)model.getSteady());
+ Vector last_steady(model.getSteady());
// calculate fix-point of the last rule for |dsigma|
/* 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);
- bool converged = fp.calcFixPoint(DecisionRule::horner, model.getSteady());
+ bool converged = fp.calcFixPoint(DecisionRule::emethod::horner, model.getSteady());
JournalRecord rec(journal);
rec << "Fix point calcs: iter=" << fp.getNumIter() << ", newton_iter="
<< fp.getNewtonTotalIter() << ", last_newton_iter=" << fp.getNewtonLastIter() << ".";
@@ -175,14 +161,14 @@ Approximation::walkStochSteady()
KORD_RAISE_X("Fix point calculation not converged", KORD_FP_NOT_CONV);
}
Vector steadyi{ss.getCol(i)};
- steadyi = (const Vector &) model.getSteady();
+ steadyi = model.getSteady();
// calculate |hh| as expectations of the last $g^{**}$
/* We form the steady state shift |dy|, which is the new steady state
minus the old steady state. Then we create |StochForwardDerivs|
object, which calculates the derivatives of $g^{**}$ expectations at
new sigma and new steady. */
- Vector dy((const Vector &)model.getSteady());
+ Vector dy(model.getSteady());
dy.add(-1.0, last_steady);
StochForwardDerivs<KOrder::fold> hh(ypart, model.nexog(), *rule_ders_ss, mom, dy,
@@ -197,9 +183,8 @@ 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<KOrder::fold>(d);
+
saveRuleDerivs(korder_stoch.getFoldDers());
check(sigma_so_far+dsigma);
@@ -207,24 +192,14 @@ Approximation::walkStochSteady()
}
// construct the resulting decision rules
- if (fdr)
- {
- delete fdr;
- fdr = nullptr;
- }
- if (udr)
- {
- delete udr;
- udr = nullptr;
- }
-
- fdr = new FoldDecisionRule(*rule_ders, ypart, model.nexog(),
- model.getSteady(), 1.0-sigma_so_far);
+ udr.reset();
+ fdr = std::make_unique<FoldDecisionRule>(*rule_ders, ypart, model.nexog(),
+ model.getSteady(), 1.0-sigma_so_far);
if (steps == 0 && dr_centralize)
{
// centralize decision rule for zero steps
DRFixPoint<KOrder::fold> fp(*rule_ders, ypart, model.getSteady(), 1.0);
- bool converged = fp.calcFixPoint(DecisionRule::horner, model.getSteady());
+ bool converged = fp.calcFixPoint(DecisionRule::emethod::horner, model.getSteady());
JournalRecord rec(journal);
rec << "Fix point calcs: iter=" << fp.getNumIter() << ", newton_iter="
<< fp.getNewtonTotalIter() << ", last_newton_iter=" << fp.getNewtonLastIter() << ".";
@@ -239,9 +214,7 @@ Approximation::walkStochSteady()
{
JournalRecordPair recp(journal);
recp << "Centralizing about fix-point." << endrec;
- FoldDecisionRule *dr_backup = fdr;
- fdr = new FoldDecisionRule(*dr_backup, model.getSteady());
- delete dr_backup;
+ fdr = std::make_unique<FoldDecisionRule>(*fdr, model.getSteady());
}
}
}
@@ -254,14 +227,9 @@ Approximation::walkStochSteady()
void
Approximation::saveRuleDerivs(const FGSContainer &g)
{
- if (rule_ders)
- {
- delete rule_ders;
- delete rule_ders_ss;
- }
- rule_ders = new FGSContainer(g);
- rule_ders_ss = new FGSContainer(4);
- for (auto & run : (*rule_ders))
+ rule_ders = std::make_unique<FGSContainer>(g);
+ rule_ders_ss = std::make_unique<FGSContainer>(4);
+ for (auto &run : *rule_ders)
{
auto ten = std::make_unique<FGSTensor>(ypart.nstat+ypart.npred, ypart.nyss(), *(run.second));
rule_ders_ss->insert(std::move(ten));
@@ -287,7 +255,7 @@ Approximation::calcStochShift(Vector &out, double at_sigma) const
"Wrong length of output vector for Approximation::calcStochShift");
out.zeros();
- ZAuxContainer zaux(rule_ders_ss, ypart.nyss(), ypart.ny(),
+ ZAuxContainer zaux(rule_ders_ss.get(), ypart.nyss(), ypart.ny(),
ypart.nys(), model.nexog());
int dfac = 1;
@@ -298,7 +266,7 @@ Approximation::calcStochShift(Vector &out, double at_sigma) const
Symmetry sym{0, d, 0, 0};
// calculate $F_{u'^d}$ via |ZAuxContainer|
- FGSTensor *ten = new FGSTensor(ypart.ny(), TensorDimens(sym, nvs));
+ auto ten = std::make_unique<FGSTensor>(ypart.ny(), TensorDimens(sym, nvs));
ten->zeros();
for (int l = 1; l <= d; l++)
{
@@ -307,13 +275,11 @@ Approximation::calcStochShift(Vector &out, double at_sigma) const
}
// multiply with shocks and add to result
- FGSTensor *tmp = new FGSTensor(ypart.ny(), TensorDimens(Symmetry{0, 0, 0, 0}, nvs));
+ auto tmp = std::make_unique<FGSTensor>(ypart.ny(), TensorDimens(Symmetry{0, 0, 0, 0}, nvs));
tmp->zeros();
ten->contractAndAdd(1, *tmp, mom.get(Symmetry{d}));
out.add(pow(at_sigma, d)/dfac, tmp->getData());
- delete ten;
- delete tmp;
}
}
}
@@ -359,7 +325,7 @@ Approximation::check(double at_sigma) const
So we invoke the Sylvester solver for the first dimension with $A=I$,
$B=-G$, $C=G^T$ and $D=g_u\Sigma g_u^T$. */
-TwoDMatrix *
+TwoDMatrix
Approximation::calcYCov() const
{
const TwoDMatrix &gy = rule_ders->get(Symmetry{1, 0, 0, 0});
@@ -367,7 +333,7 @@ Approximation::calcYCov() const
TwoDMatrix G(model.numeq(), model.numeq());
G.zeros();
G.place(gy, 0, model.nstat());
- TwoDMatrix B((const TwoDMatrix &)G);
+ TwoDMatrix B(const_cast<const TwoDMatrix &>(G));
B.mult(-1.0);
TwoDMatrix C(transpose(G));
TwoDMatrix A(model.numeq(), model.numeq());
@@ -375,10 +341,10 @@ Approximation::calcYCov() const
for (int i = 0; i < model.numeq(); i++)
A.get(i, i) = 1.0;
- TwoDMatrix *X = new TwoDMatrix((gu * model.getVcov()) * transpose(gu));
+ TwoDMatrix X((gu * model.getVcov()) * transpose(gu));
GeneralSylvester gs(1, model.numeq(), model.numeq(), 0,
- A.getData(), B.getData(), C.getData(), X->getData());
+ A.getData(), B.getData(), C.getData(), X.getData());
gs.solve();
return X;
diff --git a/dynare++/kord/approximation.hh b/dynare++/kord/approximation.hh
index 2df8641557131b4c8634ef5be5ffbba2a7818451..0cc77f9f58e65c2f760020dbf85ea83bd6a9c2b3 100644
--- a/dynare++/kord/approximation.hh
+++ b/dynare++/kord/approximation.hh
@@ -55,6 +55,8 @@
#include "korder.hh"
#include "journal.hh"
+#include <memory>
+
/* This class is used to calculate derivatives by faa Di Bruno of the
$$f(g^{**}(g^*(y^*,u,\sigma),u',\sigma),g(y^*,u,\sigma),y^*,u)$$ with
respect $u'$. In order to keep it as simple as possible, the class
@@ -108,10 +110,10 @@ class Approximation
{
DynamicModel &model;
Journal &journal;
- FGSContainer *rule_ders;
- FGSContainer *rule_ders_ss;
- FoldDecisionRule *fdr;
- UnfoldDecisionRule *udr;
+ std::unique_ptr<FGSContainer> rule_ders;
+ std::unique_ptr<FGSContainer> rule_ders_ss;
+ std::unique_ptr<FoldDecisionRule> fdr;
+ std::unique_ptr<UnfoldDecisionRule> udr;
const PartitionY ypart;
const FNormalMoments mom;
IntSequence nvs;
@@ -121,11 +123,9 @@ class Approximation
TwoDMatrix ss;
public:
Approximation(DynamicModel &m, Journal &j, int ns, bool dr_centr, double qz_crit);
- virtual
- ~Approximation();
- const FoldDecisionRule&getFoldDecisionRule() const;
- const UnfoldDecisionRule&getUnfoldDecisionRule() const;
+ const FoldDecisionRule &getFoldDecisionRule() const;
+ const UnfoldDecisionRule &getUnfoldDecisionRule() const;
const TwoDMatrix &
getSS() const
{
@@ -138,16 +138,16 @@ public:
}
void walkStochSteady();
- TwoDMatrix *calcYCov() const;
- const FGSContainer *
+ TwoDMatrix calcYCov() const;
+ const FGSContainer &
get_rule_ders() const
{
- return rule_ders;
+ return *rule_ders;
}
- const FGSContainer *
+ const FGSContainer &
get_rule_ders_ss() const
{
- return rule_ders;
+ return *rule_ders_ss;
}
protected:
void approxAtSteady();
diff --git a/dynare++/kord/decision_rule.cc b/dynare++/kord/decision_rule.cc
index 6e86cce92dbba82c9fc12d73a997e00aa0b00870..10de34f239d88b010b35aefb59f8cb8b003dfc36 100644
--- a/dynare++/kord/decision_rule.cc
+++ b/dynare++/kord/decision_rule.cc
@@ -13,29 +13,12 @@
#include <utility>
#include <memory>
-template <>
-int DRFixPoint<KOrder::fold>::max_iter = 10000;
-template <>
-int DRFixPoint<KOrder::unfold>::max_iter = 10000;
-template <>
-double DRFixPoint<KOrder::fold>::tol = 1.e-10;
-template <>
-double DRFixPoint<KOrder::unfold>::tol = 1.e-10;
-template <>
-int DRFixPoint<KOrder::fold>::max_newton_iter = 50;
-template <>
-int DRFixPoint<KOrder::unfold>::max_newton_iter = 50;
-template <>
-int DRFixPoint<KOrder::fold>::newton_pause = 100;
-template <>
-int DRFixPoint<KOrder::unfold>::newton_pause = 100;
-
// |FoldDecisionRule| conversion from |UnfoldDecisionRule|
FoldDecisionRule::FoldDecisionRule(const UnfoldDecisionRule &udr)
: DecisionRuleImpl<KOrder::fold>(ctraits<KOrder::fold>::Tpol(udr.nrows(), udr.nvars()),
udr.ypart, udr.nu, udr.ysteady)
{
- for (const auto & it : udr)
+ for (const auto &it : udr)
insert(std::make_unique<ctraits<KOrder::fold>::Ttensym>(*(it.second)));
}
@@ -44,19 +27,10 @@ UnfoldDecisionRule::UnfoldDecisionRule(const FoldDecisionRule &fdr)
: DecisionRuleImpl<KOrder::unfold>(ctraits<KOrder::unfold>::Tpol(fdr.nrows(), fdr.nvars()),
fdr.ypart, fdr.nu, fdr.ysteady)
{
- for (const auto & it : fdr)
+ for (const auto &it : fdr)
insert(std::make_unique<ctraits<KOrder::unfold>::Ttensym>(*(it.second)));
}
-SimResults::~SimResults()
-{
- for (int i = 0; i < getNumSets(); i++)
- {
- delete data[i];
- delete shocks[i];
- }
-}
-
/* This runs simulations with an output to journal file. Note that we
report how many simulations had to be thrown out due to Nan or Inf. */
@@ -92,7 +66,7 @@ SimResults::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
{
RandomShockRealization sr(vcov, seed_generator::get_new_seed());
rsrs.push_back(sr);
- gr.insert(std::make_unique<SimulationWorker>(*this, dr, DecisionRule::horner,
+ gr.insert(std::make_unique<SimulationWorker>(*this, dr, DecisionRule::emethod::horner,
num_per+num_burn, start, rsrs.back()));
}
gr.run();
@@ -103,39 +77,35 @@ SimResults::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
and shocks are thrown away. */
bool
-SimResults::addDataSet(TwoDMatrix *d, ExplicitShockRealization *sr, const ConstVector &st)
+SimResults::addDataSet(const TwoDMatrix &d, const ExplicitShockRealization &sr, const ConstVector &st)
{
- KORD_RAISE_IF(d->nrows() != num_y,
+ KORD_RAISE_IF(d.nrows() != num_y,
"Incompatible number of rows for SimResults::addDataSets");
- KORD_RAISE_IF(d->ncols() != num_per+num_burn,
+ KORD_RAISE_IF(d.ncols() != num_per+num_burn,
"Incompatible number of cols for SimResults::addDataSets");
bool ret = false;
- if (d->isFinite())
+ if (d.isFinite())
{
- data.push_back(new TwoDMatrix((const TwoDMatrix &) (*d), num_burn, num_per));
- shocks.push_back(new ExplicitShockRealization(
- ConstTwoDMatrix(sr->getShocks(), num_burn, num_per)));
+ data.emplace_back(d, num_burn, num_per);
+ shocks.emplace_back(ConstTwoDMatrix(sr.getShocks(), num_burn, num_per));
if (num_burn == 0)
start.emplace_back(st);
else
- start.emplace_back(d->getCol(num_burn-1));
+ start.emplace_back(d.getCol(num_burn-1));
ret = true;
}
- delete d;
- delete sr;
return ret;
}
void
-SimResults::writeMat(const char *base, const char *lname) const
+SimResults::writeMat(const std::string &base, const std::string &lname) const
{
- char matfile_name[100];
- sprintf(matfile_name, "%s.mat", base);
- mat_t *matfd = Mat_Create(matfile_name, nullptr);
- if (matfd != nullptr)
+ std::string matfile_name = base + ".mat";
+ mat_t *matfd = Mat_Create(matfile_name.c_str(), nullptr);
+ if (matfd)
{
- writeMat(matfd, lname);
+ writeMat(matfd, lname.c_str());
Mat_Close(matfd);
}
}
@@ -144,17 +114,14 @@ SimResults::writeMat(const char *base, const char *lname) const
appended. If there is only one matrix, the index is not appended. */
void
-SimResults::writeMat(mat_t *fd, const char *lname) const
+SimResults::writeMat(mat_t *fd, const std::string &lname) const
{
- char tmp[100];
for (int i = 0; i < getNumSets(); i++)
{
+ std::string tmp = lname + "_data";
if (getNumSets() > 1)
- sprintf(tmp, "%s_data%d", lname, i+1);
- else
- sprintf(tmp, "%s_data", lname);
- ConstTwoDMatrix m(*(data[i]));
- m.writeMat(fd, tmp);
+ tmp += std::to_string(i+1);
+ data[i].writeMat(fd, tmp);
}
}
@@ -178,14 +145,10 @@ SimResultsStats::simulate(int num_sim, const DecisionRule &dr,
/* Here we do not save the data itself, we save only mean and vcov. */
void
-SimResultsStats::writeMat(mat_t *fd, const char *lname) const
+SimResultsStats::writeMat(mat_t *fd, const std::string &lname) const
{
- char tmp[100];
- sprintf(tmp, "%s_mean", lname);
- ConstTwoDMatrix m(num_y, 1, mean);
- m.writeMat(fd, tmp);
- sprintf(tmp, "%s_vcov", lname);
- ConstTwoDMatrix(vcov).writeMat(fd, tmp);
+ ConstTwoDMatrix(num_y, 1, mean).writeMat(fd, lname + "_mean");;
+ vcov.writeMat(fd, lname + "_vcov");
}
void
@@ -195,11 +158,11 @@ SimResultsStats::calcMean()
if (data.size()*num_per > 0)
{
double mult = 1.0/data.size()/num_per;
- for (auto & i : data)
+ for (const auto &i : data)
{
for (int j = 0; j < num_per; j++)
{
- ConstVector col{i->getCol(j)};
+ ConstVector col{i.getCol(j)};
mean.add(mult, col);
}
}
@@ -213,28 +176,19 @@ SimResultsStats::calcVcov()
{
vcov.zeros();
double mult = 1.0/(data.size()*num_per - 1);
- for (auto & i : data)
- {
- const TwoDMatrix &d = *i;
- for (int j = 0; j < num_per; j++)
- {
- for (int m = 0; m < num_y; m++)
- {
- for (int n = m; n < num_y; n++)
- {
- double s = (d.get(m, j)-mean[m])*(d.get(n, j)-mean[n]);
- vcov.get(m, n) += mult*s;
- if (m != n)
- vcov.get(n, m) += mult*s;
- }
- }
- }
- }
+ for (const auto &d : data)
+ for (int j = 0; j < num_per; j++)
+ for (int m = 0; m < num_y; m++)
+ for (int n = m; n < num_y; n++)
+ {
+ double s = (d.get(m, j)-mean[m])*(d.get(n, j)-mean[n]);
+ vcov.get(m, n) += mult*s;
+ if (m != n)
+ vcov.get(n, m) += mult*s;
+ }
}
else
- {
- vcov.infs();
- }
+ vcov.infs();
}
void
@@ -256,13 +210,10 @@ SimResultsDynamicStats::simulate(int num_sim, const DecisionRule &dr,
}
void
-SimResultsDynamicStats::writeMat(mat_t *fd, const char *lname) const
+SimResultsDynamicStats::writeMat(mat_t *fd, const std::string &lname) const
{
- char tmp[100];
- sprintf(tmp, "%s_cond_mean", lname);
- ConstTwoDMatrix(mean).writeMat(fd, tmp);
- sprintf(tmp, "%s_cond_variance", lname);
- ConstTwoDMatrix(variance).writeMat(fd, tmp);
+ mean.writeMat(fd, lname + "_cond_mean");
+ variance.writeMat(fd, lname + "_cond_variance");
}
void
@@ -275,9 +226,9 @@ SimResultsDynamicStats::calcMean()
for (int j = 0; j < num_per; j++)
{
Vector meanj{mean.getCol(j)};
- for (auto & i : data)
+ for (const auto &i : data)
{
- ConstVector col{i->getCol(j)};
+ ConstVector col{i.getCol(j)};
meanj.add(mult, col);
}
}
@@ -295,9 +246,9 @@ SimResultsDynamicStats::calcVariance()
{
ConstVector meanj{mean.getCol(j)};
Vector varj{variance.getCol(j)};
- for (auto & i : data)
+ for (const auto &i : data)
{
- Vector col{i->getCol(j)};
+ Vector col{i.getCol(j)};
col.add(-1.0, meanj);
for (int k = 0; k < col.length(); k++)
col[k] = col[k]*col[k];
@@ -306,9 +257,7 @@ SimResultsDynamicStats::calcVariance()
}
}
else
- {
- variance.infs();
- }
+ variance.infs();
}
void
@@ -334,7 +283,7 @@ SimResultsIRF::simulate(const DecisionRule &dr)
{
sthread::detach_thread_group gr;
for (int idata = 0; idata < control.getNumSets(); idata++)
- gr.insert(std::make_unique<SimulationIRFWorker>(*this, dr, DecisionRule::horner,
+ gr.insert(std::make_unique<SimulationIRFWorker>(*this, dr, DecisionRule::emethod::horner,
num_per, idata, ishock, imp));
gr.run();
}
@@ -345,8 +294,8 @@ SimResultsIRF::calcMeans()
means.zeros();
if (data.size() > 0)
{
- for (auto & i : data)
- means.add(1.0, *i);
+ for (const auto &i : data)
+ means.add(1.0, i);
means.mult(1.0/data.size());
}
}
@@ -357,9 +306,9 @@ SimResultsIRF::calcVariances()
if (data.size() > 1)
{
variances.zeros();
- for (auto & i : data)
+ for (const auto &i : data)
{
- TwoDMatrix d((const TwoDMatrix &)(*i));
+ TwoDMatrix d(i);
d.add(-1.0, means);
for (int j = 0; j < d.nrows(); j++)
for (int k = 0; k < d.ncols(); k++)
@@ -368,19 +317,14 @@ SimResultsIRF::calcVariances()
}
}
else
- {
- variances.infs();
- }
+ variances.infs();
}
void
-SimResultsIRF::writeMat(mat_t *fd, const char *lname) const
+SimResultsIRF::writeMat(mat_t *fd, const std::string &lname) const
{
- char tmp[100];
- sprintf(tmp, "%s_mean", lname);
- means.writeMat(fd, tmp);
- sprintf(tmp, "%s_var", lname);
- variances.writeMat(fd, tmp);
+ means.writeMat(fd, lname + "_mean");
+ variances.writeMat(fd, lname + "_var");
}
void
@@ -416,21 +360,17 @@ RTSimResultsStats::simulate(int num_sim, const DecisionRule &dr, const Vector &s
{
RandomShockRealization sr(vcov, seed_generator::get_new_seed());
rsrs.push_back(sr);
- gr.insert(std::make_unique<RTSimulationWorker>(*this, dr, DecisionRule::horner,
+ gr.insert(std::make_unique<RTSimulationWorker>(*this, dr, DecisionRule::emethod::horner,
num_per, start, rsrs.back()));
}
gr.run();
}
void
-RTSimResultsStats::writeMat(mat_t *fd, const char *lname)
+RTSimResultsStats::writeMat(mat_t *fd, const std::string &lname)
{
- char tmp[100];
- sprintf(tmp, "%s_rt_mean", lname);
- ConstTwoDMatrix m(nc.getDim(), 1, mean);
- m.writeMat(fd, tmp);
- sprintf(tmp, "%s_rt_vcov", lname);
- ConstTwoDMatrix(vcov).writeMat(fd, tmp);
+ ConstTwoDMatrix(nc.getDim(), 1, mean).writeMat(fd, lname + "_rt_mean");
+ vcov.writeMat(fd, lname + "_rt_vcov");
}
IRFResults::IRFResults(const DynamicModel &mod, const DecisionRule &dr,
@@ -440,51 +380,42 @@ IRFResults::IRFResults(const DynamicModel &mod, const DecisionRule &dr,
{
int num_per = control.getNumPer();
JournalRecordPair pa(journal);
- pa << "Calculating IRFs against control for " << (int) irf_list_ind.size() << " shocks and for "
+ pa << "Calculating IRFs against control for " << static_cast<int>(irf_list_ind.size()) << " shocks and for "
<< num_per << " periods" << endrec;
const TwoDMatrix &vcov = mod.getVcov();
for (int ishock : irf_list_ind)
{
double stderror = sqrt(vcov.get(ishock, ishock));
- irf_res.push_back(new SimResultsIRF(control, model.numeq(), num_per,
- ishock, stderror));
- irf_res.push_back(new SimResultsIRF(control, model.numeq(), num_per,
- ishock, -stderror));
+ irf_res.emplace_back(control, model.numeq(), num_per,
+ ishock, stderror);
+ irf_res.emplace_back(control, model.numeq(), num_per,
+ ishock, -stderror);
}
for (unsigned int ii = 0; ii < irf_list_ind.size(); ii++)
{
- irf_res[2*ii]->simulate(dr, journal);
- irf_res[2*ii+1]->simulate(dr, journal);
+ irf_res[2*ii].simulate(dr, journal);
+ irf_res[2*ii+1].simulate(dr, journal);
}
}
-IRFResults::~IRFResults()
-{
- for (auto & irf_re : irf_res)
- delete irf_re;
-}
-
void
-IRFResults::writeMat(mat_t *fd, const char *prefix) const
+IRFResults::writeMat(mat_t *fd, const std::string &prefix) const
{
for (unsigned int i = 0; i < irf_list_ind.size(); i++)
{
- char tmp[100];
int ishock = irf_list_ind[i];
- const char *shockname = model.getExogNames().getName(ishock);
- sprintf(tmp, "%s_irfp_%s", prefix, shockname);
- irf_res[2*i]->writeMat(fd, tmp);
- sprintf(tmp, "%s_irfm_%s", prefix, shockname);
- irf_res[2*i+1]->writeMat(fd, tmp);
+ auto shockname = model.getExogNames().getName(ishock);
+ irf_res[2*i].writeMat(fd, prefix + "_irfp_" + shockname);
+ irf_res[2*i+1].writeMat(fd, prefix + "_irfm_" + shockname);
}
}
void
SimulationWorker::operator()(std::mutex &mut)
{
- auto *esr = new ExplicitShockRealization(sr, np);
- TwoDMatrix *m = dr.simulate(em, np, st, *esr);
+ ExplicitShockRealization esr(sr, np);
+ TwoDMatrix m{dr.simulate(em, np, st, esr)};
{
std::unique_lock<std::mutex> lk{mut};
res.addDataSet(m, esr, st);
@@ -497,11 +428,10 @@ SimulationWorker::operator()(std::mutex &mut)
void
SimulationIRFWorker::operator()(std::mutex &mut)
{
- auto *esr
- = new ExplicitShockRealization(res.control.getShocks(idata));
- esr->addToShock(ishock, 0, imp);
- TwoDMatrix *m = dr.simulate(em, np, res.control.getStart(idata), *esr);
- m->add(-1.0, res.control.getData(idata));
+ ExplicitShockRealization esr(res.control.getShocks(idata));
+ esr.addToShock(ishock, 0, imp);
+ TwoDMatrix m{dr.simulate(em, np, res.control.getStart(idata), esr)};
+ m.add(-1.0, res.control.getData(idata));
{
std::unique_lock<std::mutex> lk{mut};
res.addDataSet(m, esr, res.control.getStart(idata));
@@ -579,8 +509,7 @@ RandomShockRealization::choleskyFactor(const ConstTwoDMatrix &v)
void
RandomShockRealization::schurFactor(const ConstTwoDMatrix &v)
{
- SymSchurDecomp ssd(v);
- ssd.getFactor(factor);
+ SymSchurDecomp(v).getFactor(factor);
}
void
diff --git a/dynare++/kord/decision_rule.hh b/dynare++/kord/decision_rule.hh
index 97a92df675894677ea0d7b8444ba6636e2b7fa40..ed0165c7f827608b18d8acbc0fe51058fb8a7f8a 100644
--- a/dynare++/kord/decision_rule.hh
+++ b/dynare++/kord/decision_rule.hh
@@ -23,7 +23,9 @@
#include "korder.hh"
#include "normal_conjugate.hh"
+#include <memory>
#include <random>
+#include <string>
/* This is a general interface to a shock realizations. The interface
has only one method returning the shock realizations at the given
@@ -33,8 +35,7 @@
class ShockRealization
{
public:
- virtual ~ShockRealization()
- = default;
+ virtual ~ShockRealization() = default;
virtual void get(int n, Vector &out) = 0;
virtual int numShocks() const = 0;
};
@@ -55,19 +56,18 @@ public:
class DecisionRule
{
public:
- enum emethod { horner, trad };
- virtual ~DecisionRule()
- = default;
- virtual TwoDMatrix *simulate(emethod em, int np, const ConstVector &ystart,
- ShockRealization &sr) const = 0;
+ enum class emethod { horner, trad };
+ virtual ~DecisionRule() = default;
+ virtual TwoDMatrix simulate(emethod em, int np, const ConstVector &ystart,
+ ShockRealization &sr) const = 0;
virtual void eval(emethod em, Vector &out, const ConstVector &v) const = 0;
virtual void evaluate(emethod em, Vector &out, const ConstVector &ys,
const ConstVector &u) const = 0;
- virtual void writeMat(mat_t *fd, const char *prefix) const = 0;
- virtual DecisionRule *centralizedClone(const Vector &fixpoint) const = 0;
- virtual const Vector&getSteady() const = 0;
+ virtual void writeMat(mat_t *fd, const std::string &prefix) const = 0;
+ virtual std::unique_ptr<DecisionRule> centralizedClone(const Vector &fixpoint) const = 0;
+ virtual const Vector &getSteady() const = 0;
virtual int nexog() const = 0;
- virtual const PartitionY&getYPart() const = 0;
+ virtual const PartitionY &getYPart() const = 0;
};
/* The main purpose of this class is to implement |DecisionRule|
@@ -131,12 +131,12 @@ public:
{
return ysteady;
}
- TwoDMatrix *simulate(emethod em, int np, const ConstVector &ystart,
- ShockRealization &sr) const override;
+ TwoDMatrix simulate(emethod em, int np, const ConstVector &ystart,
+ ShockRealization &sr) const override;
void evaluate(emethod em, Vector &out, const ConstVector &ys,
const ConstVector &u) const override;
- DecisionRule *centralizedClone(const Vector &fixpoint) const override;
- void writeMat(mat_t *fd, const char *prefix) const override;
+ std::unique_ptr<DecisionRule> centralizedClone(const Vector &fixpoint) const override;
+ void writeMat(mat_t *fd, const std::string &prefix) const override;
int
nexog() const override
@@ -279,13 +279,13 @@ DecisionRuleImpl<t>::centralize(const DecisionRuleImpl &dr)
|ysteady| is added to all columns of the result. */
template <int t>
-TwoDMatrix *
+TwoDMatrix
DecisionRuleImpl<t>::simulate(emethod em, int np, const ConstVector &ystart,
ShockRealization &sr) const
{
KORD_RAISE_IF(ysteady.length() != ystart.length(),
"Start and steady lengths differ in DecisionRuleImpl::simulate");
- auto *res = new TwoDMatrix(ypart.ny(), np);
+ TwoDMatrix res(ypart.ny(), np);
// initialize vectors and subvectors for simulation
/* Here allocate the stack vector $(\Delta y^*, u)$, define the
@@ -303,7 +303,7 @@ DecisionRuleImpl<t>::simulate(emethod em, int np, const ConstVector &ystart,
dy = ystart_pred;
dy.add(-1.0, ysteady_pred);
sr.get(0, u);
- Vector out{res->getCol(0)};
+ Vector out{res.getCol(0)};
eval(em, out, dyu);
// perform all other steps of simulations
@@ -312,17 +312,17 @@ DecisionRuleImpl<t>::simulate(emethod em, int np, const ConstVector &ystart,
int i = 1;
while (i < np)
{
- ConstVector ym{res->getCol(i-1)};
+ ConstVector ym{res.getCol(i-1)};
ConstVector dym(ym, ypart.nstat, ypart.nys());
dy = dym;
sr.get(i, u);
- Vector out{res->getCol(i)};
+ Vector out{res.getCol(i)};
eval(em, out, dyu);
if (!out.isFinite())
{
if (i+1 < np)
{
- TwoDMatrix rest(*res, i+1, np-i-1);
+ TwoDMatrix rest(res, i+1, np-i-1);
rest.zeros();
}
break;
@@ -335,7 +335,7 @@ DecisionRuleImpl<t>::simulate(emethod em, int np, const ConstVector &ystart,
and leave the padded columns to zero. */
for (int j = 0; j < i; j++)
{
- Vector col{res->getCol(j)};
+ Vector col{res.getCol(j)};
col.add(1.0, ysteady);
}
@@ -371,10 +371,10 @@ DecisionRuleImpl<t>::evaluate(emethod em, Vector &out, const ConstVector &ys,
centralized constructor. */
template <int t>
-DecisionRule *
+std::unique_ptr<DecisionRule>
DecisionRuleImpl<t>::centralizedClone(const Vector &fixpoint) const
{
- return new DecisionRuleImpl<t>(*this, fixpoint);
+ return std::make_unique<DecisionRuleImpl<t>>(*this, fixpoint);
}
/* Here we only encapsulate two implementations to one, deciding
@@ -384,7 +384,7 @@ template <int t>
void
DecisionRuleImpl<t>::eval(emethod em, Vector &out, const ConstVector &v) const
{
- if (em == DecisionRule::horner)
+ if (em == emethod::horner)
_Tparent::evalHorner(out, v);
else
_Tparent::evalTrad(out, v);
@@ -394,14 +394,12 @@ DecisionRuleImpl<t>::eval(emethod em, Vector &out, const ConstVector &v) const
template <int t>
void
-DecisionRuleImpl<t>::writeMat(mat_t *fd, const char *prefix) const
+DecisionRuleImpl<t>::writeMat(mat_t *fd, const std::string &prefix) const
{
ctraits<t>::Tpol::writeMat(fd, prefix);
TwoDMatrix dum(ysteady.length(), 1);
dum.getData() = ysteady;
- char tmp[100];
- sprintf(tmp, "%s_ss", prefix);
- ConstTwoDMatrix(dum).writeMat(fd, tmp);
+ ConstTwoDMatrix(dum).writeMat(fd, prefix + "_ss");
}
/* This is exactly the same as |DecisionRuleImpl<KOrder::fold>|. The
@@ -483,20 +481,18 @@ template <int t>
class DRFixPoint : public ctraits<t>::Tpol
{
using _Tparent = typename ctraits<t>::Tpol;
- static int max_iter;
- static int max_newton_iter;
- static int newton_pause;
- static double tol;
+ 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;
- _Tparent *bigf;
- _Tparent *bigfder;
+ std::unique_ptr<_Tparent> bigf;
+ std::unique_ptr<_Tparent> bigfder;
public:
using emethod = typename DecisionRule::emethod;
DRFixPoint(const _Tg &g, const PartitionY &yp,
const Vector &ys, double sigma);
-
- ~DRFixPoint() override;
bool calcFixPoint(emethod em, Vector &out);
@@ -534,24 +530,15 @@ template <int 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), bigf(nullptr), bigfder(nullptr)
+ ysteady(ys), ypart(yp)
{
fillTensors(g, sigma);
_Tparent yspol(ypart.nstat, ypart.nys(), *this);
- bigf = new _Tparent((const _Tparent &) yspol);
+ bigf = std::make_unique<_Tparent>(const_cast<const _Tparent &>(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 = new _Tparent(*bigf, 0);
-}
-
-template <int t>
-DRFixPoint<t>::~DRFixPoint()
-{
- if (bigf)
- delete bigf;
- if (bigfder)
- delete bigfder;
+ bigfder = std::make_unique<_Tparent>(*bigf, 0);
}
/* Here we fill the tensors for the |DRFixPoint| class. We ignore the
@@ -602,7 +589,7 @@ bool
DRFixPoint<t>::solveNewton(Vector &y)
{
const double urelax_threshold = 1.e-5;
- Vector sol((const Vector &)y);
+ Vector sol(const_cast<const Vector &>(y));
Vector delta(y.length());
newton_iter_last = 0;
bool delta_finite = true;
@@ -632,7 +619,7 @@ DRFixPoint<t>::solveNewton(Vector &y)
urelax = 1.0;
while (!urelax_found && urelax > urelax_threshold)
{
- Vector soltmp((const Vector &)sol);
+ Vector soltmp(const_cast<const Vector &>(sol));
soltmp.add(-urelax, delta);
Vector f(sol.length());
bigf->evalHorner(f, soltmp);
@@ -656,7 +643,7 @@ DRFixPoint<t>::solveNewton(Vector &y)
newton_iter_total += newton_iter_last;
if (!converged)
newton_iter_last = 0;
- y = (const Vector &) sol;
+ y = const_cast<const Vector &>(sol);
return converged;
}
@@ -727,16 +714,14 @@ protected:
int num_y;
int num_per;
int num_burn;
- std::vector<TwoDMatrix *> data;
- std::vector<ExplicitShockRealization *> shocks;
+ std::vector<TwoDMatrix> data;
+ std::vector<ExplicitShockRealization> shocks;
std::vector<ConstVector> start;
public:
SimResults(int ny, int nper, int nburn = 0)
: num_y(ny), num_per(nper), num_burn(nburn)
{
}
- virtual
- ~SimResults();
void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
const TwoDMatrix &vcov, Journal &journal);
void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
@@ -754,17 +739,17 @@ public:
int
getNumSets() const
{
- return (int) data.size();
+ return static_cast<int>(data.size());
}
const TwoDMatrix &
getData(int i) const
{
- return *(data[i]);
+ return data[i];
}
const ExplicitShockRealization &
getShocks(int i) const
{
- return *(shocks[i]);
+ return shocks[i];
}
const ConstVector &
getStart(int i) const
@@ -772,9 +757,9 @@ public:
return start[i];
}
- bool addDataSet(TwoDMatrix *d, ExplicitShockRealization *sr, const ConstVector &st);
- void writeMat(const char *base, const char *lname) const;
- void writeMat(mat_t *fd, const char *lname) const;
+ bool addDataSet(const TwoDMatrix &d, const ExplicitShockRealization &sr, const ConstVector &st);
+ void writeMat(const std::string &base, const std::string &lname) const;
+ void writeMat(mat_t *fd, const std::string &lname) const;
};
/* This does the same as |SimResults| plus it calculates means and
@@ -792,7 +777,7 @@ public:
}
void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
const TwoDMatrix &vcov, Journal &journal);
- void writeMat(mat_t *fd, const char *lname) const;
+ void writeMat(mat_t *fd, const std::string &lname) const;
protected:
void calcMean();
void calcVcov();
@@ -814,7 +799,7 @@ public:
}
void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
const TwoDMatrix &vcov, Journal &journal);
- void writeMat(mat_t *fd, const char *lname) const;
+ void writeMat(mat_t *fd, const std::string &lname) const;
protected:
void calcMean();
void calcVariance();
@@ -846,7 +831,7 @@ public:
}
void simulate(const DecisionRule &dr, Journal &journal);
void simulate(const DecisionRule &dr);
- void writeMat(mat_t *fd, const char *lname) const;
+ void writeMat(mat_t *fd, const std::string &lname) const;
protected:
void calcMeans();
void calcVariances();
@@ -881,7 +866,7 @@ public:
const TwoDMatrix &vcov, Journal &journal);
void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
const TwoDMatrix &vcov);
- void writeMat(mat_t *fd, const char *lname);
+ void writeMat(mat_t *fd, const std::string &lname);
};
/* For each shock, this simulates plus and minus impulse. The class
@@ -901,15 +886,14 @@ public:
class DynamicModel;
class IRFResults
{
- std::vector<SimResultsIRF *> irf_res;
+ std::vector<SimResultsIRF> irf_res;
const DynamicModel &model;
std::vector<int> irf_list_ind;
public:
IRFResults(const DynamicModel &mod, const DecisionRule &dr,
const SimResults &control, std::vector<int> ili,
Journal &journal);
- ~IRFResults();
- void writeMat(mat_t *fd, const char *prefix) const;
+ void writeMat(mat_t *fd, const std::string &prefix) const;
};
/* This worker simulates the given decision rule and inserts the result
@@ -1003,12 +987,6 @@ public:
{
schurFactor(v);
}
- RandomShockRealization(const RandomShockRealization &sr)
- : mtwister(sr.mtwister), factor(sr.factor)
- {
- }
- ~RandomShockRealization()
- override = default;
void get(int n, Vector &out) override;
int
numShocks() const override
@@ -1027,14 +1005,10 @@ class ExplicitShockRealization : virtual public ShockRealization
{
TwoDMatrix shocks;
public:
- ExplicitShockRealization(const ConstTwoDMatrix &sh)
+ explicit ExplicitShockRealization(const ConstTwoDMatrix &sh)
: shocks(sh)
{
}
- ExplicitShockRealization(const ExplicitShockRealization &sr)
- : shocks(sr.shocks)
- {
- }
ExplicitShockRealization(ShockRealization &sr, int num_per);
void get(int n, Vector &out) override;
int
@@ -1043,7 +1017,7 @@ public:
return shocks.nrows();
}
const TwoDMatrix &
- getShocks()
+ getShocks() const
{
return shocks;
}
diff --git a/dynare++/kord/global_check.cc b/dynare++/kord/global_check.cc
index 298eababc6531b2057e2c8bd88f63ffd833e333f..43027d5e5e05eea0ea6dd7c4bb627ed42417fe75 100644
--- a/dynare++/kord/global_check.cc
+++ b/dynare++/kord/global_check.cc
@@ -69,7 +69,7 @@ ResidFunction::setYU(const ConstVector &ys, const ConstVector &xx)
ystar = new Vector(ys);
u = new Vector(xx);
yplus = new Vector(model->numeq());
- approx.getFoldDecisionRule().evaluate(DecisionRule::horner,
+ approx.getFoldDecisionRule().evaluate(DecisionRule::emethod::horner,
*yplus, *ystar, *u);
// make a tensor polynomial of in-place subtensors from decision rule
@@ -309,10 +309,9 @@ GlobalChecker::checkOnEllipseAndSave(mat_t *fd, const char *prefix,
/* Here we set |ycovfac| to the symmetric Schur decomposition factor of
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 ycov{approx.calcYCov()};
+ TwoDMatrix ycovpred((const TwoDMatrix &) ycov, model.nstat(), model.nstat(),
model.npred()+model.nboth(), model.npred()+model.nboth());
- delete ycov;
SymSchurDecomp ssd(ycovpred);
ssd.correctDefinitness(1.e-05);
TwoDMatrix ycovfac(ycovpred.nrows(), ycovpred.ncols());
@@ -400,18 +399,16 @@ GlobalChecker::checkAlongSimulationAndSave(mat_t *fd, const char *prefix,
pa << "Calculating errors at " << m
<< " simulated points" << endrec;
RandomShockRealization sr(model.getVcov(), seed_generator::get_new_seed());
- TwoDMatrix *y = approx.getFoldDecisionRule().simulate(DecisionRule::horner,
- m, model.getSteady(), sr);
+ TwoDMatrix y{approx.getFoldDecisionRule().simulate(DecisionRule::emethod::horner,
+ m, model.getSteady(), sr)};
TwoDMatrix x(model.nexog(), m);
x.zeros();
TwoDMatrix out(model.numeq(), m);
- check(max_evals, *y, x, out);
+ check(max_evals, y, x, out);
char tmp[100];
sprintf(tmp, "%s_simul_points", prefix);
- y->writeMat(fd, tmp);
+ y.writeMat(fd, tmp);
sprintf(tmp, "%s_simul_errors", prefix);
out.writeMat(fd, tmp);
-
- delete y;
}
diff --git a/dynare++/kord/journal.cc b/dynare++/kord/journal.cc
index 825f730374fc4933b23f93bd37188a21cfe484ee..5f42e3384bcc5a1fa20260e10265b033092606e4 100644
--- a/dynare++/kord/journal.cc
+++ b/dynare++/kord/journal.cc
@@ -3,110 +3,85 @@
#include "journal.hh"
#include "kord_exception.hh"
-#if !defined(__MINGW32__)
-# include <sys/resource.h>
-# include <sys/utsname.h>
-#endif
-#include <cstdlib>
-#include <unistd.h>
+#include <iomanip>
+#include <cmath>
#include <ctime>
-SystemResources _sysres;
-#if defined(__MINGW32__)
-//|sysconf| Win32 implementation
-/* Here we implement |sysconf| for MinGW. We implement only page size,
- number of physial pages, and a number of available physical pages. The
- pagesize is set to 1024 bytes, real pagesize can differ but it is not
- important. We can do this since Windows kernel32 |GlobalMemoryStatus|
- call returns number of bytes.
-
- Number of online processors is not implemented and returns -1, since
- Windows kernel32 |GetSystemInfo| call is too complicated. */
+#ifndef __MINGW32__
+# include <sys/time.h> // For getrusage()
+# include <sys/resource.h> // For getrusage()
+# include <sys/utsname.h> // For uname()
+# include <cstdlib> // For getloadavg()
+# include <unistd.h> // For sysconf()
+#else
# ifndef NOMINMAX
-# define NOMINMAX // Do not define "min" and "max" macros
+# define NOMINMAX // Do not define "min" and "max" macros
# endif
-# include <windows.h>
-
-# define _SC_PAGESIZE 1
-# define _SC_PHYS_PAGES 2
-# define _SC_AVPHYS_PAGES 3
-# define _SC_NPROCESSORS_ONLN 4
-
-long
-sysconf(int name)
-{
- switch (name)
- {
- case _SC_PAGESIZE:
- return 1024;
- case _SC_PHYS_PAGES:
- {
- MEMORYSTATUS memstat;
- GlobalMemoryStatus(&memstat);
- return memstat.dwTotalPhys/1024;
- }
- case _SC_AVPHYS_PAGES:
- {
- MEMORYSTATUS memstat;
- GlobalMemoryStatus(&memstat);
- return memstat.dwAvailPhys/1024;
- }
- case _SC_NPROCESSORS_ONLN:
- return -1;
- default:
- KORD_RAISE("Not implemented in Win32 sysconf.");
- return -1;
- }
-}
+# include <windows.h> // For GlobalMemoryStatus()
#endif
-#if defined(__APPLE__)
-# define _SC_PHYS_PAGES 2
-# define _SC_AVPHYS_PAGES 3
-#endif
+const std::chrono::time_point<std::chrono::high_resolution_clock> SystemResources::start = std::chrono::high_resolution_clock::now();
-SystemResources::SystemResources()
-{
- gettimeofday(&start, nullptr);
-}
-
-long int
+/* The pagesize is set to 1024 bytes on Windows. Real pagesize can differ but
+ it is not important. We can do this since Windows kernel32
+ GlobalMemoryStatus() call returns a number of bytes. */
+long
SystemResources::pageSize()
{
+#ifndef __MINGW32__
return sysconf(_SC_PAGESIZE);
+#else
+ return 1024;
+#endif
}
-long int
+long
SystemResources::physicalPages()
{
+#ifndef __MINGW32__
return sysconf(_SC_PHYS_PAGES);
+#else
+ MEMORYSTATUS memstat;
+ GlobalMemoryStatus(&memstat);
+ return memstat.dwTotalPhys/1024;
+#endif
+}
+
+long
+SystemResources::availablePhysicalPages()
+{
+#ifndef __MINGW32__
+ return sysconf(_SC_AVPHYS_PAGES);
+#else
+ MEMORYSTATUS memstat;
+ GlobalMemoryStatus(&memstat);
+ return memstat.dwAvailPhys/1024;
+#endif
}
-long int
+long
SystemResources::onlineProcessors()
{
+#ifndef __MINGW32__
return sysconf(_SC_NPROCESSORS_ONLN);
+#else
+ return -1;
+#endif
}
-long int
+long
SystemResources::availableMemory()
{
- return pageSize()*sysconf(_SC_AVPHYS_PAGES);
+ return pageSize()*availablePhysicalPages();
}
-/* Here we read the current values of resource usage. For MinGW, we
- implement only a number of available physical memory pages. */
-
-void
-SystemResources::getRUS(double &load_avg, long int &pg_avail,
- double &utime, double &stime, double &elapsed,
- long int &idrss, long int &majflt)
+SystemResources::SystemResources()
{
- struct timeval now;
- gettimeofday(&now, nullptr);
- elapsed = now.tv_sec-start.tv_sec + (now.tv_usec-start.tv_usec)*1.0e-6;
+ auto now = std::chrono::high_resolution_clock::now();
+ std::chrono::duration<double> duration = now - start;
+ elapsed = duration.count();
-#if !defined(__MINGW32__)
+#ifndef __MINGW32__
struct rusage rus;
getrusage(RUSAGE_SELF, &rus);
utime = rus.ru_utime.tv_sec+rus.ru_utime.tv_usec*1.0e-6;
@@ -120,26 +95,17 @@ SystemResources::getRUS(double &load_avg, long int &pg_avail,
majflt = -1;
#endif
-#define MINGCYGTMP (!defined(__MINGW32__) && !defined(__CYGWIN32__) && !defined(__CYGWIN__))
-#define MINGCYG (MINGCYGTMP && !defined(__MINGW64__) && !defined(__CYGWIN64__))
-
-#if MINGCYG
+#ifndef __MINGW32__
getloadavg(&load_avg, 1);
#else
load_avg = -1.0;
#endif
- pg_avail = sysconf(_SC_AVPHYS_PAGES);
-}
-
-SystemResourcesFlash::SystemResourcesFlash()
-{
- _sysres.getRUS(load_avg, pg_avail, utime, stime,
- elapsed, idrss, majflt);
+ pg_avail = availablePhysicalPages();
}
void
-SystemResourcesFlash::diff(const SystemResourcesFlash &pre)
+SystemResources::diff(const SystemResources &pre)
{
utime -= pre.utime;
stime -= pre.stime;
@@ -152,50 +118,63 @@ SystemResourcesFlash::diff(const SystemResourcesFlash &pre)
JournalRecord &
JournalRecord::operator<<(const IntSequence &s)
{
- operator<<("[");
+ operator<<('[');
for (int i = 0; i < s.size(); i++)
{
operator<<(s[i]);
if (i < s.size()-1)
- operator<<(",");
+ operator<<(',');
}
- operator<<("]");
+ operator<<(']');
return *this;
}
void
-JournalRecord::writePrefix(const SystemResourcesFlash &f)
+JournalRecord::writeFloatTabular(std::ostream &s, double d, int width)
{
- for (char & i : prefix)
- i = ' ';
- double mb = 1024*1024;
- sprintf(prefix, "%07.6g", f.elapsed);
- sprintf(prefix+7, ":%c%05d", recChar, ord);
- sprintf(prefix+14, ":%1.1f", f.load_avg);
- sprintf(prefix+18, ":%05.4g", f.pg_avail*_sysres.pageSize()/mb);
- sprintf(prefix+24, "%s", ": : ");
+ // Number of digits of integer part
+ int intdigits = std::max(static_cast<int>(std::floor(log10(d))+1), 1);
+
+ int prec = std::max(width - 1 - intdigits, 0);
+ s << std::fixed << std::setw(width) << std::setprecision(prec) << d;
+}
+
+void
+JournalRecord::writePrefix(const SystemResources &f)
+{
+ constexpr double mb = 1024*1024;
+ std::ostringstream s;
+ s << std::setfill('0');
+ writeFloatTabular(s, f.elapsed, 7);
+ s << ':' << recChar << std::setw(5) << ord << ':';
+ writeFloatTabular(s, f.load_avg, 3);
+ s << ':';
+ writeFloatTabular(s, f.pg_avail*SystemResources::pageSize()/mb, 5);
+ s << ": : ";
for (int i = 0; i < 2*journal.getDepth(); i++)
- prefix[i+33] = ' ';
- prefix[2*journal.getDepth()+33] = '\0';
+ s << ' ';
+ prefix = s.str();
}
void
-JournalRecordPair::writePrefixForEnd(const SystemResourcesFlash &f)
+JournalRecordPair::writePrefixForEnd(const SystemResources &f)
{
- for (char & i : prefix_end)
- i = ' ';
- double mb = 1024*1024;
- SystemResourcesFlash difnow;
+ constexpr double mb = 1024*1024;
+ SystemResources difnow;
difnow.diff(f);
- sprintf(prefix_end, "%07.6g", f.elapsed+difnow.elapsed);
- sprintf(prefix_end+7, ":E%05d", ord);
- sprintf(prefix_end+14, ":%1.1f", difnow.load_avg);
- sprintf(prefix_end+18, ":%05.4g", difnow.pg_avail*_sysres.pageSize()/mb);
- sprintf(prefix_end+24, ":%06.5g", difnow.majflt*_sysres.pageSize()/mb);
- sprintf(prefix_end+31, "%s", ": ");
+ std::ostringstream s;
+ s << std::setfill('0');
+ writeFloatTabular(s, f.elapsed+difnow.elapsed, 7);
+ s << ":E" << std::setw(5) << ord << ':';
+ writeFloatTabular(s, difnow.load_avg, 3);
+ s << ':';
+ writeFloatTabular(s, difnow.pg_avail*SystemResources::pageSize()/mb, 5);
+ s << ':';
+ writeFloatTabular(s, difnow.majflt*SystemResources::pageSize()/mb, 6);
+ s << ": ";
for (int i = 0; i < 2*journal.getDepth(); i++)
- prefix_end[i+33] = ' ';
- prefix_end[2*journal.getDepth()+33] = '\0';
+ s << ' ';
+ prefix_end = s.str();
}
JournalRecordPair::~JournalRecordPair()
@@ -222,39 +201,32 @@ endrec(JournalRecord &rec)
void
Journal::printHeader()
{
- (*this)<< "This is Dynare++, Copyright (C) 2004-2011, Ondra Kamenik\n"
- << "Dynare++ comes with ABSOLUTELY NO WARRANTY and is distributed under\n"
- << "GPL: modules integ, tl, kord, sylv, src, extern and documentation\n"
- << "LGPL: modules parser, utils\n"
- << " for GPL see http://www.gnu.org/licenses/gpl.html\n"
- << " for LGPL see http://www.gnu.org/licenses/lgpl.html\n"
- << "\n\n";
-
-#if !defined(__MINGW32__)
+ *this << "This is Dynare++, Copyright (C) 2004-2011, Ondra Kamenik\n"
+ << "Dynare++ comes with ABSOLUTELY NO WARRANTY and is distributed under\n"
+ << "GPL: modules integ, tl, kord, sylv, src, extern and documentation\n"
+ << "LGPL: modules parser, utils\n"
+ << " for GPL see http://www.gnu.org/licenses/gpl.html\n"
+ << " for LGPL see http://www.gnu.org/licenses/lgpl.html\n"
+ << "\n\n"
+ << "System info: ";
+#ifndef __MINGW32__
utsname info;
uname(&info);
- (*this)<< "System info: ";
- (*this)<< info.sysname << " " << info.release << " " << info.version << " ";
- (*this)<< info.machine << ", processors online: " << _sysres.onlineProcessors();
-
- (*this)<< "\n\nStart time: ";
- char ts[100];
- time_t curtime = time(nullptr);
- tm loctime;
- localtime_r(&curtime, &loctime);
- asctime_r(&loctime, ts);
- (*this)<< ts << "\n";
+ *this << info.sysname << " " << info.release << " " << info.version << " "
+ << info.machine << ", processors online: " << SystemResources::onlineProcessors();
#else
- (*this) << "System info: (not implemented for MINGW)\n";
- (*this) << "Start time: (not implemented for MINGW)\n\n";
+ *this << "(not implemented for MinGW)";
#endif
-
- (*this)<< " ------ elapsed time (seconds) \n";
- (*this)<< " | ------ record unique identifier \n";
- (*this)<< " | | ------ load average \n";
- (*this)<< " | | | ------ available memory (MB) \n";
- (*this)<< " | | | | ------ major faults (MB)\n";
- (*this)<< " | | | | | \n";
- (*this)<< " V V V V V \n";
- (*this)<< "\n";
+ *this << "\n\nStart time: ";
+ std::time_t t = std::time(nullptr);
+ *this << std::put_time(std::localtime(&t), "%c %Z")
+ << "\n\n"
+ << " ------ elapsed time (seconds) \n"
+ << " | ------ record unique identifier \n"
+ << " | | ------ load average \n"
+ << " | | | ------ available memory (MB) \n"
+ << " | | | | ------ major faults (MB)\n"
+ << " | | | | | \n"
+ << " V V V V V \n"
+ << "\n";
}
diff --git a/dynare++/kord/journal.hh b/dynare++/kord/journal.hh
index 70f3684244581f614f4c1a8405702190cf5d61bf..04a8afa9c56b559b4d1789602615e2c5b8c4debe 100644
--- a/dynare++/kord/journal.hh
+++ b/dynare++/kord/journal.hh
@@ -7,37 +7,35 @@
#include "int_sequence.hh"
-#include <sys/time.h>
-#include <cstdio>
-#include <cstring>
#include <iostream>
+#include <sstream>
#include <fstream>
+#include <string>
+#include <chrono>
-class SystemResources
+/* Implement static methods for accessing some system resources. An instance of
+ this class is a photograph of these resources at the time of instantiation. */
+struct SystemResources
{
- timeval start;
-public:
- SystemResources();
- static long int pageSize();
- static long int physicalPages();
- static long int onlineProcessors();
- static long int availableMemory();
- void getRUS(double &load_avg, long int &pg_avail, double &utime,
- double &stime, double &elapsed, long int &idrss,
- long int &majflt);
-};
+ // The starting time of the executable
+ static const std::chrono::time_point<std::chrono::high_resolution_clock> start;
+
+ static long pageSize();
+ static long physicalPages();
+ static long availablePhysicalPages();
+ static long onlineProcessors();
+ static long availableMemory();
-struct SystemResourcesFlash
-{
double load_avg;
- long int pg_avail;
+ long pg_avail;
double utime;
double stime;
double elapsed;
- long int idrss;
- long int majflt;
- SystemResourcesFlash();
- void diff(const SystemResourcesFlash &pre);
+ long idrss;
+ long majflt;
+
+ SystemResources();
+ void diff(const SystemResources &pre);
};
class Journal : public std::ofstream
@@ -45,7 +43,7 @@ class Journal : public std::ofstream
int ord;
int depth;
public:
- Journal(const char *fname)
+ explicit Journal(const std::string &fname)
: std::ofstream(fname), ord(0), depth(0)
{
printHeader();
@@ -82,10 +80,8 @@ public:
}
};
-#define MAXLEN 1000
-
class JournalRecord;
-JournalRecord&endrec(JournalRecord &);
+JournalRecord &endrec(JournalRecord &);
class JournalRecord
{
@@ -94,55 +90,74 @@ protected:
int ord;
public:
Journal &journal;
- char prefix[MAXLEN];
- char mes[MAXLEN];
- SystemResourcesFlash flash;
+ std::string prefix;
+ std::string mes;
+ SystemResources flash;
using _Tfunc = JournalRecord &(*)(JournalRecord &);
- JournalRecord(Journal &jr, char rc = 'M')
+ explicit JournalRecord(Journal &jr, char rc = 'M')
: recChar(rc), ord(jr.getOrd()), journal(jr)
{
- prefix[0] = '\0'; mes[0] = '\0'; writePrefix(flash);
+ writePrefix(flash);
}
- virtual ~JournalRecord()
- = default;
+ virtual ~JournalRecord() = default;
JournalRecord &operator<<(const IntSequence &s);
JournalRecord &
operator<<(_Tfunc f)
{
- (*f)(*this); return *this;
+ (*f)(*this);
+ return *this;
+ }
+ JournalRecord &
+ operator<<(char c)
+ {
+ mes += c;
+ return *this;
}
JournalRecord &
- operator<<(const char *s)
+ operator<<(const std::string &s)
{
- strcat(mes, s); return *this;
+ mes += s;
+ return *this;
}
JournalRecord &
operator<<(int i)
{
- sprintf(mes+strlen(mes), "%d", i); return *this;
+ mes += std::to_string(i);
+ return *this;
}
JournalRecord &
operator<<(double d)
{
- sprintf(mes+strlen(mes), "%f", d); return *this;
+ mes += std::to_string(d);
+ return *this;
}
protected:
- void writePrefix(const SystemResourcesFlash &f);
+ void writePrefix(const SystemResources &f);
+ /* Writes a floating point number as a field of exactly "width" characters
+ large. Note that the width will not be respected if the integer part is
+ too large. */
+ static void writeFloatTabular(std::ostream &s, double d, int width);
};
+/*
+ Constructs a "pair" of symmetric records with a RAII-like logic:
+ - when fed with "endrec", print the opening record
+ - subsequent records will have depth increased by 1
+ - when deleted, prints the symmetric closing record, and decrease depth
+*/
class JournalRecordPair : public JournalRecord
{
- char prefix_end[MAXLEN];
+ std::string prefix_end;
public:
- JournalRecordPair(Journal &jr)
+ explicit JournalRecordPair(Journal &jr)
: JournalRecord(jr, 'S')
{
- prefix_end[0] = '\0'; journal.incrementDepth();
+ journal.incrementDepth();
}
~JournalRecordPair() override;
private:
- void writePrefixForEnd(const SystemResourcesFlash &f);
+ void writePrefixForEnd(const SystemResources &f);
};
#endif
diff --git a/dynare++/src/dynare_exception.hh b/dynare++/src/dynare_exception.hh
index bbaf961e93a06e48f59ba3107c37cefb4d49dd5e..6f0beb049d075f6e613744d7d0beaf4e39ab13fb 100644
--- a/dynare++/src/dynare_exception.hh
+++ b/dynare++/src/dynare_exception.hh
@@ -6,6 +6,7 @@
#define DYNARE_EXCEPTION_H
#include <string>
+#include <cstring>
class DynareException
{
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index 5df949c55e5479fb721ef673a7439a31bc6274bf..696025194934fdbd83ee60bab2c468339a51ec8d 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -63,9 +63,9 @@ DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width, std
}
void
-copy_derivatives(mxArray *destin, const Symmetry &sym, const FGSContainer *derivs, const std::string &fieldname)
+copy_derivatives(mxArray *destin, const Symmetry &sym, const FGSContainer &derivs, const std::string &fieldname)
{
- const TwoDMatrix &x = derivs->get(sym);
+ const TwoDMatrix &x = derivs.get(sym);
int n = x.numRows();
int m = x.numCols();
mxArray *tmp = mxCreateDoubleMatrix(n, m, mxREAL);
@@ -286,7 +286,7 @@ extern "C" {
}
if (kOrder == 3 && nlhs > 5)
{
- const FGSContainer *derivs = app.get_rule_ders();
+ const FGSContainer &derivs = app.get_rule_ders();
const std::string fieldnames[] = {"gy", "gu", "gyy", "gyu", "guu", "gss",
"gyyy", "gyyu", "gyuu", "guuu", "gyss", "guss"};
// creates the char** expected by mxCreateStructMatrix()