Global reindentation of Dynare++

time-shift
Sébastien Villemot 2019-12-20 14:36:20 +01:00
parent 24a2aa2aaf
commit 22ab507c6c
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
75 changed files with 2281 additions and 2189 deletions

File diff suppressed because it is too large Load Diff

View File

@ -109,10 +109,10 @@ public:
See QuadratureImpl class declaration for details. */
template <typename _Tpit>
template<typename _Tpit>
class QuadratureImpl;
template <typename _Tpit>
template<typename _Tpit>
class IntegrationWorker : public sthread::detach_thread
{
const QuadratureImpl<_Tpit> &quad;
@ -168,7 +168,7 @@ public:
In addition, we define a method which saves all the points to a given file.
Only for debugging purposes. */
template <typename _Tpit>
template<typename _Tpit>
class QuadratureImpl : public Quadrature
{
friend class IntegrationWorker<_Tpit>;
@ -257,7 +257,7 @@ class OneDPrecalcQuadrature : public OneDQuadrature
public:
OneDPrecalcQuadrature(int nlevels, const int *npoints,
const double *wts, const double *pts)
: num_levels(nlevels), num_points(npoints),
: num_levels(nlevels), num_points(npoints),
weights(wts), points(pts), offsets(num_levels)
{
calcOffsets();

View File

@ -106,25 +106,26 @@ RadicalInverse::print() const
/* Here we have the first 170 primes. This means that we are not able to
integrate dimensions greater than 170. */
std::array<int, 170> HaltonSequence::primes = {
2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013
};
std::array<int, 170> HaltonSequence::primes =
{
2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013
};
/* This takes first dim primes and constructs dim radical inverses and
calls eval(). */
@ -145,7 +146,7 @@ HaltonSequence::HaltonSequence(int n, int mxn, int dim, const PermutationScheme
void
HaltonSequence::increase()
{
for (auto & i : ri)
for (auto &i : ri)
i.increase();
num++;
if (num <= maxn)
@ -165,7 +166,7 @@ void
HaltonSequence::print() const
{
auto ff = std::cout.flags();
for (const auto & i : ri)
for (const auto &i : ri)
i.print();
std::cout << "point=[ "
<< std::fixed << std::setprecision(6);

View File

@ -70,7 +70,6 @@ VectorFunctionSet::VectorFunctionSet(VectorFunction &f, int n)
}
}
/* Here we construct the object from the given function f and given
variance-covariance matrix Σ=vcov. The matrix A is calculated as lower
triangular and yields Σ=AA. */

View File

@ -47,7 +47,7 @@ struct QuadParams
QuadParams(int argc, char **argv);
void check_consistency() const;
private:
enum class opt {max_level, discard_weight, vcov};
enum class opt { max_level, discard_weight, vcov };
};
QuadParams::QuadParams(int argc, char **argv)
@ -62,11 +62,11 @@ QuadParams::QuadParams(int argc, char **argv)
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}
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;
@ -188,7 +188,7 @@ main(int argc, char **argv)
// Calculate weights and mass
double mass = 0.0;
std::vector<double> weights;
for (auto & point : points)
for (auto &point : points)
{
weights.push_back(std::exp(-point->dot(*point)));
mass += weights.back();

View File

@ -177,7 +177,7 @@ Function1Trans::eval(const Vector &point, const ParameterSignal &sig, Vector &ou
}
/* WallTimer class. Constructor saves the wall time, destructor cancels the
current time from the saved, and prints the message with time information */
current time from the saved, and prints the message with time information */
class WallTimer
{
std::string mes;

View File

@ -51,7 +51,7 @@ Approximation::Approximation(DynamicModel &m, Journal &j, int ns, bool dr_centr,
: model(m), journal(j),
ypart(model.nstat(), model.npred(), model.nboth(), model.nforw()),
mom(UNormalMoments(model.order(), model.getVcov())),
nvs{ypart.nys(), model.nexog(), model.nexog(), 1 },
nvs{ypart.nys(), model.nexog(), model.nexog(), 1},
steps(ns),
dr_centralize(dr_centr), qz_criterium(qz_crit), ss(ypart.ny(), steps+1)
{
@ -184,7 +184,7 @@ Approximation::walkStochSteady()
dy.add(-1.0, last_steady);
StochForwardDerivs<Storage::fold> hh(ypart, model.nexog(), *rule_ders_ss, mom, dy,
dsigma, sigma_so_far);
dsigma, sigma_so_far);
JournalRecord rec1(journal);
rec1 << "Calculation of g** expectations done" << endrec;

View File

@ -34,7 +34,7 @@
// 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)
udr.ypart, udr.nu, udr.ysteady)
{
for (const auto &it : udr)
insert(std::make_unique<ctraits<Storage::fold>::Ttensym>(*(it.second)));
@ -43,7 +43,7 @@ FoldDecisionRule::FoldDecisionRule(const UnfoldDecisionRule &udr)
// 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)
fdr.ypart, fdr.nu, fdr.ysteady)
{
for (const auto &it : fdr)
insert(std::make_unique<ctraits<Storage::unfold>::Ttensym>(*(it.second)));

View File

@ -75,7 +75,7 @@ public:
virtual void eval(emethod em, Vector &out, const ConstVector &v) const = 0;
/* makes only one step of simulation (in terms of absolute values, not
deviations) */
deviations) */
virtual void evaluate(emethod em, Vector &out, const ConstVector &ys,
const ConstVector &u) const = 0;
@ -120,7 +120,7 @@ public:
Storage::unfold. So, there are two implementations of the DecisionRule
interface. */
template <Storage t>
template<Storage t>
class DecisionRuleImpl : public ctraits<t>::Tpol, public DecisionRule
{
protected:
@ -220,7 +220,7 @@ public:
So we go through i+j=d=0q and in each loop we form the fully symmetric
tensor [g_(yu)] and insert it to the container. */
template <Storage t>
template<Storage t>
void
DecisionRuleImpl<t>::fillTensors(const _Tg &g, double sigma)
{
@ -277,7 +277,7 @@ DecisionRuleImpl<t>::fillTensors(const _Tg &g, double sigma)
where ȳ is the steady state of the original rule dr. */
template <Storage t>
template<Storage t>
void
DecisionRuleImpl<t>::centralize(const DecisionRuleImpl &dr)
{
@ -310,7 +310,7 @@ DecisionRuleImpl<t>::centralize(const DecisionRuleImpl &dr)
canceled from ystart, we simulate, and at the end ysteady is added to
all columns of the result. */
template <Storage t>
template<Storage t>
TwoDMatrix
DecisionRuleImpl<t>::simulate(emethod em, int np, const ConstVector &ystart,
ShockRealization &sr) const
@ -378,7 +378,7 @@ DecisionRuleImpl<t>::simulate(emethod em, int np, const ConstVector &ystart,
steady state (fix point) is cancelled and added once. Hence we have two
special methods. */
template <Storage t>
template<Storage t>
void
DecisionRuleImpl<t>::evaluate(emethod em, Vector &out, const ConstVector &ys,
const ConstVector &u) const
@ -401,7 +401,7 @@ DecisionRuleImpl<t>::evaluate(emethod em, Vector &out, const ConstVector &ys,
/* This is easy. We just return the newly created copy using the centralized
constructor. */
template <Storage t>
template<Storage t>
std::unique_ptr<DecisionRule>
DecisionRuleImpl<t>::centralizedClone(const Vector &fixpoint) const
{
@ -411,7 +411,7 @@ DecisionRuleImpl<t>::centralizedClone(const Vector &fixpoint) const
/* Here we only encapsulate two implementations to one, deciding according to
the parameter. */
template <Storage t>
template<Storage t>
void
DecisionRuleImpl<t>::eval(emethod em, Vector &out, const ConstVector &v) const
{
@ -423,7 +423,7 @@ DecisionRuleImpl<t>::eval(emethod em, Vector &out, const ConstVector &v) const
/* Write the decision rule and steady state to the MAT file. */
template <Storage t>
template<Storage t>
void
DecisionRuleImpl<t>::writeMat(mat_t *fd, const std::string &prefix) const
{
@ -506,7 +506,7 @@ public:
F is given by its derivatives bigf. The Jacobian of the solved system is
given by derivatives stored in bigfder. */
template <Storage t>
template<Storage t>
class DRFixPoint : public ctraits<t>::Tpol
{
using _Tpol = typename ctraits<t>::Tpol;
@ -557,7 +557,7 @@ private:
first derivative (Symmetry{1}). Then the derivative of the F polynomial is
calculated. */
template <Storage t>
template<Storage t>
DRFixPoint<t>::DRFixPoint(const _Tg &g, const PartitionY &yp,
const Vector &ys, double sigma)
: ctraits<t>::Tpol(yp.ny(), yp.nys()),
@ -577,7 +577,7 @@ DRFixPoint<t>::DRFixPoint(const _Tg &g, const PartitionY &yp,
such that d+k is between the maximum dimension and d, and add
σ/(d!k!)[g_yσ] to the tensor [g_y]. */
template <Storage t>
template<Storage t>
void
DRFixPoint<t>::fillTensors(const _Tg &g, double sigma)
{
@ -612,7 +612,7 @@ DRFixPoint<t>::fillTensors(const _Tg &g, double sigma)
underrelaxation parameter urelax, which improves the residual. If urelax
is less that urelax_threshold, we stop searching and stop the Newton. */
template <Storage t>
template<Storage t>
bool
DRFixPoint<t>::solveNewton(Vector &y)
{
@ -666,7 +666,7 @@ DRFixPoint<t>::solveNewton(Vector &y)
flastnorm = fnorm;
}
while (!converged && newton_iter_last < max_newton_iter
&&urelax > urelax_threshold);
&& urelax > urelax_threshold);
newton_iter_total += newton_iter_last;
if (!converged)
@ -688,7 +688,7 @@ DRFixPoint<t>::solveNewton(Vector &y)
The out vector is not touched if the algorithm has not convered. */
template <Storage t>
template<Storage t>
bool
DRFixPoint<t>::calcFixPoint(emethod em, Vector &out)
{

View File

@ -18,7 +18,6 @@
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
// Dynamic model abstraction
/* This file only defines a generic interface to a DSGE model. The model

View File

@ -57,7 +57,7 @@ public:
void calculate(const UnfoldedStackContainer &cont, const UGSContainer &g,
UGSTensor &out);
protected:
std::tuple<int,int,int> estimRefinement(const TensorDimens &tdims, int nr, int l);
std::tuple<int, int, int> estimRefinement(const TensorDimens &tdims, int nr, int l);
// See FaaDiBruno::calculate() folded sparse code for why we have magic_mult
constexpr static double magic_mult = 1.5;

View File

@ -178,9 +178,9 @@ FirstOrder::solve(const TwoDMatrix &fd)
std::lock_guard<std::mutex> lk{mut};
qz_criterium_global = qz_criterium;
dgges("N", "V", "S", order_eigs, &n, matE.getData().base(), &lda,
matD.getData().base(), &ldb, &sdim2, alphar.base(), alphai.base(),
beta.base(), vsl.getData().base(), &ldvsl, vsr.getData().base(), &ldvsr,
work.base(), &lwork, bwork.get(), &info);
matD.getData().base(), &ldb, &sdim2, alphar.base(), alphai.base(),
beta.base(), vsl.getData().base(), &ldvsl, vsr.getData().base(), &ldvsr,
work.base(), &lwork, bwork.get(), &info);
}
if (info)
throw KordException(__FILE__, __LINE__,

View File

@ -31,7 +31,7 @@ template<Storage>
class FirstOrderDerivs;
class FirstOrder
{
template <Storage>
template<Storage>
friend class FirstOrderDerivs;
PartitionY ypart;
int nu;
@ -89,7 +89,7 @@ protected:
/* This class only converts the derivatives g_y* and gᵤ to a folded or unfolded
container. */
template <Storage t>
template<Storage t>
class FirstOrderDerivs : public ctraits<t>::Tg
{
public:

View File

@ -90,7 +90,7 @@ protected:
public:
ResidFunction(const Approximation &app);
ResidFunction(const ResidFunction &rf);
std::unique_ptr<VectorFunction>
clone() const override
{

View File

@ -145,7 +145,7 @@ JournalRecord::writePrefix(const SystemResources &f)
std::ostringstream s;
s << std::setfill('0');
writeFloatTabular(s, f.elapsed, 7);
s << u8"" << recChar << std::setw(5) << ord << u8"";
s << u8"" << recChar << std::setw(5) << ord << u8"";
writeFloatTabular(s, f.load_avg, 3);
s << u8"";
writeFloatTabular(s, f.mem_avail/mb, 5);

View File

@ -26,18 +26,18 @@
#include <iostream>
#ifndef KORD_EXCEPTION_H
#define KORD_EXCEPTION_H
# define KORD_EXCEPTION_H
#define KORD_RAISE(mes) \
# define KORD_RAISE(mes) \
throw KordException(__FILE__, __LINE__, mes);
#define KORD_RAISE_IF(expr, mes) \
# define KORD_RAISE_IF(expr, mes) \
if (expr) throw KordException(__FILE__, __LINE__, mes);
#define KORD_RAISE_X(mes, c) \
# define KORD_RAISE_X(mes, c) \
throw KordException(__FILE__, __LINE__, mes, c);
#define KORD_RAISE_IF_X(expr, mes, c) \
# define KORD_RAISE_IF_X(expr, mes, c) \
if (expr) throw KordException(__FILE__, __LINE__, mes, c);
class KordException

View File

@ -112,117 +112,173 @@ MatrixS::MatrixS(const FSSparseTensor &f, const IntSequence &ss,
interesting here. */
template<>
ctraits<Storage::unfold>::Tg& KOrder::g<Storage::unfold>()
ctraits<Storage::unfold>::Tg &
KOrder::g<Storage::unfold>()
{
return _ug;
}
template<>
const ctraits<Storage::unfold>::Tg& KOrder::g<Storage::unfold>() const
{ return _ug;}
const ctraits<Storage::unfold>::Tg &
KOrder::g<Storage::unfold>() const
{
return _ug;
}
template<>
ctraits<Storage::fold>::Tg& KOrder::g<Storage::fold>()
ctraits<Storage::fold>::Tg &
KOrder::g<Storage::fold>()
{
return _fg;
}
template<>
const ctraits<Storage::fold>::Tg& KOrder::g<Storage::fold>() const
{ return _fg;}
const ctraits<Storage::fold>::Tg &
KOrder::g<Storage::fold>() const
{
return _fg;
}
template<>
ctraits<Storage::unfold>::Tgs& KOrder::gs<Storage::unfold>()
ctraits<Storage::unfold>::Tgs &
KOrder::gs<Storage::unfold>()
{
return _ugs;
}
template<>
const ctraits<Storage::unfold>::Tgs& KOrder::gs<Storage::unfold>() const
{ return _ugs;}
const ctraits<Storage::unfold>::Tgs &
KOrder::gs<Storage::unfold>() const
{
return _ugs;
}
template<>
ctraits<Storage::fold>::Tgs& KOrder::gs<Storage::fold>()
ctraits<Storage::fold>::Tgs &
KOrder::gs<Storage::fold>()
{
return _fgs;
}
template<>
const ctraits<Storage::fold>::Tgs& KOrder::gs<Storage::fold>() const
{ return _fgs;}
const ctraits<Storage::fold>::Tgs &
KOrder::gs<Storage::fold>() const
{
return _fgs;
}
template<>
ctraits<Storage::unfold>::Tgss& KOrder::gss<Storage::unfold>()
ctraits<Storage::unfold>::Tgss &
KOrder::gss<Storage::unfold>()
{
return _ugss;
}
template<>
const ctraits<Storage::unfold>::Tgss& KOrder::gss<Storage::unfold>() const
{ return _ugss;}
const ctraits<Storage::unfold>::Tgss &
KOrder::gss<Storage::unfold>() const
{
return _ugss;
}
template<>
ctraits<Storage::fold>::Tgss& KOrder::gss<Storage::fold>()
ctraits<Storage::fold>::Tgss &
KOrder::gss<Storage::fold>()
{
return _fgss;
}
template<>
const ctraits<Storage::fold>::Tgss& KOrder::gss<Storage::fold>() const
{ return _fgss;}
const ctraits<Storage::fold>::Tgss &
KOrder::gss<Storage::fold>() const
{
return _fgss;
}
template<>
ctraits<Storage::unfold>::TG& KOrder::G<Storage::unfold>()
ctraits<Storage::unfold>::TG &
KOrder::G<Storage::unfold>()
{
return _uG;
}
template<>
const ctraits<Storage::unfold>::TG& KOrder::G<Storage::unfold>() const
{ return _uG;}
const ctraits<Storage::unfold>::TG &
KOrder::G<Storage::unfold>() const
{
return _uG;
}
template<>
ctraits<Storage::fold>::TG& KOrder::G<Storage::fold>()
ctraits<Storage::fold>::TG &
KOrder::G<Storage::fold>()
{
return _fG;
}
template<>
const ctraits<Storage::fold>::TG& KOrder::G<Storage::fold>() const
{ return _fG;}
const ctraits<Storage::fold>::TG &
KOrder::G<Storage::fold>() const
{
return _fG;
}
template<>
ctraits<Storage::unfold>::TZstack& KOrder::Zstack<Storage::unfold>()
ctraits<Storage::unfold>::TZstack &
KOrder::Zstack<Storage::unfold>()
{
return _uZstack;
}
template<>
const ctraits<Storage::unfold>::TZstack& KOrder::Zstack<Storage::unfold>() const
{ return _uZstack;}
const ctraits<Storage::unfold>::TZstack &
KOrder::Zstack<Storage::unfold>() const
{
return _uZstack;
}
template<>
ctraits<Storage::fold>::TZstack& KOrder::Zstack<Storage::fold>()
ctraits<Storage::fold>::TZstack &
KOrder::Zstack<Storage::fold>()
{
return _fZstack;
}
template<>
const ctraits<Storage::fold>::TZstack& KOrder::Zstack<Storage::fold>() const
{ return _fZstack;}
const ctraits<Storage::fold>::TZstack &
KOrder::Zstack<Storage::fold>() const
{
return _fZstack;
}
template<>
ctraits<Storage::unfold>::TGstack& KOrder::Gstack<Storage::unfold>()
ctraits<Storage::unfold>::TGstack &
KOrder::Gstack<Storage::unfold>()
{
return _uGstack;
}
template<>
const ctraits<Storage::unfold>::TGstack& KOrder::Gstack<Storage::unfold>() const
{ return _uGstack;}
const ctraits<Storage::unfold>::TGstack &
KOrder::Gstack<Storage::unfold>() const
{
return _uGstack;
}
template<>
ctraits<Storage::fold>::TGstack& KOrder::Gstack<Storage::fold>()
ctraits<Storage::fold>::TGstack &
KOrder::Gstack<Storage::fold>()
{
return _fGstack;
}
template<>
const ctraits<Storage::fold>::TGstack& KOrder::Gstack<Storage::fold>() const
{ return _fGstack;}
const ctraits<Storage::fold>::TGstack &
KOrder::Gstack<Storage::fold>() const
{
return _fGstack;
}
template<>
ctraits<Storage::unfold>::Tm& KOrder::m<Storage::unfold>()
ctraits<Storage::unfold>::Tm &
KOrder::m<Storage::unfold>()
{
return _um;
}
template<>
const ctraits<Storage::unfold>::Tm& KOrder::m<Storage::unfold>() const
{ return _um;}
const ctraits<Storage::unfold>::Tm &
KOrder::m<Storage::unfold>() const
{
return _um;
}
template<>
ctraits<Storage::fold>::Tm& KOrder::m<Storage::fold>()
ctraits<Storage::fold>::Tm &
KOrder::m<Storage::fold>()
{
return _fm;
}
template<>
const ctraits<Storage::fold>::Tm& KOrder::m<Storage::fold>() const
{ return _fm;}
const ctraits<Storage::fold>::Tm &
KOrder::m<Storage::fold>() const
{
return _fm;
}
/* Here is the constructor of the KOrder class. We pass what we have to. The
partitioning of the y vector, a sparse container with model derivatives,

View File

@ -72,7 +72,7 @@ class UnfoldedZXContainer;
class FoldedGXContainer;
class UnfoldedGXContainer;
template <Storage type>
template<Storage type>
class ctraits
{
public:
@ -296,15 +296,15 @@ public:
/* Performs k-order step provided that k=2 or the k1-th step has been
run, this is the core method */
template <Storage t>
template<Storage t>
void performStep(int order);
/* Calculates residuals of all solved equations for k-order and reports their
sizes, it is runnable after k-order performStep() has been run */
template <Storage t>
template<Storage t>
double check(int dim) const;
template <Storage t>
template<Storage t>
Vector calcStochShift(int order, double sigma) const;
void switchToFolded();
const PartitionY &
@ -330,7 +330,7 @@ public:
protected:
/* Inserts a g derivative to the g container and also creates subtensors and
insert them to g_y* and g_y** containers */
template <Storage t>
template<Storage t>
void insertDerivative(std::unique_ptr<typename ctraits<t>::Ttensor> der);
/* Solves the sylvester equation (templated fold, and unfold) */
@ -339,27 +339,27 @@ protected:
/* Calculates derivatives of F by Faà Di Bruno for the sparse container of
system derivatives and Z stack container */
template <Storage t>
template<Storage t>
std::unique_ptr<typename ctraits<t>::Ttensor> faaDiBrunoZ(const Symmetry &sym) const;
/* Calculates derivatives of G by Faà Di Bruno for the dense container g**
and G stack */
template <Storage t>
template<Storage t>
std::unique_ptr<typename ctraits<t>::Ttensor> faaDiBrunoG(const Symmetry &sym) const;
// Recovers g_y*ⁱ
template<Storage t>
void recover_y(int i);
// Recovers g_y*ⁱuʲ
template <Storage t>
template<Storage t>
void recover_yu(int i, int j);
// Recovers g_y*ⁱσʲ
template <Storage t>
template<Storage t>
void recover_ys(int i, int j);
// Recovers g_y*ⁱuʲσ
template <Storage t>
template<Storage t>
void recover_yus(int i, int j, int k);
template <Storage t>
template<Storage t>
// Recovers g_σ
void recover_s(int i);
// Calculates specified derivatives of G and inserts them to the container
@ -367,26 +367,26 @@ protected:
void fillG(int i, int j, int k);
// Calculates Dᵢⱼₖ
template <Storage t>
template<Storage t>
typename ctraits<t>::Ttensor calcD_ijk(int i, int j, int k) const;
template <Storage t>
template<Storage t>
typename ctraits<t>::Ttensor calcD_ik(int i, int k) const;
template <Storage t>
template<Storage t>
typename ctraits<t>::Ttensor calcD_k(int k) const;
// Calculates Eᵢⱼₖ
template <Storage t>
template<Storage t>
typename ctraits<t>::Ttensor calcE_ijk(int i, int j, int k) const;
template <Storage t>
template<Storage t>
typename ctraits<t>::Ttensor calcE_ik(int i, int k) const;
template <Storage t>
template<Storage t>
typename ctraits<t>::Ttensor calcE_k(int k) const;
};
/* Here we insert the result to the container. Along the insertion, we
also create subtensors and insert as well. */
template <Storage t>
template<Storage t>
void
KOrder::insertDerivative(std::unique_ptr<typename ctraits<t>::Ttensor> der)
{
@ -405,7 +405,7 @@ KOrder::insertDerivative(std::unique_ptr<typename ctraits<t>::Ttensor> der)
where s is a given outer symmetry and k is the dimension of the symmetry. */
template <Storage t>
template<Storage t>
std::unique_ptr<typename ctraits<t>::Ttensor>
KOrder::faaDiBrunoZ(const Symmetry &sym) const
{
@ -419,7 +419,7 @@ KOrder::faaDiBrunoZ(const Symmetry &sym) const
/* The same as KOrder::faaDiBrunoZ(), but for g** and G stack. */
template <Storage t>
template<Storage t>
std::unique_ptr<typename ctraits<t>::Ttensor>
KOrder::faaDiBrunoG(const Symmetry &sym) const
{
@ -441,7 +441,7 @@ KOrder::faaDiBrunoG(const Symmetry &sym) const
Requires: everything at order i1
Provides: g_y and G_y
*/
*/
template<Storage t>
void
KOrder::recover_y(int i)
@ -475,8 +475,8 @@ KOrder::recover_y(int i)
Requires: everything at order i+j1, G_yʲ and g_yʲ.
Provides: g_yuʲ and G_yuʲ
*/
template <Storage t>
*/
template<Storage t>
void
KOrder::recover_yu(int i, int j)
{
@ -508,8 +508,8 @@ KOrder::recover_yu(int i, int j)
Provides: g_yσʲ and G_yσʲ, and finally G_yuσʲ for m=1,,j. The latter
is calculated by fillG() before the actual calculation.
*/
template <Storage t>
*/
template<Storage t>
void
KOrder::recover_ys(int i, int j)
{
@ -561,8 +561,8 @@ KOrder::recover_ys(int i, int j)
through E.
Provides: g_yuʲσ, G_yuʲσ, and G_yuʲuσ for m=1,,k
*/
template <Storage t>
*/
template<Storage t>
void
KOrder::recover_yus(int i, int j, int k)
{
@ -619,8 +619,8 @@ KOrder::recover_yus(int i, int j, int k)
F_u, and g_yuʲσ for j=1,,i1 and m+j+k=i through F_uʲσʲ.
Provides: g_σ, G_σ, and G_uσ for m=1,,i
*/
template <Storage t>
*/
template<Storage t>
void
KOrder::recover_s(int i)
{
@ -680,7 +680,7 @@ KOrder::fillG(int i, int j, int k)
So it is non zero only for even k. */
template <Storage t>
template<Storage t>
typename ctraits<t>::Ttensor
KOrder::calcD_ijk(int i, int j, int k) const
{
@ -700,7 +700,7 @@ KOrder::calcD_ijk(int i, int j, int k) const
¹
The sum can sum only for even m. */
template <Storage t>
template<Storage t>
typename ctraits<t>::Ttensor
KOrder::calcE_ijk(int i, int j, int k) const
{
@ -715,28 +715,28 @@ KOrder::calcE_ijk(int i, int j, int k) const
return res;
}
template <Storage t>
template<Storage t>
typename ctraits<t>::Ttensor
KOrder::calcD_ik(int i, int k) const
{
return calcD_ijk<t>(i, 0, k);
}
template <Storage t>
template<Storage t>
typename ctraits<t>::Ttensor
KOrder::calcD_k(int k) const
{
return calcD_ijk<t>(0, 0, k);
}
template <Storage t>
template<Storage t>
typename ctraits<t>::Ttensor
KOrder::calcE_ik(int i, int k) const
{
return calcE_ijk<t>(i, 0, k);
}
template <Storage t>
template<Storage t>
typename ctraits<t>::Ttensor
KOrder::calcE_k(int k) const
{
@ -763,7 +763,7 @@ KOrder::calcE_k(int k) const
all the recovering methods, he should find out that also all G are
provided. */
template <Storage t>
template<Storage t>
void
KOrder::performStep(int order)
{
@ -794,7 +794,7 @@ KOrder::performStep(int order)
The method returns the largest residual size. Each check simply evaluates
the equation. */
template <Storage t>
template<Storage t>
double
KOrder::check(int dim) const
{
@ -808,7 +808,7 @@ KOrder::check(int dim) const
// Check for F_yⁱuʲ=0
for (int i = 0; i <= dim; i++)
{
Symmetry sym{dim-i, i, 0, 0};
Symmetry sym{dim-i, i, 0, 0};
auto r = faaDiBrunoZ<t>(sym);
double err = r->getData().getMax();
JournalRecord(journal) << "\terror for symmetry " << sym << "\tis " << err << endrec;
@ -849,7 +849,7 @@ KOrder::check(int dim) const
return maxerror;
}
template <Storage t>
template<Storage t>
Vector
KOrder::calcStochShift(int order, double sigma) const
{

View File

@ -77,112 +77,134 @@ KOrderStoch::KOrderStoch(const PartitionY &yp, int nu,
// KOrderStoch convenience method specializations
template<>
ctraits<Storage::unfold>::Tg &KOrderStoch::g<Storage::unfold>()
ctraits<Storage::unfold>::Tg &
KOrderStoch::g<Storage::unfold>()
{
return _ug;
}
template<>
const ctraits<Storage::unfold>::Tg &KOrderStoch::g<Storage::unfold>() const
const ctraits<Storage::unfold>::Tg &
KOrderStoch::g<Storage::unfold>() const
{
return _ug;
}
template<>
ctraits<Storage::fold>::Tg &KOrderStoch::g<Storage::fold>()
ctraits<Storage::fold>::Tg &
KOrderStoch::g<Storage::fold>()
{
return _fg;
}
template<>
const ctraits<Storage::fold>::Tg &KOrderStoch::g<Storage::fold>() const
const ctraits<Storage::fold>::Tg &
KOrderStoch::g<Storage::fold>() const
{
return _fg;
}
template<>
ctraits<Storage::unfold>::Tgs &KOrderStoch::gs<Storage::unfold>()
ctraits<Storage::unfold>::Tgs &
KOrderStoch::gs<Storage::unfold>()
{
return _ugs;
}
template<>
const ctraits<Storage::unfold>::Tgs &KOrderStoch::gs<Storage::unfold>() const
const ctraits<Storage::unfold>::Tgs &
KOrderStoch::gs<Storage::unfold>() const
{
return _ugs;
}
template<>
ctraits<Storage::fold>::Tgs &KOrderStoch::gs<Storage::fold>()
ctraits<Storage::fold>::Tgs &
KOrderStoch::gs<Storage::fold>()
{
return _fgs;
}
template<>
const ctraits<Storage::fold>::Tgs &KOrderStoch::gs<Storage::fold>() const
const ctraits<Storage::fold>::Tgs &
KOrderStoch::gs<Storage::fold>() const
{
return _fgs;
}
template<>
const ctraits<Storage::unfold>::Tgss &KOrderStoch::h<Storage::unfold>() const
const ctraits<Storage::unfold>::Tgss &
KOrderStoch::h<Storage::unfold>() const
{
return *_uh;
}
template<>
const ctraits<Storage::fold>::Tgss &KOrderStoch::h<Storage::fold>() const
const ctraits<Storage::fold>::Tgss &
KOrderStoch::h<Storage::fold>() const
{
return *_fh;
}
template<>
ctraits<Storage::unfold>::TG &KOrderStoch::G<Storage::unfold>()
ctraits<Storage::unfold>::TG &
KOrderStoch::G<Storage::unfold>()
{
return _uG;
}
template<>
const ctraits<Storage::unfold>::TG &KOrderStoch::G<Storage::unfold>() const
const ctraits<Storage::unfold>::TG &
KOrderStoch::G<Storage::unfold>() const
{
return _uG;
}
template<>
ctraits<Storage::fold>::TG &KOrderStoch::G<Storage::fold>()
ctraits<Storage::fold>::TG &
KOrderStoch::G<Storage::fold>()
{
return _fG;
}
template<>
const ctraits<Storage::fold>::TG& KOrderStoch::G<Storage::fold>() const
const ctraits<Storage::fold>::TG &
KOrderStoch::G<Storage::fold>() const
{
return _fG;
}
template<>
ctraits<Storage::unfold>::TZXstack &KOrderStoch::Zstack<Storage::unfold>()
ctraits<Storage::unfold>::TZXstack &
KOrderStoch::Zstack<Storage::unfold>()
{
return _uZstack;
}
template<>
const ctraits<Storage::unfold>::TZXstack &KOrderStoch::Zstack<Storage::unfold>() const
const ctraits<Storage::unfold>::TZXstack &
KOrderStoch::Zstack<Storage::unfold>() const
{
return _uZstack;
}
template<>
ctraits<Storage::fold>::TZXstack &KOrderStoch::Zstack<Storage::fold>()
ctraits<Storage::fold>::TZXstack &
KOrderStoch::Zstack<Storage::fold>()
{
return _fZstack;
}
template<>
const ctraits<Storage::fold>::TZXstack &KOrderStoch::Zstack<Storage::fold>() const
const ctraits<Storage::fold>::TZXstack &
KOrderStoch::Zstack<Storage::fold>() const
{
return _fZstack;
}
template<>
ctraits<Storage::unfold>::TGXstack &KOrderStoch::Gstack<Storage::unfold>()
ctraits<Storage::unfold>::TGXstack &
KOrderStoch::Gstack<Storage::unfold>()
{
return _uGstack;
}
template<>
const ctraits<Storage::unfold>::TGXstack &KOrderStoch::Gstack<Storage::unfold>() const
const ctraits<Storage::unfold>::TGXstack &
KOrderStoch::Gstack<Storage::unfold>() const
{
return _uGstack;
}
template<>
ctraits<Storage::fold>::TGXstack &KOrderStoch::Gstack<Storage::fold>()
ctraits<Storage::fold>::TGXstack &
KOrderStoch::Gstack<Storage::fold>()
{
return _fGstack;
}
template<>
const ctraits<Storage::fold>::TGXstack &KOrderStoch::Gstack<Storage::fold>() const
const ctraits<Storage::fold>::TGXstack &
KOrderStoch::Gstack<Storage::fold>() const
{
return _fGstack;
}

View File

@ -36,7 +36,7 @@
containers corresponding to f(h(g*(y*,u,σ),σ),g(y,u,σ),y,u). Further, we
declare IntegDerivs and StochForwardDerivs classes which basically calculate
h as an extrapolation based on an approximation to g at lower σ.
*/
*/
#include <memory>
@ -48,7 +48,7 @@
/* This class is a container, which has a specialized constructor integrating
the policy rule at given σ. */
template <Storage t>
template<Storage t>
class IntegDerivs : public ctraits<t>::Tgss
{
public:
@ -100,8 +100,8 @@ public:
and this is exactly what the code does.
*/
template <Storage t>
*/
template<Storage t>
IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const typename ctraits<t>::Tgss &g,
const typename ctraits<t>::Tm &mom, double at_sigma)
: ctraits<t>::Tgss(4)
@ -156,7 +156,7 @@ IntegDerivs<t>::IntegDerivs(int r, const IntSequence &nvs, const typename ctrait
The class calculates derivatives of E[g(y*,u,σ)] at (ȳ*,σ¯). The derivatives
are extrapolated based on derivatives at (*,σ~). */
template <Storage t>
template<Storage t>
class StochForwardDerivs : public ctraits<t>::Tgss
{
public:
@ -180,7 +180,7 @@ public:
Centralize this polynomial about (ȳ,σ¯)
Recover general symmetry tensors from the (full symmetric) polynomial
*/
template <Storage t>
template<Storage t>
StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
const typename ctraits<t>::Tgss &g,
const typename ctraits<t>::Tm &m,
@ -260,7 +260,7 @@ StochForwardDerivs<t>::StochForwardDerivs(const PartitionY &ypart, int nu,
the σ as a second argument to h will be its fourth variable in symmetry, so
we have to do four member stack having the second and third stack dummy. */
template <class _Ttype>
template<class _Ttype>
class GXContainer : public GContainer<_Ttype>
{
public:
@ -279,8 +279,8 @@ public:
dummy
dummy
σ
*/
template <class _Ttype>
*/
template<class _Ttype>
typename GXContainer<_Ttype>::itype
GXContainer<_Ttype>::getType(int i, const Symmetry &s) const
{
@ -306,7 +306,7 @@ GXContainer<_Ttype>::getType(int i, const Symmetry &s) const
the size (number of rows) as g**. Since it is very simmilar to ZContainer,
we inherit form it and override only getType() method. */
template <class _Ttype>
template<class _Ttype>
class ZXContainer : public ZContainer<_Ttype>
{
public:
@ -326,7 +326,7 @@ public:
y
u
*/
template <class _Ttype>
template<class _Ttype>
typename ZXContainer<_Ttype>::itype
ZXContainer<_Ttype>::getType(int i, const Symmetry &s) const
{
@ -446,7 +446,7 @@ public:
const FGSContainer &hh, Journal &jr);
KOrderStoch(const PartitionY &ypart, int nu, const TensorContainer<FSSparseTensor> &fcont,
const UGSContainer &hh, Journal &jr);
template <Storage t>
template<Storage t>
void performStep(int order);
const FGSContainer &
getFoldDers() const
@ -459,9 +459,9 @@ public:
return _ug;
}
protected:
template <Storage t>
template<Storage t>
std::unique_ptr<typename ctraits<t>::Ttensor> faaDiBrunoZ(const Symmetry &sym) const;
template <Storage t>
template<Storage t>
std::unique_ptr<typename ctraits<t>::Ttensor> faaDiBrunoG(const Symmetry &sym) const;
// Convenience access methods
@ -497,7 +497,7 @@ protected:
/* This calculates a derivative of f(G(y,u,σ),g(y,u,σ),y,u) of a given
symmetry. */
template <Storage t>
template<Storage t>
std::unique_ptr<typename ctraits<t>::Ttensor>
KOrderStoch::faaDiBrunoZ(const Symmetry &sym) const
{
@ -512,7 +512,7 @@ KOrderStoch::faaDiBrunoZ(const Symmetry &sym) const
/* This calculates a derivative of G(y,u,σ)=h(g*(y,u,σ),σ) of a given
symmetry. */
template <Storage t>
template<Storage t>
std::unique_ptr<typename ctraits<t>::Ttensor>
KOrderStoch::faaDiBrunoG(const Symmetry &sym) const
{
@ -535,7 +535,7 @@ KOrderStoch::faaDiBrunoG(const Symmetry &sym) const
as g=-matA¹·RHS. Finally we have to update G by calling
Gstack<t>().multAndAdd(1, h<t>(), *G_sym_ptr). */
template <Storage t>
template<Storage t>
void
KOrderStoch::performStep(int order)
{

View File

@ -116,141 +116,147 @@ SparseGenerator::fillContainer(TensorContainer<FSSparseTensor> &c,
}
}
const double vdata [] = { // 3x3
0.1307870268, 0.1241940078, 0.1356703123,
0.1241940078, 0.1986920419, 0.2010160581,
0.1356703123, 0.2010160581, 0.2160336975
};
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 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 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 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 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
};
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
{

View File

@ -158,7 +158,7 @@ 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)
for (const auto &it : mm)
{
const string &oldname = it.first;
const AtomSubstitutions::Tshiftnameset &sset = it.second;
@ -175,7 +175,7 @@ AtomAssignings::apply_subst(const AtomSubstitutions::Toldnamemap &mm)
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)
for (const auto &itt : sset)
{
const string &newname = itt.first;
left_names.insert(newname);

View File

@ -62,7 +62,7 @@ AtomSubstitutions::substitutions_finished(VarOrdering::ord_type ot)
// add all new names derived from the old name
auto it = old2new.find(oname);
if (it != old2new.end())
for (const auto & itt : it->second)
for (const auto &itt : it->second)
na_ext.push_back(itt.first);
}
@ -77,7 +77,7 @@ AtomSubstitutions::get_new4old(const string &oldname, int tshift) const
if (it != old2new.end())
{
const Tshiftnameset &sset = it->second;
for (const auto & itt : sset)
for (const auto &itt : sset)
if (itt.second == -tshift)
return itt.first;
}
@ -88,12 +88,12 @@ void
AtomSubstitutions::print() const
{
std::cout << u8"Atom Substitutions:\nOld ⇒ New:\n";
for (const auto & it : old2new)
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)
for (const auto &it : new2old)
std::cout << " [" << it.second.first << ", " << it.second.second << "] ⇐ " << it.first << '\n';
}

View File

@ -55,7 +55,7 @@ Constants::import_constants(const Constants &c, OperationTree &otree, Tintintmap
void
Constants::setValues(EvalTree &et) const
{
for (const auto & it : cmap)
for (const auto &it : cmap)
et.set_nulary(it.first, it.second);
}
@ -250,7 +250,7 @@ vector<int>
DynamicAtoms::variables() const
{
vector<int> res;
for (const auto & var : vars)
for (const auto &var : vars)
{
const Tlagmap &lmap = var.second;
for (auto itt : lmap)
@ -372,7 +372,7 @@ DynamicAtoms::print() const
std::cout << "constants:\n";
Constants::print();
std::cout << "variables:\n";
for (const auto & var : vars)
for (const auto &var : vars)
{
const Tlagmap &lmap = var.second;
for (auto itt : lmap)
@ -501,7 +501,7 @@ VarOrdering::do_general(ord_type ordering)
// make der_atoms and positions
int off = 0;
for (auto & ord : ords)
for (auto &ord : ords)
for (unsigned int j = 0; j < ord->size(); j++, off++)
if ((*ord)[j] != -1)
{

View File

@ -347,7 +347,7 @@ namespace ogp
public:
/** This is an enum type for an ordering type implemented by
* do_general. */
enum ord_type {pbspbfbf, bfspbfpb};
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

View File

@ -339,20 +339,20 @@ namespace ogp
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;
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;
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;
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;
const vector<int> &outer2y_exo() const;
/** Return the endo_atoms_map. */
const vector<int>&get_endo_atoms_map() const;
const vector<int> &get_endo_atoms_map() const;
/** Return the exo_atoms_map. */
const vector<int>&get_exo_atoms_map() const;
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. */

View File

@ -271,7 +271,7 @@ FormulaDerivatives::derivative(const FoldMultiIndex &mi) const
void
FormulaDerivatives::print(const OperationTree &otree) const
{
for (const auto & it : ind2der)
for (const auto &it : ind2der)
{
std::cout << "derivative ";
it.first.print();

View File

@ -256,7 +256,7 @@ namespace ogp
* different order). */
void differentiate(int max_order);
/** Return i-th formula derivatives. */
const FormulaDerivatives&derivatives(int i) const;
const FormulaDerivatives &derivatives(int i) const;
/** This returns a maximum index of zero derivative formulas
* including all nulary terms. This is a mimumum length of the

View File

@ -29,12 +29,12 @@ namespace ogp
};
// set current off to the first off and add all lengths
#define YYLLOC_DEFAULT(Current, Rhs, N) \
{ \
(Current).off = (Rhs)[1].off; \
(Current).ll = 0; \
for (int i = 1; i <= N; i++) \
(Current).ll += (Rhs)[i].ll; \
#define YYLLOC_DEFAULT(Current, Rhs, N) \
{ \
(Current).off = (Rhs)[1].off; \
(Current).ll = 0; \
for (int i = 1; i <= N; i++) \
(Current).ll += (Rhs)[i].ll; \
}
#define SET_LLOC(prefix) (prefix ## lloc.off += prefix ## lloc.ll, prefix ## lloc.ll = prefix ## leng)

View File

@ -100,12 +100,12 @@ MatrixParser::end() const
}
MPIterator::MPIterator(const MatrixParser &mp)
: p(&mp), i(0), r(mp.find_first_non_empty_row())
: p(&mp), i(0), r(mp.find_first_non_empty_row())
{
}
MPIterator::MPIterator(const MatrixParser &mp, const char *dummy)
: p(&mp), i(mp.data.size()), r(mp.row_lengths.size())
: p(&mp), i(mp.data.size()), r(mp.row_lengths.size())
{
}

View File

@ -96,7 +96,7 @@ StaticFineAtoms::parsing_finished()
// go through all endo and exo insert tree indices, ignore names
// whose tree index is -1 (those which are not referenced)
for (auto & endovar : endovars)
for (auto &endovar : endovars)
{
int t = index(endovar);
if (t != -1)
@ -105,7 +105,7 @@ StaticFineAtoms::parsing_finished()
der_atoms.push_back(t);
}
}
for (auto & exovar : exovars)
for (auto &exovar : exovars)
{
int t = index(exovar);
if (t != -1)

View File

@ -287,7 +287,7 @@ OperationTree::add_derivative(int t, int v)
int tmp1 = add_derivative(terms[t].getOp1(), v);
int tmp2 = add_derivative(terms[t].getOp2(), v);
int res1 = add_binary(code_t::TIMES, terms[t].getOp1(), tmp2);
int res2 = add_binary(code_t::TIMES, tmp1, terms[t].getOp2());
int res2 = add_binary(code_t::TIMES, tmp1, terms[t].getOp2());
res = add_binary(code_t::PLUS, res1, res2);
break;
}
@ -477,7 +477,7 @@ OperationTree::select_terms_inv(int t, const opselector &sel, unordered_set<int>
void
OperationTree::forget_derivative_maps()
{
for (auto & derivative : derivatives)
for (auto &derivative : derivatives)
derivative.clear();
}

View File

@ -42,8 +42,8 @@ namespace ogp
* codes, he should update the code of #OperationTree::add_unary,
* #OperationTree::add_binary, and of course
* #OperationTree::add_derivative. */
enum class code_t {NONE, UMINUS, LOG, EXP, SIN, COS, TAN, SQRT, ERF,
ERFC, PLUS, MINUS, TIMES, DIVIDE, POWER};
enum class code_t { NONE, UMINUS, LOG, EXP, SIN, COS, TAN, SQRT, ERF,
ERFC, PLUS, MINUS, TIMES, DIVIDE, POWER };
/** Class representing a nulary, unary, or binary operation. */
class Operation
@ -64,7 +64,7 @@ namespace ogp
}
/** Constructs a unary operation. */
Operation(code_t cd, int oper1)
: code(cd), op1(oper1)
: code(cd), op1(oper1)
{
}
/** Constructs a nulary operation. */
@ -218,7 +218,7 @@ namespace ogp
* 2/pi. These will be always first four terms having indices
* zero, one and two, three. If adding anything to this
* enumeration, make sure num_constants remains the last one.*/
enum {zero, one, nan, two_over_pi, num_constants};
enum { zero, one, nan, two_over_pi, num_constants };
/** The unique constructor which initializes the object to
* contain only zero, one and nan and two_over_pi.*/

View File

@ -129,7 +129,7 @@ public:
{
return std::make_unique<Dynare>(*this);
}
~Dynare() override = default;
int
nstat() const override

View File

@ -186,7 +186,7 @@ DynareAtomValues::setValues(ogp::EvalTree &et) const
for (auto it : lmap)
{
int ll = it.first;
if (ll == 0) // this is always true because of checks
if (ll == 0) // this is always true because of checks
{
int t = it.second;
int i = atoms.outer2y_exo()[outer_i];

View File

@ -64,7 +64,7 @@ namespace ogdyn
class DynareDynamicAtoms : public ogp::SAtoms, public ogp::NularyStringConvertor
{
public:
enum class atype {endovar, exovar, param};
enum class atype { endovar, exovar, param };
protected:
using Tatypemap = map<string, atype>;
/* The map assigining a type to each name. */

View File

@ -29,8 +29,8 @@ class DynareException
std::string mes;
public:
DynareException(const std::string &m, const std::string &fname, int line, int col)
: mes{"Parse error at " + fname + ", line " + std::to_string(line) + ", column " +
std::to_string(col) + ": " + m}
: mes{"Parse error at " + fname + ", line " + std::to_string(line) + ", column "
+ std::to_string(col) + ": " + m}
{
}
DynareException(const std::string &fname, int line, const std::string &m)

View File

@ -45,7 +45,7 @@ ParsedMatrix::ParsedMatrix(const ogp::MatrixParser &mp)
}
DynareModel::DynareModel()
: atoms(), eqs(atoms)
: atoms(), eqs(atoms)
{
}

View File

@ -404,8 +404,8 @@ namespace ogdyn
in the subclass. */
void write_der0(std::ostream &os);
/* This writes the evaluation of the first order derivative of the system.
It calls pure virtual methods for writing a preamble, assignment, and
assignemnt of the resulting objects. */
It calls pure virtual methods for writing a preamble, assignment, and
assignemnt of the resulting objects. */
void write_der1(std::ostream &os);
protected:
virtual void write_der0_preamble(std::ostream &os) const = 0;

View File

@ -51,39 +51,40 @@ DynareParams::DynareParams(int argc, char **argv)
modname = argv[argc-1];
argc--;
struct option const opts [] = {
{"periods", required_argument, nullptr, static_cast<int>(opt::per)},
{"per", required_argument, nullptr, static_cast<int>(opt::per)},
{"burn", required_argument, nullptr, static_cast<int>(opt::burn)},
{"simulations", required_argument, nullptr, static_cast<int>(opt::sim)},
{"sim", required_argument, nullptr, static_cast<int>(opt::sim)},
{"rtperiods", required_argument, nullptr, static_cast<int>(opt::rtper)},
{"rtper", required_argument, nullptr, static_cast<int>(opt::rtper)},
{"rtsimulations", required_argument, nullptr, static_cast<int>(opt::rtsim)},
{"rtsim", required_argument, nullptr, static_cast<int>(opt::rtsim)},
{"condperiods", required_argument, nullptr, static_cast<int>(opt::condper)},
{"condper", required_argument, nullptr, static_cast<int>(opt::condper)},
{"condsimulations", required_argument, nullptr, static_cast<int>(opt::condsim)},
{"condsim", required_argument, nullptr, static_cast<int>(opt::condsim)},
{"prefix", required_argument, nullptr, static_cast<int>(opt::prefix)},
{"threads", required_argument, nullptr, static_cast<int>(opt::threads)},
{"steps", required_argument, nullptr, static_cast<int>(opt::steps)},
{"seed", required_argument, nullptr, static_cast<int>(opt::seed)},
{"order", required_argument, nullptr, static_cast<int>(opt::order)},
{"ss-tol", required_argument, nullptr, static_cast<int>(opt::ss_tol)},
{"check", required_argument, nullptr, static_cast<int>(opt::check)},
{"check-scale", required_argument, nullptr, static_cast<int>(opt::check_scale)},
{"check-evals", required_argument, nullptr, static_cast<int>(opt::check_evals)},
{"check-num", required_argument, nullptr, static_cast<int>(opt::check_num)},
{"qz-criterium", required_argument, nullptr, static_cast<int>(opt::qz_criterium)},
{"no-irfs", no_argument, nullptr, static_cast<int>(opt::noirfs)},
{"irfs", no_argument, nullptr, static_cast<int>(opt::irfs)},
{"centralize", no_argument, nullptr, static_cast<int>(opt::centralize)},
{"no-centralize", no_argument, nullptr, static_cast<int>(opt::no_centralize)},
{"help", no_argument, nullptr, static_cast<int>(opt::help)},
{"version", no_argument, nullptr, static_cast<int>(opt::version)},
{nullptr, 0, nullptr, 0}
};
struct option const opts[] =
{
{"periods", required_argument, nullptr, static_cast<int>(opt::per)},
{"per", required_argument, nullptr, static_cast<int>(opt::per)},
{"burn", required_argument, nullptr, static_cast<int>(opt::burn)},
{"simulations", required_argument, nullptr, static_cast<int>(opt::sim)},
{"sim", required_argument, nullptr, static_cast<int>(opt::sim)},
{"rtperiods", required_argument, nullptr, static_cast<int>(opt::rtper)},
{"rtper", required_argument, nullptr, static_cast<int>(opt::rtper)},
{"rtsimulations", required_argument, nullptr, static_cast<int>(opt::rtsim)},
{"rtsim", required_argument, nullptr, static_cast<int>(opt::rtsim)},
{"condperiods", required_argument, nullptr, static_cast<int>(opt::condper)},
{"condper", required_argument, nullptr, static_cast<int>(opt::condper)},
{"condsimulations", required_argument, nullptr, static_cast<int>(opt::condsim)},
{"condsim", required_argument, nullptr, static_cast<int>(opt::condsim)},
{"prefix", required_argument, nullptr, static_cast<int>(opt::prefix)},
{"threads", required_argument, nullptr, static_cast<int>(opt::threads)},
{"steps", required_argument, nullptr, static_cast<int>(opt::steps)},
{"seed", required_argument, nullptr, static_cast<int>(opt::seed)},
{"order", required_argument, nullptr, static_cast<int>(opt::order)},
{"ss-tol", required_argument, nullptr, static_cast<int>(opt::ss_tol)},
{"check", required_argument, nullptr, static_cast<int>(opt::check)},
{"check-scale", required_argument, nullptr, static_cast<int>(opt::check_scale)},
{"check-evals", required_argument, nullptr, static_cast<int>(opt::check_evals)},
{"check-num", required_argument, nullptr, static_cast<int>(opt::check_num)},
{"qz-criterium", required_argument, nullptr, static_cast<int>(opt::qz_criterium)},
{"no-irfs", no_argument, nullptr, static_cast<int>(opt::noirfs)},
{"irfs", no_argument, nullptr, static_cast<int>(opt::irfs)},
{"centralize", no_argument, nullptr, static_cast<int>(opt::centralize)},
{"no-centralize", no_argument, nullptr, static_cast<int>(opt::no_centralize)},
{"help", no_argument, nullptr, static_cast<int>(opt::help)},
{"version", no_argument, nullptr, static_cast<int>(opt::version)},
{nullptr, 0, nullptr, 0}
};
int ret;
int index;

View File

@ -89,11 +89,11 @@ struct DynareParams
return 10*check_num;
}
private:
enum class opt {per, burn, sim, rtper, rtsim, condper, condsim,
prefix, threads,
steps, seed, order, ss_tol, check,
check_evals, check_scale, check_num, noirfs, irfs,
help, version, centralize, no_centralize, qz_criterium};
enum class opt { per, burn, sim, rtper, rtsim, condper, condsim,
prefix, threads,
steps, seed, order, ss_tol, check,
check_evals, check_scale, check_num, noirfs, irfs,
help, version, centralize, no_centralize, qz_criterium };
void processCheckFlags(const std::string &flags);
/* This gathers strings from argv[optind] and on not starting with '-' to the
irf_list. It stops one item before the end, since this is the model

View File

@ -70,7 +70,7 @@ ForwSubstBuilder::substitute_for_term(int t, int i, int j)
// first make lagsubst be substitution setting f(x(+4)) to f(x(+1))
// this is lag = -3 (1-mlead)
map<int, int> lagsubst;
unordered_set<int> nult = model.eqs.nulary_of_term(t);// make copy of nult!
unordered_set<int> nult = model.eqs.nulary_of_term(t); // make copy of nult!
model.variable_shift_map(nult, 1-mlead, lagsubst);
int lagt = model.eqs.add_substitution(t, lagsubst);

View File

@ -87,7 +87,6 @@ PlannerBuilder::PlannerBuilder(const PlannerBuilder &pb, ogdyn::DynareModel &m)
diff_b_static(pb.diff_b_static),
diff_f_static(pb.diff_f_static),
aux_map(), static_aux_map()
{
fill_yset(m.atoms.get_name_storage(), pb.yset);
fill_aux_map(m.atoms.get_name_storage(), pb.aux_map, pb.static_aux_map);
@ -244,7 +243,7 @@ PlannerBuilder::make_static_version()
for (int yi = 0; yi < diff_b.nrows(); yi++)
diff_b_static(yi, ll-minlag)
= static_tree.add_substitution(diff_b(yi, ll-minlag),
tmap, model.eqs.getTree());
tmap, model.eqs.getTree());
// go through diff_f and fill diff_f_static
for (int ll = minlag; ll <= maxlead; ll++)
@ -384,7 +383,7 @@ MultInitSS::MultInitSS(const PlannerBuilder &pb, const Vector &pvals, Vector &yy
if (it != old2new.end())
{
const ogp::AtomSubstitutions::Tshiftnameset &sset = it->second;
for (const auto & itt : sset)
for (const auto &itt : sset)
{
const std::string &newname = itt.first;
int iouter = builder.model.atoms.name2outer_endo(newname);

View File

@ -65,7 +65,7 @@ BlockDiagonal::setZerosToRU(diag_iter edge)
of the right-most non-zero element of i-th row from the left, and
col_len[j] is distance of top-most non-zero element of j-th column
to the top. (First element has distance 1).
*/
*/
void
BlockDiagonal::setZeroBlockEdge(diag_iter edge)
{

View File

@ -42,7 +42,9 @@ class TransposedMatrix
private:
T &orig;
public:
TransposedMatrix(T &orig_arg) : orig{orig_arg} {};
TransposedMatrix(T &orig_arg) : orig{orig_arg}
{
};
};
// Syntactic sugar for representing a transposed matrix

View File

@ -31,7 +31,7 @@ IterativeSylvester::solve(SylvParams &pars, KronVector &x) const
auto kpow = matrixK->clone();
auto fpow = matrixF->clone();
while (steps < max_steps &&norm > max_norm)
while (steps < max_steps && norm > max_norm)
{
kpow->multRight(SqSylvMatrix(*kpow)); // be careful to make copy
fpow->multRight(SqSylvMatrix(*fpow)); // also here

View File

@ -148,7 +148,7 @@ void
Diagonal::changeBase(double *p)
{
int d_size = getSize();
for (auto & it : *this)
for (auto &it : *this)
{
const DiagonalBlock &b = it;
int jbar = b.getIndex();
@ -179,7 +179,7 @@ Diagonal::getEigenValues(Vector &eig) const
<< ", should be=" << 2*d_size << '.' << std::endl;
throw SYLV_MES_EXCEPTION(mes.str());
}
for (const auto & b : *this)
for (const auto &b : *this)
{
int ind = b.getIndex();
eig[2*ind] = *(b.getAlpha());
@ -303,7 +303,7 @@ Diagonal::print() const
auto ff = std::cout.flags();
std::cout << "Num real: " << getNumReal() << ", num complex: " << getNumComplex() << std::endl
<< std::fixed;
for (const auto & it : *this)
for (const auto &it : *this)
if (it.isReal())
std::cout << "real: jbar=" << it.getIndex() << ", alpha=" << *(it.getAlpha()) << std::endl;
else

View File

@ -222,7 +222,7 @@ private:
static bool isZero(double p);
};
template <class _TRef, class _TPtr>
template<class _TRef, class _TPtr>
struct _matrix_iter
{
using _Self = _matrix_iter<_TRef, _TPtr>;
@ -260,7 +260,7 @@ public:
virtual _Self &operator++() = 0;
};
template <class _TRef, class _TPtr>
template<class _TRef, class _TPtr>
class _column_iter : public _matrix_iter<_TRef, _TPtr>
{
using _Tparent = _matrix_iter<_TRef, _TPtr>;
@ -293,7 +293,7 @@ public:
}
};
template <class _TRef, class _TPtr>
template<class _TRef, class _TPtr>
class _row_iter : public _matrix_iter<_TRef, _TPtr>
{
using _Tparent = _matrix_iter<_TRef, _TPtr>;
@ -358,7 +358,7 @@ public:
explicit QuasiTriangular(const SchurDecomp &decomp);
explicit QuasiTriangular(const SchurDecompZero &decomp);
QuasiTriangular(const QuasiTriangular &t);
~QuasiTriangular() override = default;
const Diagonal &
getDiagonal() const

View File

@ -31,10 +31,10 @@ SchurDecomp::SchurDecomp(const SqSylvMatrix &m)
SqSylvMatrix auxt(m);
lapack_int lda = auxt.getLD(), ldvs = q.getLD();
lapack_int sdim;
auto wr = std::make_unique<double []>(rows);
auto wi = std::make_unique<double []>(rows);
auto wr = std::make_unique<double[]>(rows);
auto wi = std::make_unique<double[]>(rows);
lapack_int lwork = 6*rows;
auto work = std::make_unique<double []>(lwork);
auto work = std::make_unique<double[]>(lwork);
lapack_int info;
dgees("V", "N", nullptr, &rows, auxt.base(), &lda, &sdim,
wr.get(), wi.get(), q.base(), &ldvs,

View File

@ -30,7 +30,7 @@
neighbour is bubbled also in front. The method returns a new
position to, where the original block pointed by to happens to
appear at the end. from must be greater than to.
*/
*/
SchurDecompEig::diag_iter
SchurDecompEig::bubbleEigen(diag_iter from, diag_iter to)
{
@ -65,7 +65,7 @@ SchurDecompEig::bubbleEigen(diag_iter from, diag_iter to)
resolve eigenvalues from itadd to it, before the it can be
resolved.
The success is signaled by returned true.
*/
*/
bool
SchurDecompEig::tryToSwap(diag_iter &it, diag_iter &itadd)
{

View File

@ -52,7 +52,7 @@ SimilarityDecomp::getXDim(diag_iter start, diag_iter end,
/* Find solution of X for diagonal block given by start(incl.) and
end(excl.). If the solution cannot be found, or it is greater than
norm, X is not changed and flase is returned.
*/
*/
bool
SimilarityDecomp::solveX(diag_iter start, diag_iter end,
GeneralMatrix &X, double norm) const

View File

@ -133,11 +133,11 @@ mxArray *
SylvParams::IntParamItem::createMatlabArray() const
{
mxArray *res = mxCreateNumericMatrix(1, 1, mxINT32_CLASS, mxREAL);
#if MX_HAS_INTERLEAVED_COMPLEX
# if MX_HAS_INTERLEAVED_COMPLEX
*mxGetInt32s(res) = value;
#else
# else
*static_cast<int *>(mxGetData(res)) = value;
#endif
# endif
return res;
}

View File

@ -28,9 +28,9 @@
# include <dynmex.h>
#endif
enum class status {def, changed, undef};
enum class status { def, changed, undef };
template <class _Type>
template<class _Type>
struct ParamItem
{
protected:
@ -81,7 +81,7 @@ public:
class SylvParams
{
public:
enum class solve_method {iter, recurse};
enum class solve_method { iter, recurse };
protected:
class DoubleParamItem : public ParamItem<double>

View File

@ -37,7 +37,7 @@ public:
TriangularSylvester(const QuasiTriangular &k, const QuasiTriangular &f);
TriangularSylvester(const SchurDecompZero &kdecomp, const SchurDecomp &fdecomp);
TriangularSylvester(const SchurDecompZero &kdecomp, const SimilarityDecomp &fdecomp);
~TriangularSylvester() override = default;
void print() const;
void solve(SylvParams &pars, KronVector &d) const override;

View File

@ -398,7 +398,7 @@ TestRunnable::gen_sylv(const std::string &aname, const std::string &bname, const
if (m != mmc.row() || m != mmc.col()
|| n != mma.row() || n != mma.col()
|| n != mmb.row() || n < mmb.col()
|| n != mmb.row() || n < mmb.col()
|| n != mmd.row() || power(m, order) != mmd.col())
{
std::cout << " Incompatible sizes for gen_sylv.\n";
@ -447,10 +447,10 @@ TestRunnable::eig_bubble(const std::string &aname, int from, int to)
double normInf = check.getNormInf();
double onorm1 = orig.getNorm1();
double onormInf = orig.getNormInf();
std:: cout << "\tabs. error1 = " << norm1 << std::endl
<< u8"\tabs. error∞ = " << normInf << std::endl
<< "\trel. error1 = " << norm1/onorm1 << std::endl
<< u8"\trel. error∞ = " << normInf/onormInf << std::endl;
std::cout << "\tabs. error1 = " << norm1 << std::endl
<< u8"\tabs. error∞ = " << normInf << std::endl
<< "\trel. error1 = " << norm1/onorm1 << std::endl
<< u8"\trel. error∞ = " << normInf/onormInf << std::endl;
return (norm1 < eps_norm*onorm1 && normInf < eps_norm*onormInf);
}

View File

@ -90,7 +90,7 @@ public:
method, which returns a type for a given stack as the type of the
corresponding (old) stack of the former stack container. */
template <class _Ttype>
template<class _Ttype>
class FineContainer : public SizeRefinement, public StackContainer<_Ttype>
{
protected:
@ -117,8 +117,8 @@ public:
FineContainer(const _Stype &sc, int max)
: SizeRefinement(sc.getStackSizes(), sc.numConts(), max),
StackContainer<_Ttype>(numRefinements(), getNC()),
ref_conts(getNC()),
stack_cont(sc)
ref_conts(getNC()),
stack_cont(sc)
{
for (int i = 0; i < numRefinements(); i++)
_Stype::stack_sizes[i] = getRefSize(i);

View File

@ -59,7 +59,7 @@ FFSTensor::FFSTensor(const FFSTensor &t, const ConstVector &x)
symmetry. Let n be a number of variables and d the
n+d-1
dimension dim. Then the number of indices is d .
*/
*/
int
FFSTensor::calcMaxOffset(int nvar, int d)
@ -79,7 +79,7 @@ FFSTensor::FFSTensor(const FSSparseTensor &t)
nv(t.nvar())
{
zeros();
for (const auto & it : t.getMap())
for (const auto &it : t.getMap())
{
index ind(*this, it.first);
get(it.second.first, *ind) = it.second.second;

View File

@ -215,7 +215,7 @@ public:
UGSTensor(UGSTensor &&) = default;
UGSTensor(int first_row, int num, UGSTensor &t)
: UTensor(first_row, num, t), tdims(t.tdims)
: UTensor(first_row, num, t), tdims(t.tdims)
{
}

View File

@ -284,7 +284,7 @@ IntSequence::print() const
we want to compute the number of permutations of the word yyyyuuvvv.
9
This is equal to the multinomial coefficient 4,2,3.
*/
*/
int
IntSequence::noverseq()
{

View File

@ -37,7 +37,7 @@ UNormalMoments::UNormalMoments(int maxdim, const TwoDMatrix &v)
equal to 2n. See the header file for proof and details.
Here we sequentially construct the Kronecker power v and apply F.
*/
*/
void
UNormalMoments::generateMoments(int maxdim, const TwoDMatrix &v)
{

View File

@ -43,7 +43,7 @@
Also, we need PermutationSet class which contains all permutations
of an n-element set, and a bundle of permutations PermutationBundle
which contains all permutation sets up to a given number.
*/
*/
#ifndef PERMUTATION_H
#define PERMUTATION_H
@ -63,7 +63,7 @@
indices which yield identity. Also we have a constructor calculating
map, which corresponds to permutation in sort. This is, we want
(sorted s)m = s.
*/
*/
class Permutation
{

View File

@ -86,7 +86,8 @@
class PerTensorDimens : public TensorDimens
{
private:
static IntSequence sortIntSequence(IntSequence s)
static IntSequence
sortIntSequence(IntSequence s)
{
s.sort();
return s;
@ -221,7 +222,7 @@ public:
void addTo(FGSTensor &out) const;
void addTo(UGSTensor &out) const;
enum class fill_method {first, second};
enum class fill_method { first, second };
static fill_method decideFillMethod(const FSSparseTensor &t);
private:
int tailIdentitySize() const;

View File

@ -43,7 +43,7 @@ USubTensor::USubTensor(const TensorDimens &bdims,
"Tensor has not full symmetry in USubTensor()");
const EquivalenceSet &eset = TLStatic::getEquiv(bdims.dimen());
zeros();
for (const auto & it : eset)
for (const auto &it : eset)
{
if (it.numClasses() == hdims.dimen())
{

View File

@ -119,7 +119,7 @@ WorkerFoldMAASparse1::operator()(std::mutex &mut)
const Permutation &per = pset.get(iper);
IntSequence percoor(coor.size());
per.apply(coor, percoor);
for (const auto & it : eset)
for (const auto &it : eset)
{
if (it.numClasses() == t.dimen())
{
@ -163,7 +163,7 @@ FoldedStackContainer::multAndAddSparse2(const FSSparseTensor &t,
sthread::detach_thread_group gr;
FFSTensor dummy_f(0, numStacks(), t.dimen());
for (Tensor::index fi = dummy_f.begin(); fi != dummy_f.end(); ++fi)
gr.insert(std::make_unique<WorkerFoldMAASparse2>(*this, t, out, fi.getCoor()));
gr.insert(std::make_unique<WorkerFoldMAASparse2>(*this, t, out, fi.getCoor()));
gr.run();
}
@ -231,7 +231,7 @@ FoldedStackContainer::multAndAddSparse3(const FSSparseTensor &t,
Vector outcol{out.getCol(*run)};
FRSingleTensor sumcol(t.nvar(), t.dimen());
sumcol.zeros();
for (const auto & it : eset)
for (const auto &it : eset)
if (it.numClasses() == t.dimen())
{
StackProduct<FGSTensor> sp(*this, it, out.getSym());
@ -485,7 +485,7 @@ WorkerUnfoldMAASparse1::operator()(std::mutex &mut)
const Permutation &per = pset.get(iper);
IntSequence percoor(coor.size());
per.apply(coor, percoor);
for (const auto & it : eset)
for (const auto &it : eset)
if (it.numClasses() == t.dimen())
{
StackProduct<UGSTensor> sp(cont, it, out.getSym());
@ -611,7 +611,7 @@ UnfoldedStackContainer::multAndAddStacks(const IntSequence &fi,
{
Permutation sort_per(ui.getCoor());
sort_per.inverse();
for (const auto & it : eset)
for (const auto &it : eset)
if (it.numClasses() == g.dimen())
{
StackProduct<UGSTensor> sp(*this, it, sort_per, out.getSym());

View File

@ -107,12 +107,12 @@
from zero types, or unit matrices are deleted. See kron_prod.hh for an
explanation. */
template <class _Ttype>
template<class _Ttype>
class StackContainerInterface
{
public:
using _Ctype = TensorContainer<_Ttype>;
enum class itype { matrix, unit, zero};
enum class itype { matrix, unit, zero };
public:
StackContainerInterface() = default;
virtual ~StackContainerInterface() = default;
@ -145,7 +145,7 @@ public:
It does not own its tensors. */
template <class _Ttype>
template<class _Ttype>
class StackContainer : virtual public StackContainerInterface<_Ttype>
{
public:
@ -347,7 +347,7 @@ protected:
on four variables, the third being u a dummy and always returning zero if
dimension of u is positive. */
template <class _Ttype>
template<class _Ttype>
class ZContainer : public StackContainer<_Ttype>
{
public:
@ -429,7 +429,7 @@ public:
third one is dummy, and always returns zero. The first stack
corresponds to a container of g*. */
template <class _Ttype>
template<class _Ttype>
class GContainer : public StackContainer<_Ttype>
{
public:
@ -503,7 +503,7 @@ public:
equivalence can have permuted classes by some given
permutation. Nothing else is interesting. */
template <class _Ttype>
template<class _Ttype>
class StackProduct
{
public:
@ -614,7 +614,7 @@ public:
/* Here we only inherit from Kronecker product KronProdAllOptim, only to
allow for a constructor constructing from StackProduct. */
template <class _Ttype>
template<class _Ttype>
class KronProdStack : public KronProdAllOptim
{
public:

View File

@ -116,7 +116,7 @@ symiterator::operator++()
InducedSymmetries::InducedSymmetries(const Equivalence &e, const Symmetry &s)
{
for (const auto & i : e)
for (const auto &i : e)
emplace_back(s, i);
}

View File

@ -136,7 +136,8 @@ public:
insert(std::make_unique<_Ttype>(first_row, num, *(it.second)));
}
TensorContainer<_Ttype> &operator=(const TensorContainer<_Ttype> &c)
TensorContainer<_Ttype> &
operator=(const TensorContainer<_Ttype> &c)
{
n = c.n;
m.clear();
@ -182,7 +183,7 @@ public:
TL_RAISE_IF(check(t->getSym()),
"Tensor already in container in TensorContainer::insert");
if (!t->isFinite())
throw TLException(__FILE__, __LINE__, "NaN or Inf asserted in TensorContainer::insert");
throw TLException(__FILE__, __LINE__, "NaN or Inf asserted in TensorContainer::insert");
m.emplace(t->getSym(), std::move(t));
}

View File

@ -101,7 +101,8 @@ public:
struct dummy { using type = T; };
template<class T>
const T &getNext()
const T &
getNext()
{
return getNext(dummy<T>());
}
@ -128,7 +129,7 @@ private:
See documentation to TensorPolynomial::derivative() and
TensorPolynomial::evalPartially() for details. */
template <class _Ttype, class _TGStype, class _Stype>
template<class _Ttype, class _TGStype, class _Stype>
class TensorPolynomial : public TensorContainer<_Ttype>
{
int nr;
@ -138,18 +139,18 @@ class TensorPolynomial : public TensorContainer<_Ttype>
public:
TensorPolynomial(int rows, int vars)
: TensorContainer<_Ttype>(1),
nr(rows), nv(vars), maxdim(0)
nr(rows), nv(vars), maxdim(0)
{
}
TensorPolynomial(const TensorPolynomial<_Ttype, _TGStype, _Stype> &tp, int k)
: TensorContainer<_Ttype>(tp),
nr(tp.nr), nv(tp.nv), maxdim(0)
nr(tp.nr), nv(tp.nv), maxdim(0)
{
derivative(k);
}
TensorPolynomial(int first_row, int num, TensorPolynomial<_Ttype, _TGStype, _Stype> &tp)
: TensorContainer<_Ttype>(first_row, num, tp),
nr(num), nv(tp.nv), maxdim(tp.maxdim)
nr(num), nv(tp.nv), maxdim(tp.maxdim)
{
}
@ -176,7 +177,7 @@ public:
TensorPolynomial(const TensorPolynomial<_Ttype, _TGStype, _Stype> &tp, const Vector &xval)
: TensorContainer<_Ttype>(1),
nr(tp.nrows()), nv(tp.nvars() - xval.length()), maxdim(0)
nr(tp.nrows()), nv(tp.nvars() - xval.length()), maxdim(0)
{
TL_RAISE_IF(nvars() < 0,
"Length of xval too big in TensorPolynomial contract constructor");
@ -472,7 +473,7 @@ public:
tensor, with the number of variables equal to the number of variables
of the polynomial plus 1 for 1. */
template <class _Ttype, class _TGStype, class _Stype>
template<class _Ttype, class _TGStype, class _Stype>
class CompactPolynomial : public _Ttype
{
public:

View File

@ -83,7 +83,7 @@
class Tensor : public TwoDMatrix
{
public:
enum class indor {along_row, along_col};
enum class indor { along_row, along_col };
/* The index represents n-tuple α₁…αₙ. Since its movement is dependent on the
underlying tensor (with storage and symmetry), we maintain a reference to

View File

@ -59,7 +59,7 @@
# define TL_DEBUG 0
#endif
#define TL_RAISE(mes) \
#define TL_RAISE(mes) \
throw TLException(__FILE__, __LINE__, mes)
#define TL_RAISE_IF(expr, mes) \

View File

@ -142,7 +142,6 @@ public:
{
}
#endif
~TwoDMatrix() override = default;
TwoDMatrix &operator=(const TwoDMatrix &m) = default;

View File

@ -42,7 +42,7 @@ class Factory
public:
double get();
// This can be used with UGSTensor, FGSTensor
template <class _Ttype>
template<class _Ttype>
std::unique_ptr<_Ttype>
make(int r, const Symmetry &s, const IntSequence &nvs)
{
@ -53,7 +53,7 @@ public:
}
// This can be used with FFSTensor, UFSTensor, FRTensor, URTensor
template <class _Ttype>
template<class _Ttype>
std::unique_ptr<_Ttype>
make(int r, int nv, int dim)
{
@ -63,7 +63,7 @@ public:
return res;
}
template <class _Ttype, class _Ctype>
template<class _Ttype, class _Ctype>
_Ctype
makeCont(int r, const IntSequence &nvs, int maxdim)
{
@ -80,7 +80,7 @@ public:
return res;
}
template <class _Ttype, class _Ptype>
template<class _Ttype, class _Ptype>
_Ptype
makePoly(int r, int nv, int maxdim)
{

View File

@ -56,10 +56,10 @@ protected:
template<class _Ttype>
static bool index_forward(const Symmetry &s, const IntSequence &nvs);
template <class _Ttype>
template<class _Ttype>
static bool index_backward(const Symmetry &s, const IntSequence &nvs);
template <class _Ttype>
template<class _Ttype>
static bool index_offset(const Symmetry &s, const IntSequence &nvs);
static bool fold_unfold(std::unique_ptr<FTensor> folded);
@ -122,7 +122,7 @@ TestRunnable::test() const
/****************************************************/
/* definition of TestRunnable static methods */
/****************************************************/
template <class _Ttype>
template<class _Ttype>
bool
TestRunnable::index_forward(const Symmetry &s, const IntSequence &nvs)
{
@ -154,7 +154,7 @@ TestRunnable::index_forward(const Symmetry &s, const IntSequence &nvs)
return fails == 0;
}
template <class _Ttype>
template<class _Ttype>
bool
TestRunnable::index_backward(const Symmetry &s, const IntSequence &nvs)
{
@ -185,7 +185,7 @@ TestRunnable::index_backward(const Symmetry &s, const IntSequence &nvs)
return fails == 0;
}
template <class _Ttype>
template<class _Ttype>
bool
TestRunnable::index_offset(const Symmetry &s, const IntSequence &nvs)
{

View File

@ -68,7 +68,6 @@ PascalRow::print() const
std::cout << std::endl;
}
namespace PascalTriangle
{
namespace // Anonymous namespace that is a functional equivalent of “private”
@ -127,7 +126,7 @@ namespace PascalTriangle
void
print()
{
for (const auto & i : tr)
for (const auto &i : tr)
i.print();
}
}

View File

@ -49,20 +49,20 @@ namespace sthread
{
counter++;
std::thread th{[&, it] {
// The it variable is captured by value, because otherwise the iterator may move
(*it)->operator()(mut_threads);
std::unique_lock<std::mutex> lk2{mut_cv};
counter--;
/* First notify the thread waiting on the condition variable, then
unlock the mutex. We must do these two operations in that order,
otherwise there is a possibility of having the main process
destroying the condition variable before the thread tries to
notify it (if all other threads terminate at the same time and
bring the counter down to zero).
For that reason, we cannot use std::notify_all_at_thread_exit() */
cv.notify_one();
lk2.unlock();
}};
// The it variable is captured by value, because otherwise the iterator may move
(*it)->operator()(mut_threads);
std::unique_lock<std::mutex> lk2{mut_cv};
counter--;
/* First notify the thread waiting on the condition variable, then
unlock the mutex. We must do these two operations in that order,
otherwise there is a possibility of having the main process
destroying the condition variable before the thread tries to
notify it (if all other threads terminate at the same time and
bring the counter down to zero).
For that reason, we cannot use std::notify_all_at_thread_exit() */
cv.notify_one();
lk2.unlock();
}};
th.detach();
++it;
cv.wait(lk, [&] { return counter < max_parallel_threads; });