Commit df772325 authored by Sébastien Villemot's avatar Sébastien Villemot

dynare++/kord: move away from CWEB

By the way apply Dynare C++ coding style and extensions (.cc/.hh).
parent 54335242
Pipeline #545 passed with stages
in 65 minutes and 29 seconds
......@@ -126,10 +126,6 @@ mex/build/matlab/run_m2html.m
/dynare++/integ/src/quadrature-points
/dynare++/integ/src/quadrature-points.exe
/dynare++/integ/testing/tests
/dynare++/kord/*.cpp
!/dynare++/kord/tests.cpp
/dynare++/kord/*.h
/dynare++/kord/main.tex
/dynare++/kord/tests
/dynare++/kord/tests.exe
/dynare++/kord/out.txt
......
......@@ -22,7 +22,7 @@
#include "dynmex.h"
#include "mex.h"
#include "decision_rule.h"
#include "decision_rule.hh"
#include "fs_tensor.h"
#include "SylvException.h"
......
CWEBSRC = \
faa_di_bruno.cweb \
korder_stoch.cweb \
journal.cweb \
decision_rule.cweb \
dynamic_model.cweb \
random.cweb \
first_order.cweb \
normal_conjugate.cweb \
approximation.cweb \
global_check.cweb \
korder.cweb \
kord_exception.hweb \
random.hweb \
journal.hweb \
approximation.hweb \
korder_stoch.hweb \
dynamic_model.hweb \
decision_rule.hweb \
korder.hweb \
normal_conjugate.hweb \
first_order.hweb \
mersenne_twister.hweb \
global_check.hweb \
faa_di_bruno.hweb
GENERATED_FILES = \
faa_di_bruno.cpp \
korder_stoch.cpp \
journal.cpp \
decision_rule.cpp \
dynamic_model.cpp \
random.cpp \
first_order.cpp \
normal_conjugate.cpp \
approximation.cpp \
global_check.cpp \
korder.cpp \
kord_exception.h \
random.h \
journal.h \
approximation.h \
korder_stoch.h \
dynamic_model.h \
decision_rule.h \
korder.h \
normal_conjugate.h \
first_order.h \
mersenne_twister.h \
global_check.h \
faa_di_bruno.h
noinst_LIBRARIES = libkord.a
libkord_a_SOURCES = $(CWEBSRC) $(GENERATED_FILES)
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 \
mersenne_twister.hh \
normal_conjugate.cc \
normal_conjugate.hh \
random.cc \
random.hh
libkord_a_CPPFLAGS = -I../sylv/cc -I../tl/cc -I../integ/cc -I$(top_srcdir)/mex/sources $(CPPFLAGS_MATIO)
libkord_a_CXXFLAGS = $(AM_CXXFLAGS) $(PTHREAD_CFLAGS)
BUILT_SOURCES = $(GENERATED_FILES)
EXTRA_DIST = main.web dummy.ch
check_PROGRAMS = tests
tests_SOURCES = tests.cpp
tests_SOURCES = tests.cc
tests_CPPFLAGS = -I../sylv/cc -I../tl/cc -I../integ/cc -I$(top_srcdir)/mex/sources
tests_CXXFLAGS = $(AM_CXXFLAGS) $(PTHREAD_CFLAGS)
tests_LDFLAGS = $(AM_LDFLAGS) $(LDFLAGS_MATIO)
......@@ -71,23 +40,4 @@ tests_LDADD = libkord.a ../tl/cc/libtl.a ../sylv/cc/libsylv.a $(LAPACK_LIBS) $(B
check-local:
./tests
%.cpp: %.cweb dummy.ch
$(CTANGLE) -bhp $< dummy.ch $@
%.h: %.hweb dummy.ch
$(CTANGLE) -bhp $< dummy.ch $@
if HAVE_CWEAVE
if HAVE_PDFTEX
if HAVE_EPLAIN
pdf-local: kord.pdf
kord.pdf: main.web $(CWEBSRC)
$(CWEAVE) -bhp main.web
$(PDFTEX) main
mv main.pdf kord.pdf
endif
endif
endif
CLEANFILES = kord.pdf main.idx main.log main.scn main.tex main.toc out.txt
CLEANFILES = out.txt
This diff is collapsed.
This diff is collapsed.
// Copyright 2005, Ondra Kamenik
// Approximating model solution
/* The class |Approximation| in this file is a main interface to the
algorithms calculating approximations to the decision rule about
deterministic and stochastic steady states.
The approximation about a deterministic steady state is solved by
classes |@<|FirstOrder| class declaration@>| and |@<|KOrder| class
declaration@>|. The approximation about the stochastic steady state is
solved by class |@<|KOrderStoch| class declaration@>| together with a
method of |Approximation| class |@<|Approximation::walkStochSteady|
code@>|.
The approximation about the stochastic steady state is done with
explicit expression of forward derivatives of $g^{**}$. More formally,
we have to solve the decision rule $g$ from the implicit system:
$$E_t(f(g^{**}(g^*(y^*,u_t,\sigma),u_{t+1},\sigma),g(y^*,u_t,\sigma),y_t,u_t))=0$$
The term within the expectations can be Taylor expanded, and the
expectation can be driven into the formula. However, when doing this
at $\sigma\not=0$, the term $g^{**}$ at $\sigma\not=0$ is dependent on
$u_{t+1}$ and thus the integral of its approximation includes all
derivatives wrt. $u$ of $g^{**}$. Note that for $\sigma=0$, the
derivatives of $g^{**}$ in this context are constant. This is the main
difference between the approximation at deterministic steady
($\sigma=0$), and stochastic steady ($\sigma\not=0$). This means that
$k$-order derivative of the above equation at $\sigma\not=0$ depends of
all derivatives of $g^**$ (including those with order greater than
$k$).
The explicit expression of the forward $g^{**}$ means that the
derivatives of $g$ are not solved simultaneously, but that the forward
derivatives of $g^{**}$ are calculated as an extrapolation based on
the approximation at lower $\sigma$. This is exactly what does the
|@<|Approximation::walkStochSteady| code@>|. It starts at the
deterministic steady state, and in a few steps it adds to $\sigma$
explicitly expressing forward $g^{**}$ from a previous step.
Further details on the both solution methods are given in (todo: put
references here when they exist).
Very important note: all classes here used for calculation of decision
rule approximation are folded. For the time being, it seems that faa
Di Bruno formula is quicker for folded tensors, and that is why we
stick to folded tensors here. However, when the calcs are done, we
calculate also its unfolded versions, to be available for simulations
and so on. */
#ifndef APPROXIMATION_H
#define APPROXIMATION_H
#include "dynamic_model.hh"
#include "decision_rule.hh"
#include "korder.hh"
#include "journal.hh"
/* This class is used to calculate derivatives by faa Di Bruno of the
$$f(g^{**}(g^*(y^*,u,\sigma),u',\sigma),g(y^*,u,\sigma),y^*,u)$$ with
respect $u'$. In order to keep it as simple as possible, the class
represents an equivalent (with respect to $u'$) container for
$f(g^{**}(y^*,u',\sigma),0,0,0)$. The class is used only for
evaluation of approximation error in |Approximation| class, which is
calculated in |Approximation::calcStochShift| method.
Since it is a folded version, we inherit from
|StackContainer<FGSTensor>| and |FoldedStackContainer|. To construct
it, we need only the $g^{**}$ container and size of stacks. */
class ZAuxContainer : public StackContainer<FGSTensor>, public FoldedStackContainer
{
public:
typedef StackContainer<FGSTensor>::_Ctype _Ctype;
typedef StackContainer<FGSTensor>::itype itype;
ZAuxContainer(const _Ctype *gss, int ngss, int ng, int ny, int nu);
itype getType(int i, const Symmetry &s) const;
};
/* This class provides an interface to approximation algorithms. The
core method is |walkStochSteady| which calculates the approximation
about stochastic steady state in a given number of steps. The number
is given as a parameter |ns| of the constructor. If the number is
equal to zero, the resulted approximation is about the deterministic
steady state.
An object is constructed from the |DynamicModel|, and the number of
steps |ns|. Also, we pass a reference to journal. That's all. The
result of the core method |walkStochSteady| is a decision rule |dr|
and a matrix |ss| whose columns are steady states for increasing
$\sigma$ during the walk. Both can be retrived by public methods. The
first column of the matrix is the deterministic steady state, the last
is the stochastic steady state for the full size shocks.
The method |walkStochSteady| calls the following methods:
|approxAtSteady| calculates an initial approximation about the
deterministic steady, |saveRuleDerivs| saves derivatives of a rule for
the following step in |rule_ders| and |rule_ders_ss| (see
|@<|Approximation::saveRuleDerivs| code@>| for their description),
|check| reports an error of the current approximation and
|calcStochShift| (called from |check|) calculates a shift of the
system equations due to uncertainity.
dr\_centralize is a new option. dynare++ was automatically expressing
results around the fixed point instead of the deterministic steady
state. dr\_centralize controls this behavior. */
class Approximation
{
DynamicModel &model;
Journal &journal;
FGSContainer *rule_ders;
FGSContainer *rule_ders_ss;
FoldDecisionRule *fdr;
UnfoldDecisionRule *udr;
const PartitionY ypart;
const FNormalMoments mom;
IntSequence nvs;
int steps;
bool dr_centralize;
double qz_criterium;
TwoDMatrix ss;
public:
Approximation(DynamicModel &m, Journal &j, int ns, bool dr_centr, double qz_crit);
virtual
~Approximation();
const FoldDecisionRule&getFoldDecisionRule() const;
const UnfoldDecisionRule&getUnfoldDecisionRule() const;
const TwoDMatrix &
getSS() const
{
return ss;
}
const DynamicModel &
getModel() const
{
return model;
}
void walkStochSteady();
TwoDMatrix *calcYCov() const;
const FGSContainer *
get_rule_ders() const
{
return rule_ders;
}
const FGSContainer *
get_rule_ders_ss() const
{
return rule_ders;
}
protected:
void approxAtSteady();
void calcStochShift(Vector &out, double at_sigma) const;
void saveRuleDerivs(const FGSContainer &g);
void check(double at_sigma) const;
};
#endif
@q $Id: approximation.hweb 2352 2009-09-03 19:18:15Z michel $ @>
@q Copyright 2005, Ondra Kamenik @>
@*2 Approximating model solution. Start of {\tt approximation.h} file.
The class |Approximation| in this file is a main interface to the
algorithms calculating approximations to the decision rule about
deterministic and stochastic steady states.
The approximation about a deterministic steady state is solved by
classes |@<|FirstOrder| class declaration@>| and |@<|KOrder| class
declaration@>|. The approximation about the stochastic steady state is
solved by class |@<|KOrderStoch| class declaration@>| together with a
method of |Approximation| class |@<|Approximation::walkStochSteady|
code@>|.
The approximation about the stochastic steady state is done with
explicit expression of forward derivatives of $g^{**}$. More formally,
we have to solve the decision rule $g$ from the implicit system:
$$E_t(f(g^{**}(g^*(y^*,u_t,\sigma),u_{t+1},\sigma),g(y^*,u_t,\sigma),y_t,u_t))=0$$
The term within the expectations can be Taylor expanded, and the
expectation can be driven into the formula. However, when doing this
at $\sigma\not=0$, the term $g^{**}$ at $\sigma\not=0$ is dependent on
$u_{t+1}$ and thus the integral of its approximation includes all
derivatives wrt. $u$ of $g^{**}$. Note that for $\sigma=0$, the
derivatives of $g^{**}$ in this context are constant. This is the main
difference between the approximation at deterministic steady
($\sigma=0$), and stochastic steady ($\sigma\not=0$). This means that
$k$-order derivative of the above equation at $\sigma\not=0$ depends of
all derivatives of $g^**$ (including those with order greater than
$k$).
The explicit expression of the forward $g^{**}$ means that the
derivatives of $g$ are not solved simultaneously, but that the forward
derivatives of $g^{**}$ are calculated as an extrapolation based on
the approximation at lower $\sigma$. This is exactly what does the
|@<|Approximation::walkStochSteady| code@>|. It starts at the
deterministic steady state, and in a few steps it adds to $\sigma$
explicitly expressing forward $g^{**}$ from a previous step.
Further details on the both solution methods are given in (todo: put
references here when they exist).
Very important note: all classes here used for calculation of decision
rule approximation are folded. For the time being, it seems that faa
Di Bruno formula is quicker for folded tensors, and that is why we
stick to folded tensors here. However, when the calcs are done, we
calculate also its unfolded versions, to be available for simulations
and so on.
@s ZAuxContainer int
@s Approximation int
@c
#ifndef APPROXIMATION_H
#define APPROXIMATION_H
#include "dynamic_model.h"
#include "decision_rule.h"
#include "korder.h"
#include "journal.h"
@<|ZAuxContainer| class declaration@>;
@<|Approximation| class declaration@>;
#endif
@ This class is used to calculate derivatives by faa Di Bruno of the
$$f(g^{**}(g^*(y^*,u,\sigma),u',\sigma),g(y^*,u,\sigma),y^*,u)$$ with
respect $u'$. In order to keep it as simple as possible, the class
represents an equivalent (with respect to $u'$) container for
$f(g^{**}(y^*,u',\sigma),0,0,0)$. The class is used only for
evaluation of approximation error in |Approximation| class, which is
calculated in |Approximation::calcStochShift| method.
Since it is a folded version, we inherit from
|StackContainer<FGSTensor>| and |FoldedStackContainer|. To construct
it, we need only the $g^{**}$ container and size of stacks.
@<|ZAuxContainer| class declaration@>=
class ZAuxContainer : public StackContainer<FGSTensor>, public FoldedStackContainer {
public:@;
typedef StackContainer<FGSTensor>::_Ctype _Ctype;
typedef StackContainer<FGSTensor>::itype itype;
ZAuxContainer(const _Ctype* gss, int ngss, int ng, int ny, int nu);
itype getType(int i, const Symmetry& s) const;
};
@ This class provides an interface to approximation algorithms. The
core method is |walkStochSteady| which calculates the approximation
about stochastic steady state in a given number of steps. The number
is given as a parameter |ns| of the constructor. If the number is
equal to zero, the resulted approximation is about the deterministic
steady state.
An object is constructed from the |DynamicModel|, and the number of
steps |ns|. Also, we pass a reference to journal. That's all. The
result of the core method |walkStochSteady| is a decision rule |dr|
and a matrix |ss| whose columns are steady states for increasing
$\sigma$ during the walk. Both can be retrived by public methods. The
first column of the matrix is the deterministic steady state, the last
is the stochastic steady state for the full size shocks.
The method |walkStochSteady| calls the following methods:
|approxAtSteady| calculates an initial approximation about the
deterministic steady, |saveRuleDerivs| saves derivatives of a rule for
the following step in |rule_ders| and |rule_ders_ss| (see
|@<|Approximation::saveRuleDerivs| code@>| for their description),
|check| reports an error of the current approximation and
|calcStochShift| (called from |check|) calculates a shift of the
system equations due to uncertainity.
dr\_centralize is a new option. dynare++ was automatically expressing
results around the fixed point instead of the deterministic steady
state. dr\_centralize controls this behavior.
@<|Approximation| class declaration@>=
class Approximation {
DynamicModel& model;
Journal& journal;
FGSContainer* rule_ders;
FGSContainer* rule_ders_ss;
FoldDecisionRule* fdr;
UnfoldDecisionRule* udr;
const PartitionY ypart;
const FNormalMoments mom;
IntSequence nvs;
int steps;
bool dr_centralize;
double qz_criterium;
TwoDMatrix ss;
public:@;
Approximation(DynamicModel& m, Journal& j, int ns, bool dr_centr, double qz_crit);
virtual ~Approximation();
const FoldDecisionRule& getFoldDecisionRule() const;
const UnfoldDecisionRule& getUnfoldDecisionRule() const;
const TwoDMatrix& getSS() const
{@+ return ss;@+}
const DynamicModel& getModel() const
{@+ return model;@+}
void walkStochSteady();
TwoDMatrix* calcYCov() const;
const FGSContainer* get_rule_ders() const
{@+ return rule_ders;@+}
const FGSContainer* get_rule_ders_ss() const
{@+ return rule_ders;@+}
protected:@;
void approxAtSteady();
void calcStochShift(Vector& out, double at_sigma) const;
void saveRuleDerivs(const FGSContainer& g);
void check(double at_sigma) const;
};
@ End of {\tt approximation.h} file.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
// Copyright 2005, Ondra Kamenik
#include "dynamic_model.hh"
void
NameList::print() const
{
for (int i = 0; i < getNum(); i++)
printf("%s\n", getName(i));
}
void
NameList::writeMat(mat_t *fd, const char *vname) const
{
int maxlen = 0;
for (int i = 0; i < getNum(); i++)
if (maxlen < (int) strlen(getName(i)))
maxlen = (int) strlen(getName(i));
if (maxlen == 0)
return;
char *m = new char[getNum()*maxlen];
for (int i = 0; i < getNum(); i++)
for (int j = 0; j < maxlen; j++)
if (j < (int) strlen(getName(i)))
m[j*getNum()+i] = getName(i)[j];
else
m[j*getNum()+i] = ' ';
#if MATIO_MAJOR_VERSION > 1 || (MATIO_MAJOR_VERSION == 1 && MATIO_MINOR_VERSION >= 5)
size_t dims[2];
const matio_compression compression = MAT_COMPRESSION_NONE;
#else
int dims[2];
const int compression = COMPRESSION_NONE;
#endif
dims[0] = getNum();
dims[1] = maxlen;
matvar_t *v = Mat_VarCreate(vname, MAT_C_CHAR, MAT_T_UINT8, 2, dims, m, 0);
Mat_VarWrite(fd, v, compression);
Mat_VarFree(v);
delete[] m;
}
void
NameList::writeMatIndices(mat_t *fd, const char *prefix) const
{
char tmp[100];
TwoDMatrix aux(1, 1);
for (int i = 0; i < getNum(); i++)
{
sprintf(tmp, "%s_i_%s", prefix, getName(i));
aux.get(0, 0) = i+1;
aux.writeMat(fd, tmp);
}
}
@q $Id: dynamic_model.cweb 431 2005-08-16 15:41:01Z kamenik $ @>
@q Copyright 2005, Ondra Kamenik @>
@ Start of {\tt dynamic\_model.cpp} file.
@c
#include "dynamic_model.h"
@<|NameList::print| code@>;
@<|NameList::writeMat| code@>;
@<|NameList::writeMatIndices| code@>;
@
@<|NameList::print| code@>=
void NameList::print() const
{
for (int i = 0; i < getNum(); i++)
printf("%s\n", getName(i));
}
@
@<|NameList::writeMat| code@>=
void NameList::writeMat(mat_t* fd, const char* vname) const
{
int maxlen = 0;
for (int i = 0; i < getNum(); i++)
if (maxlen < (int)strlen(getName(i)))
maxlen = (int)strlen(getName(i));
if (maxlen == 0)
return;
char *m = new char[getNum()*maxlen];
for (int i = 0; i < getNum(); i++)
for (int j = 0; j < maxlen; j++)
if (j < (int)strlen(getName(i)))
m[j*getNum()+i] = getName(i)[j];
else
m[j*getNum()+i] = ' ';
# if MATIO_MAJOR_VERSION > 1 || (MATIO_MAJOR_VERSION == 1 && MATIO_MINOR_VERSION >= 5)
size_t dims[2];
const matio_compression compression = MAT_COMPRESSION_NONE;
# else
int dims[2];
const int compression = COMPRESSION_NONE;
# endif
dims[0] = getNum();
dims[1] = maxlen;
matvar_t *v = Mat_VarCreate(vname, MAT_C_CHAR, MAT_T_UINT8, 2, dims, m, 0);
Mat_VarWrite(fd, v, compression);
Mat_VarFree(v);
delete[] m;
}
@
@<|NameList::writeMatIndices| code@>=
void NameList::writeMatIndices(mat_t* fd, const char* prefix) const
{
char tmp[100];
TwoDMatrix aux(1,1);
for (int i = 0; i < getNum(); i++) {
sprintf(tmp, "%s_i_%s", prefix, getName(i));
aux.get(0,0) = i+1;
aux.writeMat(fd, tmp);
}
}
@ End of {\tt dynamic\_model.cpp} file.
// Copyright 2005, Ondra Kamenik
// Dynamic model abstraction
/* This file only defines a generic interface to an SDGE model. The model
takes the form:
$$E_t\left[f(g^{**}(g^*(y,u_t),u_{t+1}),g(y,u),y,u_t)\right]=0$$
The interface is defined via pure virtual class |DynamicModel|. */
#ifndef DYNAMIC_MODEL_H
#define DYNAMIC_MODEL_H
#include "t_container.h"
#include "sparse_tensor.h"
#include "Vector.h"
/* The class is a virtual pure class which provides an access to names
of the variables. */
class NameList
{
public:
virtual ~NameList()
{
}
virtual int getNum() const = 0;
virtual const char *getName(int i) const = 0;
void print() const;
void writeMat(mat_t *fd, const char *vname) const;
void writeMatIndices(mat_t *fd, const char *prefix) const;
};
/* This is the interface to an information on a generic SDGE
model. It is sufficient for calculations of policy rule Taylor
approximations at some (not necessarily deterministic) steady state.
We need to know a partitioning of endogenous variables $y$. We suppose
that $y$ is partitioned as
$$y=\left[\matrix{\hbox{static}\cr\hbox{pred}\cr\hbox{both}\cr\hbox{forward}}\right]$$
of which we define
$$y^*=\left[\matrix{\hbox{pred}\cr\hbox{both}}\right]\quad
y^{**}=\left[\matrix{\hbox{both}\cr\hbox{forward}}\right]$$
where ``static'' are meant those variables, which appear only at time
$t$; ``pred'' are meant those variables, which appear only at $t$ and
$t-1$; ``both'' are meant those variables, which appear at least at
$t-1$ and $t+1$; and ``forward'' are meant those variables, which
appear only at $t$ and $t+1$. This partitioning is given by methods
|nstat()|, |npred()|, |nboth()|, and |nforw()|. The number of
equations |numeq()| must be the same as a number of endogenous
variables.
In order to complete description, we need to know a number of
exogenous variables, which is a size of $u$, hence |nexog()| method.
The model contains an information about names of variables, the
variance-covariance matrix of the shocks, the derivatives of equations
of $f$ at some steady state, and the steady state. These can be
retrieved by the corresponding methods.
The derivatives of the system are calculated with respect to stacked
variables, the stack looks as:
$$\left[\matrix{y^{**}_{t+1}\cr y_t\cr y^*_{t-1}\cr u_t}\right].$$
There are only three operations. The first
|solveDeterministicSteady()| solves the deterministic steady steate
which can be retrieved by |getSteady()| later. The method
|evaluateSystem| calculates $f(y^{**},y,y^*,u)$, where $y$ and $u$ are
passed, or $f(y^{**}_{t+1}, y_t, y^*_{t-1}, u)$, where $y^{**}_{t+1}$,
$y_t$, $y^*_{t-1}$, $u$ are passed. Finally, the method
|calcDerivativesAtSteady()| calculates derivatives of $f$ at the
current steady state, and zero shocks. The derivatives can be
retrieved with |getModelDerivatives()|. All the derivatives are done
up to a given order in the model, which can be retrieved by |order()|.
The model initialization is done in a constructor of the implementing
class. The constructor usually calls a parser, which parses a given
file (usually a text file), and retrieves all necessary information
about the model, inluding variables, partitioning, variance-covariance
matrix, information helpful for calculation of the deterministic
steady state, and so on. */
class DynamicModel
{
public:
virtual DynamicModel *clone() const = 0;
virtual ~DynamicModel()
{
}
virtual int nstat() const = 0;
virtual int nboth() const = 0;
virtual int npred() const = 0;
virtual int nforw() const = 0;
virtual int nexog() const = 0;
virtual int order() const = 0;
int
numeq() const
{
return nstat()+nboth()+npred()+nforw();
}
virtual const NameList&getAllEndoNames() const = 0;
virtual const NameList&getStateNames() const = 0;
virtual const NameList&getExogNames() const = 0;
virtual const TwoDMatrix&getVcov() const = 0;
virtual const TensorContainer<FSSparseTensor>&getModelDerivatives() const = 0;
virtual const Vector&getSteady() const = 0;
virtual Vector&getSteady() = 0;
virtual void solveDeterministicSteady() = 0;
virtual void evaluateSystem(Vector &out, const Vector &yy, const Vector &xx) = 0;
virtual void evaluateSystem(Vector &out, const Vector &yym, const Vector &yy,
const Vector &yyp, const Vector &xx) = 0;
virtual void calcDerivativesAtSteady() = 0;
};
#endif