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

Remove symmetric elements in 3rd derivatives

issue#70
Sébastien Villemot 2019-06-17 15:28:33 +02:00
parent d59f9f75ff
commit 271a579808
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
3 changed files with 5 additions and 103 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

49
src/DynamicModel.cc
View File

@ -2566,7 +2566,7 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
k++;
}
// Output symetric elements, but only at order 2 and 3
// Output symetric elements at order 2
if (i == 2 && vidx[1] != vidx[2])
{
int col_idx_sym = getDynJacobianCol(vidx[2]) * dynJacobianColsNbr + getDynJacobianCol(vidx[1]);
@ -2590,41 +2590,6 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput,
k++;
}
}
if (i == 3)
{
// Use std::next_permutation() to explore all the permutations of the 3 indices
vector<int> idx3{getDynJacobianCol(vidx[1]), getDynJacobianCol(vidx[2]), getDynJacobianCol(vidx[3])};
sort(idx3.begin(), idx3.end());
int k2 = 0; // Keeps the offset of the permutation relative to k
do
{
int col_idx_sym = idx3[0]*hessianColsNbr + idx3[1]*dynJacobianColsNbr + idx3[2];
if (col_idx_sym != col_idx)
if (output_type == ExprNodeOutputType::juliaDynamicModel)
d_output[3] << " @inbounds g3[" << eq + 1 << "," << col_idx_sym + 1 << "] = "
<< "g3[" << eq + 1 << "," << col_idx + 1 << "]" << endl;
else
{
sparseHelper(3, col0_output, k+k2, 0, output_type);
col0_output << "=" << eq + 1 << ";" << endl;
sparseHelper(3, col1_output, k+k2, 1, output_type);
col1_output << "=" << col_idx_sym + 1 << ";" << endl;
sparseHelper(3, col2_output, k+k2, 2, output_type);
col2_output << "=";
sparseHelper(3, col2_output, k-1, 2, output_type);
col2_output << ";" << endl;
k2++;
}
}
while (next_permutation(idx3.begin(), idx3.end()));
if (output_type != ExprNodeOutputType::juliaDynamicModel)
k += k2;
}
}
if (output_type != ExprNodeOutputType::juliaDynamicModel)
d_output[i] << col0_output.str() << col1_output.str() << col2_output.str();
@ -6889,18 +6854,6 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c
int col_idx_sym = getDynJacobianCol(vidx[2]) * dynJacobianColsNbr + getDynJacobianCol(vidx[1]);
d_output[i] << ", " << col_idx_sym + 1;
}
if (i == 3) // Symmetric elements in 3rd derivative
{
vector<int> idx3{getDynJacobianCol(vidx[1]), getDynJacobianCol(vidx[2]), getDynJacobianCol(vidx[3])};
sort(idx3.begin(), idx3.end());
do
{
int col_idx_sym = (idx3[0]*dynJacobianColsNbr + idx3[1])*dynJacobianColsNbr + idx3[2];
if (col_idx_sym != col_idx)
d_output[i] << ", " << col_idx_sym+1;
}
while (next_permutation(idx3.begin(), idx3.end()));
}
if (i > 1)
d_output[i] << "]";

10
src/ModelTree.cc
View File

@ -1309,13 +1309,9 @@ ModelTree::computeDerivatives(int order, const set<int> &vars)
indices.push_back(var);
// At this point, indices of endogenous variables are sorted in non-decreasing order
derivatives[o][indices] = d;
// We output symmetries at order ≤ 3
if (o <= 3)
{
do
NNZDerivatives[o]++;
while (next_permutation(next(indices.begin()), indices.end()));
}
// We output symmetric elements at order = 2
if (o == 2 && indices[1] != indices[2])
NNZDerivatives[o] += 2;
else
NNZDerivatives[o]++;
}

49
src/StaticModel.cc
View File

@ -1556,7 +1556,7 @@ StaticModel::writeStaticModel(const string &basename,
k++;
}
// Output symetric elements, but only at order 2 and 3
// Output symetric elements at order 2
if (i == 2 && vidx[1] != vidx[2])
{
int col_idx_sym = getJacobCol(vidx[2]) * JacobianColsNbr + getJacobCol(vidx[1]);
@ -1580,41 +1580,6 @@ StaticModel::writeStaticModel(const string &basename,
k++;
}
}
if (i == 3)
{
// Use std::next_permutation() to explore all the permutations of the 3 indices
vector<int> idx3{getJacobCol(vidx[1]), getJacobCol(vidx[2]), getJacobCol(vidx[3])};
sort(idx3.begin(), idx3.end());
int k2 = 0; // Keeps the offset of the permutation relative to k
do
{
int col_idx_sym = idx3[0]*hessianColsNbr + idx3[1]*JacobianColsNbr + idx3[2];
if (col_idx_sym != col_idx)
if (output_type == ExprNodeOutputType::juliaStaticModel)
d_output[3] << " @inbounds g3[" << eq + 1 << "," << col_idx_sym + 1 << "] = "
<< "g3[" << eq + 1 << "," << col_idx + 1 << "]" << endl;
else
{
sparseHelper(3, col0_output, k+k2, 0, output_type);
col0_output << "=" << eq + 1 << ";" << endl;
sparseHelper(3, col1_output, k+k2, 1, output_type);
col1_output << "=" << col_idx_sym + 1 << ";" << endl;
sparseHelper(3, col2_output, k+k2, 2, output_type);
col2_output << "=";
sparseHelper(3, col2_output, k-1, 2, output_type);
col2_output << ";" << endl;
k2++;
}
}
while (next_permutation(idx3.begin(), idx3.end()));
if (output_type != ExprNodeOutputType::juliaStaticModel)
k += k2;
}
}
if (output_type != ExprNodeOutputType::juliaStaticModel)
d_output[i] << col0_output.str() << col1_output.str() << col2_output.str();
@ -2869,18 +2834,6 @@ StaticModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) co
int col_idx_sym = getJacobCol(vidx[2]) * symbol_table.endo_nbr() + getJacobCol(vidx[1]);
d_output[i] << ", " << col_idx_sym + 1;
}
if (i == 3) // Symmetric elements in 3rd derivative
{
vector<int> idx3{getJacobCol(vidx[1]), getJacobCol(vidx[2]), getJacobCol(vidx[3])};
sort(idx3.begin(), idx3.end());
do
{
int col_idx_sym = (idx3[0]*symbol_table.endo_nbr() + idx3[1])*symbol_table.endo_nbr() + idx3[2];
if (col_idx_sym != col_idx)
d_output[i] << ", " << col_idx_sym+1;
}
while (next_permutation(idx3.begin(), idx3.end()));
}
if (i > 1)
d_output[i] << "]";