From 271a57980847fb0b7c4d5c0f78e339f21ec847ad Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Mon, 17 Jun 2019 15:28:33 +0200
Subject: [PATCH] Remove symmetric elements in 3rd derivatives
---
src/DynamicModel.cc | 49 +--------------------------------------------
src/ModelTree.cc | 10 +++------
src/StaticModel.cc | 49 +--------------------------------------------
3 files changed, 5 insertions(+), 103 deletions(-)
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 67510c68..2e60f33e 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -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] << "]";
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index c9e800d4..440236eb 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -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]++;
}
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index e5f88010..6babac0a 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -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] << "]";
--
GitLab