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
Loading items

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
  • ebenetce/dynare
  • chskcau/dynare-doc-fixes
28 results
Select Git revision
Loading items
Show changes
Showing
with 0 additions and 5742 deletions
/*
* Copyright © 2008-2011 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "parser/cc/matrix_parser.hh"
#include "utils/cc/exception.hh"
#include "sylv/cc/GeneralMatrix.hh"
#include "sylv/cc/Vector.hh"
#include "sylv/cc/SymSchurDecomp.hh"
#include "sylv/cc/SylvException.hh"
#include "integ/cc/quadrature.hh"
#include "integ/cc/smolyak.hh"
#include "integ/cc/product.hh"
#include <getopt.h>
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <sstream>
#include <memory>
#include <string>
struct QuadParams
{
std::string outname;
std::string vcovname;
int max_level{3};
double discard_weight{0.0};
QuadParams(int argc, char **argv);
void check_consistency() const;
private:
enum class opt {max_level, discard_weight, vcov};
};
QuadParams::QuadParams(int argc, char **argv)
{
if (argc == 1)
{
// Print the help and exit
std::cerr << "Usage: " << argv[0] << " [--max-level INTEGER] [--discard-weight FLOAT] [--vcov FILENAME] OUTPUT_FILENAME" << std::endl;
std::exit(EXIT_FAILURE);
}
outname = argv[argc-1];
argc--;
struct option const opts [] = {
{"max-level", required_argument, nullptr, static_cast<int>(opt::max_level)},
{"discard-weight", required_argument, nullptr, static_cast<int>(opt::discard_weight)},
{"vcov", required_argument, nullptr, static_cast<int>(opt::vcov)},
{nullptr, 0, nullptr, 0}
};
int ret;
int index;
while (-1 != (ret = getopt_long(argc, argv, "", opts, &index)))
{
if (ret == '?')
{
std::cerr << "Unknown option, ignored\n";
continue;
}
switch (static_cast<opt>(ret))
{
case opt::max_level:
try
{
max_level = std::stoi(optarg);
}
catch (const std::invalid_argument &e)
{
std::cerr << "Couldn't parse integer " << optarg << ", ignored" << std::endl;
}
break;
case opt::discard_weight:
try
{
discard_weight = std::stod(optarg);
}
catch (const std::invalid_argument &e)
{
std::cerr << "Couldn't parse float " << optarg << ", ignored" << std::endl;
}
break;
case opt::vcov:
vcovname = optarg;
break;
}
}
check_consistency();
}
void
QuadParams::check_consistency() const
{
if (outname.empty())
{
std::cerr << "Error: output name not set" << std::endl;
std::exit(EXIT_FAILURE);
}
if (vcovname.empty())
{
std::cerr << "Error: vcov file name not set" << std::endl;
std::exit(EXIT_FAILURE);
}
}
int
main(int argc, char **argv)
{
QuadParams params(argc, argv);
// Open output file for writing
std::ofstream fout{params.outname, std::ios::out | std::ios::trunc};
if (fout.fail())
{
std::cerr << "Could not open " << params.outname << " for writing" << std::endl;
std::exit(EXIT_FAILURE);
}
try
{
std::ifstream f{params.vcovname};
std::ostringstream buffer;
buffer << f.rdbuf();
std::string contents{buffer.str()};
// Parse the vcov matrix
ogp::MatrixParser mp;
mp.parse(contents);
if (mp.nrows() != mp.ncols())
throw ogu::Exception(__FILE__, __LINE__,
"VCOV matrix not square");
// And put to the GeneralMatrix
GeneralMatrix vcov(mp.nrows(), mp.ncols());
vcov.zeros();
for (ogp::MPIterator it = mp.begin(); it != mp.end(); ++it)
vcov.get(it.row(), it.col()) = *it;
// Calculate the factor A of vcov, so that A·Aᵀ=VCOV
GeneralMatrix A(vcov.nrows(), vcov.nrows());
SymSchurDecomp ssd(vcov);
ssd.getFactor(A);
// Construct Gauss-Hermite quadrature
GaussHermite ghq;
// Construct Smolyak quadrature
int level = params.max_level;
SmolyakQuadrature sq(vcov.nrows(), level, ghq);
std::cout << "Dimension: " << vcov.nrows() << std::endl
<< "Maximum level: " << level << std::endl
<< "Total number of nodes: " << sq.numEvals(level) << std::endl;
// Put the points to the vector
std::vector<std::unique_ptr<Vector>> points;
for (smolpit qit = sq.start(level); qit != sq.end(level); ++qit)
points.push_back(std::make_unique<Vector>(const_cast<const Vector &>(qit.point())));
// Sort and uniq
std::sort(points.begin(), points.end(), [](auto &a, auto &b) { return a.get() < b.get(); });
auto new_end = std::unique(points.begin(), points.end());
points.erase(new_end, points.end());
std::cout << "Duplicit nodes removed: " << static_cast<unsigned long>(sq.numEvals(level)-points.size())
<< std::endl;
// Calculate weights and mass
double mass = 0.0;
std::vector<double> weights;
for (auto & point : points)
{
weights.push_back(std::exp(-point->dot(*point)));
mass += weights.back();
}
// Calculate discarded mass
double discard_mass = 0.0;
for (double weight : weights)
if (weight/mass < params.discard_weight)
discard_mass += weight;
std::cout << "Total mass discarded: " << std::fixed << discard_mass/mass << std::endl;
// Dump the results
int npoints = 0;
double upscale_weight = 1/(mass-discard_mass);
Vector x(vcov.nrows());
fout << std::setprecision(16);
for (int i = 0; i < static_cast<int>(weights.size()); i++)
if (weights[i]/mass >= params.discard_weight)
{
// Print the upscaled weight
fout << std::setw(20) << upscale_weight*weights[i];
// Multiply point with the factor A and √2
A.multVec(0.0, x, std::sqrt(2.), *(points[i]));
// Print the coordinates
for (int j = 0; j < x.length(); j++)
fout << ' ' << std::setw(20) << x[j];
fout << std::endl;
npoints++;
}
std::cout << "Final number of points: " << npoints << std::endl;
fout.close();
}
catch (const SylvException &e)
{
e.printMessage();
return EXIT_FAILURE;
}
catch (const ogu::Exception &e)
{
e.print();
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
check_PROGRAMS = tests
tests_SOURCES = tests.cc
tests_CPPFLAGS = -I../cc -I../../tl/cc -I../../sylv/cc -I../../utils/cc -I$(top_srcdir)/mex/sources
tests_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
tests_LDFLAGS = $(AM_LDFLAGS) $(LDFLAGS_MATIO)
tests_LDADD = ../../sylv/cc/libsylv.a ../cc/libinteg.a ../../tl/cc/libtl.a ../../utils/cc/libutils.a $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(LIBADD_MATIO)
check-local:
./tests
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "GeneralMatrix.hh"
#include <dynlapack.h>
#include "SylvException.hh"
#include "rfs_tensor.hh"
#include "normal_moments.hh"
#include "vector_function.hh"
#include "quadrature.hh"
#include "smolyak.hh"
#include "product.hh"
#include "quasi_mcarlo.hh"
#include <iomanip>
#include <chrono>
#include <cmath>
#include <iostream>
#include <utility>
#include <array>
#include <memory>
#include <cstdlib>
/* Evaluates unfolded (Dx)ᵏ power, where x is a vector, D is a Cholesky factor
(lower triangular) */
class MomentFunction : public VectorFunction
{
GeneralMatrix D;
int k;
public:
MomentFunction(const GeneralMatrix &inD, int kk)
: VectorFunction(inD.nrows(), UFSTensor::calcMaxOffset(inD.nrows(), kk)),
D(inD), k(kk)
{
}
MomentFunction(const MomentFunction &func) = default;
std::unique_ptr<VectorFunction>
clone() const override
{
return std::make_unique<MomentFunction>(*this);
}
void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
};
void
MomentFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{
if (point.length() != indim() || out.length() != outdim())
{
std::cerr << "Wrong length of vectors in MomentFunction::eval" << std::endl;
std::exit(EXIT_FAILURE);
}
Vector y(point);
y.zeros();
D.multaVec(y, point);
URSingleTensor ypow(y, k);
out.zeros();
out.add(1.0, ypow.getData());
}
class TensorPower : public VectorFunction
{
int k;
public:
TensorPower(int nvar, int kk)
: VectorFunction(nvar, UFSTensor::calcMaxOffset(nvar, kk)), k(kk)
{
}
TensorPower(const TensorPower &func) = default;
std::unique_ptr<VectorFunction>
clone() const override
{
return std::make_unique<TensorPower>(*this);
}
void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
};
void
TensorPower::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{
if (point.length() != indim() || out.length() != outdim())
{
std::cerr << "Wrong length of vectors in TensorPower::eval" << std::endl;
std::exit(EXIT_FAILURE);
}
URSingleTensor ypow(point, k);
out.zeros();
out.add(1.0, ypow.getData());
}
/* Evaluates (1+1/d)ᵈ(x₁·…·x_d)^(1/d), its integral over [0,1]ᵈ
is 1.0, and its variation grows exponentially */
class Function1 : public VectorFunction
{
int dim;
public:
Function1(int d)
: VectorFunction(d, 1), dim(d)
{
}
Function1(const Function1 &f)
: VectorFunction(f.indim(), f.outdim()), dim(f.dim)
{
}
std::unique_ptr<VectorFunction>
clone() const override
{
return std::make_unique<Function1>(*this);
}
void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
};
void
Function1::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{
if (point.length() != dim || out.length() != 1)
{
std::cerr << "Wrong length of vectors in Function1::eval" << std::endl;
std::exit(EXIT_FAILURE);
}
double r = 1;
for (int i = 0; i < dim; i++)
r *= point[i];
r = pow(r, 1.0/dim);
r *= pow(1.0 + 1.0/dim, static_cast<double>(dim));
out[0] = r;
}
// Evaluates Function1 but with transformation xᵢ=0.5(yᵢ+1)
// This makes the new function integrate over [−1,1]ᵈ to 1.0
class Function1Trans : public Function1
{
public:
Function1Trans(int d)
: Function1(d)
{
}
Function1Trans(const Function1Trans &func) = default;
std::unique_ptr<VectorFunction>
clone() const override
{
return std::make_unique<Function1Trans>(*this);
}
void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
};
void
Function1Trans::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
{
Vector p(point.length());
for (int i = 0; i < p.length(); i++)
p[i] = 0.5*(point[i]+1);
Function1::eval(p, sig, out);
out.mult(pow(0.5, indim()));
}
/* WallTimer class. Constructor saves the wall time, destructor cancels the
current time from the saved, and prints the message with time information */
class WallTimer
{
std::string mes;
std::chrono::time_point<std::chrono::high_resolution_clock> start;
bool new_line;
public:
WallTimer(std::string m, bool nl = true)
: mes{m}, start{std::chrono::high_resolution_clock::now()}, new_line{nl}
{
}
~WallTimer()
{
auto end = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> duration = end - start;
std::cout << mes << std::setw(8) << std::setprecision(4) << duration.count();
if (new_line)
std::cout << std::endl;
}
};
/****************************************************/
/* declaration of TestRunnable class */
/****************************************************/
class TestRunnable
{
public:
const std::string name;
int dim; // dimension of the solved problem
int nvar; // number of variables of the solved problem
TestRunnable(std::string name_arg, int d, int nv)
: name{move(name_arg)}, dim(d), nvar(nv)
{
}
virtual ~TestRunnable() = default;
bool test() const;
virtual bool run() const = 0;
protected:
static bool smolyak_normal_moments(const GeneralMatrix &m, int imom, int level);
static bool product_normal_moments(const GeneralMatrix &m, int imom, int level);
static bool qmc_normal_moments(const GeneralMatrix &m, int imom, int level);
static bool smolyak_product_cube(const VectorFunction &func, const Vector &res,
double tol, int level);
static bool qmc_cube(const VectorFunction &func, double res, double tol, int level);
};
bool
TestRunnable::test() const
{
std::cout << "Running test <" << name << ">" << std::endl;
bool passed;
{
WallTimer tim("Wall clock time ", false);
passed = run();
}
if (passed)
{
std::cout << "............................ passed" << std::endl << std::endl;
return passed;
}
else
{
std::cout << "............................ FAILED" << std::endl << std::endl;
return passed;
}
}
/****************************************************/
/* definition of TestRunnable static methods */
/****************************************************/
bool
TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level)
{
// First make m·mᵀ and then Cholesky factor
GeneralMatrix msq(m * transpose(m));
// Make vector function
int dim = m.nrows();
TensorPower tp(dim, imom);
GaussConverterFunction func(tp, msq);
// Smolyak quadrature
Vector smol_out(UFSTensor::calcMaxOffset(dim, imom));
{
WallTimer tim("\tSmolyak quadrature time: ");
GaussHermite gs;
SmolyakQuadrature quad(dim, level, gs);
quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, smol_out);
std::cout << "\tNumber of Smolyak evaluations: " << quad.numEvals(level) << std::endl;
}
// Check against theoretical moments
UNormalMoments moments(imom, msq);
smol_out.add(-1.0, moments.get(Symmetry{imom}).getData());
std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << smol_out.getMax() << std::endl;
return smol_out.getMax() < 1.e-7;
}
bool
TestRunnable::product_normal_moments(const GeneralMatrix &m, int imom, int level)
{
// First make m·mᵀ and then Cholesky factor
GeneralMatrix msq(m * transpose(m));
// Make vector function
int dim = m.nrows();
TensorPower tp(dim, imom);
GaussConverterFunction func(tp, msq);
// Product quadrature
Vector prod_out(UFSTensor::calcMaxOffset(dim, imom));
{
WallTimer tim("\tProduct quadrature time: ");
GaussHermite gs;
ProductQuadrature quad(dim, gs);
quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, prod_out);
std::cout << "\tNumber of product evaluations: " << quad.numEvals(level) << std::endl;
}
// Check against theoretical moments
UNormalMoments moments(imom, msq);
prod_out.add(-1.0, moments.get(Symmetry{imom}).getData());
std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << prod_out.getMax() << std::endl;
return prod_out.getMax() < 1.e-7;
}
bool
TestRunnable::smolyak_product_cube(const VectorFunction &func, const Vector &res,
double tol, int level)
{
if (res.length() != func.outdim())
{
std::cerr << "Incompatible dimensions of check value and function." << std::endl;
std::exit(EXIT_FAILURE);
}
GaussLegendre glq;
Vector out(func.outdim());
double smol_error;
double prod_error;
{
WallTimer tim("\tSmolyak quadrature time: ");
SmolyakQuadrature quad(func.indim(), level, glq);
quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, out);
out.add(-1.0, res);
smol_error = out.getMax();
std::cout << "\tNumber of Smolyak evaluations: " << quad.numEvals(level) << std::endl;
std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << smol_error << std::endl;
}
{
WallTimer tim("\tProduct quadrature time: ");
ProductQuadrature quad(func.indim(), glq);
quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, out);
out.add(-1.0, res);
prod_error = out.getMax();
std::cout << "\tNumber of product evaluations: " << quad.numEvals(level) << std::endl;
std::cout << "\tError: " << std::setw(16) << std::setprecision(12) << prod_error << std::endl;
}
return smol_error < tol && prod_error < tol;
}
bool
TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int level)
{
Vector r(1);
double error1;
{
WallTimer tim("\tQuasi-Monte Carlo (Warnock scrambling) time: ");
WarnockPerScheme wps;
QMCarloCubeQuadrature qmc(func.indim(), level, wps);
qmc.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, r);
error1 = std::max(res - r[0], r[0] - res);
std::cout << "\tQuasi-Monte Carlo (Warnock scrambling) error: " << std::setw(16) << std::setprecision(12) << error1 << std::endl;
}
double error2;
{
WallTimer tim("\tQuasi-Monte Carlo (reverse scrambling) time: ");
ReversePerScheme rps;
QMCarloCubeQuadrature qmc(func.indim(), level, rps);
qmc.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, r);
error2 = std::max(res - r[0], r[0] - res);
std::cout << "\tQuasi-Monte Carlo (reverse scrambling) error: " << std::setw(16) << std::setprecision(12) << error2 << std::endl;
}
double error3;
{
WallTimer tim("\tQuasi-Monte Carlo (no scrambling) time: ");
IdentityPerScheme ips;
QMCarloCubeQuadrature qmc(func.indim(), level, ips);
qmc.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, r);
error3 = std::max(res - r[0], r[0] - res);
std::cout << "\tQuasi-Monte Carlo (no scrambling) error: " << std::setw(16) << std::setprecision(12) << error3 << std::endl;
}
return error1 < tol && error2 < tol && error3 < tol;
}
/****************************************************/
/* definition of TestRunnable subclasses */
/****************************************************/
class SmolyakNormalMom1 : public TestRunnable
{
public:
SmolyakNormalMom1()
: TestRunnable("Smolyak normal moments (dim=2, level=4, order=4)", 4, 2)
{
}
bool
run() const override
{
GeneralMatrix m(2, 2);
m.zeros();
m.get(0, 0) = 1;
m.get(1, 1) = 1;
return smolyak_normal_moments(m, 4, 4);
}
};
class SmolyakNormalMom2 : public TestRunnable
{
public:
SmolyakNormalMom2()
: TestRunnable("Smolyak normal moments (dim=3, level=8, order=8)", 8, 3)
{
}
bool
run() const override
{
GeneralMatrix m(3, 3);
m.zeros();
m.get(0, 0) = 1;
m.get(0, 2) = 0.5;
m.get(1, 1) = 1;
m.get(1, 0) = 0.5;
m.get(2, 2) = 2;
m.get(2, 1) = 4;
return smolyak_normal_moments(m, 8, 8);
}
};
class ProductNormalMom1 : public TestRunnable
{
public:
ProductNormalMom1()
: TestRunnable("Product normal moments (dim=2, level=4, order=4)", 4, 2)
{
}
bool
run() const override
{
GeneralMatrix m(2, 2);
m.zeros();
m.get(0, 0) = 1;
m.get(1, 1) = 1;
return product_normal_moments(m, 4, 4);
}
};
class ProductNormalMom2 : public TestRunnable
{
public:
ProductNormalMom2()
: TestRunnable("Product normal moments (dim=3, level=8, order=8)", 8, 3)
{
}
bool
run() const override
{
GeneralMatrix m(3, 3);
m.zeros();
m.get(0, 0) = 1;
m.get(0, 2) = 0.5;
m.get(1, 1) = 1;
m.get(1, 0) = 0.5;
m.get(2, 2) = 2;
m.get(2, 1) = 4;
return product_normal_moments(m, 8, 8);
}
};
// Note that here we pass 1,1 to tls since smolyak has its own PascalTriangle
class F1GaussLegendre : public TestRunnable
{
public:
F1GaussLegendre()
: TestRunnable("Function1 Gauss-Legendre (dim=6, level=13", 1, 1)
{
}
bool
run() const override
{
Function1Trans f1(6);
Vector res(1);
res[0] = 1.0;
return smolyak_product_cube(f1, res, 1e-2, 13);
}
};
class F1QuasiMCarlo : public TestRunnable
{
public:
F1QuasiMCarlo()
: TestRunnable("Function1 Quasi-Monte Carlo (dim=6, level=1000000)", 1, 1)
{
}
bool
run() const override
{
Function1 f1(6);
return qmc_cube(f1, 1.0, 1.e-4, 1000000);
}
};
int
main()
{
std::vector<std::unique_ptr<TestRunnable>> all_tests;
// Fill in vector of all tests
all_tests.push_back(std::make_unique<SmolyakNormalMom1>());
all_tests.push_back(std::make_unique<SmolyakNormalMom2>());
all_tests.push_back(std::make_unique<ProductNormalMom1>());
all_tests.push_back(std::make_unique<ProductNormalMom2>());
all_tests.push_back(std::make_unique<F1GaussLegendre>());
all_tests.push_back(std::make_unique<F1QuasiMCarlo>());
// Find maximum dimension and maximum nvar
int dmax = 0;
int nvmax = 0;
for (const auto &test : all_tests)
{
dmax = std::max(dmax, test->dim);
nvmax = std::max(nvmax, test->nvar);
}
TLStatic::init(dmax, nvmax); // initialize library
// Launch the tests
int success = 0;
for (const auto &test : all_tests)
{
try
{
if (test->test())
success++;
}
catch (const TLException &e)
{
std::cout << "Caught TL exception in <" << test->name << ">:" << std::endl;
e.print();
}
catch (SylvException &e)
{
std::cout << "Caught Sylv exception in <" << test->name << ">:" << std::endl;
e.printMessage();
}
}
int nfailed = all_tests.size() - success;
std::cout << "There were " << nfailed << " tests that failed out of "
<< all_tests.size() << " tests run." << std::endl;
if (nfailed)
return EXIT_FAILURE;
else
return EXIT_SUCCESS;
}
noinst_LIBRARIES = libkord.a
libkord_a_SOURCES = \
approximation.cc \
approximation.hh \
decision_rule.cc \
decision_rule.hh \
dynamic_model.cc \
dynamic_model.hh \
faa_di_bruno.cc \
faa_di_bruno.hh \
first_order.cc \
first_order.hh \
global_check.cc \
global_check.hh \
kord_exception.hh \
korder.cc \
korder.hh \
korder_stoch.cc \
korder_stoch.hh \
journal.cc \
journal.hh \
normal_conjugate.cc \
normal_conjugate.hh \
seed_generator.cc \
seed_generator.hh
libkord_a_CPPFLAGS = -I../sylv/cc -I../tl/cc -I../integ/cc -I../utils/cc -I$(top_srcdir)/mex/sources $(CPPFLAGS_MATIO) -DDYNVERSION=\"$(PACKAGE_VERSION)\"
libkord_a_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
check_PROGRAMS = tests
tests_SOURCES = tests.cc
tests_CPPFLAGS = -I../sylv/cc -I../tl/cc -I../integ/cc -I../utils/cc -I$(top_srcdir)/mex/sources
tests_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
tests_LDFLAGS = $(AM_LDFLAGS) $(LDFLAGS_MATIO)
tests_LDADD = libkord.a ../tl/cc/libtl.a ../sylv/cc/libsylv.a ../utils/cc/libutils.a $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(LIBADD_MATIO)
check-local:
./tests
CLEANFILES = out.txt
/*
* Copyright © 2004 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "kord_exception.hh"
#include "decision_rule.hh"
#include "dynamic_model.hh"
#include "seed_generator.hh"
#include "SymSchurDecomp.hh"
#include <dynlapack.h>
#include <limits>
#include <utility>
#include <memory>
// FoldDecisionRule conversion from UnfoldDecisionRule
FoldDecisionRule::FoldDecisionRule(const UnfoldDecisionRule &udr)
: DecisionRuleImpl<Storage::fold>(ctraits<Storage::fold>::Tpol(udr.nrows(), udr.nvars()),
udr.ypart, udr.nu, udr.ysteady)
{
for (const auto &it : udr)
insert(std::make_unique<ctraits<Storage::fold>::Ttensym>(*(it.second)));
}
// UnfoldDecisionRule conversion from FoldDecisionRule
UnfoldDecisionRule::UnfoldDecisionRule(const FoldDecisionRule &fdr)
: DecisionRuleImpl<Storage::unfold>(ctraits<Storage::unfold>::Tpol(fdr.nrows(), fdr.nvars()),
fdr.ypart, fdr.nu, fdr.ysteady)
{
for (const auto &it : fdr)
insert(std::make_unique<ctraits<Storage::unfold>::Ttensym>(*(it.second)));
}
/* This runs simulations with an output to journal file. Note that we
report how many simulations had to be thrown out due to Nan or Inf. */
void
SimResults::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
const TwoDMatrix &vcov, Journal &journal)
{
JournalRecordPair paa(journal);
paa << "Performing " << num_sim << " stochastic simulations for "
<< num_per << " periods burning " << num_burn << " initial periods" << endrec;
simulate(num_sim, dr, start, vcov);
int thrown = num_sim - data.size();
if (thrown > 0)
{
JournalRecord rec(journal);
rec << "I had to throw " << thrown << " simulations away due to Nan or Inf" << endrec;
}
}
/* This runs a given number of simulations by creating
SimulationWorker for each simulation and inserting them to the
thread group. */
void
SimResults::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
const TwoDMatrix &vcov)
{
std::vector<RandomShockRealization> rsrs;
rsrs.reserve(num_sim);
sthread::detach_thread_group gr;
for (int i = 0; i < num_sim; i++)
{
RandomShockRealization sr(vcov, seed_generator::get_new_seed());
rsrs.push_back(sr);
gr.insert(std::make_unique<SimulationWorker>(*this, dr, DecisionRule::emethod::horner,
num_per+num_burn, start, rsrs.back()));
}
gr.run();
}
/* This adds the data with the realized shocks. It takes only periods
which are not to be burnt. If the data is not finite, the both data
and shocks are thrown away. */
bool
SimResults::addDataSet(const TwoDMatrix &d, const ExplicitShockRealization &sr, const ConstVector &st)
{
KORD_RAISE_IF(d.nrows() != num_y,
"Incompatible number of rows for SimResults::addDataSets");
KORD_RAISE_IF(d.ncols() != num_per+num_burn,
"Incompatible number of cols for SimResults::addDataSets");
bool ret = false;
if (d.isFinite())
{
data.emplace_back(d, num_burn, num_per);
shocks.emplace_back(ConstTwoDMatrix(sr.getShocks(), num_burn, num_per));
if (num_burn == 0)
start.emplace_back(st);
else
start.emplace_back(d.getCol(num_burn-1));
ret = true;
}
return ret;
}
void
SimResults::writeMat(const std::string &base, const std::string &lname) const
{
std::string matfile_name = base + ".mat";
mat_t *matfd = Mat_Create(matfile_name.c_str(), nullptr);
if (matfd)
{
writeMat(matfd, lname);
Mat_Close(matfd);
}
}
/* This save the results as matrices with given prefix and with index
appended. If there is only one matrix, the index is not appended. */
void
SimResults::writeMat(mat_t *fd, const std::string &lname) const
{
for (int i = 0; i < getNumSets(); i++)
{
std::string tmp = lname + "_data";
if (getNumSets() > 1)
tmp += std::to_string(i+1);
data[i].writeMat(fd, tmp);
}
}
void
SimResultsStats::simulate(int num_sim, const DecisionRule &dr,
const Vector &start,
const TwoDMatrix &vcov, Journal &journal)
{
SimResults::simulate(num_sim, dr, start, vcov, journal);
{
JournalRecordPair paa(journal);
paa << "Calculating means from the simulations." << endrec;
calcMean();
}
{
JournalRecordPair paa(journal);
paa << "Calculating covariances from the simulations." << endrec;
calcVcov();
}
}
/* Here we do not save the data itself, we save only mean and vcov. */
void
SimResultsStats::writeMat(mat_t *fd, const std::string &lname) const
{
ConstTwoDMatrix(num_y, 1, mean).writeMat(fd, lname + "_mean");;
vcov.writeMat(fd, lname + "_vcov");
}
void
SimResultsStats::calcMean()
{
mean.zeros();
if (data.size()*num_per > 0)
{
double mult = 1.0/data.size()/num_per;
for (const auto &i : data)
{
for (int j = 0; j < num_per; j++)
{
ConstVector col{i.getCol(j)};
mean.add(mult, col);
}
}
}
}
void
SimResultsStats::calcVcov()
{
if (data.size()*num_per > 1)
{
vcov.zeros();
double mult = 1.0/(data.size()*num_per - 1);
for (const auto &d : data)
for (int j = 0; j < num_per; j++)
for (int m = 0; m < num_y; m++)
for (int n = m; n < num_y; n++)
{
double s = (d.get(m, j)-mean[m])*(d.get(n, j)-mean[n]);
vcov.get(m, n) += mult*s;
if (m != n)
vcov.get(n, m) += mult*s;
}
}
else
vcov.infs();
}
void
SimResultsDynamicStats::simulate(int num_sim, const DecisionRule &dr,
const Vector &start,
const TwoDMatrix &vcov, Journal &journal)
{
SimResults::simulate(num_sim, dr, start, vcov, journal);
{
JournalRecordPair paa(journal);
paa << "Calculating means of the conditional simulations." << endrec;
calcMean();
}
{
JournalRecordPair paa(journal);
paa << "Calculating variances of the conditional simulations." << endrec;
calcVariance();
}
}
void
SimResultsDynamicStats::writeMat(mat_t *fd, const std::string &lname) const
{
mean.writeMat(fd, lname + "_cond_mean");
variance.writeMat(fd, lname + "_cond_variance");
}
void
SimResultsDynamicStats::calcMean()
{
mean.zeros();
if (data.size() > 0)
{
double mult = 1.0/data.size();
for (int j = 0; j < num_per; j++)
{
Vector meanj{mean.getCol(j)};
for (const auto &i : data)
{
ConstVector col{i.getCol(j)};
meanj.add(mult, col);
}
}
}
}
void
SimResultsDynamicStats::calcVariance()
{
if (data.size() > 1)
{
variance.zeros();
double mult = 1.0/(data.size()-1);
for (int j = 0; j < num_per; j++)
{
ConstVector meanj{mean.getCol(j)};
Vector varj{variance.getCol(j)};
for (const auto &i : data)
{
Vector col{i.getCol(j)};
col.add(-1.0, meanj);
for (int k = 0; k < col.length(); k++)
col[k] = col[k]*col[k];
varj.add(mult, col);
}
}
}
else
variance.infs();
}
void
SimResultsIRF::simulate(const DecisionRule &dr, Journal &journal)
{
JournalRecordPair paa(journal);
paa << "Performing " << control.getNumSets() << " IRF simulations for "
<< num_per << " periods; shock=" << ishock << ", impulse=" << imp << endrec;
simulate(dr);
int thrown = control.getNumSets() - data.size();
if (thrown > 0)
{
JournalRecord rec(journal);
rec << "I had to throw " << thrown
<< " simulations away due to Nan or Inf" << endrec;
}
calcMeans();
calcVariances();
}
void
SimResultsIRF::simulate(const DecisionRule &dr)
{
sthread::detach_thread_group gr;
for (int idata = 0; idata < control.getNumSets(); idata++)
gr.insert(std::make_unique<SimulationIRFWorker>(*this, dr, DecisionRule::emethod::horner,
num_per, idata, ishock, imp));
gr.run();
}
void
SimResultsIRF::calcMeans()
{
means.zeros();
if (data.size() > 0)
{
for (const auto &i : data)
means.add(1.0, i);
means.mult(1.0/data.size());
}
}
void
SimResultsIRF::calcVariances()
{
if (data.size() > 1)
{
variances.zeros();
for (const auto &i : data)
{
TwoDMatrix d(i);
d.add(-1.0, means);
for (int j = 0; j < d.nrows(); j++)
for (int k = 0; k < d.ncols(); k++)
variances.get(j, k) += d.get(j, k)*d.get(j, k);
d.mult(1.0/(data.size()-1));
}
}
else
variances.infs();
}
void
SimResultsIRF::writeMat(mat_t *fd, const std::string &lname) const
{
means.writeMat(fd, lname + "_mean");
variances.writeMat(fd, lname + "_var");
}
void
RTSimResultsStats::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
const TwoDMatrix &v, Journal &journal)
{
JournalRecordPair paa(journal);
paa << "Performing " << num_sim << " real-time stochastic simulations for "
<< num_per << " periods" << endrec;
simulate(num_sim, dr, start, v);
mean = nc.getMean();
mean.add(1.0, dr.getSteady());
nc.getVariance(vcov);
if (thrown_periods > 0)
{
JournalRecord rec(journal);
rec << "I had to throw " << thrown_periods << " periods away due to Nan or Inf" << endrec;
JournalRecord rec1(journal);
rec1 << "This affected " << incomplete_simulations << " out of "
<< num_sim << " simulations" << endrec;
}
}
void
RTSimResultsStats::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
const TwoDMatrix &vcov)
{
std::vector<RandomShockRealization> rsrs;
rsrs.reserve(num_sim);
sthread::detach_thread_group gr;
for (int i = 0; i < num_sim; i++)
{
RandomShockRealization sr(vcov, seed_generator::get_new_seed());
rsrs.push_back(sr);
gr.insert(std::make_unique<RTSimulationWorker>(*this, dr, DecisionRule::emethod::horner,
num_per, start, rsrs.back()));
}
gr.run();
}
void
RTSimResultsStats::writeMat(mat_t *fd, const std::string &lname)
{
ConstTwoDMatrix(nc.getDim(), 1, mean).writeMat(fd, lname + "_rt_mean");
vcov.writeMat(fd, lname + "_rt_vcov");
}
IRFResults::IRFResults(const DynamicModel &mod, const DecisionRule &dr,
const SimResults &control, std::vector<int> ili,
Journal &journal)
: model(mod), irf_list_ind(std::move(ili))
{
int num_per = control.getNumPer();
JournalRecordPair pa(journal);
pa << "Calculating IRFs against control for " << static_cast<int>(irf_list_ind.size()) << " shocks and for "
<< num_per << " periods" << endrec;
const TwoDMatrix &vcov = mod.getVcov();
for (int ishock : irf_list_ind)
{
double stderror = sqrt(vcov.get(ishock, ishock));
irf_res.emplace_back(control, model.numeq(), num_per,
ishock, stderror);
irf_res.emplace_back(control, model.numeq(), num_per,
ishock, -stderror);
}
for (unsigned int ii = 0; ii < irf_list_ind.size(); ii++)
{
irf_res[2*ii].simulate(dr, journal);
irf_res[2*ii+1].simulate(dr, journal);
}
}
void
IRFResults::writeMat(mat_t *fd, const std::string &prefix) const
{
for (unsigned int i = 0; i < irf_list_ind.size(); i++)
{
int ishock = irf_list_ind[i];
auto shockname = model.getExogNames().getName(ishock);
irf_res[2*i].writeMat(fd, prefix + "_irfp_" + shockname);
irf_res[2*i+1].writeMat(fd, prefix + "_irfm_" + shockname);
}
}
void
SimulationWorker::operator()(std::mutex &mut)
{
ExplicitShockRealization esr(sr, np);
TwoDMatrix m{dr.simulate(em, np, st, esr)};
{
std::unique_lock<std::mutex> lk{mut};
res.addDataSet(m, esr, st);
}
}
/* Here we create a new instance of ExplicitShockRealization of the
corresponding control, add the impulse, and simulate. */
void
SimulationIRFWorker::operator()(std::mutex &mut)
{
ExplicitShockRealization esr(res.control.getShocks(idata));
esr.addToShock(ishock, 0, imp);
TwoDMatrix m{dr.simulate(em, np, res.control.getStart(idata), esr)};
m.add(-1.0, res.control.getData(idata));
{
std::unique_lock<std::mutex> lk{mut};
res.addDataSet(m, esr, res.control.getStart(idata));
}
}
void
RTSimulationWorker::operator()(std::mutex &mut)
{
NormalConj nc(res.nc.getDim());
const PartitionY &ypart = dr.getYPart();
int nu = dr.nexog();
const Vector &ysteady = dr.getSteady();
// initialize vectors and subvectors for simulation
Vector dyu(ypart.nys()+nu);
ConstVector ystart_pred(ystart, ypart.nstat, ypart.nys());
ConstVector ysteady_pred(ysteady, ypart.nstat, ypart.nys());
Vector dy(dyu, 0, ypart.nys());
Vector u(dyu, ypart.nys(), nu);
Vector y(nc.getDim());
ConstVector ypred(y, ypart.nstat, ypart.nys());
// simulate the first real-time period
int ip = 0;
dy = ystart_pred;
dy.add(-1.0, ysteady_pred);
sr.get(ip, u);
dr.eval(em, y, dyu);
if (ip >= res.num_burn)
nc.update(y);
// simulate other real-time periods
while (y.isFinite() && ip < res.num_burn + res.num_per)
{
ip++;
dy = ypred;
sr.get(ip, u);
dr.eval(em, y, dyu);
if (ip >= res.num_burn)
nc.update(y);
}
{
std::unique_lock<std::mutex> lk{mut};
res.nc.update(nc);
if (res.num_per-ip > 0)
{
res.incomplete_simulations++;
res.thrown_periods += res.num_per-ip;
}
}
}
/* This calculates factorization FFᵀ=V in the Cholesky way. It does
not work for semidefinite matrices. */
void
RandomShockRealization::choleskyFactor(const ConstTwoDMatrix &v)
{
factor = v;
lapack_int rows = factor.nrows(), lda = factor.getLD();
for (int i = 0; i < rows; i++)
for (int j = i+1; j < rows; j++)
factor.get(i, j) = 0.0;
lapack_int info;
dpotrf("L", &rows, factor.base(), &lda, &info);
KORD_RAISE_IF(info != 0,
"Info!=0 in RandomShockRealization::choleskyFactor");
}
/* This calculates FFᵀ=V factorization by symmetric Schur
decomposition. It works for semidefinite matrices. */
void
RandomShockRealization::schurFactor(const ConstTwoDMatrix &v)
{
SymSchurDecomp(v).getFactor(factor);
}
void
RandomShockRealization::get(int n, Vector &out)
{
KORD_RAISE_IF(out.length() != numShocks(),
"Wrong length of out vector in RandomShockRealization::get");
Vector d(out.length());
for (int i = 0; i < d.length(); i++)
d[i] = dis(mtwister);
out.zeros();
factor.multaVec(out, ConstVector(d));
}
ExplicitShockRealization::ExplicitShockRealization(ShockRealization &sr,
int num_per)
: shocks(sr.numShocks(), num_per)
{
for (int j = 0; j < num_per; j++)
{
Vector jcol{shocks.getCol(j)};
sr.get(j, jcol);
}
}
void
ExplicitShockRealization::get(int n, Vector &out)
{
KORD_RAISE_IF(out.length() != numShocks(),
"Wrong length of out vector in ExplicitShockRealization::get");
int i = n % shocks.ncols();
ConstVector icol{shocks.getCol(i)};
out = icol;
}
void
ExplicitShockRealization::addToShock(int ishock, int iper, double val)
{
KORD_RAISE_IF(ishock < 0 || ishock > numShocks(),
"Wrong index of shock in ExplicitShockRealization::addToShock");
int j = iper % shocks.ncols();
shocks.get(ishock, j) += val;
}
void
GenShockRealization::get(int n, Vector &out)
{
KORD_RAISE_IF(out.length() != numShocks(),
"Wrong length of out vector in GenShockRealization::get");
ExplicitShockRealization::get(n, out);
Vector r(numShocks());
RandomShockRealization::get(n, r);
for (int j = 0; j < numShocks(); j++)
if (!std::isfinite(out[j]))
out[j] = r[j];
}
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "SymSchurDecomp.hh"
#include "global_check.hh"
#include "seed_generator.hh"
#include "smolyak.hh"
#include "product.hh"
#include "quasi_mcarlo.hh"
#include <utility>
#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())
{
}
ResidFunction::ResidFunction(const ResidFunction &rf)
: VectorFunction(rf), approx(rf.approx), model(rf.model->clone())
{
if (rf.yplus)
yplus = std::make_unique<Vector>(*(rf.yplus));
if (rf.ystar)
ystar = std::make_unique<Vector>(*(rf.ystar));
if (rf.u)
u = std::make_unique<Vector>(*(rf.u));
if (rf.hss)
hss = std::make_unique<FTensorPolynomial>(*(rf.hss));
}
/* This sets y* and u. We have to create ‘ystar’, ‘u’, ‘yplus’ and ‘hss’. */
void
ResidFunction::setYU(const ConstVector &ys, const ConstVector &xx)
{
ystar = std::make_unique<Vector>(ys);
u = std::make_unique<Vector>(xx);
yplus = std::make_unique<Vector>(model->numeq());
approx.getFoldDecisionRule().evaluate(DecisionRule::emethod::horner,
*yplus, *ystar, *u);
// make a tensor polynomial of in-place subtensors from decision rule
/* Note that the non-const polynomial will be used for a construction of
‘hss’ and will be used in a const context. So this const cast is safe.
Note, that there is always a folded decision rule in Approximation. */
const FoldDecisionRule &dr = approx.getFoldDecisionRule();
FTensorPolynomial dr_ss(model->nstat()+model->npred(), model->nboth()+model->nforw(),
const_cast<FoldDecisionRule &>(dr));
// make ‘ytmp_star’ be a difference of ‘yplus’ from steady
Vector ytmp_star(ConstVector(*yplus, model->nstat(), model->npred()+model->nboth()));
ConstVector ysteady_star(dr.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 = std::make_unique<FTensorPolynomial>(dr_ss, ytmp_star);
ConstVector ysteady_ss(dr.getSteady(), model->nstat()+model->npred(),
model->nboth()+model->nforw());
if (hss->check(Symmetry{0}))
hss->get(Symmetry{0}).getData().add(1.0, ysteady_ss);
else
{
auto ten = std::make_unique<FFSTensor>(hss->nrows(), hss->nvars(), 0);
ten->getData() = ysteady_ss;
hss->insert(std::move(ten));
}
}
/* Here we evaluate the residual F(y*,u,u′). We have to evaluate ‘hss’ 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 𝔼[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++)
dynamic_cast<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 which 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);
std::unique_ptr<Quadrature> quad;
int lev;
// Create the quadrature and report the decision
if (take_smolyak)
{
quad = std::make_unique<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 = std::make_unique<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 columns 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);
}
}
/* This method checks an error of the approximation by evaluating residual
𝔼[F(y*,u,u′) | y*,u] for y* equal to the steady state, and changing u. We go
through all elements of u and vary them from −mult·σ to mult·σ in ‘m’
steps. */
void
GlobalChecker::checkAlongShocksAndSave(mat_t *fd, const std::string &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 = 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)};
for (int ishock = 0; ishock < model.nexog(); ishock++)
{
TwoDMatrix err_out(model.numeq(), 2*m+1);
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);
std::string shockname{model.getExogNames().getName(ishock)};
shockname.resize(8, ' ');
rec1 << shockname << ' ' << exo_mat.get(ishock, jj)
<< "\t" << error.getMax() << endrec;
}
err_out.writeMat(fd, prefix + "_shock_" + model.getExogNames().getName(ishock) + "_errors");
}
}
/* 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ₜ to zeros. Fourth we run the check() method and save the
results. */
void
GlobalChecker::checkOnEllipseAndSave(mat_t *fd, const std::string &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_cast<const TwoDMatrix &>(ycov), model.nstat(), model.nstat(),
model.npred()+model.nboth(), model.npred()+model.nboth());
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 [0,1]ᵈ
by QMCarloCubeQuadrature and make a polar transformation to the sphere.
The polar transformation fⁱ can be written recursively w.r.t. the
dimension i as:
f⁰() = [1]
⎡cos(2πxᵢ)·fⁱ⁻¹(x₁,…,xᵢ₋₁)⎤
fⁱ(x₁,…,xᵢ) = ⎣ sin(2πxᵢ) ⎦
*/
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);
ymat.writeMat(fd, prefix + "_ellipse_points");
out.writeMat(fd, prefix + "_ellipse_errors");
}
/* Here we check the errors along a simulation. We simulate, then set ‘x’ to
zeros, check and save results. */
void
GlobalChecker::checkAlongSimulationAndSave(mat_t *fd, const std::string &prefix,
int m, int max_evals)
{
JournalRecordPair pa(journal);
pa << "Calculating errors at " << m
<< " simulated points" << endrec;
RandomShockRealization sr(model.getVcov(), seed_generator::get_new_seed());
TwoDMatrix y{approx.getFoldDecisionRule().simulate(DecisionRule::emethod::horner,
m, model.getSteady(), sr)};
TwoDMatrix x(model.nexog(), m);
x.zeros();
TwoDMatrix out(model.numeq(), m);
check(max_evals, y, x, out);
y.writeMat(fd, prefix + "_simul_points");
out.writeMat(fd, prefix + "_simul_errors");
}
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
// Global check
/* The purpose of this file is to provide classes for checking error of
approximation. If yₜ=g(y*ₜ₋₁,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 the integral:
𝔼ₜ[F(y*,u,u′)]
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 ones:
— Along shocks. The y* is set to the 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 the interval [−3σ,3σ] (where σ is
one standard error of the shock), and a number of steps. This is repeated
for each shock (element of the u vector).
— Along simulation. Some random simulation is run, and for each realization
of y* and u along the path we evaluate the residual.
— On ellipse. Let V=AAᵀ be a covariance matrix of the predetermined
variables y* based on linear approximation, then we calculate integral for
points on the ellipse { Ax | ‖x‖₂=1 }. The points are selected by means of
low discrepancy method and polar transformation. The shock u are zeros.
— Unconditional distribution.
*/
#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 the 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
the setYU() method. The object has basically two states. One is after
construction and before the call to setYU(). The second is after the call to
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;
std::unique_ptr<DynamicModel> model;
std::unique_ptr<Vector> yplus, ystar, u;
std::unique_ptr<FTensorPolynomial> hss;
public:
ResidFunction(const Approximation &app);
ResidFunction(const ResidFunction &rf);
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())
{
}
std::unique_ptr<VectorFunction>
clone() const override
{
return std::make_unique<GResidFunction>(*this);
}
void
setYU(const ConstVector &ys, const ConstVector &xx)
{
dynamic_cast<ResidFunction *>(func)->setYU(ys, xx);
}
};
/* This is a class encapsulating checking algorithms. Its core routine is
check(), which calculates integral 𝔼[F(y*,u,u′) | 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 𝔼[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 std::string &prefix,
int m, double mult, int max_evals);
void checkOnEllipseAndSave(mat_t *fd, const std::string &prefix,
int m, double mult, int max_evals);
void checkAlongSimulationAndSave(mat_t *fd, const std::string &prefix,
int m, int max_evals);
void checkUnconditionalAndSave(mat_t *fd, const std::string &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
/*
* Copyright © 2007 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#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.nrows()), kappa(ydata.ncols()), nu(ydata.ncols()-1),
lambda(ydata.nrows(), ydata.nrows())
{
mu.zeros();
for (int i = 0; i < ydata.ncols(); i++)
mu.add(1.0/ydata.ncols(), ydata.getCol(i));
lambda.zeros();
for (int i = 0; i < ydata.ncols(); i++)
{
Vector diff{ydata.getCol(i)};
diff.add(-1, mu);
lambda.addOuter(diff);
}
}
// NormalConj::update() one observation code
/* The method performs the following:
κ₀ 1
μ₁ = ──── μ₀ + ──── y
κ₀+1 κ₀+1
κ₁ = κ₀ + 1
ν₁ = ν₀ + 1
κ₀
Λ₁ = Λ₀ + ──── (y − μ₀)(y − μ₀)ᵀ
κ₀+1
*/
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 = static_cast<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/(ν−d−1)·Λ, 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_cast<const TwoDMatrix &>(lambda);
v.mult(1.0/(nu-getDim()-1));
}
else
v.nans();
}
/*
* Copyright © 2007 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
// 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 algebra can be found in Gelman, Carlin, Stern, Rubin (p.87). It goes as
follows. Prior conjugate distribution takes the following form:
Σ ↝ InvWishart_ν₀(Λ₀⁻¹)
μ|Σ ↝ 𝒩(μ₀,Σ/κ₀)
If the observations are y₁…yₙ, then the posterior distribution has the same
form with the following parameters:
κ₀ n
μₙ = ──── μ₀ + ──── ȳ
κ₀+n κ₀+n
κₙ = κ₀ + n
νₙ = ν₀ + n
κ₀·n
Λₙ = Λ₀ + S + ──── (ȳ − μ₀)(ȳ − μ₀)ᵀ
κ₀+n
where
1 ₙ
ȳ = ─ ∑ yᵢ
n ⁱ⁼¹
S = ∑ (yᵢ − ȳ)(yᵢ − ȳ)ᵀ
ⁱ⁼¹
*/
#ifndef NORMAL_CONJUGATE_H
#define NORMAL_CONJUGATE_H
#include "twod_matrix.hh"
/* The class is described by the four parameters: μ, κ, ν and Λ. */
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 κ and Λ to zeros, ν to −1 and also the mean μ
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 &) = default;
NormalConj(NormalConj &&) = default;
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
/*
* Copyright © 2004 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <chrono>
#include <random>
#include <string>
#include <utility>
#include <iostream>
#include <iomanip>
#include <vector>
#include <memory>
#include "korder.hh"
#include "SylvException.hh"
struct Rand
{
static std::mt19937 mtgen;
static std::uniform_real_distribution<> dis;
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
};
std::mt19937 Rand::mtgen;
std::uniform_real_distribution<> Rand::dis;
ConstTwoDMatrix
make_matrix(int rows, int cols, const double *p)
{
return ConstTwoDMatrix{rows, cols, ConstVector{p, rows*cols}};
}
void
Rand::init(int n1, int n2, int n3, int n4, int n5)
{
decltype(mtgen)::result_type seed = n1;
seed = 256*seed+n2;
seed = 256*seed+n3;
seed = 256*seed+n4;
seed = 256*seed+n5;
mtgen.seed(seed);
}
double
Rand::get(double m)
{
return 2*m*(dis(mtgen)-0.5);
}
int
Rand::get(int m)
{
return static_cast<int>(get(0.9999*m));
}
bool
Rand::discrete(double prob)
{
return dis(mtgen) < prob;
}
struct SparseGenerator
{
static std::unique_ptr<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);
};
std::unique_ptr<FSSparseTensor>
SparseGenerator::makeTensor(int dim, int nv, int r,
double fill, double m)
{
auto res = std::make_unique<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, static_cast<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 [] = { // 10×10 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 30×20, 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
{
public:
const std::string name;
int dim; // dimension of the solved problem
int nvar; // number of variable of the solved problem
TestRunnable(std::string n, int d, int nv)
: name{std::move(n)}, dim(d), nvar(nv)
{
}
virtual ~TestRunnable() = default;
bool test() const;
virtual bool run() const = 0;
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
{
std::cout << "Running test <" << name << ">" << std::endl;
clock_t start = clock();
auto start_real = std::chrono::steady_clock::now();
bool passed = run();
clock_t end = clock();
auto end_real = std::chrono::steady_clock::now();
std::chrono::duration<double> duration = end_real - start_real;
std::cout << "CPU time " << std::setprecision(4) << std::setw(8)
<< static_cast<double>(end-start)/CLOCKS_PER_SEC << " (CPU seconds)\n"
<< "Real time " << std::setw(8) << duration.count() << " (seconds).....................";
if (passed)
std::cout << "passed\n\n";
else
std::cout << "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++)
std::cout << "\ttensor fill for dim=" << d << " is: "
<< std::setprecision(2) << std::setw(6) << std::fixed
<< c.get(Symmetry{d}).getFillFactor()*100.0 << " %\n"
<< std::defaultfloat;
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<Storage::unfold>(d);
pertime = clock()-pertime;
std::cout << "\ttime for unfolded step dim=" << d << ": " << std::setprecision(4)
<< static_cast<double>(pertime)/CLOCKS_PER_SEC << std::endl;
clock_t checktime = clock();
double err = kord.check<Storage::unfold>(d);
checktime = clock()-checktime;
std::cout << "\ttime for step check dim=" << d << ": " << std::setprecision(4)
<< static_cast<double>(checktime)/CLOCKS_PER_SEC << '\n'
<< "\tmax error in step dim=" << d << ": " << std::setprecision(6) << err
<< std::endl;
maxerror = std::max(err, maxerror);
}
// Perform folded steps until maxdim
if (unfold_dim < maxdim)
{
clock_t swtime = clock();
kord.switchToFolded();
swtime = clock()-swtime;
std::cout << "\ttime for switching dim=" << unfold_dim << ": " << std::setprecision(4)
<< static_cast<double>(swtime)/CLOCKS_PER_SEC << std::endl;
for (int d = unfold_dim+1; d <= maxdim; d++)
{
clock_t pertime = clock();
kord.performStep<Storage::fold>(d);
pertime = clock()-pertime;
std::cout << "\ttime for folded step dim=" << d << ": " << std::setprecision(4)
<< static_cast<double>(pertime)/CLOCKS_PER_SEC << std::endl;
clock_t checktime = clock();
double err = kord.check<Storage::fold>(d);
checktime = clock()-checktime;
std::cout << "\ttime for step check dim=" << d << ": " << std::setprecision(4)
<< static_cast<double>(checktime)/CLOCKS_PER_SEC << '\n'
<< "\tmax error in step dim=" << d << ": " << std::setprecision(6) << err
<< std::endl;
maxerror = std::max(err, maxerror);
}
}
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 < 5e-7;
}
};
// 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.5;
}
};
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.5;
}
};
int
main()
{
std::vector<std::unique_ptr<TestRunnable>> all_tests;
// Fill in vector of all tests
all_tests.push_back(std::make_unique<UnfoldKOrderSmall>());
all_tests.push_back(std::make_unique<UnfoldKOrderSW>());
all_tests.push_back(std::make_unique<UnfoldFoldKOrderSW>());
// Find maximum dimension and maximum nvar
int dmax = 0;
int nvmax = 0;
for (const auto &test : all_tests)
{
if (dmax < test->dim)
dmax = test->dim;
if (nvmax < test->nvar)
nvmax = test->nvar;
}
TLStatic::init(dmax, nvmax); // initialize library
// Launch the tests
int success = 0;
for (const auto &test : all_tests)
{
try
{
if (test->test())
success++;
}
catch (const TLException &e)
{
std::cout << "Caught TL exception in <" << test->name << ">:" << std::endl;
e.print();
}
catch (SylvException &e)
{
std::cout << "Caught Sylv exception in <" << test->name << ">:" << std::endl;
e.printMessage();
}
}
int nfailed = all_tests.size() - success;
std::cout << "There were " << nfailed << " tests that failed out of "
<< all_tests.size() << " tests run." << std::endl;
if (nfailed)
return EXIT_FAILURE;
else
return EXIT_SUCCESS;
}
noinst_LIBRARIES = libparser.a
GENERATED_FILES = assign_tab.cc formula_tab.cc matrix_tab.cc assign_tab.hh formula_tab.hh matrix_tab.hh assign_ll.cc formula_ll.cc matrix_ll.cc
libparser_a_SOURCES = \
location.hh \
atom_assignings.cc \
atom_assignings.hh \
atom_substitutions.cc \
atom_substitutions.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 formula.yy matrix.yy assign.ll formula.ll matrix.ll
%_tab.cc %_tab.hh: %.yy
$(YACC) -W -o$*_tab.cc $<
%_tab.$(OBJEXT): CXXFLAGS += -Wno-old-style-cast
%_ll.cc: %.ll
$(LEX) -i -o$@ $<
%_ll.$(OBJEXT): CXXFLAGS += -Wno-old-style-cast
// -*- C++ -*-
/*
* Copyright © 2006-2011 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
%code requires
{
#include "location.hh"
#define ASGN_LTYPE ogp::location_type
}
%code
{
#include "atom_assignings.hh"
#include <string>
void asgn_error(std::string);
int asgn_lex();
extern ogp::AtomAssignings* aparser;
}
%union
{
int integer;
char *string;
char character;
}
%token EQUAL_SIGN SEMICOLON CHARACTER BLANK
%token <string> NAME;
%define api.prefix {asgn_}
%locations
%defines
%define parse.error verbose
%%
root : assignments | %empty;
assignments : assignments BLANK | assignments assignment | assignment | BLANK;
assignment : NAME EQUAL_SIGN material SEMICOLON {
aparser->add_assignment(@1.off, $1, @1.ll, @3.off-@1.off, @3.ll + @4.ll);}
| NAME space EQUAL_SIGN material SEMICOLON {
aparser->add_assignment(@1.off, $1, @1.ll, @4.off-@1.off, @4.ll + @5.ll);}
;
material : material CHARACTER | material NAME | material BLANK | NAME | CHARACTER | BLANK;
space : space BLANK | BLANK;
%%
void
asgn_error(std::string mes)
{
aparser->error(std::move(mes));
}
/*
* Copyright © 2006 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "atom_assignings.hh"
#include "location.hh"
#include "parser_exception.hh"
#include "utils/cc/exception.hh"
#include <limits>
#include <iostream>
#include <sstream>
#include <iomanip>
using namespace ogp;
AtomAssignings::AtomAssignings(const AtomAssignings &aa, ogp::StaticAtoms &a)
: atoms(a), expr(aa.expr, atoms), left_names(aa.left_names),
lname2expr(aa.lname2expr), order(aa.order)
{
}
/** A global symbol for passing info to the AtomAssignings from
* asgn_parse(). */
AtomAssignings *aparser;
/** The declaration of functions defined in asgn_ll.cc and asgn_tab.cc
* generated from assign.lex assign.y */
void *asgn__scan_string(const char *);
void asgn__destroy_buffer(void *);
void asgn_parse();
extern location_type asgn_lloc;
void
AtomAssignings::parse(const string &stream)
{
asgn_lloc.off = 0;
asgn_lloc.ll = 0;
void *p = asgn__scan_string(stream.c_str());
aparser = this;
asgn_parse();
asgn__destroy_buffer(p);
}
void
AtomAssignings::error(string mes)
{
throw ParserException(std::move(mes), asgn_lloc.off);
}
void
AtomAssignings::add_assignment_to_double(string name, double val)
{
// if left hand side is a registered atom, insert it to tree
int t;
try
{
if (atoms.check(name))
t = expr.add_nulary(name);
else
t = -1;
}
catch (const ParserException &e)
{
t = -1;
}
// register left hand side in order
order.push_back(t);
// add the double to the tree
std::ostringstream buf;
buf << std::setprecision(std::numeric_limits<double>::max_digits10)
<< val;
try
{
expr.parse(buf.str());
}
catch (const ParserException &e)
{
// should never happen
throw ParserException(string("Error parsing double ")+buf.str()+": "+e.message(), 0);
}
// register name of the left hand side and put to lname2expr
left_names.insert(name);
lname2expr.emplace(std::move(name), order.size()-1);
}
void
AtomAssignings::add_assignment(int asgn_off, const string &str, int name_len,
int right_off, int right_len)
{
// the order of doing things here is important: since the
// FormulaParser requires that all references from the i-th tree
// refere to trees with index lass than i, so to capture also a
// nulary term for the left hand side, it must be inserted to the
// expression tree before the expression is parsed.
// find the name in the atoms
string name = str.substr(0, name_len);
// if left hand side is a registered atom, insert it to tree
int t;
try
{
t = atoms.check(name);
if (t == -1)
t = expr.add_nulary(name);
}
catch (const ParserException &e)
{
atoms.register_name(name);
t = expr.add_nulary(name);
}
// register left hand side in order
order.push_back(t);
// parse expression on the right
try
{
expr.parse(str.substr(right_off, right_len));
}
catch (const ParserException &e)
{
throw ParserException(e, asgn_off+right_off);
}
// register name of the left hand side and put to lname2expr
left_names.insert(name);
if (lname2expr.find(name) != lname2expr.end())
{
// Prevent the occurrence of #415
std::cerr << "Changing the value of " << name << " is not supported. Aborting." << std::endl;
exit(EXIT_FAILURE);
}
lname2expr[name] = order.size()-1;
}
void
AtomAssignings::apply_subst(const AtomSubstitutions::Toldnamemap &mm)
{
// go through all old variables and see what are their derived new
// variables
for (const auto & it : mm)
{
const string &oldname = it.first;
const AtomSubstitutions::Tshiftnameset &sset = it.second;
if (!sset.empty())
{
int told = atoms.index(oldname);
if (told < 0 && !atoms.get_name_storage().query(oldname))
atoms.register_name(oldname);
if (told == -1)
told = expr.add_nulary(oldname);
// at least one substitution here, so make an expression
expr.add_formula(told);
// say that this expression is not assigned to any atom
order.push_back(-1);
// now go through all new names derived from the old name and
// reference to the newly added formula
for (const auto & itt : sset)
{
const string &newname = itt.first;
left_names.insert(newname);
lname2expr.emplace(newname, expr.nformulas()-1);
}
}
}
}
void
AtomAssignings::print() const
{
std::cout << "Atom Assignings\nExpressions:\n";
expr.print();
std::cout << "Left names:\n";
for (auto it : lname2expr)
std::cout << it.first << u8" ⇒ " << expr.formula(it.second) << " (t=" << order[it.second] << ")\n";
}
void
AtomAsgnEvaluator::setValues(EvalTree &et) const
{
// set values of constants
aa.atoms.setValues(et);
// set values of variables to NaN or to user set values
double nan = std::numeric_limits<double>::quiet_NaN();
for (int i = 0; i < aa.atoms.nvar(); i++)
{
const string &ss = aa.atoms.name(i);
int t = aa.atoms.index(ss);
if (t >= 0)
{
auto it = user_values.find(t);
if (it == user_values.end())
et.set_nulary(t, nan);
else
et.set_nulary(t, it->second);
}
}
}
void
AtomAsgnEvaluator::set_user_value(const string &name, double val)
{
int t = aa.atoms.index(name);
if (t >= 0)
{
auto it = user_values.find(t);
if (it == user_values.end())
user_values.emplace(t, val);
else
it->second = val;
}
}
void
AtomAsgnEvaluator::load(int i, double res)
{
// set the value
operator[](i) = res;
// if i-th expression is atom, set its value to this EvalTree
int t = aa.order[i];
if (t >= 0)
etree.set_nulary(t, res);
}
double
AtomAsgnEvaluator::get_value(const string &name) const
{
auto it = aa.lname2expr.find(name);
if (it == aa.lname2expr.end())
return std::numeric_limits<double>::quiet_NaN();
else
return operator[](it->second);
}
/*
* Copyright © 2006 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OGP_ATOM_ASSIGNINGS_H
#define OGP_ATOM_ASSIGNINGS_H
#include "static_atoms.hh"
#include "formula_parser.hh"
#include "atom_substitutions.hh"
#include <vector>
#include <map>
namespace ogp
{
class AtomAsgnEvaluator;
/** This class represents atom assignments used in parameters
* settings and initval initialization. It maintains atoms of the
* all expressions on the right hand side, the parsed formulas of
* the right hand sides, and the information about the left hand
* sides. See documentation to the order member below. */
class AtomAssignings
{
friend class AtomAsgnEvaluator;
protected:
using Tvarintmap = std::map<string, int>;
/** All atoms which should be sufficient for formulas at the
* right hand sides. The atoms should be filled with names
* (preregistered). This is a responsibility of the caller. */
StaticAtoms &atoms;
/** The formulas of right hand sides. */
FormulaParser expr;
/** Name storage of the names from left hand sides. */
NameStorage left_names;
/** Information on left hand sides. This maps a name to the
* index of its assigned expression in expr. More than one
* name may reference to the same expression. */
Tvarintmap lname2expr;
/** Information on left hand sides. If order[i] >= 0, then it
* says that i-th expression in expr is assigned to atom with
* order[i] tree index. */
std::vector<int> order;
public:
/** Construct the object using the provided static atoms. */
AtomAssignings(StaticAtoms &a) : atoms(a), expr(atoms)
{
}
/** Make a copy with provided reference to (posibly different)
* static atoms. */
AtomAssignings(const AtomAssignings &aa, StaticAtoms &a);
virtual ~AtomAssignings() = default;
/** Parse the assignments from the given string. */
void parse(const string &stream);
/** Process a syntax error from bison. */
void error(string mes);
/** Add an assignment of the given name to the given
* double. Can be called by a user, anytime. */
void add_assignment_to_double(string name, double val);
/** Add an assignment. Called from assign.y. */
void add_assignment(int asgn_off, const string &str, int name_len,
int right_off, int right_len);
/** This applies old2new map (possibly from atom
* substitutions) to this object. It registers new variables
* in the atoms, and adds the expressions to expr, and left
* names to lname2expr. The information about dynamical part
* of substitutions is ignored, since we are now in the static
* world. */
void apply_subst(const AtomSubstitutions::Toldnamemap &mm);
/** Debug print. */
void print() const;
};
/** This class basically evaluates the atom assignments
* AtomAssignings, so it inherits from ogp::FormulaEvaluator. It
* is also a storage for the results of the evaluation stored as a
* vector, so the class inherits from std::vector<double> and
* ogp::FormulaEvalLoader. As the expressions for atoms are
* evaluated, the results are values for atoms which will be
* used in subsequent evaluations. For this reason, the class
* inherits also from AtomValues. */
class AtomAsgnEvaluator : public FormulaEvalLoader,
public AtomValues,
protected FormulaEvaluator,
public std::vector<double>
{
protected:
using Tusrvalmap = std::map<int, double>;
Tusrvalmap user_values;
const AtomAssignings &aa;
public:
AtomAsgnEvaluator(const AtomAssignings &a)
: FormulaEvaluator(a.expr),
std::vector<double>(a.expr.nformulas()), aa(a)
{
}
~AtomAsgnEvaluator() override = default;
/** This sets all initial values to NaNs, all constants and
* all values set by user by call set_value. This is called by
* FormulaEvaluator::eval() method, which is called by eval()
* method passing this argument as AtomValues. So the
* ogp::EvalTree will be always this->etree. */
void setValues(EvalTree &et) const override;
/** User setting of the values. For example in initval,
* parameters are known and should be set to their values. In
* constrast endogenous variables are set initially to NaNs by
* AtomValues::setValues. */
void set_user_value(const string &name, double val);
/** This sets the result of i-th expression in aa to res, and
* also checks whether the i-th expression is an atom. If so,
* it sets the value of the atom in ogp::EvalTree
* this->etree. */
void load(int i, double res) override;
/** After the user values have been set, the assignments can
* be evaluated. For this purpose we have eval() method. The
* result is that this object as std::vector<double> will
* contain the values. It is ordered given by formulas in
* expr. */
void
eval()
{
FormulaEvaluator::eval(*this, *this);
}
/** This returns a value for a given name. If the name is not
* found among atoms, or there is no assignment for the atom,
* NaN is returned. */
double get_value(const string &name) const;
};
};
#endif
/*
* Copyright © 2006 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "atom_substitutions.hh"
#include "utils/cc/exception.hh"
using namespace ogp;
AtomSubstitutions::AtomSubstitutions(const AtomSubstitutions &as, const FineAtoms &oa,
FineAtoms &na)
: new2old(as.new2old), old2new(as.old2new), old_atoms(oa), new_atoms(na)
{
}
void
AtomSubstitutions::add_substitution(string newname, string oldname, int tshift)
{
// insert to new2old map
new2old.emplace(newname, Tshiftname(oldname, tshift));
// insert to old2new map
auto it = old2new.find(oldname);
if (it != old2new.end())
it->second.emplace(std::move(newname), -tshift);
else
{
Tshiftnameset snset;
snset.emplace(std::move(newname), -tshift);
old2new.emplace(std::move(oldname), snset);
}
// put to info
info.num_substs++;
}
void
AtomSubstitutions::substitutions_finished(VarOrdering::ord_type ot)
{
// create an external ordering of new_atoms from old_atoms
const vector<string> &oa_ext = old_atoms.get_allvar();
vector<string> na_ext;
for (const auto &oname : oa_ext)
{
// add the old name itself
na_ext.push_back(oname);
// add all new names derived from the old name
auto it = old2new.find(oname);
if (it != old2new.end())
for (const auto & itt : it->second)
na_ext.push_back(itt.first);
}
// call parsing finished for the new_atoms
new_atoms.parsing_finished(ot, na_ext);
}
string
AtomSubstitutions::get_new4old(const string &oldname, int tshift) const
{
auto it = old2new.find(oldname);
if (it != old2new.end())
{
const Tshiftnameset &sset = it->second;
for (const auto & itt : sset)
if (itt.second == -tshift)
return itt.first;
}
return "";
}
void
AtomSubstitutions::print() const
{
std::cout << u8"Atom Substitutions:\nOld ⇒ New:\n";
for (const auto & it : old2new)
for (const auto &itt : it.second)
std::cout << " " << it.first << u8" ⇒ [" << itt.first << ", " << itt.second << "]\n";
std::cout << u8"Old ⇐ New:\n";
for (const auto & it : new2old)
std::cout << " [" << it.second.first << ", " << it.second.second << "] ⇐ " << it.first << '\n';
}
void
SAtoms::substituteAllLagsAndLeads(FormulaParser &fp, AtomSubstitutions &as)
{
string name;
int mlead, mlag;
endovarspan(mlead, mlag);
// substitute all endo lagged more than 1
while (!(name = findEndoWithLeadInInterval(mlag, -2)).empty())
makeAuxVariables(name, -1, -2, mlag, fp, as);
// substitute all endo leaded more than 1
while (!(name = findEndoWithLeadInInterval(2, mlead)).empty())
makeAuxVariables(name, 1, 2, mlead, fp, as);
exovarspan(mlead, mlag);
// substitute all lagged exo
while (!(name = findExoWithLeadInInterval(mlag, -1)).empty())
makeAuxVariables(name, -1, -1, mlag, fp, as);
// substitute all leaded exo
while (!(name = findExoWithLeadInInterval(1, mlead)).empty())
makeAuxVariables(name, 1, 1, mlead, fp, as);
// notify that substitution have been finished
as.substitutions_finished(order_type);
}
void
SAtoms::substituteAllLagsAndExo1Leads(FormulaParser &fp, AtomSubstitutions &as)
{
string name;
int mlead, mlag;
endovarspan(mlead, mlag);
// substitute all endo lagged more than 1
while (!(name = findEndoWithLeadInInterval(mlag, -2)).empty())
makeAuxVariables(name, -1, -2, mlag, fp, as);
exovarspan(mlead, mlag);
// substitute all lagged exo
while (!(name = findExoWithLeadInInterval(mlag, -1)).empty())
makeAuxVariables(name, -1, -1, mlag, fp, as);
// substitute all leaded exo by 1
while (!(name = findExoWithLeadInInterval(1, 1)).empty())
makeAuxVariables(name, 1, 1, 1, fp, as);
// notify that substitution have been finished
as.substitutions_finished(order_type);
}
string
SAtoms::findNameWithLeadInInterval(const vector<string> &names,
int ll1, int ll2) const
{
for (auto name : names)
{
auto it = vars.find(name);
if (it != vars.end())
{
const DynamicAtoms::Tlagmap &lmap = it->second;
for (auto itt : lmap)
if (itt.first >= ll1 && itt.first <= ll2)
return name;
}
}
// nothing found
return "";
}
void
SAtoms::attemptAuxName(const string &str, int ll, string &out) const
{
char c = (ll >= 0) ? ((ll == 0) ? 'e' : 'p') : 'm';
string absll = std::to_string(std::abs(ll));
int iter = 1;
do
{
out = str + '_';
for (int i = 0; i < iter; i++)
out += c;
if (ll != 0)
out += absll;
iter++;
}
while (varnames.query(out));
}
void
SAtoms::makeAuxVariables(const string &name, int step, int start, int limit_lead,
FormulaParser &fp, AtomSubstitutions &as)
{
if (!(step == 1 || step == -1))
throw ogu::Exception(__FILE__, __LINE__,
"Wrong value of step in SAtoms::makeAuxVariables");
if (step*start > step*limit_lead)
throw ogu::Exception(__FILE__, __LINE__,
"Wrong value of start in SAtoms::makeAuxVariables");
// make sure that we do not go further than necessary, this is
// that the limit lead is not behind maxlead or minlag
int mlead, mlag;
varspan(name, mlead, mlag);
if (step == -1)
limit_lead = std::max(limit_lead, mlag);
else
limit_lead = std::min(limit_lead, mlead);
// Comment to comments: name="a"; start=-3; step=-1;
// recover tree index of a previous atom, i.e. set tprev to a tree
// index of atom "a(-2)"
int tprev = index(name, start-step);
if (tprev == -1)
tprev = fp.add_nulary(name + '(' + std::to_string(start-step) + ')');
int ll = start;
do
{
// either create atom "a_m2(0)" with tree index taux and add
// equation "a_m2(0)=a(-2)"
// or
// check if "a_m2(0)" has not been already created (with
// different step), in this case do not add equation "a_m2(0)
// = a(-2)"
string newname, newname_str;
int taux;
if ((newname = as.get_new4old(name, ll-step)).empty())
{
attemptAuxName(name, ll-step, newname_str);
newname = newname_str;
register_uniq_endo(newname);
taux = fp.add_nulary(newname + "(0)");
// add to substitutions
as.add_substitution(newname, name, ll-step);
// add equation "a_m2(0) = a(-2)", this is taux = tprev
fp.add_formula(fp.add_binary(code_t::MINUS, taux, tprev));
}
else
{
// example: exogenous EPS and occurrence at both EPS(-1)
// EPS(+1)
// first call makeAuxVariables("EPS",1,1,...) will make endo EPS_p0 = EPS
// second call makeAuxVariables("EPS",-1,-1,...) will use this EPS_p0
// to substitute for EPS(-1)
taux = index(newname, 0);
if (taux < 0)
throw ogu::Exception(__FILE__, __LINE__,
"Couldn't find tree index of previously substituted variable");
}
// create atom "a_m2(-1)" or turn "a(-3)" if any to "a_m2(-1)"; tree index t
int t = index(name, ll);
if (t == -1)
{
// no "a(-3)", make t <-> a_m2(-1)
t = fp.add_nulary(newname + '(' + std::to_string(step) + ')');
}
else
{
// turn a(-3) to a_m2(-1)
unassign_variable(name, ll, t);
assign_variable(newname, step, t);
}
// next iteration starts with tprev <-> "a_m2(-1)" (this will be made equal to "a_m3(0)")
tprev = t;
ll += step;
}
while (step*ll <= step*limit_lead);
}
/*
* Copyright © 2006 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OGP_ATOM_SUBSTITUTIONS_H
#define OGP_ATOM_SUBSTITUTIONS_H
#include "fine_atoms.hh"
#include <string>
namespace ogp
{
using std::string;
using std::map;
using std::pair;
/** This class tracts an information about the performed
* substitutions. In fact, there is only one number to keep track
* about, this is a number of substitutions. */
struct SubstInfo
{
int num_substs{0};
SubstInfo() = default;
};
/** This class tracks all atom substitutions during the job and
* then builds structures when all substitutions are finished. */
class AtomSubstitutions
{
public:
using Tshiftname = pair<string, int>;
using Tshiftmap = map<string, Tshiftname>;
using Tshiftnameset = set<Tshiftname>;
using Toldnamemap = map<string, Tshiftnameset>;
protected:
/** This maps a new name to a shifted old name. This is, one
* entry looks as "a_m3 ==> a(-3)", saying that a variable
* "a_m3" corresponds to a variable "a" lagged by 3. */
Tshiftmap new2old;
/** This is inverse to new2old, which is not unique. For old
* name, say "a", it says what new names are derived with what
* shifts from the "a". For example, it can map "a" to a two
* element set {["a_m3", +3], ["a_p2", -2]}. This says that
* leading "a_m3" by 3 one gets old "a" and lagging "a_p2" by
* 2 one gets also old "a". */
Toldnamemap old2new;
/** This is a reference to old atoms with multiple leads and
* lags. They are supposed to be used with parsing finished
* being had called, so that the external ordering is
* available. */
const FineAtoms &old_atoms;
/** This is a reference to new atoms. All name pointers point
* to storage of these atoms. */
FineAtoms &new_atoms;
/** Substitutions information. */
SubstInfo info;
public:
/** Create the object with reference to the old and new
* atoms. In the beginning, old atoms are supposed to be with
* parsing_finished() called, and new atoms a simple copy of
* old atoms. The new atoms will be an instance of SAtoms. All
* substitution job is done by a substitution method of the
* new atoms. */
AtomSubstitutions(const FineAtoms &oa, FineAtoms &na)
: old_atoms(oa), new_atoms(na)
{
}
/** Construct a copy of the object using a different instances
* of old atoms and new atoms, which are supposed to be
* semantically same as the atoms from as. */
AtomSubstitutions(const AtomSubstitutions &as, const FineAtoms &oa, FineAtoms &na);
virtual ~AtomSubstitutions() = default;
/** This is called during the substitution job from the
* substitution method of the new atoms. This says that the
* new name, say "a_m3" is a substitution of old name "a"
* shifted by -3. */
void add_substitution(string newname, string oldname, int tshift);
/** This is called when all substitutions are finished. This
* forms the new external ordering of the new atoms and calls
* parsing_finished() for the new atoms with the given ordering type. */
void substitutions_finished(VarOrdering::ord_type ot);
/** Returns a new name for old name and given tshift. For "a"
* and tshift=-3, it returns "a_m3". If there is no such
* substitution, it returns an empty string. */
string get_new4old(const string &oldname, int tshift) const;
/** Return new2old. */
const Tshiftmap &
get_new2old() const
{
return new2old;
}
/** Return old2new. */
const Toldnamemap &
get_old2new() const
{
return old2new;
}
/** Return substitution info. */
const SubstInfo &
get_info() const
{
return info;
}
/** Return old atoms. */
const FineAtoms &
get_old_atoms() const
{
return old_atoms;
}
/** Return new atoms. */
const FineAtoms &
get_new_atoms() const
{
return new_atoms;
}
/** Debug print. */
void print() const;
};
class SAtoms : public FineAtoms
{
public:
SAtoms()
: FineAtoms()
{
}
SAtoms(const SAtoms &sa) = default;
/** This substitutes all lags and leads for all exogenous and
* all lags and leads greater than 1 for all endogenous
* variables. This is useful for perfect foresight problems
* where we can do that. */
void substituteAllLagsAndLeads(FormulaParser &fp, AtomSubstitutions &as);
/** This substitutes all lags of all endo and exo and one step
* leads of all exo variables. This is useful for stochastic
* models where we cannot solve leads more than 1. */
void substituteAllLagsAndExo1Leads(FormulaParser &fp, AtomSubstitutions &as);
protected:
/** This finds an endogenous variable name which occurs between
* ll1 and ll2 included. */
string
findEndoWithLeadInInterval(int ll1, int ll2) const
{
return findNameWithLeadInInterval(get_endovars(), ll1, ll2);
}
/** This finds an exogenous variable name which occurs between
* ll1 and ll2 included. */
string
findExoWithLeadInInterval(int ll1, int ll2) const
{
return findNameWithLeadInInterval(get_exovars(), ll1, ll2);
}
/** This attempts to find a non registered name of the form
* <str>_m<abs(ll)> or <str>_p<abs(ll)>. A letter 'p' is
* chosen if ll is positive, 'm' if negative. If a name of
* such form is already registered, one more character (either
* 'p' or 'm') is added and the test is performed again. The
* resulting name is returned in a string out. */
void attemptAuxName(const string &str, int ll, string &out) const;
/** This makes auxiliary variables to eliminate all leads/lags
* greater/less than or equal to start up to the limit_lead
* for a variable with the given name. If the limit_lead is
* greater/less than the maxlead/minlag of the variable, than
* maxlead/minlag is used. This process is recorded in
* AtomSubstitutions. The new auxiliary variables and their
* atoms are created in this object. The auxiliary equations
* are created in the given FormulaParser. The value of step
* is allowed to be either -1 (lags) or +1 (leads). */
void makeAuxVariables(const string &name, int step, int start, int limit_lead,
FormulaParser &fp, AtomSubstitutions &as);
private:
/** This is a worker routine for findEndoWithLeadInInterval
* and findExoWithLeadInInterval. */
string findNameWithLeadInInterval(const vector<string> &names,
int ll1, int ll2) const;
};
};
#endif
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "utils/cc/exception.hh"
#include "dynamic_atoms.hh"
using namespace ogp;
void
NameStorage::insert(string name)
{
if (!query(name))
{
name_store.push_back(name);
name_set.insert(std::move(name));
}
}
void
NameStorage::print() const
{
for (auto i : name_store)
std::cout << i << '\n';
}
void
Constants::import_constants(const Constants &c, OperationTree &otree, Tintintmap &tmap)
{
for (auto it : c.cmap)
{
int told = it.first;
int tnew = otree.add_nulary();
tmap.emplace(told, tnew);
add_constant(tnew, it.second);
}
}
void
Constants::setValues(EvalTree &et) const
{
for (const auto & it : cmap)
et.set_nulary(it.first, it.second);
}
void
Constants::add_constant(int t, double val)
{
cmap.emplace(t, val);
cinvmap.emplace(val, t);
}
bool
Constants::is_constant(int t) const
{
if (t < OperationTree::num_constants)
return true;
auto it = cmap.find(t);
return (it != cmap.end());
}
double
Constants::get_constant_value(int t) const
{
auto it = cmap.find(t);
if (it != cmap.end())
return it->second;
else
throw ogu::Exception(__FILE__, __LINE__,
"Tree index is not constant in Constants::get_constant_value");
}
int
Constants::check(const string &str) const
{
double d = std::stod(str);
auto it = cinvmap.find(d);
if (it != cinvmap.end())
return it->second;
else
return -1;
}
void
Constants::print() const
{
for (const auto &it : cmap)
std::cout << "$" << it.first << ": " << it.second << "\n";
}
int
DynamicAtoms::check(const string &name) const
{
if (is_string_constant(name))
return Constants::check(name);
return check_variable(name);
}
int
DynamicAtoms::check_variable(const string &name) const
{
string str;
int ll;
parse_variable(name, str, ll);
auto it = vars.find(str);
if (it != vars.end())
{
const Tlagmap &lmap = it->second;
auto itt = lmap.find(ll);
if (itt != lmap.end())
return itt->second;
}
return -1;
}
void
DynamicAtoms::assign(const string &name, int t)
{
if (is_string_constant(name))
assign_constant(name, t);
else
assign_variable(name, t);
}
void
DynamicAtoms::assign_constant(const string &name, int t)
{
double val = std::stod(name);
add_constant(t, val);
}
// parse the name and then call assing_variable(varname, ll, t)
void
DynamicAtoms::assign_variable(const string &name, int t)
{
int ll;
string str;
parse_variable(name, str, ll);
// here str is just name without lead/lag
varnames.insert(str);
assign_variable(str, ll, t);
}
void
DynamicAtoms::assign_variable(const string &varname, int ll, int t)
{
if (indices.end() != indices.find(t))
throw ogu::Exception(__FILE__, __LINE__,
"Attempt to assign already allocated tree index");
auto it = vars.find(varname);
if (it != vars.end())
{
Tlagmap &lmap = it->second;
if (lmap.end() != lmap.find(ll))
throw ogu::Exception(__FILE__, __LINE__,
"Attempt to assign already allocated variable");
lmap.emplace(ll, t);
}
else
{
Tlagmap lmap;
lmap.emplace(ll, t);
vars.emplace(varname, lmap);
}
indices.emplace(t, varname);
nv++;
minlag = std::min(ll, minlag);
maxlead = std::max(ll, maxlead);
}
void
DynamicAtoms::unassign_variable(const string &varname, int ll, int t)
{
auto it = vars.find(varname);
if (it != vars.end())
{
Tlagmap &lmap = it->second;
auto itt = lmap.find(ll);
if (itt != lmap.end())
{
if (itt->second == t)
{
// erase it from the lagmap; if it becomes empty,
// erase the lagmap from varmap
lmap.erase(itt);
if (lmap.size() == 0)
vars.erase(it);
// erase it from the indices
auto ittt = indices.find(t);
if (ittt != indices.end())
indices.erase(ittt);
nv--;
if (ll == minlag || ll == maxlead)
update_minmaxll();
}
else
throw ogu::Exception(__FILE__, __LINE__,
"Tree index inconsistent in DynamicAtoms::unassign_variable");
}
else
throw ogu::Exception(__FILE__, __LINE__,
"Lead/lag of the variable not found in DynamicAtoms::unassign_variable");
}
else
throw ogu::Exception(__FILE__, __LINE__,
"Variable not found in DynamicAtoms::unassign_variable");
}
void
DynamicAtoms::update_minmaxll()
{
minlag = std::numeric_limits<int>::max();
maxlead = std::numeric_limits<int>::min();
for (const auto &it : vars)
{
const Tlagmap &lmap = it.second;
for (auto itt : lmap)
{
int ll = itt.first;
minlag = std::min(ll, minlag);
maxlead = std::max(ll, maxlead);
}
}
}
vector<int>
DynamicAtoms::variables() const
{
vector<int> res;
for (const auto & var : vars)
{
const Tlagmap &lmap = var.second;
for (auto itt : lmap)
res.push_back(itt.second);
}
return res;
}
void
DynamicAtoms::varspan(int t, int &mlead, int &mlag) const
{
auto it = indices.find(t);
if (indices.end() == it)
{
mlead = std::numeric_limits<int>::min();
mlag = std::numeric_limits<int>::max();
return;
}
varspan(it->second, mlead, mlag);
}
void
DynamicAtoms::varspan(const string &name, int &mlead, int &mlag) const
{
auto it = vars.find(name);
if (vars.end() == it)
{
mlead = std::numeric_limits<int>::min();
mlag = std::numeric_limits<int>::max();
return;
}
const Tlagmap &lmap = it->second;
auto beg = lmap.begin();
auto end = lmap.rbegin();
mlag = beg->first;
mlead = end->first;
}
void
DynamicAtoms::varspan(const vector<string> &names, int &mlead, int &mlag) const
{
mlead = std::numeric_limits<int>::min();
mlag = std::numeric_limits<int>::max();
for (const auto &name : names)
{
int lag, lead;
varspan(name, lead, lag);
mlead = std::max(lead, mlead);
mlag = std::min(lag, mlag);
}
}
bool
DynamicAtoms::is_named_atom(int t) const
{
return indices.end() != indices.find(t);
}
int
DynamicAtoms::index(const string &name, int ll) const
{
auto it = vars.find(name);
if (vars.end() != it)
{
const Tlagmap &lmap = it->second;
auto itt = lmap.find(ll);
if (lmap.end() != itt)
return itt->second;
}
return -1;
}
bool
DynamicAtoms::is_referenced(const string &name) const
{
return vars.find(name) != vars.end();
}
const DynamicAtoms::Tlagmap &
DynamicAtoms::lagmap(const string &name) const
{
auto it = vars.find(name);
if (vars.end() == it)
throw ogu::Exception(__FILE__, __LINE__,
std::string("Couldn't find the name ")
+ name + " in DynamicAtoms::lagmap");
return it->second;
}
const string &
DynamicAtoms::name(int t) const
{
auto it = indices.find(t);
if (indices.end() == it)
throw ogu::Exception(__FILE__, __LINE__,
"Couldn't find tree index in DynamicAtoms::name");
return it->second;
}
int
DynamicAtoms::lead(int t) const
{
const string &nam = name(t);
const Tlagmap &lmap = lagmap(nam);
auto it = lmap.begin();
while (it != lmap.end() && it->second != t)
++it;
if (lmap.end() == it)
throw ogu::Exception(__FILE__, __LINE__,
"Couldn't find the three index in DynamicAtoms::lead");
return it->first;
}
void
DynamicAtoms::print() const
{
std::cout << "names:\n";
varnames.print();
std::cout << "constants:\n";
Constants::print();
std::cout << "variables:\n";
for (const auto & var : vars)
{
const Tlagmap &lmap = var.second;
for (auto itt : lmap)
std::cout << "$" << itt.second << ": " << var.first << "(" << itt.first << ")\n";
}
std::cout << "indices:\n";
for (auto indice : indices)
std::cout << "t=" << indice.first << u8" ⇒ " << indice.second << "\n";
}
/** Note that the str has been parsed by the lexicographic
* analyzer. It can be either a variable or a double. So it is easy to
* recognize it by the first character. */
bool
DynamicAtoms::is_string_constant(const string &str)
{
return str[0] == '.' || str[0] == '-' || (str[0] >= '0' && str[0] <= '9');
}
VarOrdering::VarOrdering(const VarOrdering &vo, const vector<string> &vnames,
const DynamicAtoms &a)
: n_stat(vo.n_stat), n_pred(vo.n_pred), n_both(vo.n_both), n_forw(vo.n_forw),
der_atoms(vo.der_atoms), positions(vo.positions),
outer2y(vo.outer2y), y2outer(vo.y2outer), varnames(vnames), atoms(a)
{
}
bool
VarOrdering::check(int t) const
{
return positions.find(t) != positions.end();
}
int
VarOrdering::get_pos_of(int t) const
{
auto it = positions.find(t);
if (it != positions.end())
return it->second;
else
throw ogu::Exception(__FILE__, __LINE__,
"Couldn't find the tree index in VarOrdering::get_pos_of");
}
void
VarOrdering::do_general(ord_type ordering)
{
// auxiliary vectors for setting der_atoms and map
vector<int> pred_minus;
vector<int> both_minus;
vector<int> stat;
vector<int> pred_pad;
vector<int> both_pad;
vector<int> forw_pad;
vector<int> both_plus;
vector<int> forw_plus;
// auxiliary vectors for setting y2outer and outer2y
vector<int> y2o_stat;
vector<int> y2o_pred;
vector<int> y2o_both;
vector<int> y2o_forw;
for (unsigned int i = 0; i < varnames.size(); i++)
{
const string &ss = varnames[i];
int lead;
int lag;
atoms.varspan(ss, lead, lag);
if (lag == 0 && lead == 0)
{
stat.push_back(atoms.index(ss, 0));
y2o_stat.push_back(i);
}
else if (lag == -1 && lead < 1)
{
pred_minus.push_back(atoms.index(ss, -1));
pred_pad.push_back(atoms.index(ss, 0));
y2o_pred.push_back(i);
}
else if (lag > -1 && lead == 1)
{
forw_pad.push_back(atoms.index(ss, 0));
forw_plus.push_back(atoms.index(ss, 1));
y2o_forw.push_back(i);
}
else if (lag == -1 && lead == 1)
{
both_minus.push_back(atoms.index(ss, -1));
both_pad.push_back(atoms.index(ss, 0));
both_plus.push_back(atoms.index(ss, 1));
y2o_both.push_back(i);
}
else
throw ogu::Exception(__FILE__, __LINE__,
"A wrong lag/lead of a variable in VarOrdering::do_pbspbfbf");
}
// here we fill ords according to ordering
vector<int> *ords[8];
if (ordering == pbspbfbf)
{
ords[0] = &pred_minus;
ords[1] = &both_minus;
ords[2] = &stat;
ords[3] = &pred_pad;
ords[4] = &both_pad;
ords[5] = &forw_pad;
ords[6] = &both_plus;
ords[7] = &forw_plus;
}
else if (ordering == bfspbfpb)
{
ords[0] = &both_plus;
ords[1] = &forw_plus;
ords[2] = &stat;
ords[3] = &pred_pad;
ords[4] = &both_pad;
ords[5] = &forw_pad;
ords[6] = &pred_minus;
ords[7] = &both_minus;
}
else // BEWARE: when implementing a new ordering, check also the code below setting y2outer
throw ogu::Exception(__FILE__, __LINE__,
"Ordering not implemented in VarOrdering::do_general");
// make der_atoms and positions
int off = 0;
for (auto & ord : ords)
for (unsigned int j = 0; j < ord->size(); j++, off++)
if ((*ord)[j] != -1)
{
der_atoms.push_back((*ord)[j]);
positions.emplace((*ord)[j], off);
}
// set integer constants
n_stat = stat.size();
n_pred = pred_pad.size();
n_both = both_pad.size();
n_forw = forw_pad.size();
// make y2outer mapping
y2outer.insert(y2outer.end(), y2o_stat.begin(), y2o_stat.end());
y2outer.insert(y2outer.end(), y2o_pred.begin(), y2o_pred.end());
y2outer.insert(y2outer.end(), y2o_both.begin(), y2o_both.end());
y2outer.insert(y2outer.end(), y2o_forw.begin(), y2o_forw.end());
// make outer2y mapping
outer2y.resize(y2outer.size(), -1);
for (unsigned int i = 0; i < y2outer.size(); i++)
outer2y[y2outer[i]] = i;
}
void
VarOrdering::do_increasing_time()
{
// get maxlead and minlag of the variables
int mlag, mlead;
atoms.varspan(varnames, mlead, mlag);
// setup the matrix of tree indices, if there is no occurrence,
// the index is set to -1
vector<int> ll_init(varnames.size(), -1);
vector<vector<int>> tree_ind(mlead-mlag+1, ll_init);
for (unsigned int iv = 0; iv < varnames.size(); iv++)
{
try
{
const DynamicAtoms::Tlagmap &lmap = atoms.lagmap(varnames[iv]);
for (auto it : lmap)
{
int ll = it.first;
int t = it.second;
tree_ind[ll-mlag][iv] = t;
}
}
catch (const ogu::Exception &e)
{
// ignore the error of not found variable in the tree
}
}
// setup der_atoms and positions
for (int ll = mlag; ll <= mlead; ll++)
for (unsigned int iv = 0; iv < varnames.size(); iv++)
{
int t = tree_ind[ll-mlag][iv];
if (t != -1)
{
der_atoms.push_back(t);
int pos = (ll-mlag)*varnames.size() + iv;
positions.emplace(t, pos);
}
}
// set outer2y and y2outer to identities
for (unsigned int iv = 0; iv < varnames.size(); iv++)
{
outer2y.push_back(iv);
y2outer.push_back(iv);
}
// set n_stat, n_pred, n_both, and n_forw
for (auto varname : varnames)
{
int mmlag, mmlead;
atoms.varspan(varname, mmlead, mmlag);
if (mmlead == 0 && mmlag == 0)
n_stat++;
else if (mmlead <= 0 && mmlag < 0)
n_pred++;
else if (mmlead > 0 && mmlag >= 0)
n_forw++;
else if (mmlead > 0 && mmlag < 0)
n_both++;
else if (mmlead < mmlag)
// variable does not occur in the tree, cound as static
n_stat++;
else
throw ogu::Exception(__FILE__, __LINE__,
"A wrong lag/lead of a variable in VarOrdering::do_increasing_time");
}
}
void
VarOrdering::print() const
{
std::cout << "nstat=" << n_stat << ", npred=" << n_pred << ", nboth=" << n_both
<< ", nforw=" << n_forw << "\n"
<< "der_atoms:\n";
for (int der_atom : der_atoms)
std::cout << " " << der_atom;
std::cout << "\nmap:\n";
for (auto position : positions)
std::cout << " [" << position.first << u8"→" << position.second << "]";
std::cout << "\ny2outer:\n";
for (int i : y2outer)
std::cout << " " << i;
std::cout << "\nouter2y:\n";
for (int i : outer2y)
std::cout << " " << i;
std::cout << "\n";
}
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OGP_DYNAMIC_ATOMS_H
#define OGP_DYNAMIC_ATOMS_H
#include "formula_parser.hh"
#include <vector>
#include <map>
#include <set>
#include <string>
#include <limits>
#include <memory>
namespace ogp
{
using std::vector;
using std::map;
using std::set;
using std::string;
/** Class storing names. We will keep names of variables in
* various places, and all these pointers will point to one
* storage, which will be responsible for allocation and
* deallocation. The main function of the class is to allocate
* space for names, and return a pointer of the stored name if
* required. */
class NameStorage
{
protected:
/** Vector of names allocated, this is the storage. */
vector<string> name_store;
/** Map useful to quickly decide if the name is already
* allocated or not. */
set<string> name_set;
public:
/** Query for the name. If the name has been stored, it
* true, otherwise false. */
bool
query(const string &name) const
{
return name_set.find(name) != name_set.end();
}
/** Insert the name if it has not been inserted yet. */
void insert(string name);
int
num() const
{
return static_cast<int>(name_store.size());
}
const string &
get_name(int i) const
{
return name_store[i];
}
/** Debug print. */
void print() const;
};
class Constants : public AtomValues
{
public:
/** Type for a map mapping tree indices to double values. */
using Tconstantmap = map<int, double>;
using Tintintmap = map<int, int>;
protected:
/** Map mapping a tree index of a constant to its double value. */
Tconstantmap cmap;
public:
Constants() = default;
/** Copy constructor. */
Constants(const Constants &c)
: cmap(c.cmap), cinvmap(c.cinvmap)
{
}
/** Copy constructor registering the constants in the given
* tree. The mapping from old tree indices to new ones is
* traced in tmap. */
Constants(const Constants &c, OperationTree &otree, Tintintmap &tmap)
{
import_constants(c, otree, tmap);
}
/** Import constants registering their tree indices in the
* given tree. The mapping form old tree indices to new ones
* is traced in tmap. */
void import_constants(const Constants &c, OperationTree &otree, Tintintmap &tmap);
/** Implements AtomValues interface. This sets the values to
* the evaluation tree EvalTree. */
void setValues(EvalTree &et) const override;
/** This adds a constant with the given tree index. The
* constant must be checked previously and asserted that it
* does not exist. */
void add_constant(int t, double val);
/** Returns true if the tree index is either an hardwired
* constant (initial number OperationTree:num_constants in
* OperationTree) or the tree index is a registered constant
* by add_constant method. */
bool is_constant(int t) const;
double get_constant_value(int t) const;
/** Return -1 if the given string representation of a constant
* is not among the constants (double represenations). If it
* is, its tree index is returned. */
int check(const string &str) const;
/** Debug print. */
void print() const;
const Tconstantmap &
get_constantmap() const
{
return cmap;
}
private:
/** Inverse map to Tconstantmap. */
using Tconstantinvmap = map<double, int>;
/** This is an inverse map to cmap. This is only used for fast
* queries for the existing double constants in check
* method and add_constant. */
Tconstantinvmap cinvmap;
};
/** This class is a parent to Atoms classes which distinguish between
* constants (numerical literals), and variables with lags and
* leads. This abstraction does not distinguish between a parameter
* and a variable without lag or lead. In this sense, everything is a
* variable.*/
class DynamicAtoms : public Atoms, public Constants
{
public:
/** Definition of a type mapping lags to the indices of the variables. */
using Tlagmap = map<int, int>;
protected:
/** Definition of a type mapping names of the atoms to Tlagmap. */
using Tvarmap = map<string, Tlagmap>;
/** Definition of a type mapping indices of variables to the variable names. */
using Tindexmap = map<int, string>;
/** This is just a storage for variable names, since all other
* instances of a variable name just point to the memory
* allocated by this object. */
NameStorage varnames;
/** This is the map for variables. Each variable name is
* mapped to the Tlagmap, which maps lags/leads to the nulary
* term indices in the tree. */
Tvarmap vars;
/** This is almost inverse map to the vars. It maps variable
* indices to the names. A returned name can be in turn used
* as a key in vars. */
Tindexmap indices;
/** Number of variables. */
int nv{0};
/** Minimum lag, if there is at least one lag, than this is a negative number. */
int minlag{std::numeric_limits<int>::max()};
/** Maximum lead, if there is at least one lead, than this is a positive number. */
int maxlead{std::numeric_limits<int>::min()};
public:
/** Construct empty DynamicAtoms. */
DynamicAtoms() = default;
/** Check the nulary term identified by its string
* representation. The nulary term can be either a constant or
* a variable. If constant, -1 is returned so that it could be
* assigned regardless if the same constant has already
* appeared or not. If variable, then -1 is returned only if
* the variable has not been assigned an index, otherwise the
* assigned index is returned. */
int check(const string &name) const override;
/** Assign the nulary term identified by its string
* representation. This method should be called when check()
* returns -1. */
void assign(const string &name, int t) override;
/** Return a number of all variables. */
int
nvar() const override
{
return nv;
}
/** Return the vector of variable indices. */
vector<int> variables() const override;
/** Return max lead and min lag for a variable given by the
* index. If a variable cannot be found, the method retursn
* the smallest integer as maxlead and the largest integer as
* minlag. */
void varspan(int t, int &mlead, int &mlag) const;
/** Return max lead and min lag for a variable given by the
* name (without lead, lag). The same is valid if the variable
* name cannot be found. */
void varspan(const string &name, int &mlead, int &mlag) const;
/** Return max lead and min lag for a vector of variables given by the names. */
void varspan(const vector<string> &names, int &mlead, int &mlag) const;
/** Return true for all tree indices corresponding to a
* variable in the sense of this class. (This is parameters,
* exo and endo). Since the semantics of 'variable' will be
* changed in subclasses, we use name 'named atom'. These are
* all atoms but constants. */
bool is_named_atom(int t) const;
/** Return index of the variable described by the variable
* name and lag/lead. If it doesn't exist, return -1. */
int index(const string &name, int ll) const;
/** Return true if a variable is referenced, i.e. it has lag
* map. */
bool is_referenced(const string &name) const;
/** Return the lag map for the variable name. */
const Tlagmap &lagmap(const string &name) const;
/** Return the variable name for the tree index. It throws an
* exception if the tree index t is not a named atom. */
const string &name(int t) const;
/** Return the lead/lag for the tree index. It throws an
* exception if the tree index t is not a named atom. */
int lead(int t) const;
/** Return maximum lead. */
int
get_maxlead() const
{
return maxlead;
}
/** Return minimum lag. */
int
get_minlag() const
{
return minlag;
}
/** Return the name storage to allow querying to other
* classes. */
const NameStorage &
get_name_storage() const
{
return varnames;
}
/** Assign the variable with a given lead. The varname must be
* from the varnames storage. The method checks if the
* variable iwht the given lead/lag is not assigned. If so, an
* exception is thrown. */
void assign_variable(const string &varname, int ll, int t);
/** Unassign the variable with a given lead and given tree
* index. The tree index is only provided as a check. An
* exception is thrown if the name, ll, and the tree index t
* are not consistent. The method also updates nv, indices,
* maxlead and minlag. The varname must be from the varnames
* storage. */
void unassign_variable(const string &varname, int ll, int t);
/** Debug print. */
void print() const override;
protected:
/** Do the check for the variable. A subclass may need to
* reimplement this so that it could raise an error if the
* variable is not among a given list. */
virtual int check_variable(const string &name) const;
/** Assign the constant. */
void assign_constant(const string &name, int t);
/** Assign the variable. */
void assign_variable(const string &name, int t);
/** The method just updates minlag or/and maxlead. Note that
* when assigning variables, the update is done when inserting
* to the maps, however, if removing a variable, we need to
* call this method. */
void update_minmaxll();
/** The method parses the string to recover a variable name
* and lag/lead ll. The variable name doesn't contain a lead/lag. */
virtual void parse_variable(const string &in, string &out, int &ll) const = 0;
public:
/** Return true if the str represents a double.*/
static bool is_string_constant(const string &str);
};
/** This class is a parent of all orderings of the dynamic atoms
* of variables which can appear before t, at t, or after t. It
* encapsulates the ordering, and the information about the number
* of static (appearing only at time t) predetermined (appearing
* before t and possibly at t), both (appearing before t and after
* t and possibly at t) and forward looking (appearing after t and
* possibly at t).
*
* The constructor takes a list of variable names. The class also
* provides mapping from the ordering of the variables in the list
* (outer) to the new ordering (at time t) and back.
*
* The user of the subclass must call do_ordering() after
* initialization.
*
* The class contains a few preimplemented methods for
* ordering. The class is used in this way: Make a subclass, and
* implement pure virtual do_ordering() by just plugging a
* preimplemented method, or plugging your own implementation. The
* method do_ordering() is called by the user after the constructor.
*/
class VarOrdering
{
protected:
/** Number of static variables. */
int n_stat;
/** Number of predetermined variables. */
int n_pred;
/** Number of both variables. */
int n_both;
/** Number of forward looking variables. */
int n_forw;
/** This is a set of tree indices corresponding to the
* variables at all times as they occur in the formulas. In
* fact, since this is used only for derivatives, the ordering
* of this vector is only important for ordering of the
* derivatives, in other contexts the ordering is not
* important, so it is rather a set of indices.*/
vector<int> der_atoms;
/** This maps tree index of the variable to the position in
* the row of the ordering. One should be careful with making
* space in the positions for variables not appearing at time
* t. For instance in the pred(t-1), both(t-1), stat(t),
* pred(t), both(t), forw(t), both(t+1), forw(t+1) ordering,
* the variables x(t-1), y(t-1), x(t+1), z(t-1), z(t), and
* z(t+1) having tree indices 6,5,4,3,2,1 will be ordered as
* follows: y(t-1), x(t-1), z(t-1), [y(t)], [x(t)], z(t),
* x(t+1), where a bracketed expresion means non-existent by
* occupying a space. The map thus will look as follows:
* {5→0, 6→1, 3→2, 2→5, 3→6}. Note that nothing is mapped
* to positions 3 and 4. */
map<int, int> positions;
/** This maps an ordering of the list of variables in
* constructor to the new ordering (at time t). The length is
* the number of variables. */
vector<int> outer2y;
/** This maps a new ordering to the ordering of the list of
* variables in constructor (at time t). The length is the
* number of variables. */
vector<int> y2outer;
/** This is just a reference for variable names to keep it
* from constructor to do_ordering() implementations. */
const vector<string> &varnames;
/** This is just a reference to atoms to keep it from
* constructor to do_ordering() implementations. */
const DynamicAtoms &atoms;
public:
/** This is an enum type for an ordering type implemented by
* do_general. */
enum ord_type {pbspbfbf, bfspbfpb};
/** Construct the ordering of the variables given by the names
* with their dynamic occurrences defined by the atoms. It
* calls the virtual method do_ordering which can be
* reimplemented. */
VarOrdering(const vector<string> &vnames, const DynamicAtoms &a)
: n_stat(0), n_pred(0), n_both(0), n_forw(0), varnames(vnames), atoms(a)
{
}
VarOrdering(const VarOrdering &vo, const vector<string> &vnames,
const DynamicAtoms &a);
VarOrdering(const VarOrdering &vo) = delete;
virtual std::unique_ptr<VarOrdering> clone(const vector<string> &vnames,
const DynamicAtoms &a) const = 0;
/** Destructor does nothing here. */
virtual ~VarOrdering() = default;
/** This is the method setting the ordering and the map. A
* subclass must reimplement it, possibly using a
* preimplemented ordering. This method must be called by the
* user after the class has been created. */
virtual void do_ordering() = 0;
/** Return number of static. */
int
nstat() const
{
return n_stat;
}
/** Return number of predetermined. */
int
npred() const
{
return n_pred;
}
/** Return number of both. */
int
nboth() const
{
return n_both;
}
/** Return number of forward looking. */
int
nforw() const
{
return n_forw;
}
/** Return the set of tree indices for derivatives. */
const vector<int> &
get_der_atoms() const
{
return der_atoms;
}
/** Return the y2outer. */
const vector<int> &
get_y2outer() const
{
return y2outer;
}
/** Return the outer2y. */
const vector<int> &
get_outer2y() const
{
return outer2y;
}
/** Query the atom given by the tree index. True is returned
* if the atom is one of the variables in the object. */
bool check(int t) const;
/** Return the position of the atom (nulary term) given by a
* tree index. It is a lookup to the map. If the atom cannot
* be found, the exception is raised. */
int get_pos_of(int t) const;
/** This returns a length of ordered row of atoms. In all
* cases so far, it does not depend on the ordering and it is
* as follows. */
int
length() const
{
return n_stat+2*n_pred+3*n_both+2*n_forw;
}
/** Debug print. */
void print() const;
protected:
/** This is a general ordering method which orders the
* variables by the given ordering ord_type. See documentation
* for respective do_ methods. */
void do_general(ord_type ordering);
/** This is a preimplemented ordering for do_ordering()
* method. It assumes that the variables appear only at time
* t-1, t, t+1. It orders the atoms as pred(t-1), both(t-1),
* stat(t), pred(t), both(t), forw(t), both(t+1),
* forw(t+1). It builds the der_atoms, the map of positions,
* as well as y2outer and outer2y. */
void
do_pbspbfbf()
{
do_general(pbspbfbf);
}
/** This is a preimplemented ordering for do_ordering()
* method. It assumes that the variables appear only at time
* t-1, t, t+1. It orders the atoms as both(t+1), forw(t+1),
* stat(t), pred(t), both(t), forw(t), pred(t-1),
* both(t-1). It builds the der_atoms, the map of positions,
* as well as y2outer and outer2y. */
void
do_bfspbfpb()
{
do_general(bfspbfpb);
}
/** This is a preimplemented ordering for do_ordering()
* method. It makes no assumptions about occurences of
* variables at different times. It orders the atoms with
* increasing time keeping the given ordering within one
* time. This implies that y2outer and outer2y will be
* identities. The der_atoms will be just a sequence of atoms
* from the least to the most time preserving the order of atoms
* within one time. */
void do_increasing_time();
};
};
#endif
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "utils/cc/exception.hh"
#include "parser_exception.hh"
#include "fine_atoms.hh"
using namespace ogp;
AllvarOuterOrdering::AllvarOuterOrdering(const vector<string> &allvar_outer,
const FineAtoms &a)
: atoms(a), allvar(),
endo2all(a.get_endovars().size(), -1),
exo2all(a.get_exovars().size(), -1)
{
// fill in the allvar from allvar_outer
for (const auto &s : allvar_outer)
{
if (atoms.varnames.query(s))
allvar.push_back(s);
else
throw ogu::Exception(__FILE__, __LINE__,
string("Variable ") + s + " is not a declared symbol in AllvarOuterOrdering constructor");
}
// fill in endo2all and exo2all
for (unsigned int i = 0; i < allvar.size(); i++)
{
auto it = atoms.endo_outer_map.find(allvar[i]);
if (it != atoms.endo_outer_map.end())
endo2all[it->second] = i;
else
{
it = atoms.exo_outer_map.find(allvar[i]);
if (it != atoms.exo_outer_map.end())
exo2all[it->second] = i;
else
throw ogu::Exception(__FILE__, __LINE__,
string("Name ") + allvar[i] + " is neither endogenous nor exogenous variable in AllvarOuterOrdering constructor");
}
}
// check whether everything has been filled
unsigned int iendo = 0;
while (iendo < endo2all.size() && endo2all[iendo] != -1)
iendo++;
unsigned int iexo = 0;
while (iexo < exo2all.size() && exo2all[iexo] != -1)
iexo++;
if (iendo < endo2all.size())
throw ogu::Exception(__FILE__, __LINE__,
string("Endogenous variable ") + atoms.get_endovars()[iendo]
+" not found in outer all ordering in AllvarOuterOrdering constructor");
if (iexo < exo2all.size())
throw ogu::Exception(__FILE__, __LINE__,
string("Exogenous variable ") + atoms.get_exovars()[iexo]
+" not found in outer all ordering in AllvarOuterOrdering constructor");
}
AllvarOuterOrdering::AllvarOuterOrdering(const AllvarOuterOrdering &avo,
const FineAtoms &a)
: atoms(a), allvar(avo.allvar),
endo2all(avo.endo2all),
exo2all(avo.exo2all)
{
}
FineAtoms::FineAtoms(const FineAtoms &fa)
: DynamicAtoms(fa), params(), endovars(), exovars(),
der_atoms(fa.der_atoms),
endo_atoms_map(fa.endo_atoms_map),
exo_atoms_map(fa.exo_atoms_map)
{
// fill in params
for (const auto &param : fa.params)
{
if (!varnames.query(param))
throw ogu::Exception(__FILE__, __LINE__,
string("Parameter ") + param + " does not exist in FineAtoms copy cosntructor");
params.push_back(param);
param_outer_map.emplace(param, params.size()-1);
}
// fill in endovars
for (const auto &endovar : fa.endovars)
{
if (!varnames.query(endovar))
throw ogu::Exception(__FILE__, __LINE__,
string("Endo variable ") + endovar + " does not exist in FineAtoms copy constructor");
endovars.push_back(endovar);
endo_outer_map.emplace(endovar, endovars.size()-1);
}
// fill in exovars
for (const auto &exovar : fa.exovars)
{
if (!varnames.query(exovar))
throw ogu::Exception(__FILE__, __LINE__,
string("Exo variable ") + exovar + " does not exist in FineAtoms copy cosntructor");
exovars.push_back(exovar);
exo_outer_map.emplace(exovar, exovars.size()-1);
}
if (fa.endo_order)
endo_order = fa.endo_order->clone(endovars, *this);
if (fa.exo_order)
exo_order = fa.exo_order->clone(exovars, *this);
if (fa.allvar_order)
allvar_order = std::make_unique<AllvarOuterOrdering>(*(fa.allvar_order), *this);
}
int
FineAtoms::check_variable(const string &name) const
{
string str;
int ll;
parse_variable(name, str, ll);
if (varnames.query(str))
return DynamicAtoms::check_variable(name);
else
{
throw ParserException(string("Variable <")+str+"> not declared.", 0);
return -1;
}
}
int
FineAtoms::num_exo_periods() const
{
int mlead, mlag;
exovarspan(mlead, mlag);
return mlead-mlag+1;
}
void
FineAtoms::parsing_finished(VarOrdering::ord_type ot)
{
make_internal_orderings(ot);
// by default, concatenate outer endo and outer exo and make it as
// allvar outer:
vector<string> allvar_tmp;
allvar_tmp.insert(allvar_tmp.end(), endovars.begin(), endovars.end());
allvar_tmp.insert(allvar_tmp.end(), exovars.begin(), exovars.end());
allvar_order = std::make_unique<AllvarOuterOrdering>(allvar_tmp, *this);
}
void
FineAtoms::parsing_finished(VarOrdering::ord_type ot,
const vector<string> &allvar)
{
make_internal_orderings(ot);
allvar_order = std::make_unique<AllvarOuterOrdering>(allvar, *this);
}
const vector<string> &
FineAtoms::get_allvar() const
{
if (!allvar_order)
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::get_allvars called before parsing_finished");
return allvar_order->get_allvar();
}
const vector<int> &
FineAtoms::outer_endo2all() const
{
if (!allvar_order)
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::outer_endo2all called before parsing_finished");
return allvar_order->get_endo2all();
}
const vector<int> &
FineAtoms::outer_exo2all() const
{
if (!allvar_order)
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::outer_exo2all called before parsing_finished");
return allvar_order->get_exo2all();
}
vector<int>
FineAtoms::variables() const
{
if (endo_order)
return der_atoms;
else
{
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::variables called before parsing_finished");
return {};
}
}
int
FineAtoms::nstat() const
{
if (endo_order)
return endo_order->nstat();
else
{
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::nstat called before parsing_finished");
return -1;
}
}
int
FineAtoms::npred() const
{
if (endo_order)
return endo_order->npred();
else
{
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::npred called before parsing_finished");
return -1;
}
}
int
FineAtoms::nboth() const
{
if (endo_order)
return endo_order->nboth();
else
{
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::nboth called before parsing_finished");
return -1;
}
}
int
FineAtoms::nforw() const
{
if (endo_order)
return endo_order->nforw();
else
{
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::nforw called before parsing_finished");
return -1;
}
}
int
FineAtoms::get_pos_of_endo(int t) const
{
if (endo_order)
return endo_order->get_pos_of(t);
else
{
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::get_pos_of_endo called before parsing_finished");
return -1;
}
}
int
FineAtoms::get_pos_of_exo(int t) const
{
if (exo_order)
return exo_order->get_pos_of(t);
else
{
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::get_pos_of_exo called before parsing_finished");
return -1;
}
}
int
FineAtoms::get_pos_of_all(int t) const
{
if (endo_order && exo_order)
{
if (endo_order->check(t))
return endo_order->get_pos_of(t);
else if (exo_order->check(t))
return endo_order->length() + exo_order->get_pos_of(t);
else
{
throw ogu::Exception(__FILE__, __LINE__,
"Atom is not endo nor exo in FineAtoms::get_pos_of_all");
return -1;
}
}
else
{
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::get_pos_of_exo called before parsing_finished");
return -1;
}
}
const vector<int> &
FineAtoms::y2outer_endo() const
{
if (!endo_order)
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::y2outer_endo called before parsing_finished");
return endo_order->get_y2outer();
}
const vector<int> &
FineAtoms::outer2y_endo() const
{
if (!endo_order)
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::outer2y_endo called before parsing_finished");
return endo_order->get_outer2y();
}
const vector<int> &
FineAtoms::y2outer_exo() const
{
if (!exo_order)
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::y2outer_endo called before parsing_finished");
return exo_order->get_y2outer();
}
const vector<int> &
FineAtoms::outer2y_exo() const
{
if (!exo_order)
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::outer2y_exo called before parsing_finished");
return exo_order->get_outer2y();
}
const vector<int> &
FineAtoms::get_endo_atoms_map() const
{
if (!endo_order)
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::get_endo_atoms_map called before parsing_finished");
return endo_atoms_map;
}
const vector<int> &
FineAtoms::get_exo_atoms_map() const
{
if (!exo_order)
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::get_exo_atoms_map called before parsing_finished");
return exo_atoms_map;
}
int
FineAtoms::name2outer_param(const string &name) const
{
auto it = param_outer_map.find(name);
if (it == param_outer_map.end())
throw ogu::Exception(__FILE__, __LINE__,
"Name is not a parameter in FineAtoms::name2outer_param");
return it->second;
}
int
FineAtoms::name2outer_endo(const string &name) const
{
auto it = endo_outer_map.find(name);
if (it == endo_outer_map.end())
throw ogu::Exception(__FILE__, __LINE__,
"Name is not an endogenous variable in FineAtoms::name2outer_endo");
return it->second;
}
int
FineAtoms::name2outer_exo(const string &name) const
{
auto it = exo_outer_map.find(name);
if (it == exo_outer_map.end())
throw ogu::Exception(__FILE__, __LINE__,
"Name is not an exogenous variable in FineAtoms::name2outer_exo");
return it->second;
}
int
FineAtoms::name2outer_allvar(const string &name) const
{
if (!allvar_order)
throw ogu::Exception(__FILE__, __LINE__,
"FineAtoms::name2outer_allvar called beore parsing_finished");
auto it = endo_outer_map.find(name);
if (it != endo_outer_map.end())
return allvar_order->get_endo2all()[it->second];
else
{
it = exo_outer_map.find(name);
if (it != exo_outer_map.end())
return allvar_order->get_exo2all()[it->second];
}
throw ogu::Exception(__FILE__, __LINE__,
string("Name ") + name + " is neither endo nor exo variable in FineAtoms::name2outer_allvar");
return -1;
}
void
FineAtoms::register_uniq_endo(string name)
{
if (varnames.query(name))
throw ogp::ParserException(string("Endogenous variable <")+name+"> is not unique.", 0);
varnames.insert(name);
endovars.push_back(name);
endo_outer_map.emplace(std::move(name), endovars.size()-1);
}
void
FineAtoms::register_uniq_exo(string name)
{
if (varnames.query(name))
throw ogp::ParserException(string("Exogenous variable <")+name+"> is not unique.", 0);
varnames.insert(name);
exovars.push_back(name);
exo_outer_map.emplace(std::move(name), exovars.size()-1);
}
void
FineAtoms::register_uniq_param(string name)
{
if (varnames.query(name))
throw ogp::ParserException(string("Parameter <")+name+"> is not unique.", 0);
varnames.insert(name);
params.push_back(name);
param_outer_map.emplace(std::move(name), params.size()-1);
}
void
FineAtoms::make_internal_orderings(VarOrdering::ord_type ot)
{
bool endo_ordering_done = false;
bool exo_ordering_done = false;
order_type = ot;
int mlead, mlag;
endovarspan(mlead, mlag);
if (mlag >= -1 && mlead <= 1)
{
// make endo ordering
if (ot == VarOrdering::pbspbfbf)
endo_order = std::make_unique<EndoVarOrdering1>(endovars, *this);
else
endo_order = std::make_unique<EndoVarOrdering2>(endovars, *this);
endo_order->do_ordering();
endo_ordering_done = true;
}
exovarspan(mlead, mlag);
if (mlag == 0 && mlead == 0)
{
// make exo ordering
exo_order = std::make_unique<ExoVarOrdering>(exovars, *this);
exo_order->do_ordering();
exo_ordering_done = true;
}
if (endo_ordering_done && exo_ordering_done)
{
// concatenate der atoms from endo_order and exo_order
der_atoms.clear();
der_atoms.insert(der_atoms.end(),
endo_order->get_der_atoms().begin(),
endo_order->get_der_atoms().end());
der_atoms.insert(der_atoms.end(),
exo_order->get_der_atoms().begin(),
exo_order->get_der_atoms().end());
// create endo_atoms_map; der_atoms is a concatenation, so it is easy
int endo_atoms = endo_order->get_der_atoms().size();
endo_atoms_map.clear();
for (int i = 0; i < endo_atoms; i++)
endo_atoms_map.push_back(i);
// create exo_atoms_map
int exo_atoms = exo_order->get_der_atoms().size();
exo_atoms_map.clear();
for (int i = 0; i < exo_atoms; i++)
exo_atoms_map.push_back(endo_atoms + i);
}
}
void
FineAtoms::print() const
{
DynamicAtoms::print();
if (endo_order)
{
std::cout << "Endo ordering:\n";
endo_order->print();
}
else
std::cout << "Endo ordering not created.\n";
if (exo_order)
{
std::cout << "Exo ordering:\n";
exo_order->print();
}
else
std::cout << "Exo ordering not created.\n";
std::cout << "endo atoms map:\n";
for (unsigned int i = 0; i < endo_atoms_map.size(); i++)
std::cout << i << u8" → " << endo_atoms_map[i] << "\n";
std::cout << "exo atoms map:\n";
for (unsigned int i = 0; i < exo_atoms_map.size(); i++)
std::cout << i << u8" → " << exo_atoms_map[i] << "\n";
}
/*
* Copyright © 2005 Ondra Kamenik
* Copyright © 2019 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OGP_FINE_ATOMS_H
#define OGP_FINE_ATOMS_H
#include "dynamic_atoms.hh"
#include <vector>
#include <string>
#include <memory>
namespace ogp
{
using std::vector;
using std::string;
/** This is just ordering used for endogenous variables. It
* assumes that we have only time t-1, t, and t+1, orders them as
* pred(t-1), both(t-1), stat(t), pred(t), both(t), forw(t),
* both(t+1), forw(t+1). */
class EndoVarOrdering1 : public VarOrdering
{
public:
EndoVarOrdering1(const vector<string> &vnames, const DynamicAtoms &a)
: VarOrdering(vnames, a)
{
}
EndoVarOrdering1(const EndoVarOrdering1 &vo, const vector<string> &vnames,
const DynamicAtoms &a)
: VarOrdering(vo, vnames, a)
{
}
std::unique_ptr<VarOrdering>
clone(const vector<string> &vnames, const DynamicAtoms &a) const override
{
return std::make_unique<EndoVarOrdering1>(*this, vnames, a);
}
void
do_ordering() override
{
do_pbspbfbf();
}
};
/** This is just another ordering used for endogenous
* variables. It assumes that we have only time t-1, t, and t+1,
* orders them as both(t+1), forw(t+1), pred(t-1), both(t-1),
* stat(t), pred(t), both(t), forw(t). */
class EndoVarOrdering2 : public VarOrdering
{
public:
EndoVarOrdering2(const vector<string> &vnames, const DynamicAtoms &a)
: VarOrdering(vnames, a)
{
}
EndoVarOrdering2(const EndoVarOrdering2 &vo, const vector<string> &vnames,
const DynamicAtoms &a)
: VarOrdering(vo, vnames, a)
{
}
std::unique_ptr<VarOrdering>
clone(const vector<string> &vnames, const DynamicAtoms &a) const override
{
return std::make_unique<EndoVarOrdering2>(*this, vnames, a);
}
void
do_ordering() override
{
do_bfspbfpb();
}
};
/** This is just ordering used for exogenous variables. It makes
* no assumptions about their timing. It orders them from the
* least time to the latest time. */
class ExoVarOrdering : public VarOrdering
{
public:
ExoVarOrdering(const vector<string> &vnames, const DynamicAtoms &a)
: VarOrdering(vnames, a)
{
}
ExoVarOrdering(const ExoVarOrdering &vo, const vector<string> &vnames,
const DynamicAtoms &a)
: VarOrdering(vo, vnames, a)
{
}
std::unique_ptr<VarOrdering>
clone(const vector<string> &vnames, const DynamicAtoms &a) const override
{
return std::make_unique<ExoVarOrdering>(*this, vnames, a);
}
void
do_ordering() override
{
do_increasing_time();
}
};
class FineAtoms;
/** This class provides an outer ordering of all variables (endo
* and exo). It maps the ordering to the particular outer
* orderings of endo and exo. It works tightly with the FineAtoms
* class. */
class AllvarOuterOrdering
{
protected:
/** Type for a map mapping a variable name to an integer. */
using Tvarintmap = map<string, int>;
/** Reference to atoms. */
const FineAtoms &atoms;
/** The vector of all endo and exo variables in outer
* ordering. The pointers point to storage in atoms. */
vector<string> allvar;
/** The mapping from outer endogenous to outer all. For
* example endo2all[0] is the order of the first outer
* endogenous variable in the allvar ordering. */
vector<int> endo2all;
/** The mapping from outer exogenous to outer all. For example
* exo2all[0] is the order of the first outer exogenous
* variables in the allvar ordering. */
vector<int> exo2all;
public:
/** Construct the allvar outer ordering from the provided
* sequence of endo and exo names. The names can have an
* arbitrary storage, the storage is transformed to the atoms
* storage. An exception is thrown if either the list is not
* exhaustive, or some string is not a variable. */
AllvarOuterOrdering(const vector<string> &allvar_outer, const FineAtoms &a);
/** Copy constructor using the storage of provided atoms. */
AllvarOuterOrdering(const AllvarOuterOrdering &allvar_outer, const FineAtoms &a);
/** Return endo2all mapping. */
const vector<int> &
get_endo2all() const
{
return endo2all;
}
/** Return exo2all mapping. */
const vector<int> &
get_exo2all() const
{
return exo2all;
}
/** Return the allvar ordering. */
const vector<string> &
get_allvar() const
{
return allvar;
}
};
/** This class refines the DynamicAtoms by distinguishing among
* parameters (no lag and leads) and endogenous and exogenous
* variables (with lags and leads). For parameters, endogenous and
* exogenous, it defines outer orderings and internal
* orderings. The internal orderings are created by
* parsing_finished() method when it is sure that no new variables
* would be registered. The outer orderings are given by the order
* of calls of registering methods.
*
* In addition, the class also defines outer ordering of
* endogenous and exogenous variables. This is input as a
* parameter to parsing_finished(). By default, this whole outer
* ordering is just a concatenation of outer ordering of
* endogenous and exogenous variables.
*
* The internal ordering of all endo and exo variables is just a
* concatenation of endo and exo variables in their internal
* orderings. This is the ordering with respect to which all
* derivatives are taken. */
class FineAtoms : public DynamicAtoms
{
friend class AllvarOuterOrdering;
protected:
using Tvarintmap = map<string, int>;
private:
/** The vector of parameters names. The order gives the order
* the data is communicated with outside world. */
vector<string> params;
/** A map mapping a name of a parameter to an index in the outer
* ordering. */
Tvarintmap param_outer_map;
/** The vector of endogenous variables. This defines the order
* like parameters. */
vector<string> endovars;
/** A map mapping a name of an endogenous variable to an index
* in the outer ordering. */
Tvarintmap endo_outer_map;
/** The vector of exogenous variables. Also defines the order
* like parameters and endovars. */
vector<string> exovars;
/** A map mapping a name of an exogenous variable to an index
* in the outer ordering. */
Tvarintmap exo_outer_map;
protected:
/** This is the internal ordering of all atoms corresponding
* to endogenous variables. It is constructed by
* parsing_finished() method, which should be called after all
* parsing jobs have been finished. */
std::unique_ptr<VarOrdering> endo_order;
/** This is the internal ordering of all atoms corresponding
* to exogenous variables. It has the same handling as
* endo_order. */
std::unique_ptr<VarOrdering> exo_order;
/** This is the all variables outer ordering. It is
* constructed by parsing finished. */
std::unique_ptr<AllvarOuterOrdering> allvar_order;
/** This vector defines a set of atoms as tree indices used
* for differentiation. The order of the atoms in this vector
* defines ordering of the derivative tensors. The ordering is
* a concatenation of atoms from endo_order and then
* exo_order. This vector is setup by parsing_finished() and
* is returned by variables(). */
vector<int> der_atoms;
/** This is a mapping from endogenous atoms to all atoms in
* der_atoms member. The mapping maps index in endogenous atom
* ordering to index (not value) in der_atoms. It is useful if
* one wants to evaluate derivatives wrt only endogenous
* variables. It is set by parsing_finished(). By definition,
* it is monotone. */
vector<int> endo_atoms_map;
/** This is a mapping from exogenous atoms to all atoms in
* der_atoms member. It is the same as endo_atoms_map for
* atoms of exogenous variables. */
vector<int> exo_atoms_map;
public:
FineAtoms() = default;
FineAtoms(const FineAtoms &fa);
/** Deletes endo_order and exo_order. */
~FineAtoms() override = default;
/** Overrides DynamicAtoms::check_variable so that the error
* would be raised if the variable name is not declared. A
* variable is declared by inserting it to
* DynamicAtoms::varnames. This is a responsibility of a
* subclass. */
int check_variable(const string &name) const override;
/** This calculates min lag and max lead of endogenous variables. */
void
endovarspan(int &mlead, int &mlag) const
{
varspan(endovars, mlead, mlag);
}
/** This calculates mim lag and max lead of exogenous variables. */
void
exovarspan(int &mlead, int &mlag) const
{
varspan(exovars, mlead, mlag);
}
/** This calculates the number of periods in which at least
* one exogenous variable occurs. */
int num_exo_periods() const;
/** Return an (external) ordering of parameters. */
const vector<string> &
get_params() const
{
return params;
}
/** Return an external ordering of endogenous variables. */
const vector<string> &
get_endovars() const
{
return endovars;
}
/** Return an external ordering of exogenous variables. */
const vector<string> &
get_exovars() const
{
return exovars;
}
/** This constructs internal orderings and makes the indices
* returned by variables method available. Further it
* constructs outer ordering of all variables by a simple
* concatenation of outer endogenous and outer exogenous. In
* addition, it makes nstat, npred, nboth, nforw available. */
void parsing_finished(VarOrdering::ord_type ot);
/** This does the same thing as
* parsing_finished(VarOrdering::ord_type) plus it allows for
* inputing a different outer ordering of all variables. The
* ordering is input as a list of strings, their storage can
* be arbitrary. */
void parsing_finished(VarOrdering::ord_type ot, const vector<string> &avo);
/** Return the external ordering of all variables (endo and
* exo). This is either the second argument to
* parsing_finished or the default external ordering. This
* must be called only after parsing_finished. */
const vector<string> &get_allvar() const;
/** Return the map from outer ordering of endo variables to
* the allvar ordering. This must be called only after
* parsing_finished. */
const vector<int> &outer_endo2all() const;
/** Return the map from outer ordering of exo variables to
* the allvar ordering. This must be called only after
* parsing_finished. */
const vector<int> &outer_exo2all() const;
/** Return the atoms with respect to which we are going to
* differentiate. This must be called after
* parsing_finished. */
vector<int> variables() const override;
/** Return the number of static. */
int nstat() const;
/** Return the number of predetermined. */
int npred() const;
/** Return the number of both. */
int nboth() const;
/** Return the number of forward looking. */
int nforw() const;
/** Return the index of an endogenous atom given by tree index in
* the endo ordering. This must be also called only after
* parsing_finished(). */
int get_pos_of_endo(int t) const;
/** Return the index of an exogenous atom given by tree index in
* the exo ordering. This must be also called only after
* parsing_finished(). */
int get_pos_of_exo(int t) const;
/** Return the index of either endogenous or exogenous atom
* given by tree index in the concatenated ordering of
* endogenous and exogenous atoms. This must be also called
* only after parsing_finished(). */
int get_pos_of_all(int t) const;
/** Return the mapping from endogenous at time t to outer
* ordering of endogenous. */
const vector<int>&y2outer_endo() const;
/** Return the mapping from the outer ordering of endogenous to endogenous
* at time t. */
const vector<int>&outer2y_endo() const;
/** Return the mapping from exogenous at time t to outer
* ordering of exogenous. */
const vector<int>&y2outer_exo() const;
/** Return the mapping from the outer ordering of exogenous to exogenous
* at time t. */
const vector<int>&outer2y_exo() const;
/** Return the endo_atoms_map. */
const vector<int>&get_endo_atoms_map() const;
/** Return the exo_atoms_map. */
const vector<int>&get_exo_atoms_map() const;
/** Return an index in the outer ordering of a given
* parameter. An exception is thrown if the name is not a
* parameter. */
int name2outer_param(const string &name) const;
/** Return an index in the outer ordering of a given
* endogenous variable. An exception is thrown if the name is not a
* and endogenous variable. */
int name2outer_endo(const string &name) const;
/** Return an index in the outer ordering of a given
* exogenous variable. An exception is thrown if the name is not a
* and exogenous variable. */
int name2outer_exo(const string &name) const;
/** Return an index in the outer ordering of all variables
* (endo and exo) for a given name. An exception is thrown if
* the name is not a variable. This must be called only after
* parsing_finished(). */
int name2outer_allvar(const string &name) const;
/** Return the number of endogenous variables at time t-1, these are state
* variables. */
int
nys() const
{
return npred()+nboth();
}
/** Return the number of endogenous variables at time t+1. */
int
nyss() const
{
return nboth()+nforw();
}
/** Return the number of endogenous variables. */
int
ny() const
{
return endovars.size();
}
/** Return the number of exogenous variables. */
int
nexo() const
{
return static_cast<int>(exovars.size());
}
/** Return the number of parameters. */
int
np() const
{
return static_cast<int>(params.size());
}
/** Register unique endogenous variable name. The order of
* calls defines the endo outer ordering. The method is
* virtual, since a superclass may want to do some additional
* action. */
virtual void register_uniq_endo(string name);
/** Register unique exogenous variable name. The order of
* calls defines the exo outer ordering. The method is
* virtual, since a superclass may want to do somem additional
* action. */
virtual void register_uniq_exo(string name);
/** Register unique parameter name. The order of calls defines
* the param outer ordering. The method is
* virtual, since a superclass may want to do somem additional
* action. */
virtual void register_uniq_param(string name);
/** Debug print. */
void print() const override;
private:
/** This performs the common part of parsing_finished(), which
* is a construction of internal orderings. */
void make_internal_orderings(VarOrdering::ord_type ot);
protected:
/** This remembers the ordering type of the last call make_internal_ordering. */
VarOrdering::ord_type order_type;
};
};
#endif