stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

k-order DLL: convert model derivatives from Dynare to Dynare++ format at an arbitrary order 

Ref #217
time-shift
Sébastien Villemot 2019-04-12 18:12:21 +02:00
parent 03ac8c8182
commit d7dd7214c7
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
1 changed files with 22 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
mex/sources/k_order_perturbation/k_ord_dynare.cc
View File

@ -116,48 +116,28 @@ KordpDynare::populateDerivativesContainer(int ord)
}
s[0]++;
}
else if (ord == 2)
{
for (int i = 0; i < g.nrows(); i++)
{
int j = static_cast<int>(g.get(i, 0))-1; // hessian indices start with 1
int i1 = static_cast<int>(g.get(i, 1))-1;
if (j < 0 || i1 < 0)
continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix())
int s0 = i1 / nJcols;
int s1 = i1 % nJcols;
s[0] = dynToDynpp[s0];
s[1] = dynToDynpp[s1];
if (s[1] >= s[0])
{
double x = g.get(i, 2);
mdTi->insert(s, j, x);
}
}
}
else if (ord == 3)
{
int nJcols2 = nJcols*nJcols;
for (int i = 0; i < g.nrows(); i++)
{
int j = static_cast<int>(g.get(i, 0))-1;
int i1 = static_cast<int>(g.get(i, 1))-1;
if (j < 0 || i1 < 0)
continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix())
int s0 = i1 / nJcols2;
int i2 = i1 % nJcols2;
int s1 = i2 / nJcols;
int s2 = i2 % nJcols;
s[0] = dynToDynpp[s0];
s[1] = dynToDynpp[s1];
s[2] = dynToDynpp[s2];
if (s.isSorted())
{
double x = g.get(i, 2);
mdTi->insert(s, j, x);
}
}
}
else // ord ≥ 2
for (int i = 0; i < g.nrows(); i++)
{
int j = static_cast<int>(g.get(i, 0))-1;
int i1 = static_cast<int>(g.get(i, 1))-1;
if (j < 0 || i1 < 0)
continue; // Discard empty entries (see comment in DynamicModelAC::unpackSparseMatrix())
for (int k = 0; k < ord; k++)
{
s[k] = dynToDynpp[i1 % nJcols];
i1 /= nJcols;
}
if ((ord == 2 || ord == 3) && !s.isSorted())
continue; // Skip symmetric elements (only needed at order 2 and 3)
else if (ord > 3)
s.sort(); // For higher order, canonicalize the multi-index
double x = g.get(i, 2);
mdTi->insert(s, j, x);
}
md.insert(std::move(mdTi));
}