Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • chskcau/dynare-doc-fixes
27 results
Select Git revision
Show changes
Show whitespace changes
Inline Side-by-side
Showing
with 0 additions and 4419 deletions
dynare++/kord/faa_di_bruno.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2005, Ondra Kamenik
// Faa Di Bruno evaluator
/* This defines a class which implements Faa Di Bruno Formula
$$\left[B_{s^k}\right]_{\alpha_1\ldots\alpha_l}=\left[f_{z^l}\right]_{\beta_1\ldots\beta_l}
\sum_{c\in M_{l,k}}\prod_{m=1}^l\left[z_{s^k(c_m)}\right]^{\beta_m}_{c_m(\alpha)}$$
where $s^k$ is a general symmetry of dimension $k$ and $z$ is a stack of
functions. */
#ifndef FAA_DI_BRUNO_H
#define FAA_DI_BRUNO_H
#include "journal.hh"
#include "stack_container.hh"
#include "t_container.hh"
#include "sparse_tensor.hh"
#include "gs_tensor.hh"
/* Nothing special here. See |@<|FaaDiBruno::calculate| folded sparse
code@>| for reason of having |magic_mult|. */
class FaaDiBruno
{
Journal &journal;
public:
FaaDiBruno(Journal &jr)
: journal(jr)
{
}
void calculate(const StackContainer<FGSTensor> &cont, const TensorContainer<FSSparseTensor> &f,
FGSTensor &out);
void calculate(const FoldedStackContainer &cont, const FGSContainer &g,
FGSTensor &out);
void calculate(const StackContainer<UGSTensor> &cont, const TensorContainer<FSSparseTensor> &f,
UGSTensor &out);
void calculate(const UnfoldedStackContainer &cont, const UGSContainer &g,
UGSTensor &out);
protected:
int estimRefinment(const TensorDimens &tdims, int nr, int l, int &avmem_mb, int &tmpmem_mb);
static double magic_mult;
};
#endif
dynare++/kord/first_order.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2004, Ondra Kamenik
#include "kord_exception.hh"
#include "first_order.hh"
#include <dynlapack.h>
double qz_criterium = 1.000001;
/* 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)
{
return (*alphar **alphar + *alphai **alphai < *beta **beta * qz_criterium * qz_criterium);
}
/* Here we solve the linear approximation. The result are the matrices
$g_{y^*}$ and $g_u$. The method solves the first derivatives of $g$ so
that the following equation would be true:
$$E_t[F(y^*_{t-1},u_t,u_{t+1},\sigma)] =
E_t[f(g^{**}(g^*(y_{t-1}^*,u_t,\sigma), u_{t+1}, \sigma), g(y_{t-1}^*,u_t,\sigma),
y^*_{t-1},u_t)]=0$$
where $f$ is a given system of equations.
It is known that $g_{y^*}$ is given by $F_{y^*}=0$, $g_u$ is given by
$F_u=0$, and $g_\sigma$ is zero. The only input to the method are the
derivatives |fd| of the system $f$, and partitioning of the vector $y$
(from object). */
void
FirstOrder::solve(const TwoDMatrix &fd)
{
JournalRecordPair pa(journal);
pa << "Recovering first order derivatives " << endrec;
::qz_criterium = FirstOrder::qz_criterium;
// **********************
// solve derivatives |gy|
// **********************
/* The derivatives $g_{y^*}$ are retrieved from the equation
$F_{y^*}=0$. The calculation proceeds as follows:
\orderedlist
\li For each variable appearing at both $t-1$ and $t-1$ we add a dummy
variable, so that the predetermined variables and forward looking would
be disjoint. This is, the matrix of the first derivatives of the
system written as:
$$\left[\matrix{f_{y^{**}_+}&f_{ys}&f_{yp}&f_{yb}&f_{yf}&f_{y^*_-}}\right],$$
where $f_{ys}$, $f_{yp}$, $f_{yb}$, and $f_{yf}$ are derivatives wrt
static, predetermined, both, forward looking at time $t$, is rewritten
to the matrix:
$$\left[
\matrix{f_{y^{**}_+}&f_{ys}&f_{yp}&f_{yb}&0&f_{yf}&f_{y^*_-}\cr
0 &0 &0 &I &-I&0 &0}
\right],$$
where the second line has number of rows equal to the number of both variables.
\li Next, provided that forward looking and predetermined are
disjoint, the equation $F_{y^*}=0$ is written as:
$$\left[f_+{y^{**}_+}\right]\left[g_{y^*}^{**}\right]\left[g_{y^*}^*\right]
+\left[f_{ys}\right]\left[g^s_{y^*}\right]
+\left[f_{y^*}\right]\left[g^*_{y^*}\right]
+\left[f_{y^{**}}\right]\left[g^{**}_{y^*}\right]+\left[f_{y^*_-}\right]=0$$
This is rewritten as
$$\left[\matrix{f_{y^*}&0&f_{y^{**}_+}}\right]
\left[\matrix{I\cr g^s_{y^*}\cr g^{**}_{y^*}}\right]\left[g_{y^*}^*\right]+
\left[\matrix{f_{y^*_-}&f_{ys}&f_{y^{**}}}\right]
\left[\matrix{I\cr g^s_{y^*}\cr g^{**}_{y^*}}\right]=0
$$
Now, in the above equation, there are the auxiliary variables standing
for copies of both variables at time $t+1$. This equation is then
rewritten as:
$$
\left[\matrix{f_{yp}&f_{yb}&0&f_{y^{**}_+}\cr 0&I&0&0}\right]
\left[\matrix{I\cr g^s_{y^*}\cr g^{**}_{y^*}}\right]\left[g_{y^*}^*\right]+
\left[\matrix{f_{y^*_-}&f_{ys}&0&f_{yf}\cr 0&0&-I&0}\right]
\left[\matrix{I\cr g^s_{y^*}\cr g^{**}_{y^*}}\right]=0
$$
The two matrices are denoted as $D$ and $-E$, so the equation takes the form:
$$D\left[\matrix{I\cr g^s_{y^*}\cr g^{**}_{y^*}}\right]\left[g_{y^*}^*\right]=
E\left[\matrix{I\cr g^s_{y^*}\cr g^{**}_{y^*}}\right]$$
\li Next we solve the equation by Generalized Schur decomposition:
$$
\left[\matrix{T_{11}&T_{12}\cr 0&T_{22}}\right]
\left[\matrix{Z_{11}^T&Z_{21}^T\cr Z_{12}^T&Z_{22}^T}\right]
\left[\matrix{I\cr X}\right]\left[g_{y^*}^*\right]=
\left[\matrix{S_{11}&S_{12}\cr 0&S_{22}}\right]
\left[\matrix{Z_{11}^T&Z_{21}^T\cr Z_{12}^T&Z_{22}^T}\right]
\left[\matrix{I\cr X}\right]
$$
We reorder the eigenvalue pair so that $S_{ii}/T_{ii}$ with modulus
less than one would be in the left-upper part.
\li The Blanchard--Kahn stability argument implies that the pairs
with modulus less that one will be in and only in $S_{11}/T_{11}$.
The exploding paths will be then eliminated when
$$
\left[\matrix{Z_{11}^T&Z_{21}^T\cr Z_{12}^T&Z_{22}^T}\right]
\left[\matrix{I\cr X}\right]=
\left[\matrix{Y\cr 0}\right]
$$
From this we have, $Y=Z_{11}^{-1}$, and $X=Z_{21}Y$, or equivalently
$X=-Z_{22}^{-T}Z_{12}^T$. From the equation, we get
$\left[g_{y^*}^*\right]=Y^{-1}T_{11}^{-1}S_{11}Y$, which is
$Z_{11}T_{11}^{-1}S_{11}Z_{11}^{-1}$.
\li We then copy the derivatives to storage |gy|. Note that the
derivatives of both variables are in $X$ and in
$\left[g_{y^*}^*\right]$, so we check whether the two submatrices are
the same. The difference is only numerical error.
\endorderedlist */
// setup submatrices of |f|
/* Here we setup submatrices of the derivatives |fd|. */
int off = 0;
ConstTwoDMatrix fyplus(fd, off, ypart.nyss());
off += ypart.nyss();
ConstTwoDMatrix fyszero(fd, off, ypart.nstat);
off += ypart.nstat;
ConstTwoDMatrix fypzero(fd, off, ypart.npred);
off += ypart.npred;
ConstTwoDMatrix fybzero(fd, off, ypart.nboth);
off += ypart.nboth;
ConstTwoDMatrix fyfzero(fd, off, ypart.nforw);
off += ypart.nforw;
ConstTwoDMatrix fymins(fd, off, ypart.nys());
off += ypart.nys();
ConstTwoDMatrix fuzero(fd, off, nu);
off += nu;
// form matrix $D$
lapack_int n = ypart.ny()+ypart.nboth;
TwoDMatrix matD(n, n);
matD.zeros();
matD.place(fypzero, 0, 0);
matD.place(fybzero, 0, ypart.npred);
matD.place(fyplus, 0, ypart.nys()+ypart.nstat);
for (int i = 0; i < ypart.nboth; i++)
matD.get(ypart.ny()+i, ypart.npred+i) = 1.0;
lapack_int ldb = matD.getLD();
// form matrix $E$
TwoDMatrix matE(n, n);
matE.zeros();
matE.place(fymins, 0, 0);
matE.place(fyszero, 0, ypart.nys());
matE.place(fyfzero, 0, ypart.nys()+ypart.nstat+ypart.nboth);
for (int i = 0; i < ypart.nboth; i++)
matE.get(ypart.ny()+i, ypart.nys()+ypart.nstat+i) = -1.0;
matE.mult(-1.0);
lapack_int lda = matE.getLD();
// solve generalized Schur
TwoDMatrix vsl(n, n);
TwoDMatrix vsr(n, n);
lapack_int ldvsl = vsl.getLD(), ldvsr = vsr.getLD();
lapack_int lwork = 100*n+16;
Vector work(lwork);
auto *bwork = new lapack_int[n];
lapack_int info;
lapack_int sdim2 = sdim;
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);
if (info)
{
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$. */
ConstGeneralMatrix z11(vsr, 0, 0, ypart.nys(), ypart.nys());
ConstGeneralMatrix z12(vsr, 0, ypart.nys(), ypart.nys(), n-ypart.nys());
ConstGeneralMatrix z21(vsr, ypart.nys(), 0, n-ypart.nys(), ypart.nys());
ConstGeneralMatrix z22(vsr, ypart.nys(), ypart.nys(), n-ypart.nys(), n-ypart.nys());
// calculate derivatives of static and forward
/* Here we calculate $X=-Z_{22}^{-T}Z_{12}^T$, where $X$ is |sfder| in the
code. */
GeneralMatrix sfder(z12, "transpose");
z22.multInvLeftTrans(sfder);
sfder.mult(-1);
// calculate derivatives of predetermined
/* Here we calculate
$g_{y^*}^*=Z_{11}T^{-1}_{11}S_{11}Z_{11}^{-1}
=Z_{11}T^{-1}_{11}(Z_{11}^{-T}S^T_{11})^T$. */
ConstGeneralMatrix s11(matE, 0, 0, ypart.nys(), ypart.nys());
ConstGeneralMatrix t11(matD, 0, 0, ypart.nys(), ypart.nys());
GeneralMatrix dumm(s11, "transpose");
z11.multInvLeftTrans(dumm);
GeneralMatrix preder(dumm, "transpose");
t11.multInvLeft(preder);
preder.multLeft(z11);
// copy derivatives to |gy|
gy.place(preder, ypart.nstat, 0);
GeneralMatrix sder(sfder, 0, 0, ypart.nstat, ypart.nys());
gy.place(sder, 0, 0);
GeneralMatrix fder(sfder, ypart.nstat+ypart.nboth, 0, ypart.nforw, ypart.nys());
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 bder2(preder, ypart.npred, 0, ypart.nboth, ypart.nys());
bder.add(-1, bder2);
b_error = bder.getData().getMax();
// **********************
// solve derivatives |gu|
// **********************
/* The equation $F_u=0$ can be written as
$$
\left[f_{y^{**}_+}\right]\left[g^{**}_{y^*}\right]\left[g_u^*\right]+
\left[f_y\right]\left[g_u\right]+\left[f_u\right]=0
$$
and rewritten as
$$
\left[f_y +
\left[\matrix{0&f_{y^{**}_+}g^{**}_{y^*}&0}\right]\right]g_u=f_u
$$
This is exactly done here. The matrix
$\left[f_y +\left[\matrix{0&f_{y^{**}_+}g^{**}_{y^*}&0}\right]\right]$ is |matA|
in the code. */
GeneralMatrix matA(ypart.ny(), ypart.ny());
matA.zeros();
ConstGeneralMatrix gss(gy, ypart.nstat+ypart.npred, 0, ypart.nyss(), ypart.nys());
GeneralMatrix aux(fyplus, gss);
matA.place(aux, 0, ypart.nstat);
ConstGeneralMatrix fyzero(fd, 0, ypart.nyss(), ypart.ny(), ypart.ny());
matA.add(1.0, fyzero);
gu.zeros();
gu.add(-1.0, fuzero);
ConstGeneralMatrix(matA).multInvLeft(gu);
journalEigs();
if (!gy.isFinite() || !gu.isFinite())
{
throw KordException(__FILE__, __LINE__,
"NaN or Inf asserted in first order derivatives in FirstOrder::solve");
}
}
void
FirstOrder::journalEigs()
{
if (bk_cond)
{
JournalRecord jr(journal);
jr << "Blanchard-Kahn conditition satisfied, model stable" << endrec;
}
else
{
JournalRecord jr(journal);
jr << "Blanchard-Kahn condition not satisfied, model not stable: sdim=" << sdim
<< " " << "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;
}
}
}
dynare++/kord/first_order.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2004, Ondra Kamenik
// First order at deterministic steady
#ifndef FIRST_ORDER_H
#define FIRST_ORDER_H
#include "korder.hh"
template<int>
class FirstOrderDerivs;
class FirstOrder
{
template <int>
friend class FirstOrderDerivs;
PartitionY ypart;
int nu;
TwoDMatrix gy;
TwoDMatrix gu;
bool bk_cond;
double b_error;
int sdim;
Vector alphar;
Vector alphai;
Vector beta;
double qz_criterium;
Journal &journal;
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)
: ypart(num_stat, num_pred, num_both, num_forw),
nu(num_u),
gy(ypart.ny(), ypart.nys()),
gu(ypart.ny(), nu),
alphar(ypart.ny()+ypart.nboth),
alphai(ypart.ny()+ypart.nboth),
beta(ypart.ny()+ypart.nboth),
qz_criterium(qz_crit),
journal(jr)
{
solve(FFSTensor(f));
}
bool
isStable() const
{
return bk_cond;
}
const TwoDMatrix &
getGy() const
{
return gy;
}
const TwoDMatrix &
getGu() const
{
return gu;
}
protected:
void solve(const TwoDMatrix &f);
void journalEigs();
};
/* This class only converts the derivatives $g_{y^*}$ and $g_u$ to a
folded or unfolded container. */
template <int t>
class FirstOrderDerivs : public ctraits<t>::Tg
{
public:
FirstOrderDerivs(const FirstOrder &fo)
: ctraits<t>::Tg(4)
{
IntSequence nvs(4);
nvs[0] = fo.ypart.nys(); nvs[1] = fo.nu; nvs[2] = fo.nu; nvs[3] = 1;
_Ttensor *ten = new _Ttensor(fo.ypart.ny(), TensorDimens(Symmetry(1, 0, 0, 0), nvs));
ten->zeros(); ten->add(1.0, fo.gy);
this->insert(ten);
ten = new _Ttensor(fo.ypart.ny(), TensorDimens(Symmetry(0, 1, 0, 0), nvs));
ten->zeros(); ten->add(1.0, fo.gu);
this->insert(ten);
}
};
#endif
dynare++/kord/global_check.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2005, Ondra Kamenik
#include "SymSchurDecomp.hh"
#include "global_check.hh"
#include "smolyak.hh"
#include "product.hh"
#include "quasi_mcarlo.hh"
#include <cmath>
/* Here we just set a reference to the approximation, and create a new
|DynamicModel|. */
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)
{
}
ResidFunction::ResidFunction(const ResidFunction &rf)
: VectorFunction(rf), approx(rf.approx), model(rf.model->clone()),
yplus(nullptr), ystar(nullptr), u(nullptr), hss(nullptr)
{
if (rf.yplus)
yplus = new Vector(*(rf.yplus));
if (rf.ystar)
ystar = new Vector(*(rf.ystar));
if (rf.u)
u = new 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;
}
/* This sets $y^*$ and $u$. We have to create |ystar|, |u|, |yplus| and
|hss|. */
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());
approx.getFoldDecisionRule().evaluate(DecisionRule::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
construction of |hss| and will be used in |const| context. So this
dirty thing is safe.
Note, that there is always a folded decision rule in |Approximation|. */
union {const FoldDecisionRule *c; FoldDecisionRule *n;} dr;
dr.c = &(approx.getFoldDecisionRule());
FTensorPolynomial dr_ss(model->nstat()+model->npred(), model->nboth()+model->nforw(),
*(dr.n));
// 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(),
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(),
model->nboth()+model->nforw());
if (hss->check(Symmetry(0)))
{
hss->get(Symmetry(0))->getData().add(1.0, ysteady_ss);
}
else
{
auto *ten = new FFSTensor(hss->nrows(), hss->nvars(), 0);
ten->getData() = ysteady_ss;
hss->insert(ten);
}
}
/* Here we evaluate the residual $F(y^*,u,u')$. We have to evaluate |hss|
for $u'=$|point| and then we evaluate the system $f$. */
void
ResidFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{
KORD_RAISE_IF(point.length() != hss->nvars(),
"Wrong dimension of input vector in ResidFunction::eval");
KORD_RAISE_IF(out.length() != model->numeq(),
"Wrong dimension of output vector in ResidFunction::eval");
Vector yss(hss->nrows());
hss->evalHorner(yss, point);
model->evaluateSystem(out, *ystar, *yplus, yss, *u);
}
/* This checks the $E[F(y^*,u,u')]$ for a given $y^*$ and $u$ by
integrating with a given quadrature. Note that the input |ys| is $y^*$
not whole $y$. */
void
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);
quad.integrate(vfs, level, out);
}
/* This method is a bulk version of |@<|GlobalChecker::check| vector
code@>|. It decides between Smolyak and product quadrature according
to |max_evals| constraint.
Note that |y| can be either full (all endogenous variables including
static and forward looking), or just $y^*$ (state variables). The
method is able to recognize it. */
void
GlobalChecker::check(int max_evals, const ConstTwoDMatrix &y,
const ConstTwoDMatrix &x, TwoDMatrix &out)
{
JournalRecordPair pa(journal);
pa << "Checking approximation error for " << y.ncols()
<< " states with at most " << max_evals << " evaluations" << endrec;
// decide about type of quadrature
GaussHermite gh;
SmolyakQuadrature dummy_sq(model.nexog(), 1, gh);
int smol_evals;
int smol_level;
dummy_sq.designLevelForEvals(max_evals, smol_level, smol_evals);
ProductQuadrature dummy_pq(model.nexog(), gh);
int prod_evals;
int prod_level;
dummy_pq.designLevelForEvals(max_evals, prod_level, prod_evals);
bool take_smolyak = (smol_evals < prod_evals) && (smol_level >= prod_level-1);
Quadrature *quad;
int lev;
// create the quadrature and report the decision
if (take_smolyak)
{
quad = new SmolyakQuadrature(model.nexog(), smol_level, gh);
lev = smol_level;
JournalRecord rec(journal);
rec << "Selected Smolyak (level,evals)=(" << smol_level << ","
<< smol_evals << ") over product (" << prod_level << ","
<< prod_evals << ")" << endrec;
}
else
{
quad = new ProductQuadrature(model.nexog(), gh);
lev = prod_level;
JournalRecord rec(journal);
rec << "Selected product (level,evals)=(" << prod_level << ","
<< prod_evals << ") over Smolyak (" << smol_level << ","
<< smol_evals << ")" << endrec;
}
// check all column of |y| and |x|
int first_row = (y.nrows() == model.numeq()) ? model.nstat() : 0;
ConstTwoDMatrix ysmat(y, first_row, 0, model.npred()+model.nboth(), y.ncols());
for (int j = 0; j < y.ncols(); j++)
{
ConstVector yj{ysmat.getCol(j)};
ConstVector xj{x.getCol(j)};
Vector outj{out.getCol(j)};
check(*quad, lev, yj, xj, outj);
}
delete quad;
}
/* This method checks an error of the approximation by evaluating
residual $E[F(y^*,u,u')\vert y^*,u]$ for $y^*$ being the steady state, and
changing $u$. We go through all elements of $u$ and vary them from
$-mult\cdot\sigma$ to $mult\cdot\sigma$ in |m| steps. */
void
GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const char *prefix,
int m, double mult, int max_evals)
{
JournalRecordPair pa(journal);
pa << "Calculating errors along shocks +/- "
<< mult << " std errors, granularity " << m << endrec;
// setup |y_mat| of steady states for checking
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();
}
// setup |exo_mat| for checking
TwoDMatrix exo_mat(model.nexog(), 2*m *model.nexog()+1);
exo_mat.zeros();
for (int ishock = 0; ishock < model.nexog(); ishock++)
{
double max_sigma = sqrt(model.getVcov().get(ishock, ishock));
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;
}
}
TwoDMatrix errors(model.numeq(), 2*m *model.nexog()+1);
check(max_evals, y_mat, exo_mat, errors);
// report errors along shock and save them
TwoDMatrix res(model.nexog(), 2*m+1);
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;
Vector error{err_out.getCol(j)};
if (j != m)
{
if (j < m)
jj = 1 + 2*m*ishock+j;
else
jj = 1 + 2*m*ishock+j-1;
ConstVector coljj{errors.getCol(jj)};
error = coljj;
}
else
{
jj = 0;
error = err0;
}
JournalRecord rec1(journal);
sprintf(erbuf, "%12.7g ", error.getMax());
rec1 << shock << " " << exo_mat.get(ishock, jj)
<< "\t" << erbuf << endrec;
}
char tmp[100];
sprintf(tmp, "%s_shock_%s_errors", prefix, model.getExogNames().getName(ishock));
err_out.writeMat(fd, tmp);
}
}
/* This method checks errors on ellipse of endogenous states
(predetermined variables). The ellipse is shaped according to
covariance matrix of endogenous variables based on the first order
approximation and scaled by |mult|. The points on the
ellipse are chosen as polar images of the low discrepancy grid in a
cube.
The method works as follows: First we calculate symmetric Schur factor of
covariance matrix of the states. Second we generate low discrepancy
points on the unit sphere. Third we transform the sphere with the
variance-covariance matrix factor and multiplier |mult| and initialize
matrix of $u_t$ to zeros. Fourth we run the |check| method and save
the results. */
void
GlobalChecker::checkOnEllipseAndSave(mat_t *fd, const char *prefix,
int m, double mult, int max_evals)
{
JournalRecordPair pa(journal);
pa << "Calculating errors at " << m
<< " ellipse points scaled by " << mult << endrec;
// make factor of covariance of variables
/* 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(),
model.npred()+model.nboth(), model.npred()+model.nboth());
delete ycov;
SymSchurDecomp ssd(ycovpred);
ssd.correctDefinitness(1.e-05);
TwoDMatrix ycovfac(ycovpred.nrows(), ycovpred.ncols());
KORD_RAISE_IF(!ssd.isPositiveSemidefinite(),
"Covariance matrix of the states not positive \
semidefinite in GlobalChecker::checkOnEllipseAndSave");
ssd.getFactor(ycovfac);
// put low discrepancy sphere points to |ymat|
/* Here we first calculate dimension |d| of the sphere, which is a
number of state variables minus one. We go through the |d|-dimensional
cube $\langle 0,1\rangle^d$ by |QMCarloCubeQuadrature| and make a
polar transformation to the sphere. The polar transformation $f^i$ can
be written recursively wrt. the dimension $i$ as:
$$\eqalign{
f^0() &= \left[1\right]\cr
f^i(x_1,\ldots,x_i) &=
\left[\matrix{cos(2\pi x_i)\cdot f^{i-1}(x_1,\ldots,x_{i-1})\cr sin(2\pi x_i)}\right]
}$$ */
int d = model.npred()+model.nboth()-1;
TwoDMatrix ymat(model.npred()+model.nboth(), (d == 0) ? 2 : m);
if (d == 0)
{
ymat.get(0, 0) = 1;
ymat.get(0, 1) = -1;
}
else
{
int icol = 0;
ReversePerScheme ps;
QMCarloCubeQuadrature qmc(d, m, ps);
qmcpit beg = qmc.start(m);
qmcpit end = qmc.end(m);
for (qmcpit run = beg; run != end; ++run, icol++)
{
Vector ycol{ymat.getCol(icol)};
Vector x(run.point());
x.mult(2*M_PI);
ycol[0] = 1;
for (int i = 0; i < d; i++)
{
Vector subsphere(ycol, 0, i+1);
subsphere.mult(cos(x[i]));
ycol[i+1] = sin(x[i]);
}
}
}
// transform sphere |ymat| and prepare |umat| for checking
/* Here we multiply the sphere points in |ymat| with the Cholesky
factor to obtain the ellipse, scale the ellipse by the given |mult|,
and initialize matrix of shocks |umat| to zero. */
TwoDMatrix umat(model.nexog(), ymat.ncols());
umat.zeros();
ymat.mult(mult);
ymat.multLeft(ycovfac);
ConstVector ys(model.getSteady(), model.nstat(),
model.npred()+model.nboth());
for (int icol = 0; icol < ymat.ncols(); icol++)
{
Vector ycol{ymat.getCol(icol)};
ycol.add(1.0, ys);
}
// check on ellipse and save
/* Here we check the points and save the results to MAT-4 file. */
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);
}
/* 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,
int m, int max_evals)
{
JournalRecordPair pa(journal);
pa << "Calculating errors at " << m
<< " simulated points" << endrec;
RandomShockRealization sr(model.getVcov(), system_random_generator.int_uniform());
TwoDMatrix *y = approx.getFoldDecisionRule().simulate(DecisionRule::horner,
m, model.getSteady(), sr);
TwoDMatrix x(model.nexog(), m);
x.zeros();
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);
delete y;
}
dynare++/kord/global_check.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2005, Ondra Kamenik
// Global check
/* The purpose of this file is to provide classes for checking error of
approximation. If $y_t=g(y^*_{t-1},u)$ is an approximate solution,
then we check for the error of residuals of the system equations. Let
$F(y^*,u,u')=f(g^{**}(g^*(y^*,u'),u),g(y^*,u),y^*,u)$, then we
calculate integral
$$E[F(y^*,u,u')]$$@q'@>
which we want to be zero for all $y^*$, and $u$.
There are a few possibilities how and where the integral is
evaluated. Currently we offer the following:
\numberedlist
\li Along shocks. The $y^*$ is set to steady state, and $u$ is set to
zero but one element is going from minus through plus shocks in few
steps. The user gives the scaling factor, for instance interval
$\langle<-3\sigma,3\sigma\rangle$ (where $sigma$ is a standard error
of the shock), and a number of steps. This is repeated for each shock
(element of the $u$ vector).
\li Along simulation. Some random simulation is run, and for each
realization of $y^*$ and $u$ along the path we evaluate the residual.
\li On ellipse. Let $V=AA^T$ be a covariance matrix of the
predetermined variables $y^*$ based on linear approximation, then we
calculate integral for points on the ellipse $\{Ax\vert\, \Vert
x\Vert_2=1\}$. The points are selected by means of low discrepancy
method and polar transformation. The shock $u$ are zeros.
\li Unconditional distribution.
\endnumberedlist */
#ifndef GLOBAL_CHECK_H
#define GLOBAL_CHECK_H
#include <matio.h>
#include <memory>
#include "vector_function.hh"
#include "quadrature.hh"
#include "dynamic_model.hh"
#include "journal.hh"
#include "approximation.hh"
/* This is a class for implementing |VectorFunction| interface
evaluating the residual of equations, this is
$$F(y^*,u,u')=f(g^{**}(g^*(y^*,u),u'),y^*,u)$$
is written as a function of $u'$.
When the object is constructed, one has to specify $(y^*,u)$, this is
done by |setYU| method. The object has basically two states. One is
after construction and before call to |setYU|. The second is after
call |setYU|. We distinguish between the two states, an object in the
second state contains |yplus|, |ystar|, |u|, and |hss|.
The vector |yplus| is $g^*(y^*,u)$. |ystar| is $y^*$, and polynomial
|hss| is partially evaluated $g^**(yplus, u)$.
The pointer to |DynamicModel| is important, since the |DynamicModel|
evaluates the function $f$. When copying the object, we have to make
also a copy of |DynamicModel|. */
class ResidFunction : public VectorFunction
{
protected:
const Approximation &approx;
DynamicModel *model;
Vector *yplus;
Vector *ystar;
Vector *u;
FTensorPolynomial *hss;
public:
ResidFunction(const Approximation &app);
ResidFunction(const ResidFunction &rf);
~ResidFunction() override;
std::unique_ptr<VectorFunction>
clone() const override
{
return std::make_unique<ResidFunction>(*this);
}
void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
void setYU(const ConstVector &ys, const ConstVector &xx);
};
/* This is a |ResidFunction| wrapped with |GaussConverterFunction|. */
class GResidFunction : public GaussConverterFunction
{
public:
GResidFunction(const Approximation &app)
: GaussConverterFunction(std::make_unique<ResidFunction>(app), app.getModel().getVcov())
{
}
GResidFunction(const GResidFunction &rf) = default;
~GResidFunction() override = default;
std::unique_ptr<VectorFunction>
clone() const override
{
return std::make_unique<GResidFunction>(*this);
}
void
setYU(const ConstVector &ys, const ConstVector &xx)
{
((ResidFunction *) func)->setYU(ys, xx);
}
};
/* This is a class encapsulating checking algorithms. Its core routine
is |check|, which calculates integral $E[F(y^*,u,u')\vert y^*,u]$ for
given realizations of $y^*$ and $u$. The both are given in
matrices. The methods checking along shocks, on ellipse and anlong a
simulation path, just fill the matrices and call the core |check|.
The method |checkUnconditionalAndSave| evaluates unconditional
$E[F(y,u,u')]$.
The object also maintains a set of |GResidFunction| functions |vfs| in
order to save (possibly expensive) copying of |DynamicModel|s. */
class GlobalChecker
{
const Approximation &approx;
const DynamicModel &model;
Journal &journal;
GResidFunction rf;
VectorFunctionSet vfs;
public:
GlobalChecker(const Approximation &app, int n, Journal &jr)
: approx(app), model(approx.getModel()), journal(jr),
rf(approx), vfs(rf, n)
{
}
void check(int max_evals, const ConstTwoDMatrix &y,
const ConstTwoDMatrix &x, TwoDMatrix &out);
void checkAlongShocksAndSave(mat_t *fd, const char *prefix,
int m, double mult, int max_evals);
void checkOnEllipseAndSave(mat_t *fd, const char *prefix,
int m, double mult, int max_evals);
void checkAlongSimulationAndSave(mat_t *fd, const char *prefix,
int m, int max_evals);
void checkUnconditionalAndSave(mat_t *fd, const char *prefix,
int m, int max_evals);
protected:
void check(const Quadrature &quad, int level,
const ConstVector &y, const ConstVector &x, Vector &out);
};
/* Signalled resid function. Not implemented yet. todo: */
class ResidFunctionSig : public ResidFunction
{
public:
ResidFunctionSig(const Approximation &app, const Vector &ys, const Vector &xx);
};
#endif
dynare++/kord/journal.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright (C) 2004-2011, Ondra Kamenik
#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 <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 NOMINMAX
# 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;
}
}
#endif
#if defined(__APPLE__)
# define _SC_PHYS_PAGES 2
# define _SC_AVPHYS_PAGES 3
#endif
SystemResources::SystemResources()
{
gettimeofday(&start, nullptr);
}
long int
SystemResources::pageSize()
{
return sysconf(_SC_PAGESIZE);
}
long int
SystemResources::physicalPages()
{
return sysconf(_SC_PHYS_PAGES);
}
long int
SystemResources::onlineProcessors()
{
return sysconf(_SC_NPROCESSORS_ONLN);
}
long int
SystemResources::availableMemory()
{
return pageSize()*sysconf(_SC_AVPHYS_PAGES);
}
/* 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)
{
struct timeval now;
gettimeofday(&now, nullptr);
elapsed = now.tv_sec-start.tv_sec + (now.tv_usec-start.tv_usec)*1.0e-6;
#if !defined(__MINGW32__)
struct rusage rus;
getrusage(RUSAGE_SELF, &rus);
utime = rus.ru_utime.tv_sec+rus.ru_utime.tv_usec*1.0e-6;
stime = rus.ru_stime.tv_sec+rus.ru_stime.tv_usec*1.0e-6;
idrss = rus.ru_idrss;
majflt = rus.ru_majflt;
#else
utime = -1.0;
stime = -1.0;
idrss = -1;
majflt = -1;
#endif
#define MINGCYGTMP (!defined(__MINGW32__) && !defined(__CYGWIN32__) && !defined(__CYGWIN__))
#define MINGCYG (MINGCYGTMP && !defined(__MINGW64__) && !defined(__CYGWIN64__))
#if MINGCYG
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);
}
void
SystemResourcesFlash::diff(const SystemResourcesFlash &pre)
{
utime -= pre.utime;
stime -= pre.stime;
elapsed -= pre.elapsed;
idrss -= pre.idrss;
majflt -= pre.majflt;
}
// |JournalRecord::operator<<| symmetry code
JournalRecord &
JournalRecord::operator<<(const IntSequence &s)
{
operator<<("[");
for (int i = 0; i < s.size(); i++)
{
operator<<(s[i]);
if (i < s.size()-1)
operator<<(",");
}
operator<<("]");
return *this;
}
void
JournalRecord::writePrefix(const SystemResourcesFlash &f)
{
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", ": : ");
for (int i = 0; i < 2*journal.getDepth(); i++)
prefix[i+33] = ' ';
prefix[2*journal.getDepth()+33] = '\0';
}
void
JournalRecordPair::writePrefixForEnd(const SystemResourcesFlash &f)
{
for (char & i : prefix_end)
i = ' ';
double mb = 1024*1024;
SystemResourcesFlash 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", ": ");
for (int i = 0; i < 2*journal.getDepth(); i++)
prefix_end[i+33] = ' ';
prefix_end[2*journal.getDepth()+33] = '\0';
}
JournalRecordPair::~JournalRecordPair()
{
journal.decrementDepth();
writePrefixForEnd(flash);
journal << prefix_end;
journal << mes;
journal << std::endl;
journal.flush();
}
JournalRecord &
endrec(JournalRecord &rec)
{
rec.journal << rec.prefix;
rec.journal << rec.mes;
rec.journal << std::endl;
rec.journal.flush();
rec.journal.incrementOrd();
return 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__)
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";
#else
(*this) << "System info: (not implemented for MINGW)\n";
(*this) << "Start time: (not implemented for MINGW)\n\n";
#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";
}
dynare++/kord/journal.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2004, Ondra Kamenik
// Resource usage journal
#ifndef JOURNAL_H
#define JOURNAL_H
#include "int_sequence.hh"
#include <sys/time.h>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <fstream>
class 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);
};
struct SystemResourcesFlash
{
double load_avg;
long int pg_avail;
double utime;
double stime;
double elapsed;
long int idrss;
long int majflt;
SystemResourcesFlash();
void diff(const SystemResourcesFlash &pre);
};
class Journal : public std::ofstream
{
int ord;
int depth;
public:
Journal(const char *fname)
: std::ofstream(fname), ord(0), depth(0)
{
printHeader();
}
~Journal() override
{
flush();
}
void printHeader();
void
incrementOrd()
{
ord++;
}
int
getOrd() const
{
return ord;
}
void
incrementDepth()
{
depth++;
}
void
decrementDepth()
{
depth--;
}
int
getDepth() const
{
return depth;
}
};
#define MAXLEN 1000
class JournalRecord;
JournalRecord&endrec(JournalRecord &);
class JournalRecord
{
protected:
char recChar;
int ord;
public:
Journal &journal;
char prefix[MAXLEN];
char mes[MAXLEN];
SystemResourcesFlash flash;
using _Tfunc = JournalRecord &(*)(JournalRecord &);
JournalRecord(Journal &jr, char rc = 'M')
: recChar(rc), ord(jr.getOrd()), journal(jr)
{
prefix[0] = '\0'; mes[0] = '\0'; writePrefix(flash);
}
virtual ~JournalRecord()
= default;
JournalRecord &operator<<(const IntSequence &s);
JournalRecord &
operator<<(_Tfunc f)
{
(*f)(*this); return *this;
}
JournalRecord &
operator<<(const char *s)
{
strcat(mes, s); return *this;
}
JournalRecord &
operator<<(int i)
{
sprintf(mes+strlen(mes), "%d", i); return *this;
}
JournalRecord &
operator<<(double d)
{
sprintf(mes+strlen(mes), "%f", d); return *this;
}
protected:
void writePrefix(const SystemResourcesFlash &f);
};
class JournalRecordPair : public JournalRecord
{
char prefix_end[MAXLEN];
public:
JournalRecordPair(Journal &jr)
: JournalRecord(jr, 'S')
{
prefix_end[0] = '\0'; journal.incrementDepth();
}
~JournalRecordPair() override;
private:
void writePrefixForEnd(const SystemResourcesFlash &f);
};
#endif
dynare++/kord/kord_exception.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2005, Ondra Kamenik
// Exception
/* This is a simple code defining an exception and two convenience macros. */
#ifndef KORD_EXCEPTION_H
#define KORD_EXCEPTION_H
#include <cstring>
#include <cstdio>
#define KORD_RAISE(mes) \
throw KordException(__FILE__, __LINE__, mes);
#define KORD_RAISE_IF(expr, mes) \
if (expr) throw KordException(__FILE__, __LINE__, mes);
#define KORD_RAISE_X(mes, c) \
throw KordException(__FILE__, __LINE__, mes, c);
#define KORD_RAISE_IF_X(expr, mes, c) \
if (expr) throw KordException(__FILE__, __LINE__, mes, c);
class KordException
{
protected:
char fname[50];
int lnum;
char message[500];
int cd;
public:
KordException(const char *f, int l, const char *mes, int c = 255)
{
strncpy(fname, f, 50); fname[49] = '\0';
strncpy(message, mes, 500); message[499] = '\0';
lnum = l;
cd = c;
}
virtual ~KordException()
= default;
virtual void
print() const
{
printf("At %s:%d:(%d):%s\n", fname, lnum, cd, message);
}
virtual int
code() const
{
return cd;
}
const char *
get_message() const
{
return message;
}
};
// |KordException| error code definitions
#define KORD_FP_NOT_CONV 254
#define KORD_FP_NOT_FINITE 253
#define KORD_MD_NOT_STABLE 252
#endif
dynare++/kord/korder.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2004, Ondra Kamenik
#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. */
void
PLUMatrix::calcPLU()
{
lapack_int info;
lapack_int rows = nrows(), lda = ld;
inv = (const Vector &) getData();
dgetrf(&rows, &rows, inv.base(), &lda, ipiv, &info);
}
/* Here we just call the LAPACK machinery to multiply by the inverse. */
void
PLUMatrix::multInv(TwoDMatrix &m) const
{
KORD_RAISE_IF(m.nrows() != ncols(),
"The matrix is not square in PLUMatrix::multInv");
lapack_int info;
lapack_int lda = ld;
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);
KORD_RAISE_IF(info != 0,
"Info!=0 in PLUMatrix::multInv");
}
/* Here we construct the matrix $A$. Its dimension is |ny|, and it is
$$A=\left[f_{y}\right]+
\left[0 \left[f_{y^{**}_+}\right]\cdot\left[g^{**}_{y^*}\right] 0\right]$$,
where the first zero spans |nstat| columns, and last zero spans
|nforw| columns. */
MatrixA::MatrixA(const FSSparseTensor &f, const IntSequence &ss,
const TwoDMatrix &gy, const PartitionY &ypart)
: PLUMatrix(ypart.ny())
{
zeros();
IntSequence c(1); c[0] = 1;
FGSTensor f_y(f, ss, c, TensorDimens(ss, c));
add(1.0, f_y);
ConstTwoDMatrix gss_ys(ypart.nstat+ypart.npred, ypart.nyss(), gy);
c[0] = 0;
FGSTensor f_yss(f, ss, c, TensorDimens(ss, c));
TwoDMatrix sub(*this, ypart.nstat, ypart.nys());
sub.multAndAdd(ConstTwoDMatrix(f_yss), gss_ys);
calcPLU();
}
/* Here we construct the matrix $S$. Its dimension is |ny|, and it is
$$S=\left[f_{y}\right]+
\left[0\quad\left[f_{y^{**}_+}\right]\cdot\left[g^{**}_{y^*}\right]\quad
0\right]+ \left[0\quad 0\quad\left[f_{y^{**}_+}\right]\right]$$
It is, in fact, the matrix $A$ plus the third summand. The first zero
in the summand spans |nstat| columns, the second zero spans |npred|
columns. */
MatrixS::MatrixS(const FSSparseTensor &f, const IntSequence &ss,
const TwoDMatrix &gy, const PartitionY &ypart)
: PLUMatrix(ypart.ny())
{
zeros();
IntSequence c(1); c[0] = 1;
FGSTensor f_y(f, ss, c, TensorDimens(ss, c));
add(1.0, f_y);
ConstTwoDMatrix gss_ys(ypart.nstat+ypart.npred, ypart.nyss(), gy);
c[0] = 0;
FGSTensor f_yss(f, ss, c, TensorDimens(ss, c));
TwoDMatrix sub(*this, ypart.nstat, ypart.nys());
sub.multAndAdd(ConstTwoDMatrix(f_yss), gss_ys);
TwoDMatrix sub2(*this, ypart.nstat+ypart.npred, ypart.nyss());
sub2.add(1.0, f_yss);
calcPLU();
}
// |KOrder| member access method specializations
/* These are the specializations of container access methods. Nothing
interesting here. */
template<>
ctraits<KOrder::unfold>::Tg& KOrder::g<KOrder::unfold>()
{
return _ug;
}
template<>
const ctraits<KOrder::unfold>::Tg& KOrder::g<KOrder::unfold>() const
{ return _ug;}
template<>
ctraits<KOrder::fold>::Tg& KOrder::g<KOrder::fold>()
{
return _fg;
}
template<>
const ctraits<KOrder::fold>::Tg& KOrder::g<KOrder::fold>() const
{ return _fg;}
template<>
ctraits<KOrder::unfold>::Tgs& KOrder::gs<KOrder::unfold>()
{
return _ugs;
}
template<>
const ctraits<KOrder::unfold>::Tgs& KOrder::gs<KOrder::unfold>() const
{ return _ugs;}
template<>
ctraits<KOrder::fold>::Tgs& KOrder::gs<KOrder::fold>()
{
return _fgs;
}
template<>
const ctraits<KOrder::fold>::Tgs& KOrder::gs<KOrder::fold>() const
{ return _fgs;}
template<>
ctraits<KOrder::unfold>::Tgss& KOrder::gss<KOrder::unfold>()
{
return _ugss;
}
template<>
const ctraits<KOrder::unfold>::Tgss& KOrder::gss<KOrder::unfold>() const
{ return _ugss;}
template<>
ctraits<KOrder::fold>::Tgss& KOrder::gss<KOrder::fold>()
{
return _fgss;
}
template<>
const ctraits<KOrder::fold>::Tgss& KOrder::gss<KOrder::fold>() const
{ return _fgss;}
template<>
ctraits<KOrder::unfold>::TG& KOrder::G<KOrder::unfold>()
{
return _uG;
}
template<>
const ctraits<KOrder::unfold>::TG& KOrder::G<KOrder::unfold>() const
{ return _uG;}
template<>
ctraits<KOrder::fold>::TG& KOrder::G<KOrder::fold>()
{
return _fG;
}
template<>
const ctraits<KOrder::fold>::TG& KOrder::G<KOrder::fold>() const
{ return _fG;}
template<>
ctraits<KOrder::unfold>::TZstack& KOrder::Zstack<KOrder::unfold>()
{
return _uZstack;
}
template<>
const ctraits<KOrder::unfold>::TZstack& KOrder::Zstack<KOrder::unfold>() const
{ return _uZstack;}
template<>
ctraits<KOrder::fold>::TZstack& KOrder::Zstack<KOrder::fold>()
{
return _fZstack;
}
template<>
const ctraits<KOrder::fold>::TZstack& KOrder::Zstack<KOrder::fold>() const
{ return _fZstack;}
template<>
ctraits<KOrder::unfold>::TGstack& KOrder::Gstack<KOrder::unfold>()
{
return _uGstack;
}
template<>
const ctraits<KOrder::unfold>::TGstack& KOrder::Gstack<KOrder::unfold>() const
{ return _uGstack;}
template<>
ctraits<KOrder::fold>::TGstack& KOrder::Gstack<KOrder::fold>()
{
return _fGstack;
}
template<>
const ctraits<KOrder::fold>::TGstack& KOrder::Gstack<KOrder::fold>() const
{ return _fGstack;}
template<>
ctraits<KOrder::unfold>::Tm& KOrder::m<KOrder::unfold>()
{
return _um;
}
template<>
const ctraits<KOrder::unfold>::Tm& KOrder::m<KOrder::unfold>() const
{ return _um;}
template<>
ctraits<KOrder::fold>::Tm& KOrder::m<KOrder::fold>()
{
return _fm;
}
template<>
const ctraits<KOrder::fold>::Tm& KOrder::m<KOrder::fold>() const
{ return _fm;}
/* Here is the constructor of the |KOrder| class. We pass what we have
to. The partitioning of the $y$ vector, a sparse container with model
derivatives, then the first order approximation, these are $g_y$ and
$g_u$ matrices, and covariance matrix of exogenous shocks |v|.
We build the members, it is nothing difficult. Note that we do not make
a physical copy of sparse tensors, so during running the class, the
outer world must not change them.
In the body, we have to set |nvs| array, and initialize $g$ and $G$
containers to comply to preconditions of |performStep|. */
KOrder::KOrder(int num_stat, int num_pred, int num_both, int num_forw,
const TensorContainer<FSSparseTensor> &fcont,
const TwoDMatrix &gy, const TwoDMatrix &gu, const TwoDMatrix &v,
Journal &jr)
: ypart(num_stat, num_pred, num_both, num_forw),
ny(ypart.ny()), nu(gu.ncols()), maxk(fcont.getMaxDim()),
nvs(4),
_ug(4), _fg(4), _ugs(4), _fgs(4), _ugss(4), _fgss(4),
_uG(4), _fG(4),
_uZstack(&_uG, ypart.nyss(), &_ug, ny, ypart.nys(), nu),
_fZstack(&_fG, ypart.nyss(), &_fg, ny, ypart.nys(), nu),
_uGstack(&_ugs, ypart.nys(), nu),
_fGstack(&_fgs, ypart.nys(), nu),
_um(maxk, v), _fm(_um), f(fcont),
matA(*(f.get(Symmetry(1))), _uZstack.getStackSizes(), gy, ypart),
matS(*(f.get(Symmetry(1))), _uZstack.getStackSizes(), gy, ypart),
matB(*(f.get(Symmetry(1))), _uZstack.getStackSizes()),
journal(jr)
{
KORD_RAISE_IF(gy.ncols() != ypart.nys(),
"Wrong number of columns in gy in KOrder constructor");
KORD_RAISE_IF(v.ncols() != nu,
"Wrong number of columns of Vcov in KOrder constructor");
KORD_RAISE_IF(nu != v.nrows(),
"Wrong number of rows of Vcov in KOrder constructor");
KORD_RAISE_IF(maxk < 2,
"Order of approximation must be at least 2 in KOrder constructor");
KORD_RAISE_IF(gy.nrows() != ypart.ny(),
"Wrong number of rows in gy in KOrder constructor");
KORD_RAISE_IF(gu.nrows() != ypart.ny(),
"Wrong number of rows in gu in KOrder constructor");
KORD_RAISE_IF(gu.ncols() != nu,
"Wrong number of columns in gu in KOrder constructor");
// set nvs:
nvs[0] = ypart.nys(); nvs[1] = nu; nvs[2] = nu; nvs[3] = 1;
// put $g_y$ and $g_u$ to the container
/* Note that $g_\sigma$ is zero by the nature and we do not insert it to
the container. We insert a new physical copies. */
UGSTensor *tgy = new UGSTensor(ny, TensorDimens(Symmetry(1, 0, 0, 0), nvs));
tgy->getData() = gy.getData();
insertDerivative<unfold>(tgy);
UGSTensor *tgu = new UGSTensor(ny, TensorDimens(Symmetry(0, 1, 0, 0), nvs));
tgu->getData() = gu.getData();
insertDerivative<unfold>(tgu);
// put $G_y$, $G_u$ and $G_{u'}$ to the container
/* Also note that since $g_\sigma$ is zero, so $G_\sigma$. */
UGSTensor *tGy = faaDiBrunoG<unfold>(Symmetry(1, 0, 0, 0));
G<unfold>().insert(tGy);
UGSTensor *tGu = faaDiBrunoG<unfold>(Symmetry(0, 1, 0, 0));
G<unfold>().insert(tGu);
UGSTensor *tGup = faaDiBrunoG<unfold>(Symmetry(0, 0, 1, 0));
G<unfold>().insert(tGup);
}
// |KOrder::sylvesterSolve| unfolded specialization
/* Here we have an unfolded specialization of |sylvesterSolve|. We
simply create the sylvester object and solve it. Note that the $g^*_y$
is not continuous in memory as assumed by the sylvester code, so we
make a temporary copy and pass it as matrix $C$.
If the $B$ matrix is empty, in other words there are now forward
looking variables, then the system becomes $AX=D$ which is solved by
simple |matA.multInv()|.
If one wants to display the diagnostic messages from the Sylvester
module, then after the |sylv.solve()| one needs to call
|sylv.getParams().print("")|. */
template<>
void
KOrder::sylvesterSolve<KOrder::unfold>(ctraits<unfold>::Ttensor &der) const
{
JournalRecordPair pa(journal);
pa << "Sylvester equation for dimension = " << der.getSym()[0] << endrec;
if (ypart.nys() > 0 && ypart.nyss() > 0)
{
KORD_RAISE_IF(!der.isFinite(),
"RHS of Sylverster is not finite");
TwoDMatrix gs_y(*(gs<unfold>().get(Symmetry(1, 0, 0, 0))));
GeneralSylvester sylv(der.getSym()[0], ny, ypart.nys(),
ypart.nstat+ypart.npred,
matA.getData(), matB.getData(),
gs_y.getData(), der.getData());
sylv.solve();
}
else if (ypart.nys() > 0 && ypart.nyss() == 0)
{
matA.multInv(der);
}
}
// |KOrder::sylvesterSolve| folded specialization
/* Here is the folded specialization of sylvester. We unfold the right
hand side. Then we solve it by the unfolded version of
|sylvesterSolve|, and fold it back and copy to output vector. */
template<>
void
KOrder::sylvesterSolve<KOrder::fold>(ctraits<fold>::Ttensor &der) const
{
ctraits<unfold>::Ttensor tmp(der);
sylvesterSolve<unfold>(tmp);
ctraits<fold>::Ttensor ftmp(tmp);
der.getData() = (const Vector &) (ftmp.getData());
}
void
KOrder::switchToFolded()
{
JournalRecordPair pa(journal);
pa << "Switching from unfolded to folded" << endrec;
int maxdim = g<unfold>().getMaxDim();
for (int dim = 1; dim <= maxdim; dim++)
{
SymmetrySet ss(dim, 4);
for (symiterator si(ss); !si.isEnd(); ++si)
{
if ((*si)[2] == 0 && g<unfold>().check(*si))
{
auto *ft = new FGSTensor(*(g<unfold>().get(*si)));
insertDerivative<fold>(ft);
if (dim > 1)
{
gss<unfold>().remove(*si);
gs<unfold>().remove(*si);
g<unfold>().remove(*si);
}
}
if (G<unfold>().check(*si))
{
auto *ft = new FGSTensor(*(G<unfold>().get(*si)));
G<fold>().insert(ft);
if (dim > 1)
{
G<fold>().remove(*si);
}
}
}
}
}
dynare++/kord/korder.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2004, Ondra Kamenik
// Higher order at deterministic steady
/*
The main purpose of this file is to implement a perturbation method
algorithm for an SDGE model for higher order approximations. The input
of the algorithm are sparse tensors as derivatives of the dynamic
system, then dimensions of vector variables, then the first order
approximation to the decision rule and finally a covariance matrix of
exogenous shocks. The output are higher order derivatives of decision
rule $y_t=g(y^*_{t-1},u_t,\sigma)$. The class provides also a method
for checking a size of residuals of the solved equations.
The algorithm is implemented in |KOrder| class. The class contains
both unfolded and folded containers to allow for switching (usually
from unfold to fold) during the calculations. The algorithm is
implemented in a few templated methods. To do this, we need some
container type traits, which are in |ctraits| struct. Also, the
|KOrder| class contains some information encapsulated in other
classes, which are defined here. These include: |PartitionY|,
|MatrixA|, |MatrixS| and |MatrixB|.
*/
#ifndef KORDER_H
#define KORDER_H
#include "int_sequence.hh"
#include "fs_tensor.hh"
#include "gs_tensor.hh"
#include "t_container.hh"
#include "stack_container.hh"
#include "normal_moments.hh"
#include "t_polynomial.hh"
#include "faa_di_bruno.hh"
#include "journal.hh"
#include "kord_exception.hh"
#include "GeneralSylvester.hh"
#include <dynlapack.h>
#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|. */
class FoldedZXContainer;
class UnfoldedZXContainer;
class FoldedGXContainer;
class UnfoldedGXContainer;
template <int 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>;
};
/* The |PartitionY| class defines the partitioning of state variables
$y$. The vector $y$, and subvector $y^*$, and $y^{**}$ are defined.
$$y=\left[\matrix{\hbox{static}\cr\hbox{predeter}\cr\hbox{both}\cr
\hbox{forward}}\right],\quad
y^*=\left[\matrix{\hbox{predeter}\cr\hbox{both}}\right],\quad
y^{**}=\left[\matrix{\hbox{both}\cr\hbox{forward}}\right],$$
where ``static'' means variables appearing only at time $t$,
``predeter'' means variables appearing at time $t-1$, but not at
$t+1$, ``both'' means variables appearing both at $t-1$ and $t+1$
(regardless appearance at $t$), and ``forward'' means variables
appearing at $t+1$, but not at $t-1$.
The class maintains the four lengths, and returns the whole length,
length of $y^s$, and length of $y^{**}$.
*/
struct PartitionY
{
const int nstat;
const int npred;
const int nboth;
const int nforw;
PartitionY(int num_stat, int num_pred,
int num_both, int num_forw)
: nstat(num_stat), npred(num_pred),
nboth(num_both), nforw(num_forw)
{
}
int
ny() const
{
return nstat+npred+nboth+nforw;
}
int
nys() const
{
return npred+nboth;
}
int
nyss() const
{
return nboth+nforw;
}
};
/* This is an abstraction for a square matrix with attached PLU
factorization. It can calculate the PLU factorization and apply the
inverse with some given matrix.
We use LAPACK $PLU$ decomposition for the inverse. We store the $L$
and $U$ in the |inv| array and |ipiv| is the permutation $P$. */
class PLUMatrix : public TwoDMatrix
{
public:
PLUMatrix(int n)
: TwoDMatrix(n, n),
inv(nrows()*ncols()),
ipiv(new lapack_int[nrows()])
{
}
PLUMatrix(const PLUMatrix &plu);
~PLUMatrix() override
{
delete [] ipiv;
}
void multInv(TwoDMatrix &m) const;
private:
Vector inv;
lapack_int *ipiv;
protected:
void calcPLU();
};
/* The class |MatrixA| is used for matrix $\left[f_{y}\right]+ \left[0
\left[f_{y^{**}_+}\right]\cdot\left[g^{**}_{y^*}\right] 0\right]$,
which is central for the perturbation method step. */
class MatrixA : public PLUMatrix
{
public:
MatrixA(const FSSparseTensor &f, const IntSequence &ss,
const TwoDMatrix &gy, const PartitionY &ypart);
};
/* The class |MatrixS| slightly differs from |MatrixA|. It is used for
matrix $$\left[f_{y}\right]+ \left[0
\quad\left[f_{y^{**}_+}\right]\cdot\left[g^{**}_{y^*}\right]\quad
0\right]+\left[0\quad 0\quad\left[f_{y^{**}_+}\right]\right]$$, which is
needed when recovering $g_{\sigma^k}$. */
class MatrixS : public PLUMatrix
{
public:
MatrixS(const FSSparseTensor &f, const IntSequence &ss,
const TwoDMatrix &gy, const PartitionY &ypart);
};
/* The $B$ matrix is equal to $\left[f_{y^{**}_+}\right]$. We have just
a constructor. */
class MatrixB : public TwoDMatrix
{
public:
MatrixB(const FSSparseTensor &f, const IntSequence &ss)
: TwoDMatrix(FGSTensor(f, ss, IntSequence(1, 0),
TensorDimens(ss, IntSequence(1, 0))))
{
}
};
/* Here we have the class for the higher order approximations. It
contains the following data:
\halign{\kern\parindent\vrule height12pt width0pt
\vtop{\hsize=4cm\noindent\raggedright #}&\kern0.5cm\vtop{\hsize=10cm\noindent #}\cr
variable sizes ypart& |PartitionY| struct maintaining partitions of
$y$, see |@<|PartitionY| struct declaration@>|\cr
tensor variable dimension |nvs|& variable sizes of all tensors in
containers, sizes of $y^*$, $u$, $u'$@q'@> and $\sigma$\cr
tensor containers & folded and unfolded containers for $g$, $g_{y^*}$,
$g_{y^**}$ (the latter two collect appropriate subtensors of $g$, they
do not allocate any new space), $G$, $G$ stack, $Z$ stack\cr
dynamic model derivatives & just a reference to the container of
sparse tensors of the system derivatives, lives outside the class\cr
moments & both folded and unfolded normal moment containers, both are
calculated at initialization\cr
matrices & matrix $A$, matrix $S$, and matrix $B$, see |@<|MatrixA| class
declaration@>| and |@<|MatrixB| class declaration@>|\cr
}
\kern 0.4cm
The methods are the following:
\halign{\kern\parindent\vrule height12pt width0pt
\vtop{\hsize=4cm\noindent\raggedright #}&\kern0.5cm\vtop{\hsize=10cm\noindent #}\cr
member access & we declare template methods for accessing containers
depending on |fold| and |unfold| flag, we implement their
specializations\cr
|performStep| & this performs $k$-order step provided that $k=2$ or
the $k-1$-th step has been run, this is the core method\cr
|check| & this calculates residuals of all solved equations for
$k$-order and reports their sizes, it is runnable after $k$-order
|performStep| has been run\cr
|insertDerivative| & inserts a $g$ derivative to the $g$ container and
also creates subtensors and insert them to $g_{y^*}$ and $g_{y^{**}}$
containers\cr
|sylvesterSolve| & solve the sylvester equation (templated fold, and
unfold)\cr
|faaDiBrunoZ| & calculates derivatives of $F$ by Faa Di Bruno for the
sparse container of system derivatives and $Z$ stack container\cr
|faaDiBrunoG| & calculates derivatives of $G$ by Faa Di Bruno for the
dense container $g^{**}$ and $G$ stack\cr
|recover_y| & recovers $g_{y^{*i}}$\cr
|recover_yu| & recovers $g_{y^{*i}u^j}$\cr
|recover_ys| & recovers $g_{y^{*i}\sigma^j}$\cr
|recover_yus| & recovers $g_{y^{*i}u^j\sigma^k}$\cr
|recover_s| & recovers $g_{\sigma^i}$\cr
|fillG| & calculates specified derivatives of $G$ and inserts them to
the container\cr
|calcE_ijk|& calculates $E_{ijk}$\cr
|calcD_ijk|& calculates $D_{ijk}$\cr
}
\kern 0.3cm
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.
*/
class KOrder
{
protected:
const PartitionY ypart;
const int ny;
const int nu;
const int maxk;
IntSequence nvs;
/* These are containers. The names are not important because they do
not appear anywhere else since we access them by template functions. */
UGSContainer _ug;
FGSContainer _fg;
UGSContainer _ugs;
FGSContainer _fgs;
UGSContainer _ugss;
FGSContainer _fgss;
UGSContainer _uG;
FGSContainer _fG;
UnfoldedZContainer _uZstack;
FoldedZContainer _fZstack;
UnfoldedGContainer _uGstack;
FoldedGContainer _fGstack;
UNormalMoments _um;
FNormalMoments _fm;
const TensorContainer<FSSparseTensor> &f;
const MatrixA matA;
const MatrixS matS;
const MatrixB matB;
/* 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;
Journal &journal;
public:
KOrder(int num_stat, int num_pred, int num_both, int num_forw,
const TensorContainer<FSSparseTensor> &fcont,
const TwoDMatrix &gy, const TwoDMatrix &gu, const TwoDMatrix &v,
Journal &jr);
enum { fold, unfold };
template <int t>
void performStep(int order);
template <int t>
double check(int dim) const;
template <int t>
Vector *calcStochShift(int order, double sigma) const;
void switchToFolded();
const PartitionY &
getPartY() const
{
return ypart;
}
const FGSContainer &
getFoldDers() const
{
return _fg;
}
const UGSContainer &
getUnfoldDers() const
{
return _ug;
}
static bool
is_even(int i)
{
return (i/2)*2 == i;
}
protected:
template <int t>
void insertDerivative(_Ttensor *der);
template<int t>
void sylvesterSolve(_Ttensor &der) const;
template <int t>
_Ttensor *faaDiBrunoZ(const Symmetry &sym) const;
template <int t>
_Ttensor *faaDiBrunoG(const Symmetry &sym) const;
template<int t>
void recover_y(int i);
template <int t>
void recover_yu(int i, int j);
template <int t>
void recover_ys(int i, int j);
template <int t>
void recover_yus(int i, int j, int k);
template <int t>
void recover_s(int i);
template<int 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;
};
/* Here we insert the result to the container. Along the insertion, we
also create subtensors and insert as well. */
template <int t>
void
KOrder::insertDerivative(_Ttensor *der)
{
g<t>().insert(der);
gs<t>().insert(new _Ttensor(ypart.nstat, ypart.nys(), *der));
gss<t>().insert(new _Ttensor(ypart.nstat+ypart.npred,
ypart.nyss(), *der));
}
/* Here we implement Faa Di Bruno formula
$$\sum_{l=1}^k\left[f_{z^l}\right]_{\gamma_1\ldots\gamma_l}
\sum_{c\in M_{l,k}}\prod_{m=1}^l\left[z_{s(c_m)}\right]^{\gamma_m},
$$
where $s$ is a given outer symmetry and $k$ is the dimension of the
symmetry. */
template <int t>
_Ttensor *
KOrder::faaDiBrunoZ(const Symmetry &sym) const
{
JournalRecordPair pa(journal);
pa << "Faa Di Bruno Z container for " << sym << endrec;
_Ttensor *res = new _Ttensor(ny, TensorDimens(sym, nvs));
FaaDiBruno bruno(journal);
bruno.calculate(Zstack<t>(), f, *res);
return res;
}
/* The same as |@<|KOrder::faaDiBrunoZ| templated code@>|, but for
$g^{**}$ and $G$ stack. */
template <int 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 = new _Ttensor(ypart.nyss(), tdims);
FaaDiBruno bruno(journal);
bruno.calculate(Gstack<t>(), gss<t>(), *res);
return res;
}
/* Here we solve $\left[F_{y^i}\right]=0$. First we calculate
conditional $G_{y^i}$ (it misses $l=1$ and $l=i$ since $g_{y^i}$ does
not exist yet). Then calculate conditional $F_{y^i}$ and we have the
right hand side of equation. Since we miss two orders, we solve by
Sylvester, and insert the solution as the derivative $g_{y^i}$. Then
we need to update $G_{y^i}$ running |multAndAdd| for both dimensions
$1$ and $i$.
{\bf Requires:} everything at order $\leq i-1$
{\bf Provides:} $g_{y^i}$, and $G_{y^i}$ */
template<int t>
void
KOrder::recover_y(int i)
{
Symmetry sym(i, 0, 0, 0);
JournalRecordPair pa(journal);
pa << "Recovering symmetry " << sym << endrec;
_Ttensor *G_yi = faaDiBrunoG<t>(sym);
G<t>().insert(G_yi);
_Ttensor *g_yi = faaDiBrunoZ<t>(sym);
g_yi->mult(-1.0);
sylvesterSolve<t>(*g_yi);
insertDerivative<t>(g_yi);
_Ttensor *gss_y = gss<t>().get(Symmetry(1, 0, 0, 0));
gs<t>().multAndAdd(*gss_y, *G_yi);
_Ttensor *gss_yi = gss<t>().get(sym);
gs<t>().multAndAdd(*gss_yi, *G_yi);
}
/* Here we solve $\left[F_{y^iu^j}\right]=0$ to obtain $g_{y^iu^j}$ for
$j>0$. We calculate conditional $G_{y^iu^j}$ (this misses only $l=1$)
and calculate conditional $F_{y^iu^j}$ and we have the right hand
side. It is solved by multiplication of inversion of $A$. Then we insert
the result, and update $G_{y^iu^j}$ by |multAndAdd| for $l=1$.
{\bf Requires:} everything at order $\leq i+j-1$, $G_{y^{i+j}}$, and
$g_{y^{i+j}}$.
{\bf Provides:} $g_{y^iu^j}$, and $G_{y^iu^j}$ */
template <int t>
void
KOrder::recover_yu(int i, int j)
{
Symmetry sym(i, j, 0, 0);
JournalRecordPair pa(journal);
pa << "Recovering symmetry " << sym << endrec;
_Ttensor *G_yiuj = faaDiBrunoG<t>(sym);
G<t>().insert(G_yiuj);
_Ttensor *g_yiuj = faaDiBrunoZ<t>(sym);
g_yiuj->mult(-1.0);
matA.multInv(*g_yiuj);
insertDerivative<t>(g_yiuj);
gs<t>().multAndAdd(*(gss<t>().get(Symmetry(1, 0, 0, 0))), *G_yiuj);
}
/* Here we solve
$\left[F_{y^i\sigma^j}\right]+\left[D_{ij}\right]+\left[E_{ij}\right]=0$
to obtain $g_{y^i\sigma^j}$. We calculate conditional
$G_{y^i\sigma^j}$ (missing dimensions $1$ and $i+j$), calculate
conditional $F_{y^i\sigma^j}$. Before we can calculate $D_{ij}$ and
$E_{ij}$, we have to calculate $G_{y^iu'^m\sigma^{j-m}}$ for
$m=1,\ldots,j$. Then we add the $D_{ij}$ and $E_{ij}$ to obtain the
right hand side. Then we solve the sylvester to obtain
$g_{y^i\sigma^j}$. Then we update $G_{y^i\sigma^j}$ for $l=1$ and
$l=i+j$.
{\bf Requires:} everything at order $\leq i+j-1$, $g_{y^{i+j}}$,
$G_{y^iu'^j}$ and $g_{y^iu^j}$ through $D_{ij}$,
$g_{y^iu^m\sigma^{j-m}}$ for
$m=1,\ldots,j-1$ through $E_{ij}$.
{\bf Provides:} $g_{y^i\sigma^j}$ and $G_{y^i\sigma^j}$, and finally
$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>
void
KOrder::recover_ys(int i, int j)
{
Symmetry sym(i, 0, 0, j);
JournalRecordPair pa(journal);
pa << "Recovering symmetry " << sym << endrec;
fillG<t>(i, 0, j);
if (is_even(j))
{
_Ttensor *G_yisj = faaDiBrunoG<t>(sym);
G<t>().insert(G_yisj);
_Ttensor *g_yisj = faaDiBrunoZ<t>(sym);
{
_Ttensor *D_ij = calcD_ik<t>(i, j);
g_yisj->add(1.0, *D_ij);
delete D_ij;
}
if (j >= 3)
{
_Ttensor *E_ij = calcE_ik<t>(i, j);
g_yisj->add(1.0, *E_ij);
delete E_ij;
}
g_yisj->mult(-1.0);
sylvesterSolve<t>(*g_yisj);
insertDerivative<t>(g_yisj);
Gstack<t>().multAndAdd(1, gss<t>(), *G_yisj);
Gstack<t>().multAndAdd(i+j, gss<t>(), *G_yisj);
}
}
/* Here we solve
$\left[F_{y^iu^j\sigma^k}\right]+\left[D_{ijk}\right]+\left[E_{ijk}\right]=0$
to obtain $g_{y^iu^j\sigma^k}$. First we calculate conditional
$G_{y^iu^j\sigma^k}$ (missing only for dimension $l=1$), then we
evaluate conditional $F_{y^iu^j\sigma^k}$. Before we can calculate
$D_{ijk}$, and $E_{ijk}$, we need to insert
$G_{y^iu^ju'^m\sigma^{k-m}}$ for $m=1,\ldots, k$. This is done by
|fillG|. Then we have right hand side and we multiply by $A^{-1}$ to
obtain $g_{y^iu^j\sigma^k}$. Finally we have to update
$G_{y^iu^j\sigma^k}$ by |multAndAdd| for dimension $l=1$.
{\bf Requires:} everything at order $\leq i+j+k$, $g_{y^{i+j}\sigma^k}$
through $G_{y^iu^j\sigma^k}$ involved in right hand side, then
$g_{y^iu^{j+k}}$ through $D_{ijk}$, and $g_{y^iu^{j+m}\sigma^{k-m}}$
for $m=1,\ldots,k-1$ through $E_{ijk}$.
{\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>
void
KOrder::recover_yus(int i, int j, int k)
{
Symmetry sym(i, j, 0, k);
JournalRecordPair pa(journal);
pa << "Recovering symmetry " << sym << endrec;
fillG<t>(i, j, k);
if (is_even(k))
{
_Ttensor *G_yiujsk = faaDiBrunoG<t>(sym);
G<t>().insert(G_yiujsk);
_Ttensor *g_yiujsk = faaDiBrunoZ<t>(sym);
{
_Ttensor *D_ijk = calcD_ijk<t>(i, j, k);
g_yiujsk->add(1.0, *D_ijk);
delete D_ijk;
}
if (k >= 3)
{
_Ttensor *E_ijk = calcE_ijk<t>(i, j, k);
g_yiujsk->add(1.0, *E_ijk);
delete E_ijk;
}
g_yiujsk->mult(-1.0);
matA.multInv(*g_yiujsk);
insertDerivative<t>(g_yiujsk);
Gstack<t>().multAndAdd(1, gss<t>(), *G_yiujsk);
}
}
/* Here we solve
$\left[F_{\sigma^i}\right]+\left[D_i\right]+\left[E_i\right]=0$ to
recover $g_{\sigma^i}$. First we calculate conditional $G_{\sigma^i}$
(missing dimension $l=1$ and $l=i$), then we calculate conditional
$F_{\sigma^i}$. Before we can calculate $D_i$ and $E_i$, we have to
obtain $G_{u'm\sigma^{i-m}}$ for $m=1,\ldots,i$. Than
adding $D_i$ and $E_i$ we have the right hand side. We solve by
$S^{-1}$ multiplication and update $G_{\sigma^i}$ by calling
|multAndAdd| for dimension $l=1$.
Recall that the solved equation here is:
$$
\left[f_y\right]\left[g_{\sigma^k}\right]+
\left[f_{y^{**}_+}\right]\left[g^{**}_{y^*}\right]\left[g^*_{\sigma^k}\right]+
\left[f_{y^{**}_+}\right]\left[g^{**}_{\sigma^k}\right]=\hbox{RHS}
$$
This is a sort of deficient sylvester equation (sylvester equation for
dimension=0), we solve it by $S^{-1}$. See |@<|MatrixS| constructor
code@>| to see how $S$ looks like.
{\bf Requires:} everything at order $\leq i-1$, $g_{y^i}$ and
$g_{y^{i-j}\sigma^j}$, then $g_{u^k}$ through $F_{u'^k}$, and
$g_{y^mu^j\sigma^k}$ for $j=1,\ldots,i-1$ and $m+j+k=i$ through
$F_{u'j\sigma^{i-j}}$.
{\bf Provides:} $g_{\sigma^i}$, $G_{\sigma^i}$, and
$G_{u'^m\sigma^{i-m}}$ for $m=1,\ldots,i$ */
template <int t>
void
KOrder::recover_s(int i)
{
Symmetry sym(0, 0, 0, i);
JournalRecordPair pa(journal);
pa << "Recovering symmetry " << sym << endrec;
fillG<t>(0, 0, i);
if (is_even(i))
{
_Ttensor *G_si = faaDiBrunoG<t>(sym);
G<t>().insert(G_si);
_Ttensor *g_si = faaDiBrunoZ<t>(sym);
{
_Ttensor *D_i = calcD_k<t>(i);
g_si->add(1.0, *D_i);
delete D_i;
}
if (i >= 3)
{
_Ttensor *E_i = calcE_k<t>(i);
g_si->add(1.0, *E_i);
delete E_i;
}
g_si->mult(-1.0);
matS.multInv(*g_si);
insertDerivative<t>(g_si);
Gstack<t>().multAndAdd(1, gss<t>(), *G_si);
Gstack<t>().multAndAdd(i, gss<t>(), *G_si);
}
}
/* Here we calculate and insert $G_{y^iu^ju'^m\sigma^{k-m}}$ for
$m=1,\ldots, k$. The derivatives are inserted only for $k-m$ being
even. */
template<int t>
void
KOrder::fillG(int i, int j, int k)
{
for (int m = 1; m <= k; m++)
{
if (is_even(k-m))
{
_Ttensor *G_yiujupms = faaDiBrunoG<t>(Symmetry(i, j, m, k-m));
G<t>().insert(G_yiujupms);
}
}
}
/* Here we calculate
$$\left[D_{ijk}\right]_{\alpha_1\ldots\alpha_i\beta_1\ldots\beta_j}=
\left[F_{y^iu^ju'^k}\right]
_{\alpha_1\ldots\alpha_i\beta_1\ldots\beta_j\gamma_1\ldots\gamma_k}
\left[\Sigma\right]^{\gamma_1\ldots\gamma_k}$$
So it is non zero only for even $k$. */
template <int 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();
if (is_even(k))
{
_Ttensor *tmp = faaDiBrunoZ<t>(Symmetry(i, j, k, 0));
tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry(k))));
delete tmp;
}
return res;
}
/* Here we calculate
$$\left[E_{ijk}\right]_{\alpha_1\ldots\alpha_i\beta_1\ldots\beta_j}=
\sum_{m=1}^{k-1}\left(\matrix{k\cr m}\right)\left[F_{y^iu^ju'^m\sigma^{k-m}}\right]
_{\alpha_1\ldots\alpha_i\beta_1\ldots\beta_j\gamma_1\ldots\gamma_m}
\left[\Sigma\right]^{\gamma_1\ldots\gamma_m}$$
The sum can sum only for even $m$. */
template <int 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();
for (int n = 2; n <= k-1; n += 2)
{
_Ttensor *tmp = faaDiBrunoZ<t>(Symmetry(i, j, n, k-n));
tmp->mult((double) (Tensor::noverk(k, n)));
tmp->contractAndAdd(2, *res, *(m<t>().get(Symmetry(n))));
delete tmp;
}
return res;
}
template <int t>
_Ttensor *
KOrder::calcD_ik(int i, int k) const
{
return calcD_ijk<t>(i, 0, k);
}
template <int t>
_Ttensor *
KOrder::calcD_k(int k) const
{
return calcD_ijk<t>(0, 0, k);
}
template <int t>
_Ttensor *
KOrder::calcE_ik(int i, int k) const
{
return calcE_ijk<t>(i, 0, k);
}
template <int t>
_Ttensor *
KOrder::calcE_k(int k) const
{
return calcE_ijk<t>(0, 0, k);
}
/* Here is the core routine. It calls methods recovering derivatives in
the right order. Recall, that the code, namely Faa Di Bruno's formula,
is implemented as to be run conditionally on the current contents of
containers. So, if some call of Faa Di Bruno evaluates derivatives,
and some derivatives are not present in the container, then it is
considered to be zero. So, we have to be very careful to put
everything in the right order. The order here can be derived from
dependencies, or it is in the paper.
The method recovers all the derivatives of the given |order|.
The precondition of the method is that all tensors of order |order-1|,
which are not zero, exist (including $G$). The postcondition of of the
method is derivatives of $g$ and $G$ of order |order| are calculated
and stored in the containers. Responsibility of precondition lays upon
the constructor (for |order==2|), or upon the previous call of
|performStep|.
From the code, it is clear, that all $g$ are calculated. If one goes
through all the recovering methods, he should find out that also all
$G$ are provided. */
template <int t>
void
KOrder::performStep(int order)
{
KORD_RAISE_IF(order-1 != g<t>().getMaxDim(),
"Wrong order for KOrder::performStep");
JournalRecordPair pa(journal);
pa << "Performing step for order = " << order << endrec;
recover_y<t>(order);
for (int i = 0; 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_ys<t>(order-j, j);
}
for (int i = order-1; i >= 1; i--)
{
recover_yus<t>(0, i, order-i);
}
recover_s<t>(order);
}
/* Here we check for residuals of all the solved equations at the given
order. The method returns the largest residual size. Each check simply
evaluates the equation. */
template <int t>
double
KOrder::check(int dim) const
{
KORD_RAISE_IF(dim > g<t>().getMaxDim(),
"Wrong dimension for KOrder::check");
JournalRecordPair pa(journal);
pa << "Checking residuals for order = " << dim << endrec;
double maxerror = 0.0;
// check for $F_{y^iu^j}=0
for (int i = 0; i <= dim; i++)
{
Symmetry sym(dim-i, i, 0, 0);
_Ttensor *r = faaDiBrunoZ<t>(sym);
double err = r->getData().getMax();
JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
if (err > maxerror)
maxerror = err;
delete r;
}
// check for $F_{y^iu^ju'^k}+D_{ijk}+E_{ijk}=0$
SymmetrySet ss(dim, 3);
for (symiterator si(ss); !si.isEnd(); ++si)
{
int i = (*si)[0];
int j = (*si)[1];
int k = (*si)[2];
if (i+j > 0 && k > 0)
{
Symmetry sym(i, j, 0, k);
_Ttensor *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;
double err = r->getData().getMax();
JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
delete r;
}
}
// check for $F_{\sigma^i}+D_i+E_i=0
_Ttensor *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;
double err = r->getData().getMax();
Symmetry sym(0, 0, 0, dim);
JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
if (err > maxerror)
maxerror = err;
delete r;
return maxerror;
}
template <int t>
Vector *
KOrder::calcStochShift(int order, double sigma) const
{
auto *res = new Vector(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;
}
return res;
}
#endif
dynare++/kord/korder_stoch.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2005, Ondra Kamenik
#include "korder_stoch.hh"
/* Same as |@<|MatrixA| constructor code@>|, but the submatrix |gss_ys| is
passed directly. */
MatrixAA::MatrixAA(const FSSparseTensor &f, const IntSequence &ss,
const TwoDMatrix &gss_ys, const PartitionY &ypart)
: PLUMatrix(ypart.ny())
{
zeros();
IntSequence c(1); c[0] = 1;
FGSTensor f_y(f, ss, c, TensorDimens(ss, c));
add(1.0, f_y);
c[0] = 0;
FGSTensor f_yss(f, ss, c, TensorDimens(ss, c));
TwoDMatrix sub(*this, ypart.nstat, ypart.nys());
sub.multAndAdd(f_yss, gss_ys);
calcPLU();
}
// |KOrderStoch| folded constructor code
KOrderStoch::KOrderStoch(const PartitionY &yp, int nu,
const TensorContainer<FSSparseTensor> &fcont,
const FGSContainer &hh, Journal &jr)
: nvs(4), ypart(yp), journal(jr),
_ug(4), _fg(4), _ugs(4), _fgs(4), _uG(4), _fG(4),
_uh(nullptr), _fh(&hh),
_uZstack(&_uG, ypart.nyss(), &_ug, ypart.ny(), ypart.nys(), nu),
_fZstack(&_fG, ypart.nyss(), &_fg, ypart.ny(), ypart.nys(), nu),
_uGstack(&_ugs, ypart.nys(), nu),
_fGstack(&_fgs, ypart.nys(), nu),
f(fcont),
matA(*(fcont.get(Symmetry(1))), _uZstack.getStackSizes(), *(hh.get(Symmetry(1, 0, 0, 0))),
ypart)
{
nvs[0] = ypart.nys();
nvs[1] = nu;
nvs[2] = nu;
nvs[3] = 1;
}
// |KOrderStoch| unfolded constructor code
KOrderStoch::KOrderStoch(const PartitionY &yp, int nu,
const TensorContainer<FSSparseTensor> &fcont,
const UGSContainer &hh, Journal &jr)
: nvs(4), ypart(yp), journal(jr),
_ug(4), _fg(4), _ugs(4), _fgs(4), _uG(4), _fG(4),
_uh(&hh), _fh(nullptr),
_uZstack(&_uG, ypart.nyss(), &_ug, ypart.ny(), ypart.nys(), nu),
_fZstack(&_fG, ypart.nyss(), &_fg, ypart.ny(), ypart.nys(), nu),
_uGstack(&_ugs, ypart.nys(), nu),
_fGstack(&_fgs, ypart.nys(), nu),
f(fcont),
matA(*(fcont.get(Symmetry(1))), _uZstack.getStackSizes(), *(hh.get(Symmetry(1, 0, 0, 0))),
ypart)
{
nvs[0] = ypart.nys();
nvs[1] = nu;
nvs[2] = nu;
nvs[3] = 1;
}
// |KOrderStoch| convenience method specializations
template<>
ctraits<KOrder::unfold>::Tg& KOrderStoch::g<KOrder::unfold>()
{
return _ug;
}
template<>
const ctraits<KOrder::unfold>::Tg& KOrderStoch::g<KOrder::unfold>() const
{ return _ug;}
template<>
ctraits<KOrder::fold>::Tg& KOrderStoch::g<KOrder::fold>()
{
return _fg;
}
template<>
const ctraits<KOrder::fold>::Tg& KOrderStoch::g<KOrder::fold>() const
{ return _fg;}
template<>
ctraits<KOrder::unfold>::Tgs& KOrderStoch::gs<KOrder::unfold>()
{
return _ugs;
}
template<>
const ctraits<KOrder::unfold>::Tgs& KOrderStoch::gs<KOrder::unfold>() const
{ return _ugs;}
template<>
ctraits<KOrder::fold>::Tgs& KOrderStoch::gs<KOrder::fold>()
{
return _fgs;
}
template<>
const ctraits<KOrder::fold>::Tgs& KOrderStoch::gs<KOrder::fold>() const
{ return _fgs;}
template<>
const ctraits<KOrder::unfold>::Tgss& KOrderStoch::h<KOrder::unfold>() const
{ return *_uh;}
template<>
const ctraits<KOrder::fold>::Tgss& KOrderStoch::h<KOrder::fold>() const
{ return *_fh;}
template<>
ctraits<KOrder::unfold>::TG& KOrderStoch::G<KOrder::unfold>()
{
return _uG;
}
template<>
const ctraits<KOrder::unfold>::TG& KOrderStoch::G<KOrder::unfold>() const
{ return _uG;}
template<>
ctraits<KOrder::fold>::TG& KOrderStoch::G<KOrder::fold>()
{
return _fG;
}
template<>
const ctraits<KOrder::fold>::TG& KOrderStoch::G<KOrder::fold>() const
{ return _fG;}
template<>
ctraits<KOrder::unfold>::TZXstack& KOrderStoch::Zstack<KOrder::unfold>()
{
return _uZstack;
}
template<>
const ctraits<KOrder::unfold>::TZXstack& KOrderStoch::Zstack<KOrder::unfold>() const
{ return _uZstack;}
template<>
ctraits<KOrder::fold>::TZXstack& KOrderStoch::Zstack<KOrder::fold>()
{
return _fZstack;
}
template<>
const ctraits<KOrder::fold>::TZXstack& KOrderStoch::Zstack<KOrder::fold>() const
{ return _fZstack;}
template<>
ctraits<KOrder::unfold>::TGXstack& KOrderStoch::Gstack<KOrder::unfold>()
{
return _uGstack;
}
template<>
const ctraits<KOrder::unfold>::TGXstack& KOrderStoch::Gstack<KOrder::unfold>() const
{ return _uGstack;}
template<>
ctraits<KOrder::fold>::TGXstack& KOrderStoch::Gstack<KOrder::fold>()
{
return _fGstack;
}
template<>
const ctraits<KOrder::fold>::TGXstack& KOrderStoch::Gstack<KOrder::fold>() const
{ return _fGstack;}
dynare++/kord/korder_stoch.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2005, Ondra Kamenik
// Higher order at stochastic steady
/* This file defines a number of classes of which |KOrderStoch| is the
main purpose. Basically, |KOrderStoch| calculates first and higher
order Taylor expansion of a policy rule at $\sigma>0$ with explicit
forward $g^{**}$. More formally, we have to solve a policy rule $g$
from
$$E_t\left[f(g^{**}(g^*(y^*_t,u_t,\sigma),u_{t+1},\sigma),g(y^*,u_t,\sigma),y^*,u_t)\right]$$
As an introduction in {\tt approximation.hweb} argues, $g^{**}$ at
tine $t+1$ must be given from outside. Let the explicit
$E_t(g^{**}(y^*,u_{t+1},\sigma)$ be equal to $h(y^*,\sigma)$. Then we
have to solve
$$f(h(g^*(y^*,u,\sigma),\sigma),g(y,u,\sigma),y,u),$$
which is much easier than fully implicit system for $\sigma=0$.
Besides the class |KOrderStoch|, we declare here also classes for the
new containers corresponding to
$f(h(g^*(y^*,u,\sigma),\sigma),g(y,u,\sigma),y,u)$. Further, we
declare |IntegDerivs| and |StochForwardDerivs| classes which basically
calculate $h$ as an extrapolation based on an approximation to $g$ at
lower $\sigma$. */
#include "korder.hh"
#include "faa_di_bruno.hh"
#include "journal.hh"
/* This class is a container, which has a specialized constructor
integrating the policy rule at given $\sigma$. */
template <int t>
class IntegDerivs : public ctraits<t>::Tgss
{
public:
IntegDerivs(int r, const IntSequence &nvs, const _Tgss &g, const __Tm &mom,
double at_sigma);
};
/* This constructor integrates a rule (namely its $g^{**}$ part) with
respect to $u=\tilde\sigma\eta$, and stores to the object the
derivatives of this integral $h$ at $(y^*,u,\sigma)=(\tilde
y^*,0,\tilde\sigma)$. The original container of $g^{**}$, the moments of
the stochastic shocks |mom| and the $\tilde\sigma$ are input.
The code follows the following derivation
\def\lims{\vbox{\baselineskip=0pt\lineskip=1pt
\setbox0=\hbox{$\scriptstyle n+k=p$}\hbox to\wd0{\hss$\scriptstyle m=0$\hss}\box0}}
$$
\eqalign{h(y,\sigma)&=E_t\left[g(y,u',\sigma)\right]=\cr
&=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+j+k=d}\pmatrix{d\cr i,j,k}\left[g_{y^iu^j\sigma^k}\right]
(y^*-\tilde y^*)^i\sigma^j\Sigma^j(\sigma-\tilde\sigma)^k\cr
&=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+m+n+k=d}\pmatrix{d\cr i,m+n,k}
\left[g_{y^iu^{m+n}\sigma^k}\right]
\hat y^{*i}\Sigma^{m+n}\pmatrix{m+n\cr m,n}{\tilde\sigma}^m\hat\sigma^{k+n}\cr
&=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+m+n+k=d}\pmatrix{d\cr i,m,n,k}
\left[g_{y^iu^{m+n}\sigma^k}\right]
\Sigma^{m+n}{\tilde\sigma}^m\hat y^{*i}\hat\sigma^{k+n}\cr
&=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+p=d}\sum_{\lims}\pmatrix{d\cr i,m,n,k}
\left[g_{y^iu^{m+n}\sigma^k}\right]
\Sigma^{m+n}{\tilde\sigma}^m\hat y^{*i}\hat\sigma^{k+n}\cr
&=\tilde y+\sum_{d=1}{1\over d!}\sum_{i+p=d}\pmatrix{d\cr i,p}
\left[\sum_{\lims}\pmatrix{p\cr n,k}{1\over m!}
\left[g_{y^iu^{m+n}\sigma^k}\right]
\Sigma^{m+n}{\tilde\sigma}^m\right]\hat y^{*i}\hat\sigma^{k+n},
}
$$
where $\pmatrix{a\cr b_1,\ldots, b_n}$ is a generalized combination
number, $p=k+n$, $\hat\sigma=\sigma-\tilde\sigma$, $\hat
y^*=y^*-\tilde y$, and we dropped writing the multidimensional indexes
in Einstein summation.
This implies that
$$h_{y^i\sigma^p}=\sum_{\lims}\pmatrix{p\cr n,k}{1\over m!}
\left[g_{y^iu^{m+n}\sigma^k}\right]
\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)
: ctraits<t>::Tgss(4)
{
int maxd = g.getMaxDim();
for (int d = 1; d <= maxd; d++)
{
for (int i = 0; i <= d; i++)
{
int p = d-i;
Symmetry sym(i, 0, 0, p);
_Ttensor *ten = new _Ttensor(r, TensorDimens(sym, nvs));
// calculate derivative $h_{y^i\sigma^p}$
/* This code calculates
$$h_{y^i\sigma^p}=\sum_{\lims}\pmatrix{p\cr n,k}{1\over m!}
\left[g_{y^iu^{m+n}\sigma^k}\right]
\Sigma^{m+n}{\tilde\sigma}^m$$
and stores it in |ten|. */
ten->zeros();
for (int n = 0; n <= p; n++)
{
int k = p-n;
int povern = Tensor::noverk(p, n);
int mfac = 1;
for (int m = 0; i+m+n+k <= maxd; m++, mfac *= m)
{
double mult = (pow(at_sigma, m)*povern)/mfac;
Symmetry sym_mn(i, m+n, 0, k);
if (m+n == 0 && g.check(sym_mn))
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)));
gtmp.mult(mult);
gtmp.contractAndAdd(1, *ten, *(mom.get(Symmetry(m+n))));
}
}
}
this->insert(ten);
}
}
}
/* This class calculates an extrapolation of expectation of forward
derivatives. It is a container, all calculations are done in a
constructor.
The class calculates derivatives of $E[g(y*,u,\sigma)]$ at $(\bar
y^*,\bar\sigma)$. The derivatives are extrapolated based on
derivatives at $(\tilde y^*,\tilde\sigma)$. */
template <int t>
class StochForwardDerivs : public ctraits<t>::Tgss
{
public:
StochForwardDerivs(const PartitionY &ypart, int nu,
const _Tgss &g, const __Tm &m,
const Vector &ydelta, double sdelta,
double at_sigma);
};
/* This is the constructor which performs the integration and the
extrapolation. Its parameters are: |g| is the container of derivatives
at $(\tilde y,\tilde\sigma)$; |m| are the moments of stochastic
shocks; |ydelta| is a difference of the steady states $\bar y-\tilde
y$; |sdelta| is a difference between new sigma and old sigma
$\bar\sigma-\tilde\sigma$, and |at_sigma| is $\tilde\sigma$. There is
no need of inputing the $\tilde y$.
We do the operation in four steps:
\orderedlist
\li Integrate $g^{**}$, the derivatives are at $(\tilde y,\tilde\sigma)$
\li Form the (full symmetric) polynomial from the derivatives stacking
$\left[\matrix{y^*\cr\sigma}\right]$
\li Centralize this polynomial about $(\bar y,\bar\sigma)$
\li Recover general symmetry tensors from the (full symmetric) polynomial
\endorderedlist */
template <int t>
StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
const _Tgss &g, const __Tm &m,
const Vector &ydelta, double sdelta,
double at_sigma)
: ctraits<t>::Tgss(4)
{
int maxd = g.getMaxDim();
int r = ypart.nyss();
// make |g_int| be integral of $g^{**}$ at $(\tilde y,\tilde\sigma)$
/* This simply constructs |IntegDerivs| class. Note that the |nvs| of
the tensors has zero dimensions for shocks, this is because we need to
make easily stacks of the form $\left[\matrix{y^*\cr\sigma}\right]$ in
the next step. */
IntSequence nvs(4);
nvs[0] = ypart.nys(); nvs[1] = 0; nvs[2] = 0; nvs[3] = 1;
IntegDerivs<t> g_int(r, nvs, g, m, at_sigma);
// 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);
for (int d = 1; d <= maxd; d++)
{
auto *ten = new _Ttensym(r, ypart.nys()+1, d);
ten->zeros();
for (int i = 0; i <= d; i++)
{
int k = d-i;
if (g_int.check(Symmetry(i, 0, 0, k)))
ten->addSubTensor(*(g_int.get(Symmetry(i, 0, 0, k))));
}
g_int_sym.insert(ten);
}
// make |g_int_cent| the centralized polynomial about $(\bar y,\bar\sigma)$
/* Here we centralize the polynomial to $(\bar y,\bar\sigma)$ knowing
that the polynomial was centralized about $(\tilde
y,\tilde\sigma)$. This is done by derivating and evaluating the
derivated polynomial at $(\bar y-\tilde
y,\bar\sigma-\tilde\sigma)$. The stack of this vector is |delta| in
the code. */
Vector delta(ypart.nys()+1);
Vector dy(delta, 0, ypart.nys());
ConstVector dy_in(ydelta, ypart.nstat, ypart.nys());
dy = dy_in;
delta[ypart.nys()] = sdelta;
_Tpol g_int_cent(r, ypart.nys()+1);
for (int d = 1; d <= maxd; d++)
{
g_int_sym.derivative(d-1);
_Ttensym *der = g_int_sym.evalPartially(d, delta);
g_int_cent.insert(der);
}
// pull out general symmetry tensors from |g_int_cent|
/* Here we only recover the general symmetry derivatives from the full
symmetric polynomial. Note that the derivative get the true |nvs|. */
IntSequence ss(4);
ss[0] = ypart.nys(); ss[1] = 0; ss[2] = 0; ss[3] = 1;
IntSequence pp(4);
pp[0] = 0; pp[1] = 1; pp[2] = 2; pp[3] = 3;
IntSequence true_nvs(nvs);
true_nvs[1] = nu; true_nvs[2] = nu;
for (int d = 1; d <= maxd; d++)
{
if (g_int_cent.check(Symmetry(d)))
{
for (int i = 0; i <= d; i++)
{
Symmetry sym(i, 0, 0, d-i);
IntSequence coor(sym, pp);
_Ttensor *ten = new _Ttensor(*(g_int_cent.get(Symmetry(d))), ss, coor,
TensorDimens(sym, true_nvs));
this->insert(ten);
}
}
}
}
/* This container corresponds to $h(g^*(y,u,\sigma),\sigma)$. Note that
in our application, the $\sigma$ as a second argument to $h$ will be
its fourth variable in symmetry, so we have to do four member stack
having the second and third stack dummy. */
template <class _Ttype>
class GXContainer : public GContainer<_Ttype>
{
public:
using _Stype = StackContainerInterface<_Ttype>;
using _Ctype = typename StackContainer<_Ttype>::_Ctype;
using itype = typename StackContainer<_Ttype>::itype;
GXContainer(const _Ctype *gs, int ngs, int nu)
: GContainer<_Ttype>(gs, ngs, nu)
{
}
itype getType(int i, const Symmetry &s) const override;
};
/* This routine corresponds to this stack:
$$\left[\matrix{g^*(y,u,\sigma)\cr dummy\cr dummy\cr\sigma}\right]$$ */
template <class _Ttype>
typename GXContainer<_Ttype>::itype
GXContainer<_Ttype>::getType(int i, const Symmetry &s) const
{
if (i == 0)
if (s[2] > 0)
return _Stype::zero;
else
return _Stype::matrix;
if (i == 1)
return _Stype::zero;
if (i == 2)
return _Stype::zero;
if (i == 3)
if (s == Symmetry(0, 0, 0, 1))
return _Stype::unit;
else
return _Stype::zero;
KORD_RAISE("Wrong stack index in GXContainer::getType");
}
/* This container corresponds to $f(H(y,u,\sigma),g(y,u,sigma),y,u)$,
where the $H$ has the size (number of rows) as $g^{**}$. Since it is
very simmilar to |ZContainer|, we inherit form it and override only
|getType| method. */
template <class _Ttype>
class ZXContainer : public ZContainer<_Ttype>
{
public:
using _Stype = StackContainerInterface<_Ttype>;
using _Ctype = typename StackContainer<_Ttype>::_Ctype;
using itype = typename StackContainer<_Ttype>::itype;
ZXContainer(const _Ctype *gss, int ngss, const _Ctype *g, int ng, int ny, int nu)
: ZContainer<_Ttype>(gss, ngss, g, ng, ny, nu)
{
}
itype getType(int i, const Symmetry &s) const override;
};
/* This |getType| method corresponds to this stack:
$$\left[\matrix{H(y,u,\sigma)\cr g(y,u,\sigma)\cr y\cr u}\right]$$ */
template <class _Ttype>
typename ZXContainer<_Ttype>::itype
ZXContainer<_Ttype>::getType(int i, const Symmetry &s) const
{
if (i == 0)
if (s[2] > 0)
return _Stype::zero;
else
return _Stype::matrix;
if (i == 1)
if (s[2] > 0)
return _Stype::zero;
else
return _Stype::matrix;
if (i == 2)
if (s == Symmetry(1, 0, 0, 0))
return _Stype::unit;
else
return _Stype::zero;
if (i == 3)
if (s == Symmetry(0, 1, 0, 0))
return _Stype::unit;
else
return _Stype::zero;
KORD_RAISE("Wrong stack index in ZXContainer::getType");
}
class UnfoldedGXContainer : public GXContainer<UGSTensor>, public UnfoldedStackContainer
{
public:
using _Ctype = TensorContainer<UGSTensor>;
UnfoldedGXContainer(const _Ctype *gs, int ngs, int nu)
: GXContainer<UGSTensor>(gs, ngs, nu)
{
}
};
class FoldedGXContainer : public GXContainer<FGSTensor>, public FoldedStackContainer
{
public:
using _Ctype = TensorContainer<FGSTensor>;
FoldedGXContainer(const _Ctype *gs, int ngs, int nu)
: GXContainer<FGSTensor>(gs, ngs, nu)
{
}
};
class UnfoldedZXContainer : public ZXContainer<UGSTensor>, public UnfoldedStackContainer
{
public:
using _Ctype = TensorContainer<UGSTensor>;
UnfoldedZXContainer(const _Ctype *gss, int ngss, const _Ctype *g, int ng, int ny, int nu)
: ZXContainer<UGSTensor>(gss, ngss, g, ng, ny, nu)
{
}
};
class FoldedZXContainer : public ZXContainer<FGSTensor>, public FoldedStackContainer
{
public:
using _Ctype = TensorContainer<FGSTensor>;
FoldedZXContainer(const _Ctype *gss, int ngss, const _Ctype *g, int ng, int ny, int nu)
: ZXContainer<FGSTensor>(gss, ngss, g, ng, ny, nu)
{
}
};
/* This matrix corresponds to
$$\left[f_{y}\right]+ \left[0
\left[f_{y^{**}_+}\right]\cdot\left[h^{**}_{y^*}\right] 0\right]$$
This is very the same as |MatrixA|, the only difference that the
|MatrixA| is constructed from whole $h_{y^*}$, not only from
$h^{**}_{y^*}$, hence the new abstraction. */
class MatrixAA : public PLUMatrix
{
public:
MatrixAA(const FSSparseTensor &f, const IntSequence &ss,
const TwoDMatrix &gyss, const PartitionY &ypart);
};
/* This class calculates derivatives of $g$ given implicitly by
$f(h(g^*(y,u,\sigma),\sigma),g(y,u,\sigma),y,u)$, where $h(y,\sigma)$
is given from outside.
Structurally, the class is very similar to |KOrder|, but calculations
are much easier. The two constructors construct an object from sparse
derivatives of $f$, and derivatives of $h$. The caller must ensure
that the both derivatives are done at the same point.
The calculation for order $k$ (including $k=1$) is done by a call
|performStep(k)|. The derivatives can be retrived by |getFoldDers()|
or |getUnfoldDers()|. */
class KOrderStoch
{
protected:
IntSequence nvs;
PartitionY ypart;
Journal &journal;
UGSContainer _ug;
FGSContainer _fg;
UGSContainer _ugs;
FGSContainer _fgs;
UGSContainer _uG;
FGSContainer _fG;
const UGSContainer *_uh;
const FGSContainer *_fh;
UnfoldedZXContainer _uZstack;
FoldedZXContainer _fZstack;
UnfoldedGXContainer _uGstack;
FoldedGXContainer _fGstack;
const TensorContainer<FSSparseTensor> &f;
MatrixAA matA;
public:
KOrderStoch(const PartitionY &ypart, int nu, const TensorContainer<FSSparseTensor> &fcont,
const FGSContainer &hh, Journal &jr);
KOrderStoch(const PartitionY &ypart, int nu, const TensorContainer<FSSparseTensor> &fcont,
const UGSContainer &hh, Journal &jr);
template <int t>
void performStep(int order);
const FGSContainer &
getFoldDers() const
{
return _fg;
}
const UGSContainer &
getUnfoldDers() const
{
return _ug;
}
protected:
template <int t>
_Ttensor *faaDiBrunoZ(const Symmetry &sym) const;
template <int 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;
};
/* This calculates a derivative of $f(G(y,u,\sigma),g(y,u,\sigma),y,u)$
of a given symmetry. */
template <int t>
_Ttensor *
KOrderStoch::faaDiBrunoZ(const Symmetry &sym) const
{
JournalRecordPair pa(journal);
pa << "Faa Di Bruno ZX container for " << sym << endrec;
_Ttensor *res = new _Ttensor(ypart.ny(), TensorDimens(sym, nvs));
FaaDiBruno bruno(journal);
bruno.calculate(Zstack<t>(), f, *res);
return res;
}
/* This calculates a derivative of
$G(y,u,\sigma)=h(g^*(y,u,\sigma),\sigma)$ of a given symmetry. */
template <int 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 = new _Ttensor(ypart.nyss(), tdims);
FaaDiBruno bruno(journal);
bruno.calculate(Gstack<t>(), h<t>(), *res);
return res;
}
/* This retrieves all $g$ derivatives of a given dimension from implicit
$f(h(g^*(y,u,\sigma),\sigma),g(y,u,\sigma),y,u)$. It suppose that all
derivatives of smaller dimensions have been retrieved.
So, we go through all symmetries $s$, calculate $G_s$ conditional on
$g_s=0$, insert the derivative to the $G$ container, then calculate
$F_s$ conditional on $g_s=0$. This is a righthand side. The left hand
side is $matA\cdot g_s$. The $g_s$ is retrieved as
$$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>
void
KOrderStoch::performStep(int order)
{
int maxd = g<t>().getMaxDim();
KORD_RAISE_IF(order-1 != maxd && (order != 1 || maxd != -1),
"Wrong order for KOrderStoch::performStep");
SymmetrySet ss(order, 4);
for (symiterator si(ss); !si.isEnd(); ++si)
{
if ((*si)[2] == 0)
{
JournalRecordPair pa(journal);
pa << "Recovering symmetry " << *si << endrec;
_Ttensor *G_sym = faaDiBrunoG<t>(*si);
G<t>().insert(G_sym);
_Ttensor *g_sym = faaDiBrunoZ<t>(*si);
g_sym->mult(-1.0);
matA.multInv(*g_sym);
g<t>().insert(g_sym);
gs<t>().insert(new _Ttensor(ypart.nstat, ypart.nys(), *g_sym));
Gstack<t>().multAndAdd(1, h<t>(), *G_sym);
}
}
}
dynare++/kord/mersenne_twister.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2007, Ondra Kamenik
// Mersenne Twister PRNG
/* This file provides a class for generating random numbers with
encapsulated state. It is based on the work of Makoto Matsumoto and
Takuji Nishimura, implementation inspired by code of Richard Wagner
and Geoff Kuenning. */
#ifndef MERSENNE_TWISTER_H
#define MERSENNE_TWISTER_H
#include "random.hh"
#include <cstring>
class MersenneTwister : public RandomGenerator
{
protected:
using uint32 = unsigned int;
enum {STATE_SIZE = 624};
enum {RECUR_OFFSET = 397};
uint32 statevec[STATE_SIZE];
int stateptr;
public:
MersenneTwister(uint32 iseed);
MersenneTwister(const MersenneTwister &mt);
virtual ~MersenneTwister()
= default;
uint32 lrand();
double drand();
double
uniform() override
{
return drand();
}
protected:
void seed(uint32 iseed);
void refresh();
private:
static uint32
hibit(uint32 u)
{
return u & 0x80000000UL;
}
static uint32
lobit(uint32 u)
{
return u & 0x00000001UL;
}
static uint32
lobits(uint32 u)
{
return u & 0x7fffffffUL;
}
static uint32
mixbits(uint32 u, uint32 v)
{
return hibit(u) | lobits(v);
}
static uint32
twist(uint32 m, uint32 s0, uint32 s1)
{
return m ^ (mixbits(s0, s1)>>1) ^ (-lobit(s1) & 0x9908b0dfUL);
}
};
inline MersenneTwister::MersenneTwister(uint32 iseed)
{
seed(iseed);
}
inline MersenneTwister::MersenneTwister(const MersenneTwister &mt)
: stateptr(mt.stateptr)
{
memcpy(statevec, mt.statevec, sizeof(uint32)*STATE_SIZE);
}
inline MersenneTwister::uint32
MersenneTwister::lrand()
{
if (stateptr >= STATE_SIZE)
refresh();
register uint32 v = statevec[stateptr++];
v ^= v >> 11;
v ^= (v << 7) & 0x9d2c5680;
v ^= (v << 15) & 0xefc60000;
return (v ^ (v >> 18));
}
inline double
MersenneTwister::drand()
{
uint32 a = lrand() >> 5;
uint32 b = lrand() >> 6;
return (a*67108864.0+b) * (1.0/9007199254740992.0);
}
// PRNG of D. Knuth
inline void
MersenneTwister::seed(uint32 iseed)
{
statevec[0] = iseed & 0xffffffffUL;
for (int i = 1; i < STATE_SIZE; i++)
{
register uint32 val = statevec[i-1] >> 30;
val ^= statevec[i-1];
val *= 1812433253ul;
val += i;
statevec[i] = val & 0xffffffffUL;
}
refresh();
}
inline void
MersenneTwister::refresh()
{
register uint32 *p = statevec;
for (int i = STATE_SIZE-RECUR_OFFSET; i--; ++p)
*p = twist(p[RECUR_OFFSET], p[0], p[1]);
for (int i = RECUR_OFFSET; --i; ++p)
*p = twist(p[RECUR_OFFSET-STATE_SIZE], p[0], p[1]);
*p = twist(p[RECUR_OFFSET-STATE_SIZE], p[0], statevec[0]);
stateptr = 0;
}
#endif
dynare++/kord/normal_conjugate.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2007, Ondra Kamenik
#include "normal_conjugate.hh"
#include "kord_exception.hh"
// |NormalConj| diffuse prior constructor
NormalConj::NormalConj(int d)
: mu(d), kappa(0), nu(-1), lambda(d, d)
{
mu.zeros();
lambda.zeros();
}
// |NormalConj| data update constructor
NormalConj::NormalConj(const ConstTwoDMatrix &ydata)
: mu(ydata.numRows()), kappa(ydata.numCols()), nu(ydata.numCols()-1),
lambda(ydata.numRows(), ydata.numRows())
{
mu.zeros();
for (int i = 0; i < ydata.numCols(); i++)
mu.add(1.0/ydata.numCols(), ydata.getCol(i));
lambda.zeros();
for (int i = 0; i < ydata.numCols(); i++)
{
Vector diff{ydata.getCol(i)};
diff.add(-1, mu);
lambda.addOuter(diff);
}
}
NormalConj::NormalConj(const NormalConj &nc)
= default;
// |NormalConj::update| one observation code
/* The method performs the following:
$$\eqalign{
\mu_1 = &\; {\kappa_0\over \kappa_0+1}\mu_0 + {1\over \kappa_0+1}y\cr
\kappa_1 = &\; \kappa_0 + 1\cr
\nu_1 = &\; \nu_0 + 1\cr
\Lambda_1 = &\; \Lambda_0 + {\kappa_0\over\kappa_0+1}(y-\mu_0)(y-\mu_0)^T,
}$$ */
void
NormalConj::update(const ConstVector &y)
{
KORD_RAISE_IF(y.length() != mu.length(),
"Wrong length of a vector in NormalConj::update");
mu.mult(kappa/(1.0+kappa));
mu.add(1.0/(1.0+kappa), y);
Vector diff(y);
diff.add(-1, mu);
lambda.addOuter(diff, kappa/(1.0+kappa));
kappa++;
nu++;
}
// |NormalConj::update| multiple observations code
/* The method evaluates the formula in the header file. */
void
NormalConj::update(const ConstTwoDMatrix &ydata)
{
NormalConj nc(ydata);
update(nc);
}
// |NormalConj::update| with |NormalConj| code
void
NormalConj::update(const NormalConj &nc)
{
double wold = ((double) kappa)/(kappa+nc.kappa);
double wnew = 1-wold;
mu.mult(wold);
mu.add(wnew, nc.mu);
Vector diff(nc.mu);
diff.add(-1, mu);
lambda.add(1.0, nc.lambda);
lambda.addOuter(diff);
kappa = kappa + nc.kappa;
nu = nu + nc.kappa;
}
/* This returns ${1\over \nu-d-1}\Lambda$, which is the mean of the
variance in the posterior distribution. If the number of degrees of
freedom is less than $d$, then NaNs are returned. */
void
NormalConj::getVariance(TwoDMatrix &v) const
{
if (nu > getDim()+1)
{
v = (const TwoDMatrix &) lambda;
v.mult(1.0/(nu-getDim()-1));
}
else
v.nans();
}
dynare++/kord/normal_conjugate.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2007, Ondra Kamenik
// Conjugate family for normal distribution
/* The main purpose here is to implement a class representing conjugate
distributions for mean and variance of the normal distribution. The
class has two main methods: the first one is to update itself with
respect to one observation, the second one is to update itself with
respect to anothe object of the class. In the both methods, the
previous state of the class corresponds to the prior distribution, and
the final state corresponds to the posterior distribution.
The algrebra can be found in Gelman, Carlin, Stern, Rubin (p.87). It
goes as follows: Prior conjugate distribution takes the following form:
$$\eqalign{
\Sigma \sim& {\rm InvWishart}_{\nu_0}(\Lambda_0^{-1}) \cr
\mu\vert\Sigma \sim& N(\mu_0,\Sigma/\kappa_0)
}$$
If the observations are $y_1\ldots y_n$, then the posterior distribution has the
same form with the following parameters:
$$\eqalign{
\mu_n = &\; {\kappa_0\over \kappa_0+n}\mu_0 + {n\over \kappa_0+n}\bar y\cr
\kappa_n = &\; \kappa_0 + n\cr
\nu_n = &\; \nu_0 + n\cr
\Lambda_n = &\; \Lambda_0 + S + {\kappa_0 n\over\kappa_0+n}(\bar y-\mu_0)(\bar y-\mu_0)^T,
}$$
where
$$\eqalign{
\bar y = &\; {1\over n}\sum_{i=1}^ny_i\cr
S = &\; \sum_{i=1}^n(y_i-\bar y)(y_i-\bar y)^T
}$$ */
#ifndef NORMAL_CONJUGATE_H
#define NORMAL_CONJUGATE_H
#include "twod_matrix.hh"
/* The class is described by the four parameters: $\mu$, $\kappa$, $\nu$ and
$\Lambda$. */
class NormalConj
{
protected:
Vector mu;
int kappa;
int nu;
TwoDMatrix lambda;
public:
/* We provide the following constructors: The first constructs diffuse
(Jeffrey's) prior. It sets $\kappa$, and $\Lambda$ to zeros, $nu$ to
$-1$ and also the mean $\mu$ to zero (it should not be
referenced). The second constructs the posterior using the diffuse
prior and the observed data (columnwise). The third is a copy
constructor. */
NormalConj(int d);
NormalConj(const ConstTwoDMatrix &ydata);
NormalConj(const NormalConj &nc);
virtual ~NormalConj()
= default;
void update(const ConstVector &y);
void update(const ConstTwoDMatrix &ydata);
void update(const NormalConj &nc);
int
getDim() const
{
return mu.length();
}
const Vector &
getMean() const
{
return mu;
}
void getVariance(TwoDMatrix &v) const;
};
#endif
dynare++/kord/random.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2007, Ondra Kamenik
#include "random.hh"
#include <cstdlib>
#include <limits>
#include <cmath>
SystemRandomGenerator system_random_generator;
int
RandomGenerator::int_uniform()
{
double s = std::numeric_limits<int>::max()*uniform();
return (int) s;
}
/* This implements Marsaglia Polar Method. */
double
RandomGenerator::normal()
{
double x1, x2;
double w;
do
{
x1 = 2*uniform()-1;
x2 = 2*uniform()-1;
w = x1*x1 + x2*x2;
}
while (w >= 1.0 || w < 1.0e-30);
return x1*std::sqrt((-2.0*std::log(w))/w);
}
double
SystemRandomGenerator::uniform()
{
#if !defined(__MINGW32__)
return drand48();
#else
return ((double) rand())/RAND_MAX;
#endif
}
void
SystemRandomGenerator::initSeed(int seed)
{
#if !defined(__MINGW32__)
srand48(seed);
#else
srand(seed);
#endif
}
dynare++/kord/random.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2007, Ondra Kamenik
// Random number generation
#ifndef RANDOM_H
#define RANDOM_H
/* This is a general interface to an object able to generate random
numbers. Subclass needs to implement |uniform| method, other is, by
default, implemented here. */
class RandomGenerator
{
public:
virtual double uniform() = 0;
int int_uniform();
double normal();
};
/* This implements |RandomGenerator| interface with system |drand| or
|rand|. It is not thread aware. */
class SystemRandomGenerator : public RandomGenerator
{
public:
double uniform() override;
void initSeed(int seed);
};
extern SystemRandomGenerator system_random_generator;
#endif
dynare++/kord/tests.cc deleted 100644 → 0
View file @ 219d2bb3
/* $Id: tests.cpp 148 2005-04-19 15:12:26Z kamenik $ */
/* Copyright 2004, Ondra Kamenik */
#include <cstdlib>
#include "korder.hh"
#include "SylvException.hh"
struct Rand
{
static void init(int n1, int n2, int n3, int n4, int n5);
static double get(double m);
static int get(int m);
static bool discrete(double prob); // answers true with given probability
};
ConstTwoDMatrix
make_matrix(int rows, int cols, const double *p)
{
return ConstTwoDMatrix{rows, cols,
ConstVector{std::shared_ptr<const double>{p, [](const double *arr) {}}, rows*cols}};
}
void
Rand::init(int n1, int n2, int n3, int n4, int n5)
{
long int seed = n1;
seed = 256*seed+n2;
seed = 256*seed+n3;
seed = 256*seed+n4;
seed = 256*seed+n5;
srand48(seed);
}
double
Rand::get(double m)
{
return 2*m*(drand48()-0.5);
}
int
Rand::get(int m)
{
return (int) (Rand::get(0.9999*m));
}
bool
Rand::discrete(double prob)
{
return drand48() < prob;
}
struct SparseGenerator
{
static FSSparseTensor *makeTensor(int dim, int nv, int r,
double fill, double m);
static void fillContainer(TensorContainer<FSSparseTensor> &c,
int maxdim, int nv, int r, double m);
};
FSSparseTensor *
SparseGenerator::makeTensor(int dim, int nv, int r,
double fill, double m)
{
auto *res = new FSSparseTensor(dim, nv, r);
FFSTensor dummy(0, nv, dim);
for (Tensor::index fi = dummy.begin(); fi != dummy.end(); ++fi)
{
for (int i = 0; i < r; i++)
{
if (Rand::discrete(fill))
{
double x = Rand::get(m);
res->insert(fi.getCoor(), i, x);
}
}
}
return res;
}
void
SparseGenerator::fillContainer(TensorContainer<FSSparseTensor> &c,
int maxdim, int nv, int r,
double m)
{
Rand::init(maxdim, nv, r, (int) (5*m), 0);
double fill = 0.5;
for (int d = 1; d <= maxdim; d++)
{
c.insert(makeTensor(d, nv, r, fill, m));
fill *= 0.3;
}
}
const double vdata [] = { // 3x3
0.1307870268, 0.1241940078, 0.1356703123,
0.1241940078, 0.1986920419, 0.2010160581,
0.1356703123, 0.2010160581, 0.2160336975
};
const double gy_data [] = { // 8x4
0.3985178619, -0.5688233582, 0.9572900437, -0.6606847776, 0.1453004017,
0.3025310675, -0.8627437750, -0.6903410191, 0.4751910580, -0.7270018589,
-0.0939612498, -0.1463831989, 0.6742110220, 0.6046671043, 0.5215893126,
-1.0412969986, -0.3524898417, -1.0986703430, 0.8006531522, 0.8879776376,
-0.1037608317, -0.5587378073, -0.1010366945, 0.9462411248, -0.2439199881,
1.3420621236, -0.7820285935, 0.3205293447, 0.3606124791, 0.2975422208,
-0.5452861965, 1.6320340279
};
const double gu_data [] = { // just some numbers, no structure
1.8415286914, -0.2638743845, 1.7690713274, 0.9668585956, 0.2303143646,
-0.2229624279, -0.4381991822, 1.0082401405, -0.3186555860, -0.0624691529,
-0.5189085756, 1.4269672156, 0.1163282969, 1.4020183445, -0.0952660426,
0.2099097124, 0.6912400502, -0.5180935114, 0.5288316624, 0.2188053448,
0.5715516767, 0.7813893410, -0.6385073106, 0.8335131513, 0.3605202168,
-1.1167944865, -1.2263750934, 0.6113636081, 0.6964915482, -0.6451217688,
0.4062810500, -2.0552251116, -1.6383406284, 0.0198915095, 0.0111014458,
-1.2421792262, -1.0724161722, -0.4276904972, 0.1801494950, -2.0716473264
};
const double vdata2 [] = { // 10x10 positive definite
0.79666, -0.15536, 0.05667, -0.21026, 0.20262, 0.28505, 0.60341, -0.09703, 0.32363, 0.13299,
-0.15536, 0.64380, -0.01131, 0.00980, 0.03755, 0.43791, 0.21784, -0.31755, -0.55911, -0.29655,
0.05667, -0.01131, 0.56165, -0.34357, -0.40584, 0.20990, 0.28348, 0.20398, -0.19856, 0.35820,
-0.21026, 0.00980, -0.34357, 0.56147, 0.10972, -0.34146, -0.49906, -0.19685, 0.21088, -0.31560,
0.20262, 0.03755, -0.40584, 0.10972, 0.72278, 0.02155, 0.04089, -0.19696, 0.03446, -0.12919,
0.28505, 0.43791, 0.20990, -0.34146, 0.02155, 0.75867, 0.77699, -0.31125, -0.55141, -0.02155,
0.60341, 0.21784, 0.28348, -0.49906, 0.04089, 0.77699, 1.34553, -0.18613, -0.25811, -0.19016,
-0.09703, -0.31755, 0.20398, -0.19685, -0.19696, -0.31125, -0.18613, 0.59470, 0.08386, 0.41750,
0.32363, -0.55911, -0.19856, 0.21088, 0.03446, -0.55141, -0.25811, 0.08386, 0.98917, -0.12992,
0.13299, -0.29655, 0.35820, -0.31560, -0.12919, -0.02155, -0.19016, 0.41750, -0.12992, 0.89608
};
const double gy_data2 [] = { // 600 items make gy 30x20, whose gy(6:25,:) has spectrum within unit
0.39414, -0.29766, 0.08948, -0.19204, -0.00750, 0.21159, 0.05494, 0.06225, 0.01771, 0.21913,
-0.01373, 0.20086, -0.06086, -0.10955, 0.14424, -0.08390, 0.03948, -0.14713, 0.11674, 0.05091,
0.24039, 0.28307, -0.11835, 0.13030, 0.11682, -0.27444, -0.19311, -0.16654, 0.12867, 0.25116,
-0.19781, 0.45242, -0.15862, 0.24428, -0.11966, 0.11483, -0.32279, 0.29727, 0.20934, -0.18190,
-0.15080, -0.09477, -0.30551, -0.02672, -0.26919, 0.11165, -0.06390, 0.03449, -0.26622, 0.22197,
0.45141, -0.41683, 0.09760, 0.31094, -0.01652, 0.05809, -0.04514, -0.05645, 0.00554, 0.47980,
0.11726, 0.42459, -0.13136, -0.30902, -0.14648, 0.11455, 0.02947, -0.03835, -0.04044, 0.03559,
-0.26575, -0.01783, 0.31243, -0.14412, -0.13218, -0.05080, 0.18576, 0.13840, -0.05560, 0.35530,
-0.25573, -0.11560, 0.15187, -0.18431, 0.08193, -0.32278, 0.17560, -0.05529, -0.10020, -0.23088,
-0.20979, -0.49245, 0.09915, -0.16909, -0.03443, 0.19497, 0.18473, 0.25662, 0.29605, -0.20531,
-0.39244, -0.43369, 0.05588, 0.24823, -0.14236, -0.08311, 0.16371, -0.19975, 0.30605, -0.17087,
-0.01270, 0.00123, -0.22426, -0.13810, 0.05079, 0.06971, 0.01922, -0.09952, -0.23177, -0.41962,
-0.41991, 0.41430, -0.04247, -0.13706, -0.12048, -0.28906, -0.22813, -0.25057, -0.18579, -0.20642,
-0.47976, 0.25490, -0.05138, -0.30794, 0.31651, 0.02034, 0.12954, -0.20110, 0.13336, -0.40775,
-0.30195, -0.13704, 0.12396, 0.28152, 0.02986, 0.27669, 0.24623, 0.08635, -0.11956, -0.02949,
0.37401, 0.20838, 0.24801, -0.26872, 0.11195, 0.00315, -0.19069, 0.12839, -0.23036, -0.48228,
0.08434, -0.39872, -0.28896, -0.28754, 0.24668, 0.23285, 0.25437, 0.10456, -0.14124, 0.20483,
-0.19117, -0.33836, -0.24875, 0.08207, -0.03930, 0.20364, 0.15384, -0.15270, 0.24372, -0.11199,
-0.46591, 0.30319, 0.05745, 0.09084, 0.06058, 0.31884, 0.05071, -0.28899, -0.30793, -0.03566,
0.02286, 0.28178, 0.00736, -0.31378, -0.18144, -0.22346, -0.27239, 0.31043, -0.26228, 0.22181,
-0.15096, -0.36953, -0.06032, 0.21496, 0.29545, -0.13112, 0.16420, -0.07573, -0.43111, -0.43057,
0.26716, -0.31209, -0.05866, -0.29101, -0.27437, -0.18727, 0.28732, -0.19014, 0.08837, 0.30405,
0.06103, -0.35612, 0.00173, 0.25134, -0.08987, -0.22766, -0.03254, -0.18662, -0.08491, 0.49401,
-0.12145, -0.02961, -0.03668, -0.30043, -0.08555, 0.01701, -0.12544, 0.10969, -0.48202, 0.07245,
0.20673, 0.11408, 0.04343, -0.01815, -0.31594, -0.23632, -0.06258, -0.27474, 0.12180, 0.16613,
-0.37931, 0.30219, 0.15765, 0.25489, 0.17529, -0.17020, -0.30060, 0.22058, -0.02450, -0.42143,
0.49642, 0.46899, -0.28552, -0.22549, -0.01333, 0.21567, 0.22251, 0.21639, -0.19194, -0.19140,
-0.24106, 0.10952, -0.11019, 0.29763, -0.02039, -0.25748, 0.23169, 0.01357, 0.09802, -0.19022,
0.37604, -0.40777, 0.18131, -0.10258, 0.29573, -0.31773, 0.09069, -0.02198, -0.26594, 0.48302,
-0.10041, 0.20210, -0.05609, -0.01169, -0.17339, 0.17862, -0.22502, 0.29009, -0.45160, 0.19771,
0.27634, 0.31695, -0.09993, 0.17167, 0.12394, 0.28088, -0.12502, -0.16967, -0.06296, -0.17036,
0.27320, 0.01595, 0.16955, 0.30146, -0.15173, -0.29807, 0.08178, -0.06811, 0.21655, 0.26348,
0.06316, 0.45661, -0.29756, -0.05742, -0.14715, -0.03037, -0.16656, -0.08768, 0.38078, 0.40679,
-0.32779, -0.09106, 0.16107, -0.07301, 0.07700, -0.22694, -0.15692, -0.02548, 0.38749, -0.12203,
-0.02980, -0.22067, 0.00680, -0.23058, -0.29112, 0.23032, -0.16026, 0.23392, -0.09990, 0.03628,
-0.42592, -0.33474, -0.09499, -0.17442, -0.20110, 0.24618, -0.06418, -0.06715, 0.40754, 0.29377,
0.29543, -0.16832, -0.08468, 0.06491, -0.01410, 0.19988, 0.24950, 0.14626, -0.27851, 0.06079,
0.48134, -0.13475, 0.25398, 0.11738, 0.23369, -0.00661, -0.16811, -0.04557, -0.12030, -0.39527,
-0.35760, 0.01840, -0.15941, 0.03290, 0.09988, -0.08307, 0.06644, -0.24637, 0.34112, -0.08026,
0.00951, 0.27656, 0.16247, 0.28217, 0.17198, -0.16389, -0.03835, -0.02675, -0.08032, -0.21045,
-0.38946, 0.23207, 0.10987, -0.31674, -0.28653, -0.27430, -0.29109, -0.00648, 0.38431, -0.38478,
-0.41195, -0.19364, -0.20977, -0.05524, 0.05558, -0.20109, 0.11803, -0.19884, 0.43318, -0.39255,
0.26612, -0.21771, 0.12471, 0.12856, -0.15104, -0.11676, 0.17582, -0.25330, 0.00298, -0.31712,
0.21532, -0.20319, 0.14507, -0.04588, -0.22995, -0.06470, 0.18849, -0.13444, 0.37107, 0.07387,
-0.14008, 0.09896, 0.13727, -0.28417, -0.09461, -0.18703, 0.04080, 0.02343, -0.49988, 0.17993,
0.23189, -0.30581, -0.18334, -0.09667, -0.27699, -0.05998, 0.09118, -0.32453, 0.46251, 0.41500,
-0.45314, -0.00544, 0.08529, 0.29099, -0.00937, -0.31650, 0.26163, 0.14506, 0.37498, -0.16454,
0.35215, 0.31642, -0.09161, -0.31452, -0.04792, -0.04677, -0.19523, 0.27998, 0.05491, 0.44461,
-0.01258, -0.27887, 0.18361, -0.04539, -0.02977, 0.30821, 0.29454, -0.17932, 0.16193, 0.23934,
0.47923, 0.25373, 0.23258, 0.31484, -0.17958, -0.01136, 0.17681, 0.12869, 0.03235, 0.43762,
0.13734, -0.09433, -0.03735, 0.17949, 0.14122, -0.17814, 0.06359, 0.16044, 0.12249, -0.22314,
0.40775, 0.05147, 0.12389, 0.04290, -0.01642, 0.00082, -0.18056, 0.02875, 0.32690, 0.17712,
0.34001, -0.21581, -0.01086, -0.18180, 0.17480, -0.17774, -0.07503, 0.28438, -0.19747, 0.29595,
-0.28002, -0.02073, -0.16522, -0.18234, -0.20565, 0.29620, 0.07502, 0.01429, -0.31418, 0.43693,
-0.12212, 0.11178, -0.28503, 0.04683, 0.00072, 0.05566, 0.18857, 0.26101, -0.38891, -0.21216,
-0.21850, -0.15147, -0.30749, -0.23762, 0.14984, 0.03535, -0.02862, -0.00105, -0.39907, -0.06909,
-0.36094, 0.21717, 0.15930, -0.18924, 0.13741, 0.01039, 0.13613, 0.00659, 0.07676, -0.13711,
0.24285, -0.07564, -0.28349, -0.15658, 0.03135, -0.30909, -0.22534, 0.17363, -0.19376, 0.26038,
0.05546, -0.22607, 0.32420, -0.02552, -0.05400, 0.13388, 0.04643, -0.31535, -0.06181, 0.30237,
-0.04680, -0.29441, 0.12231, 0.03960, -0.01188, 0.01406, 0.25402, 0.03315, 0.25026, -0.10922
};
const double gu_data2 [] = { // raw data 300 items
0.26599, 0.41329, 0.31846, 0.92590, 0.43050, 0.17466, 0.02322, 0.72621, 0.37921, 0.70597,
0.97098, 0.14023, 0.57619, 0.09938, 0.02281, 0.92341, 0.72654, 0.71000, 0.76687, 0.70182,
0.88752, 0.49524, 0.42549, 0.42806, 0.57615, 0.76051, 0.15341, 0.47457, 0.60066, 0.40880,
0.20668, 0.41949, 0.97620, 0.94318, 0.71491, 0.56402, 0.23553, 0.94387, 0.78567, 0.06362,
0.85252, 0.86262, 0.25190, 0.03274, 0.93216, 0.37971, 0.08797, 0.14596, 0.73871, 0.06574,
0.67447, 0.28575, 0.43911, 0.92133, 0.12327, 0.87762, 0.71060, 0.07141, 0.55443, 0.53310,
0.91529, 0.25121, 0.07593, 0.94490, 0.28656, 0.82174, 0.68887, 0.67337, 0.99291, 0.03316,
0.02849, 0.33891, 0.25594, 0.90071, 0.01248, 0.67871, 0.65953, 0.65369, 0.97574, 0.31578,
0.23678, 0.39220, 0.06706, 0.80943, 0.57694, 0.08220, 0.18151, 0.19969, 0.37096, 0.37858,
0.70153, 0.46816, 0.76511, 0.02520, 0.39387, 0.25527, 0.39050, 0.60141, 0.30322, 0.46195,
0.12025, 0.33616, 0.04174, 0.00196, 0.68886, 0.74445, 0.15869, 0.18994, 0.95195, 0.62874,
0.82874, 0.53369, 0.34383, 0.50752, 0.97023, 0.22695, 0.62407, 0.25840, 0.71279, 0.28785,
0.31611, 0.20391, 0.19702, 0.40760, 0.85158, 0.68369, 0.63760, 0.09879, 0.11924, 0.32920,
0.53052, 0.15900, 0.21229, 0.84080, 0.33933, 0.93651, 0.42705, 0.06199, 0.50092, 0.47192,
0.57152, 0.01818, 0.31404, 0.50173, 0.87725, 0.50530, 0.10717, 0.04035, 0.32901, 0.33538,
0.04780, 0.40984, 0.78216, 0.91288, 0.11314, 0.25248, 0.23823, 0.74001, 0.48089, 0.55531,
0.82486, 0.01058, 0.05409, 0.44357, 0.52641, 0.68188, 0.94629, 0.61627, 0.33037, 0.11961,
0.57988, 0.19653, 0.91902, 0.59838, 0.52974, 0.28364, 0.45767, 0.65836, 0.63045, 0.76140,
0.27918, 0.27256, 0.46035, 0.77418, 0.92918, 0.14095, 0.89645, 0.25146, 0.21172, 0.47910,
0.95451, 0.34377, 0.29927, 0.79220, 0.97654, 0.67591, 0.44385, 0.38434, 0.44860, 0.28170,
0.90712, 0.20337, 0.00292, 0.55046, 0.62255, 0.45127, 0.80896, 0.43965, 0.59145, 0.23801,
0.33601, 0.30119, 0.89935, 0.40850, 0.98226, 0.75430, 0.68318, 0.65407, 0.68067, 0.32942,
0.11756, 0.27626, 0.83879, 0.72174, 0.75430, 0.13702, 0.03402, 0.58781, 0.07393, 0.23067,
0.92537, 0.29445, 0.43437, 0.47685, 0.54548, 0.66082, 0.23805, 0.60208, 0.94337, 0.21363,
0.72637, 0.57181, 0.77679, 0.63931, 0.72860, 0.38901, 0.94920, 0.04535, 0.12863, 0.40550,
0.90095, 0.21418, 0.13953, 0.99639, 0.02526, 0.70018, 0.21828, 0.20294, 0.20191, 0.30954,
0.39490, 0.68955, 0.11506, 0.15748, 0.40252, 0.91680, 0.61547, 0.78443, 0.19693, 0.67630,
0.56552, 0.58556, 0.53554, 0.53507, 0.09831, 0.21229, 0.83135, 0.26375, 0.89287, 0.97069,
0.70615, 0.42041, 0.43117, 0.21291, 0.26086, 0.26978, 0.77340, 0.43833, 0.46179, 0.54418,
0.67878, 0.42776, 0.61454, 0.55915, 0.36363, 0.31999, 0.42442, 0.86649, 0.62513, 0.02047
};
class TestRunnable
{
char name[100];
public:
int dim; // dimension of the solved problem
int nvar; // number of variable of the solved problem
TestRunnable(const char *n, int d, int nv)
: dim(d), nvar(nv)
{
strncpy(name, n, 100);
}
bool test() const;
virtual bool run() const = 0;
const char *
getName() const
{
return name;
}
protected:
static double korder_unfold_fold(int maxdim, int unfold_dim,
int nstat, int npred, int nboth, int forw,
const TwoDMatrix &gy, const TwoDMatrix &gu,
const TwoDMatrix &v);
};
bool
TestRunnable::test() const
{
printf("Running test <%s>\n", name);
clock_t start = clock();
bool passed = run();
clock_t end = clock();
printf("CPU time %8.4g (CPU seconds)..................",
((double) (end-start))/CLOCKS_PER_SEC);
if (passed)
{
printf("passed\n\n");
return passed;
}
else
{
printf("FAILED\n\n");
return passed;
}
}
double
TestRunnable::korder_unfold_fold(int maxdim, int unfold_dim,
int nstat, int npred, int nboth, int nforw,
const TwoDMatrix &gy, const TwoDMatrix &gu,
const TwoDMatrix &v)
{
TensorContainer<FSSparseTensor> c(1);
int ny = nstat+npred+nboth+nforw;
int nu = v.nrows();
int nz = nboth+nforw+ny+nboth+npred+nu;
SparseGenerator::fillContainer(c, maxdim, nz, ny, 5.0);
for (int d = 1; d <= maxdim; d++)
{
printf("\ttensor fill for dim=%d is: %3.2f %%\n",
d, c.get(Symmetry(d))->getFillFactor()*100.0);
}
Journal jr("out.txt");
KOrder kord(nstat, npred, nboth, nforw, c, gy, gu, v, jr);
// perform unfolded steps until unfold_dim
double maxerror = 0.0;
for (int d = 2; d <= unfold_dim; d++)
{
clock_t pertime = clock();
kord.performStep<KOrder::unfold>(d);
pertime = clock()-pertime;
printf("\ttime for unfolded step dim=%d: %8.4g\n",
d, ((double) (pertime))/CLOCKS_PER_SEC);
clock_t checktime = clock();
double err = kord.check<KOrder::unfold>(d);
checktime = clock()-checktime;
printf("\ttime for step check dim=%d: %8.4g\n",
d, ((double) (checktime))/CLOCKS_PER_SEC);
printf("\tmax error in step dim=%d: %10.6g\n",
d, err);
if (maxerror < err)
maxerror = err;
}
// perform folded steps until maxdim
if (unfold_dim < maxdim)
{
clock_t swtime = clock();
kord.switchToFolded();
swtime = clock()-swtime;
printf("\ttime for switching dim=%d: %8.4g\n",
unfold_dim, ((double) (swtime))/CLOCKS_PER_SEC);
for (int d = unfold_dim+1; d <= maxdim; d++)
{
clock_t pertime = clock();
kord.performStep<KOrder::fold>(d);
pertime = clock()-pertime;
printf("\ttime for folded step dim=%d: %8.4g\n",
d, ((double) (pertime))/CLOCKS_PER_SEC);
clock_t checktime = clock();
double err = kord.check<KOrder::fold>(d);
checktime = clock()-checktime;
printf("\ttime for step check dim=%d: %8.4g\n",
d, ((double) (checktime))/CLOCKS_PER_SEC);
printf("\tmax error in step dim=%d: %10.6g\n",
d, err);
if (maxerror < err)
maxerror = err;
}
}
return maxerror;
}
class UnfoldKOrderSmall : public TestRunnable
{
public:
UnfoldKOrderSmall()
: TestRunnable("unfold-3 fold-4 korder (stat=2,pred=3,both=1,forw=2,u=3,dim=4)",
4, 18)
{
}
bool
run() const override
{
TwoDMatrix gy{make_matrix(8, 4, gy_data)};
TwoDMatrix gu{make_matrix(8, 3, gu_data)};
TwoDMatrix v{make_matrix(3, 3, vdata)};
double err = korder_unfold_fold(4, 3, 2, 3, 1, 2,
gy, gu, v);
return err < 0.08;
}
};
// same dimension as Smets & Wouters
class UnfoldKOrderSW : public TestRunnable
{
public:
UnfoldKOrderSW()
: TestRunnable("unfold S&W korder (stat=5,pred=12,both=8,forw=5,u=10,dim=4)",
4, 73)
{
}
bool
run() const override
{
TwoDMatrix gy{make_matrix(30, 20, gy_data2)};
TwoDMatrix gu{make_matrix(30, 10, gu_data2)};
TwoDMatrix v{make_matrix(10, 10, vdata2)};
v.mult(0.001);
gu.mult(.01);
double err = korder_unfold_fold(4, 4, 5, 12, 8, 5,
gy, gu, v);
return err < 0.08;
}
};
class UnfoldFoldKOrderSW : public TestRunnable
{
public:
UnfoldFoldKOrderSW()
: TestRunnable("unfold-2 fold-3 S&W korder (stat=5,pred=12,both=8,forw=5,u=10,dim=3)",
4, 73)
{
}
bool
run() const override
{
TwoDMatrix gy{make_matrix(30, 20, gy_data2)};
TwoDMatrix gu{make_matrix(30, 10, gu_data2)};
TwoDMatrix v{make_matrix(10, 10, vdata2)};
v.mult(0.001);
gu.mult(.01);
double err = korder_unfold_fold(4, 3, 5, 12, 8, 5,
gy, gu, v);
return err < 0.08;
}
};
int
main()
{
TestRunnable *all_tests[50];
// fill in vector of all tests
int num_tests = 0;
all_tests[num_tests++] = new UnfoldKOrderSmall();
all_tests[num_tests++] = new UnfoldKOrderSW();
all_tests[num_tests++] = new UnfoldFoldKOrderSW();
// find maximum dimension and maximum nvar
int dmax = 0;
int nvmax = 0;
for (int i = 0; i < num_tests; i++)
{
if (dmax < all_tests[i]->dim)
dmax = all_tests[i]->dim;
if (nvmax < all_tests[i]->nvar)
nvmax = all_tests[i]->nvar;
}
tls.init(dmax, nvmax); // initialize library
// launch the tests
int success = 0;
for (int i = 0; i < num_tests; i++)
{
try
{
if (all_tests[i]->test())
success++;
}
catch (const TLException &e)
{
printf("Caugth TL exception in <%s>:\n", all_tests[i]->getName());
e.print();
}
catch (SylvException &e)
{
printf("Caught Sylv exception in <%s>:\n", all_tests[i]->getName());
e.printMessage();
}
}
printf("There were %d tests that failed out of %d tests run.\n",
num_tests - success, num_tests);
// destroy
for (int i = 0; i < num_tests; i++)
{
delete all_tests[i];
}
return 0;
}
dynare++/parser/cc/Makefile.am deleted 100644 → 0
View file @ 219d2bb3
noinst_LIBRARIES = libparser.a
GENERATED_FILES = assign_tab.cc csv_tab.cc formula_tab.cc matrix_tab.cc namelist_tab.cc assign_tab.hh csv_tab.hh formula_tab.hh matrix_tab.hh namelist_tab.hh assign_ll.cc csv_ll.cc formula_ll.cc matrix_ll.cc namelist_ll.cc
libparser_a_SOURCES = \
location.hh \
namelist.hh \
atom_assignings.cc \
atom_assignings.hh \
atom_substitutions.cc \
atom_substitutions.hh \
csv_parser.cc \
csv_parser.hh \
dynamic_atoms.cc \
dynamic_atoms.hh \
fine_atoms.cc \
fine_atoms.hh \
formula_parser.cc \
formula_parser.hh \
matrix_parser.cc \
matrix_parser.hh \
parser_exception.cc \
parser_exception.hh \
static_atoms.cc \
static_atoms.hh \
static_fine_atoms.cc \
static_fine_atoms.hh \
tree.cc \
tree.hh \
$(GENERATED_FILES)
libparser_a_CPPFLAGS = -I../.. $(BOOST_CPPFLAGS)
BUILT_SOURCES = $(GENERATED_FILES)
EXTRA_DIST = assign.yy csv.yy formula.yy matrix.yy namelist.yy assign.ll csv.ll formula.ll matrix.ll namelist.ll
%_tab.cc %_tab.hh: %.yy
$(YACC) -d -o$*_tab.cc $<
%_ll.cc: %.ll
$(LEX) -i -o$@ $<
dynare++/parser/cc/assign.ll deleted 100644 → 0
View file @ 219d2bb3
/* -*- C++ -*- */
%{
#include "location.hh"
#include "assign_tab.hh"
extern YYLTYPE asgn_lloc;
#define YY_USER_ACTION SET_LLOC(asgn_);
%}
%option nounput
%option noyy_top_state
%option stack
%option yylineno
%option prefix="asgn_"
%option never-interactive
%x CMT
%%
/* comments */
<*>"/*" {yy_push_state(CMT);}
<CMT>[^*\n]*
<CMT>"*"+[^*/\n]*
<CMT>"*"+"/" {yy_pop_state();}
<CMT>[\n]
"//".*\n
/* spaces */
[ \t\r\n] {return BLANK;}
/* names */
[A-Za-z_][A-Za-z0-9_]* {
asgn_lval.string = asgn_text;
return NAME;
}
; {return SEMICOLON;}
= {return EQUAL_SIGN;}
. {
asgn_lval.character = asgn_text[0];
return CHARACTER;
}
%%
int asgn_wrap()
{
return 1;
}
void asgn__destroy_buffer(void* p)
{
asgn__delete_buffer((YY_BUFFER_STATE)p);
}