diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 2acd9d1160688514804afcc4a358fc49abc44c9c..3c6b1e7d4ae5fc9d4f5cc485666ebae1f23a98b1 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -63,7 +63,7 @@ DynamicModel::copyHelper(const DynamicModel &m)
     {
       block_derivatives_equation_variable_laglead_nodeid_t v;
       for (const auto &it2 : it)
-        v.push_back(make_pair(it2.first, make_pair(it2.second.first, f(it2.second.second))));
+        v.emplace_back(get<0>(it2), get<1>(it2), get<2>(it2), f(get<3>(it2)));
       blocks_derivatives.push_back(v);
     }
 
@@ -258,7 +258,7 @@ DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_n
 void
 DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, const map_idx_t &map_idx) const
 {
-  auto it = first_chain_rule_derivatives.find({ eqr, { varr, lag } });
+  auto it = first_chain_rule_derivatives.find({ eqr, varr, lag });
   if (it != first_chain_rule_derivatives.end())
     (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
   else
@@ -304,18 +304,16 @@ DynamicModel::computeTemporaryTermsOrdered()
                   eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
                 }
             }
-          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+          for (const auto &it : blocks_derivatives[block])
             {
-              expr_t id = it->second.second;
+              expr_t id = get<3>(it);
               id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
             }
-          for (derivative_t::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
-            it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
-          for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
-            it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
-          set<int> temporary_terms_in_use;
-          temporary_terms_in_use.clear();
-          v_temporary_terms_inuse[block] = temporary_terms_in_use;
+          for (const auto &it : derivative_endo[block])
+            it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
+          for (const auto &it : derivative_other_endo[block])
+            it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
+          v_temporary_terms_inuse[block] = {};
         }
     }
   else
@@ -337,15 +335,15 @@ DynamicModel::computeTemporaryTermsOrdered()
                   eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
                 }
             }
-          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+          for (const auto &it : blocks_derivatives[block])
             {
-              expr_t id = it->second.second;
+              expr_t id = get<3>(it);
               id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
             }
-          for (derivative_t::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
-            it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
-          for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
-            it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
+          for (const auto &it : derivative_endo[block])
+            it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
+          for (const auto &it : derivative_other_endo[block])
+            it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
         }
       for (unsigned int block = 0; block < nb_blocks; block++)
         {
@@ -364,19 +362,19 @@ DynamicModel::computeTemporaryTermsOrdered()
                   eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
                 }
             }
-          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+          for (const auto &it : blocks_derivatives[block])
             {
-              expr_t id = it->second.second;
+              expr_t id = get<3>(it);
               id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
             }
-          for (derivative_t::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
-            it->second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-          for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
-            it->second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-          for (derivative_t::const_iterator it = derivative_exo[block].begin(); it != derivative_exo[block].end(); it++)
-            it->second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
-          for (derivative_t::const_iterator it = derivative_exo_det[block].begin(); it != derivative_exo_det[block].end(); it++)
-            it->second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+          for (const auto &it : derivative_endo[block])
+            it.second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+          for (const auto &it : derivative_other_endo[block])
+            it.second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+          for (const auto &it : derivative_exo[block])
+            it.second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+          for (const auto &it : derivative_exo_det[block])
+            it.second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
           v_temporary_terms_inuse[block] = temporary_terms_in_use;
         }
       computeTemporaryTermsMapping();
@@ -436,16 +434,16 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
 
       int prev_lag;
       unsigned int prev_var, count_col, count_col_endo, count_col_exo, count_col_exo_det, count_col_other_endo;
-      map<pair<int, pair<int, int>>, expr_t> tmp_block_endo_derivative;
-      for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
-        tmp_block_endo_derivative[{ it->second.first, { it->first.second, it->first.first } }] = it->second.second;
+      map<tuple<int, int, int>, expr_t> tmp_block_endo_derivative;
+      for (const auto &it : blocks_derivatives[block])
+        tmp_block_endo_derivative[{ get<2>(it), get<1>(it), get<0>(it) }] = get<3>(it);
       prev_var = 999999999;
       prev_lag = -9999999;
       count_col_endo = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_endo_derivative.begin(); it != tmp_block_endo_derivative.end(); it++)
+      for (const auto &it : tmp_block_endo_derivative)
         {
-          int lag = it->first.first;
-          unsigned int var = it->first.second.first;
+          int lag = get<0>(it.first);
+          unsigned int var = get<1>(it.first);
           if (var != prev_var || lag != prev_lag)
             {
               prev_var = var;
@@ -453,16 +451,16 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
               count_col_endo++;
             }
         }
-      map<pair<int, pair<int, int>>, expr_t> tmp_block_exo_derivative;
-      for (auto it = derivative_exo[block].begin(); it != (derivative_exo[block]).end(); it++)
-        tmp_block_exo_derivative[{ it->first.first, { it->first.second.second, it->first.second.first } }] = it->second;
+      map<tuple<int, int, int>, expr_t> tmp_block_exo_derivative;
+      for (const auto &it : derivative_exo[block])
+        tmp_block_exo_derivative[{ get<0>(it.first), get<2>(it.first), get<1>(it.first) }] = it.second;
       prev_var = 999999999;
       prev_lag = -9999999;
       count_col_exo = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_exo_derivative.begin(); it != tmp_block_exo_derivative.end(); it++)
+      for (const auto &it : tmp_block_exo_derivative)
         {
-          int lag = it->first.first;
-          unsigned int var = it->first.second.first;
+          int lag = get<0>(it.first);
+          unsigned int var = get<1>(it.first);
           if (var != prev_var || lag != prev_lag)
             {
               prev_var = var;
@@ -470,16 +468,16 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
               count_col_exo++;
             }
         }
-      map<pair<int, pair<int, int>>, expr_t> tmp_block_exo_det_derivative;
-      for (auto it = derivative_exo_det[block].begin(); it != (derivative_exo_det[block]).end(); it++)
-        tmp_block_exo_det_derivative[{ it->first.first, { it->first.second.second, it->first.second.first } }] = it->second;
+      map<tuple<int, int, int>, expr_t> tmp_block_exo_det_derivative;
+      for (const auto &it : derivative_exo_det[block])
+        tmp_block_exo_det_derivative[{ get<0>(it.first), get<2>(it.first), get<1>(it.first) }] = it.second;
       prev_var = 999999999;
       prev_lag = -9999999;
       count_col_exo_det = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_exo_derivative.begin(); it != tmp_block_exo_derivative.end(); it++)
+      for (const auto &it : tmp_block_exo_derivative)
         {
-          int lag = it->first.first;
-          unsigned int var = it->first.second.first;
+          int lag = get<0>(it.first);
+          unsigned int var = get<1>(it.first);
           if (var != prev_var || lag != prev_lag)
             {
               prev_var = var;
@@ -487,16 +485,16 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
               count_col_exo_det++;
             }
         }
-      map<pair<int, pair<int, int>>, expr_t> tmp_block_other_endo_derivative;
-      for (auto it = derivative_other_endo[block].begin(); it != (derivative_other_endo[block]).end(); it++)
-        tmp_block_other_endo_derivative[{ it->first.first, { it->first.second.second, it->first.second.first } }] = it->second;
+      map<tuple<int, int, int>, expr_t> tmp_block_other_endo_derivative;
+      for (const auto &it : derivative_other_endo[block])
+        tmp_block_other_endo_derivative[{ get<0>(it.first), get<2>(it.first), get<1>(it.first) }] = it.second;
       prev_var = 999999999;
       prev_lag = -9999999;
       count_col_other_endo = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_other_endo_derivative.begin(); it != tmp_block_other_endo_derivative.end(); it++)
+      for (const auto &it : tmp_block_other_endo_derivative)
         {
-          int lag = it->first.first;
-          unsigned int var = it->first.second.first;
+          int lag = get<0>(it.first);
+          unsigned int var = get<1>(it.first);
           if (var != prev_var || lag != prev_lag)
             {
               prev_var = var;
@@ -508,22 +506,22 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
       tmp1_output.str("");
       tmp1_output << packageDir(basename + ".block") << "/dynamic_" << block+1 << ".m";
       output.open(tmp1_output.str(), ios::out | ios::binary);
-      output << "%\n";
-      output << "% " << tmp1_output.str() << " : Computes dynamic model for Dynare\n";
-      output << "%\n";
-      output << "% Warning : this file is generated automatically by Dynare\n";
-      output << "%           from model file (.mod)\n\n";
-      output << "%/\n";
+      output << "%" << endl
+             << "% " << tmp1_output.str() << " : Computes dynamic model for Dynare" << endl
+             << "%" << endl
+             << "% Warning : this file is generated automatically by Dynare" << endl
+             << "%           from model file (.mod)" << endl << endl
+             << "%/" << endl;
       if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
         {
-          output << "function [y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, jacobian_eval, y_kmin, periods)\n";
+          output << "function [y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, jacobian_eval, y_kmin, periods)" << endl;
         }
       else if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
-        output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)\n";
+        output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl;
       else if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE)
-        output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)\n";
+        output << "function [residual, y, g1, g2, g3, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, it_, jacobian_eval)" << endl;
       else
-        output << "function [residual, y, g1, g2, g3, b, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, periods, jacobian_eval, y_kmin, y_size, Periods)\n";
+        output << "function [residual, y, g1, g2, g3, b, varargout] = dynamic_" << block+1 << "(y, x, params, steady_state, periods, jacobian_eval, y_kmin, y_size, Periods)" << endl;
       BlockType block_type;
       if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
         block_type = SIMULTAN;
@@ -548,47 +546,42 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
       //The Temporary terms
       if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
         {
-          output << "  if(jacobian_eval)\n";
-          output << "    g1 = spalloc(" << block_mfs  << ", " << count_col_endo << ", " << nze << ");\n";
-          output << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");\n";
-          output << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");\n";
-          output << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");\n";
-          output << "  end;\n";
+          output << "  if(jacobian_eval)" << endl
+                 << "    g1 = spalloc(" << block_mfs  << ", " << count_col_endo << ", " << nze << ");" << endl
+                 << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");" << endl
+                 << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");" << endl
+                 << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");" << endl
+                 << "  end;" << endl;
         }
       else
         {
-          output << "  if(jacobian_eval)\n";
-          output << "    g1 = spalloc(" << block_size << ", " << count_col_endo << ", " << nze << ");\n";
-          output << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");\n";
-          output << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");\n";
-          output << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");\n";
-          output << "  else\n";
+          output << "  if(jacobian_eval)" << endl
+                 << "    g1 = spalloc(" << block_size << ", " << count_col_endo << ", " << nze << ");" << endl
+                 << "    g1_x=spalloc(" << block_size << ", " << count_col_exo  << ", " << nze_exo << ");" << endl
+                 << "    g1_xd=spalloc(" << block_size << ", " << count_col_exo_det  << ", " << nze_exo_det << ");" << endl
+                 << "    g1_o=spalloc(" << block_size << ", " << count_col_other_endo << ", " << nze_other_endo << ");" << endl
+                 << "  else" << endl;
           if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
-            {
-              output << "    g1 = spalloc(" << block_mfs << "*Periods, "
-                     << block_mfs << "*(Periods+" << max_leadlag_block[block].first+max_leadlag_block[block].second+1 << ")"
-                     << ", " << nze << "*Periods);\n";
-            }
+            output << "    g1 = spalloc(" << block_mfs << "*Periods, "
+                   << block_mfs << "*(Periods+" << max_leadlag_block[block].first+max_leadlag_block[block].second+1 << ")"
+                   << ", " << nze << "*Periods);" << endl;
           else
-            {
-              output << "    g1 = spalloc(" << block_mfs
-                     << ", " << block_mfs << ", " << nze << ");\n";
-            }
-          output << "  end;\n";
+            output << "    g1 = spalloc(" << block_mfs
+                   << ", " << block_mfs << ", " << nze << ");" << endl;
+          output << "  end;" << endl;
         }
 
-      output << "  g2=0;g3=0;\n";
+      output << "  g2=0;g3=0;" << endl;
       if (v_temporary_terms_inuse[block].size())
         {
           tmp_output.str("");
           for (int it : v_temporary_terms_inuse[block])
             tmp_output << " T" << it;
-          output << "  global" << tmp_output.str() << ";\n";
+          output << "  global" << tmp_output.str() << ";" << endl;
         }
       if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
         {
           temporary_terms_t tt2;
-          tt2.clear();
           for (int i = 0; i < (int) block_size; i++)
             {
               if (v_temporary_terms[block][i].size() && global_temporary_terms)
@@ -606,13 +599,13 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
             }
         }
       if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
-        output << "  residual=zeros(" << block_mfs << ",1);\n";
+        output << "  residual=zeros(" << block_mfs << ",1);" << endl;
       else if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
-        output << "  residual=zeros(" << block_mfs << ",y_kmin+periods);\n";
+        output << "  residual=zeros(" << block_mfs << ",y_kmin+periods);" << endl;
       if (simulation_type == EVALUATE_BACKWARD)
-        output << "  for it_ = (y_kmin+periods):y_kmin+1\n";
+        output << "  for it_ = (y_kmin+periods):y_kmin+1" << endl;
       if (simulation_type == EVALUATE_FORWARD)
-        output << "  for it_ = y_kmin+1:(y_kmin+periods)\n";
+        output << "  for it_ = y_kmin+1:(y_kmin+periods)" << endl;
 
       if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
         {
@@ -633,7 +626,6 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
       for (unsigned int i = 0; i < block_size; i++)
         {
           temporary_terms_t tt2;
-          tt2.clear();
           if (v_temporary_terms[block].size())
             {
               output << "  " << "% //Temporary variables" << endl;
@@ -682,7 +674,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
                   if (isBlockEquationRenormalized(block, i))
                     {
                       rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
-                      output << "\n    ";
+                      output << endl << "    ";
                       tmp_output.str("");
                       eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
                       lhs = eq_node->get_arg1();
@@ -694,10 +686,10 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
                 }
               else
                 {
-                  cerr << "Type mismatch for equation " << equation_ID+1  << "\n";
+                  cerr << "Type mismatch for equation " << equation_ID+1  << endl;
                   exit(EXIT_FAILURE);
                 }
-              output << ";\n";
+              output << ";" << endl;
               break;
             case SOLVE_BACKWARD_SIMPLE:
             case SOLVE_FORWARD_SIMPLE:
@@ -727,10 +719,10 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
               output << tmp_output.str();
               output << ") - (";
               rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
-              output << ");\n";
+              output << ");" << endl;
 #ifdef CONDITION
               if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
-                output << "  condition(" << i+1 << ")=0;\n";
+                output << "  condition(" << i+1 << ")=0;" << endl;
 #endif
             }
         }
@@ -746,11 +738,11 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
       prev_var = 999999999;
       prev_lag = -9999999;
       count_col = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_endo_derivative.begin(); it != tmp_block_endo_derivative.end(); it++)
+      for (const auto &it : tmp_block_endo_derivative)
         {
-          int lag = it->first.first;
-          unsigned int var = it->first.second.first;
-          unsigned int eq = it->first.second.second;
+          int lag;
+          unsigned int var, eq;
+          tie(lag, var, eq) = it.first;
           int eqr = getBlockEquationID(block, eq);
           int varr = getBlockVariableID(block, var);
           if (var != prev_var || lag != prev_lag)
@@ -760,7 +752,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
               count_col++;
             }
 
-          expr_t id = it->second;
+          expr_t id = it.second;
 
           output << "      g1(" << eq+1 << ", " << count_col << ") = ";
           id->writeOutput(output, local_output_type, local_temporary_terms, {});
@@ -772,11 +764,11 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
       prev_var = 999999999;
       prev_lag = -9999999;
       count_col = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_exo_derivative.begin(); it != tmp_block_exo_derivative.end(); it++)
+      for (const auto &it : tmp_block_exo_derivative)
         {
-          int lag = it->first.first;
-          unsigned int var = it->first.second.first;
-          unsigned int eq = it->first.second.second;
+          int lag;
+          unsigned int var, eq;
+          tie(lag, var, eq) = it.first;
           int eqr = getBlockInitialEquationID(block, eq);
           if (var != prev_var || lag != prev_lag)
             {
@@ -784,7 +776,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
               prev_lag = lag;
               count_col++;
             }
-          expr_t id = it->second;
+          expr_t id = it.second;
           output << "      g1_x(" << eqr+1 << ", " << count_col << ") = ";
           id->writeOutput(output, local_output_type, local_temporary_terms, {});
           output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::exogenous, var))
@@ -795,11 +787,11 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
       prev_var = 999999999;
       prev_lag = -9999999;
       count_col = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_exo_det_derivative.begin(); it != tmp_block_exo_det_derivative.end(); it++)
+      for (const auto &it : tmp_block_exo_det_derivative)
         {
-          int lag = it->first.first;
-          unsigned int var = it->first.second.first;
-          unsigned int eq = it->first.second.second;
+          int lag;
+          unsigned int var, eq;
+          tie(lag, var, eq) = it.first;
           int eqr = getBlockInitialEquationID(block, eq);
           if (var != prev_var || lag != prev_lag)
             {
@@ -807,7 +799,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
               prev_lag = lag;
               count_col++;
             }
-          expr_t id = it->second;
+          expr_t id = it.second;
           output << "      g1_xd(" << eqr+1 << ", " << count_col << ") = ";
           id->writeOutput(output, local_output_type, local_temporary_terms, {});
           output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::exogenous, var))
@@ -818,11 +810,11 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
       prev_var = 999999999;
       prev_lag = -9999999;
       count_col = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_other_endo_derivative.begin(); it != tmp_block_other_endo_derivative.end(); it++)
+      for (const auto &it : tmp_block_other_endo_derivative)
         {
-          int lag = it->first.first;
-          unsigned int var = it->first.second.first;
-          unsigned int eq = it->first.second.second;
+          int lag;
+          unsigned int var, eq;
+          tie(lag, var, eq) = it.first;
           int eqr = getBlockInitialEquationID(block, eq);
           if (var != prev_var || lag != prev_lag)
             {
@@ -830,7 +822,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
               prev_lag = lag;
               count_col++;
             }
-          expr_t id = it->second;
+          expr_t id = it.second;
 
           output << "      g1_o(" << eqr+1 << ", " << /*var+1+(lag+block_max_lag)*block_size*/ count_col << ") = ";
           id->writeOutput(output, local_output_type, local_temporary_terms, {});
@@ -839,30 +831,30 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
                  << ") " << var+1
                  << ", equation=" << eq+1 << endl;
         }
-      output << "      varargout{1}=g1_x;\n";
-      output << "      varargout{2}=g1_xd;\n";
-      output << "      varargout{3}=g1_o;\n";
+      output << "      varargout{1}=g1_x;" << endl
+             << "      varargout{2}=g1_xd;" << endl
+             << "      varargout{3}=g1_o;" << endl;
 
       switch (simulation_type)
         {
         case EVALUATE_FORWARD:
         case EVALUATE_BACKWARD:
-          output << "    end;" << endl;
-          output << "  end;" << endl;
+          output << "    end;" << endl
+                 << "  end;" << endl;
           break;
         case SOLVE_BACKWARD_SIMPLE:
         case SOLVE_FORWARD_SIMPLE:
         case SOLVE_BACKWARD_COMPLETE:
         case SOLVE_FORWARD_COMPLETE:
           output << "  else" << endl;
-          for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+          for (const auto &it : blocks_derivatives[block])
             {
-              unsigned int eq = it->first.first;
-              unsigned int var = it->first.second;
+              unsigned int eq, var;
+              expr_t id;
+              int lag;
+              tie(eq, var, lag, id) = it;
               unsigned int eqr = getBlockEquationID(block, eq);
               unsigned int varr = getBlockVariableID(block, var);
-              expr_t id = it->second.second;
-              int lag = it->second.first;
               if (lag == 0)
                 {
                   output << "    g1(" << eq+1 << ", " << var+1-block_recursive << ") = ";
@@ -874,20 +866,20 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
                 }
 
             }
-          output << "  end;\n";
+          output << "  end;" << endl;
           break;
         case SOLVE_TWO_BOUNDARIES_SIMPLE:
         case SOLVE_TWO_BOUNDARIES_COMPLETE:
           output << "    else" << endl;
-          for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+          for (const auto &it : blocks_derivatives[block])
             {
-              unsigned int eq = it->first.first;
-              unsigned int var = it->first.second;
+              unsigned int eq, var;
+              int lag;
+              expr_t id;
+              tie(eq, var, lag, id) = it;
               unsigned int eqr = getBlockEquationID(block, eq);
               unsigned int varr = getBlockVariableID(block, var);
               ostringstream tmp_output;
-              expr_t id = it->second.second;
-              int lag = it->second.first;
               if (eq >= block_recursive && var >= block_recursive)
                 {
                   if (lag == 0)
@@ -930,17 +922,17 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
                 }
 
 #ifdef CONDITION
-              output << "  if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
-              output << "    condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
+              output << "  if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))" << endl
+                     << "    condition(" << eqr << ")=u(" << u << "+Per_u_);" << endl;
 #endif
             }
           for (unsigned int i = 0; i < block_size; i++)
             {
               if (i >= block_recursive)
-                output << "  " << Uf[getBlockEquationID(block, i)] << ";\n";
+                output << "  " << Uf[getBlockEquationID(block, i)] << ";" << endl;
 #ifdef CONDITION
-              output << "  if (fabs(condition(" << i+1 << "))<fabs(u(" << i << "+Per_u_)))\n";
-              output << "    condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);\n";
+              output << "  if (fabs(condition(" << i+1 << "))<fabs(u(" << i << "+Per_u_)))" << endl
+                     << "    condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);" << endl;
 #endif
             }
 #ifdef CONDITION
@@ -953,14 +945,14 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
                   unsigned int var = ModelBlock->Block_List[block].IM_lead_lag[m].Var_Index[i];
                   unsigned int u = ModelBlock->Block_List[block].IM_lead_lag[m].u[i];
                   unsigned int eqr = ModelBlock->Block_List[block].IM_lead_lag[m].Equ[i];
-                  output << "  u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n";
+                  output << "  u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");" << endl;
                 }
             }
           for (i = 0; i < ModelBlock->Block_List[block].Size; i++)
-            output << "  u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n";
+            output << "  u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");" << endl;
 #endif
-          output << "    end;" << endl;
-          output << "  end;" << endl;
+          output << "    end;" << endl
+                 << "  end;" << endl;
           break;
         default:
           break;
@@ -1013,8 +1005,8 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
   for (int i = 0; i < symbol_table.exo_nbr(); i++)
     exo.push_back(i);
 
-  map<pair< int, pair<int, int>>, expr_t> first_derivatives_reordered_endo;
-  map<pair< pair<int, SymbolType>, pair<int, int>>, expr_t>  first_derivatives_reordered_exo;
+  map<tuple<int, int, int>, expr_t> first_derivatives_reordered_endo;
+  map<tuple<int, SymbolType, int, int>, expr_t> first_derivatives_reordered_exo;
   for (const auto & first_derivative : derivatives[1])
     {
       int deriv_id = first_derivative.first[1];
@@ -1023,18 +1015,17 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
       unsigned int var = symbol_table.getTypeSpecificID(symb);
       int lag = getLagByDerivID(deriv_id);
       if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
-        first_derivatives_reordered_endo[{ lag, make_pair(var, eq) }] = first_derivative.second;
+        first_derivatives_reordered_endo[{ lag, var, eq }] = first_derivative.second;
       else if (getTypeByDerivID(deriv_id) == SymbolType::exogenous || getTypeByDerivID(deriv_id) == SymbolType::exogenousDet)
-        first_derivatives_reordered_exo[{ { lag, getTypeByDerivID(deriv_id) }, { var, eq } }] = first_derivative.second;
+        first_derivatives_reordered_exo[{ lag, getTypeByDerivID(deriv_id), var, eq }] = first_derivative.second;
     }
   int prev_var = -1;
   int prev_lag = -999999999;
   int count_col_endo = 0;
-  for (map<pair< int, pair<int, int>>, expr_t>::const_iterator it = first_derivatives_reordered_endo.begin();
-       it != first_derivatives_reordered_endo.end(); it++)
+  for (const auto &it : first_derivatives_reordered_endo)
     {
-      int var = it->first.second.first;
-      int lag = it->first.first;
+      int var, lag;
+      tie(lag, var, ignore) = it.first;
       if (prev_var != var || prev_lag != lag)
         {
           prev_var = var;
@@ -1050,9 +1041,9 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
 
   for (const auto & it : first_derivatives_reordered_exo)
     {
-      int var{it.first.second.first};
-      int lag{it.first.first.first};
-      SymbolType type{it.first.first.second};
+      int var, lag;
+      SymbolType type;
+      tie(lag, type, var, ignore) = it.first;
       if (prev_var != var || prev_lag != lag || prev_type != type)
         {
           prev_var = var;
@@ -1102,8 +1093,7 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
   fjmp_if_eval.write(code_file, instruction_number);
   int prev_instruction_number = instruction_number;
 
-  vector<vector<pair<pair<int, int>, int >>> my_derivatives;
-  my_derivatives.resize(symbol_table.endo_nbr());
+  vector<vector<tuple<int, int, int>>> my_derivatives(symbol_table.endo_nbr());;
   count_u = symbol_table.endo_nbr();
   for (const auto & first_derivative : derivatives[1])
     {
@@ -1119,7 +1109,7 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
           fnumexpr.write(code_file, instruction_number);
           if (!my_derivatives[eq].size())
             my_derivatives[eq].clear();
-          my_derivatives[eq].emplace_back(make_pair(var, lag), count_u);
+          my_derivatives[eq].emplace_back(var, lag, count_u);
           d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
 
           FSTPU_ fstpu(count_u);
@@ -1133,12 +1123,11 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
       fldr.write(code_file, instruction_number);
       if (my_derivatives[i].size())
         {
-          for (vector<pair<pair<int, int>, int>>::const_iterator it = my_derivatives[i].begin();
-               it != my_derivatives[i].end(); it++)
+          for (auto it = my_derivatives[i].begin(); it != my_derivatives[i].end(); it++)
             {
-              FLDU_ fldu(it->second);
+              FLDU_ fldu(get<2>(*it));
               fldu.write(code_file, instruction_number);
-              FLDV_ fldv{static_cast<int>(SymbolType::endogenous), static_cast<unsigned int>(it->first.first), it->first.second};
+              FLDV_ fldv{static_cast<int>(SymbolType::endogenous), static_cast<unsigned int>(get<0>(*it)), get<1>(*it)};
               fldv.write(code_file, instruction_number);
               FBINARY_ fbinary{static_cast<int>(BinaryOpcode::times)};
               fbinary.write(code_file, instruction_number);
@@ -1171,13 +1160,12 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
   prev_var = -1;
   prev_lag = -999999999;
   count_col_endo = 0;
-  for (map<pair< int, pair<int, int>>, expr_t>::const_iterator it = first_derivatives_reordered_endo.begin();
-       it != first_derivatives_reordered_endo.end(); it++)
+  for (const auto &it : first_derivatives_reordered_endo)
     {
-      unsigned int eq = it->first.second.second;
-      int var = it->first.second.first;
-      int lag = it->first.first;
-      expr_t d1 = it->second;
+      unsigned int eq;
+      int var, lag;
+      tie(lag, var, eq) = it.first;
+      expr_t d1 = it.second;
       FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var, lag);
       fnumexpr.write(code_file, instruction_number);
       if (prev_var != var || prev_lag != lag)
@@ -1195,10 +1183,9 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m
   count_col_exo = 0;
   for (const auto & it : first_derivatives_reordered_exo)
     {
-      int eq{it.first.second.second};
-      int var{it.first.second.first};
-      int lag{it.first.first.first};
-      expr_t d1{it.second};
+      int lag, var, eq;
+      tie(lag, ignore, var, eq) = it.first;
+      expr_t d1 = it.second;
       FNUMEXPR_ fnumexpr(FirstExoDerivative, eq, var, lag);
       fnumexpr.write(code_file, instruction_number);
       if (prev_var != var || prev_lag != lag)
@@ -1292,25 +1279,25 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
                                       simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE, linear_decomposition);
           file_open = true;
         }
-      map<pair<int, pair<int, int>>, expr_t> tmp_block_endo_derivative;
-      for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
-        tmp_block_endo_derivative[{ it->second.first, { it->first.second, it->first.first } }] = it->second.second;
-      map<pair<int, pair<int, int>>, expr_t> tmp_exo_derivative;
-      for (auto it = derivative_exo[block].begin(); it != (derivative_exo[block]).end(); it++)
-        tmp_exo_derivative[{ it->first.first, { it->first.second.second, it->first.second.first } }] = it->second;
-      map<pair<int, pair<int, int>>, expr_t> tmp_exo_det_derivative;
-      for (auto it = derivative_exo_det[block].begin(); it != (derivative_exo_det[block]).end(); it++)
-        tmp_exo_det_derivative[{ it->first.first, { it->first.second.second, it->first.second.first } }] = it->second;
-      map<pair<int, pair<int, int>>, expr_t> tmp_other_endo_derivative;
-      for (auto it = derivative_other_endo[block].begin(); it != (derivative_other_endo[block]).end(); it++)
-        tmp_other_endo_derivative[{ it->first.first, { it->first.second.second, it->first.second.first } }] = it->second;
+      map<tuple<int, int, int>, expr_t> tmp_block_endo_derivative;
+      for (const auto &it : blocks_derivatives[block])
+        tmp_block_endo_derivative[{ get<2>(it), get<1>(it), get<0>(it) }] = get<3>(it);
+      map<tuple<int, int, int>, expr_t> tmp_exo_derivative;
+      for (const auto &it : derivative_exo[block])
+        tmp_exo_derivative[{ get<0>(it.first), get<2>(it.first), get<1>(it.first) }] = it.second;
+      map<tuple<int, int, int>, expr_t> tmp_exo_det_derivative;
+      for (const auto &it : derivative_exo_det[block])
+        tmp_exo_det_derivative[{ get<0>(it.first), get<2>(it.first), get<1>(it.first) }] = it.second;
+      map<tuple<int, int, int>, expr_t> tmp_other_endo_derivative;
+      for (const auto &it : derivative_other_endo[block])
+        tmp_other_endo_derivative[{ get<0>(it.first), get<2>(it.first), get<1>(it.first) }] = it.second;
       int prev_var = -1;
       int prev_lag = -999999999;
       int count_col_endo = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_endo_derivative.begin(); it != tmp_block_endo_derivative.end(); it++)
+      for (const auto &it : tmp_block_endo_derivative)
         {
-          int lag = it->first.first;
-          int var = it->first.second.first;
+          int lag, var;
+          tie(lag, var, ignore) = it.first;
           if (prev_var != var || prev_lag != lag)
             {
               prev_var = var;
@@ -1321,31 +1308,31 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
       unsigned int count_col_det_exo = 0;
       vector<unsigned int> exo_det;
       for (const auto & it : exo_det_block[block])
-        for (auto it1 = it.second.begin(); it1 != it.second.end(); it1++)
+        for (const auto &it1 : it.second)
           {
             count_col_det_exo++;
-            if (find(exo_det.begin(), exo_det.end(), *it1) == exo_det.end())
-              exo_det.push_back(*it1);
+            if (find(exo_det.begin(), exo_det.end(), it1) == exo_det.end())
+              exo_det.push_back(it1);
           }
 
       unsigned int count_col_exo = 0;
       vector<unsigned int> exo;
       for (const auto & it : exo_block[block])
-        for (auto it1 = it.second.begin(); it1 != it.second.end(); it1++)
+        for (const auto &it1 : it.second)
           {
             count_col_exo++;
-            if (find(exo.begin(), exo.end(), *it1) == exo.end())
-              exo.push_back(*it1);
+            if (find(exo.begin(), exo.end(), it1) == exo.end())
+              exo.push_back(it1);
           }
 
       vector<unsigned int> other_endo;
       unsigned int count_col_other_endo = 0;
       for (const auto & it : other_endo_block[block])
-        for (auto it1 = it.second.begin(); it1 != it.second.end(); it1++)
+        for (const auto &it1 : it.second)
           {
             count_col_other_endo++;
-            if (find(other_endo.begin(), other_endo.end(), *it1) == other_endo.end())
-              other_endo.push_back(*it1);
+            if (find(other_endo.begin(), other_endo.end(), it1) == other_endo.end())
+              other_endo.push_back(it1);
           }
 
       FBEGINBLOCK_ fbeginblock(block_mfs,
@@ -1380,7 +1367,6 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
         {
           //The Temporary terms
           temporary_terms_t tt2;
-          tt2.clear();
           if (v_temporary_terms[block][i].size() && !linear_decomposition)
             {
               for (auto it : v_temporary_terms[block][i])
@@ -1396,7 +1382,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
                   // Insert current node into tt2
                   tt2.insert(it);
 #ifdef DEBUGC
-                  cout << "FSTPT " << v << "\n";
+                  cout << "FSTPT " << v << endl;
                   instruction_number++;
                   code_file.write(&FOK, sizeof(FOK));
                   code_file.write(reinterpret_cast<char *>(&k), sizeof(k));
@@ -1406,11 +1392,10 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
                 }
             }
 #ifdef DEBUGC
-          for (temporary_terms_t::const_iterator it = v_temporary_terms[block][i].begin();
-               it != v_temporary_terms[block][i].end(); it++)
+          for (const auto &it : v_temporary_terms[block][i])
             {
-              map_idx_t::const_iterator ii = map_idx.find((*it)->idx);
-              cout << "map_idx[" << (*it)->idx <<"]=" << ii->second << "\n";
+              auto ii = map_idx.find(it->idx);
+              cout << "map_idx[" << it->idx <<"]=" << ii->second << endl;
             }
 #endif
 
@@ -1503,11 +1488,11 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
             case SOLVE_TWO_BOUNDARIES_COMPLETE:
             case SOLVE_TWO_BOUNDARIES_SIMPLE:
               count_u = feedback_variables.size();
-              for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+              for (const auto &it : blocks_derivatives[block])
                 {
-                  int lag = it->second.first;
-                  unsigned int eq = it->first.first;
-                  unsigned int var = it->first.second;
+                  unsigned int eq, var;
+                  int lag;
+                  tie(eq, var, lag, ignore) = it;
                   unsigned int eqr = getBlockEquationID(block, eq);
                   unsigned int varr = getBlockVariableID(block, var);
                   if (eq >= block_recursive and var >= block_recursive)
@@ -1595,11 +1580,11 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
       prev_var = -1;
       prev_lag = -999999999;
       count_col_endo = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_block_endo_derivative.begin(); it != tmp_block_endo_derivative.end(); it++)
+      for (const auto &it : tmp_block_endo_derivative)
         {
-          int lag = it->first.first;
-          unsigned int eq = it->first.second.second;
-          int var = it->first.second.first;
+          int lag, var;
+          unsigned int eq;
+          tie(lag, var, eq) = it.first;
           unsigned int eqr = getBlockEquationID(block, eq);
           unsigned int varr = getBlockVariableID(block, var);
           if (prev_var != var || prev_lag != lag)
@@ -1617,11 +1602,10 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
       prev_var = -1;
       prev_lag = -999999999;
       count_col_exo = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_exo_derivative.begin(); it != tmp_exo_derivative.end(); it++)
+      for (const auto &it : tmp_exo_derivative)
         {
-          int lag = it->first.first;
-          int eq = it->first.second.second;
-          int var = it->first.second.first;
+          int lag, eq, var;
+          tie(lag, var, eq) = it.first;
           int eqr = getBlockInitialEquationID(block, eq);
           int varr = getBlockInitialExogenousID(block, var);
           if (prev_var != var || prev_lag != lag)
@@ -1630,7 +1614,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
               prev_lag = lag;
               count_col_exo++;
             }
-          expr_t id = it->second;
+          expr_t id = it.second;
 
           FNUMEXPR_ fnumexpr(FirstExoDerivative, eqr, varr, lag);
           fnumexpr.write(code_file, instruction_number);
@@ -1641,11 +1625,10 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
       prev_var = -1;
       prev_lag = -999999999;
       int count_col_exo_det = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_exo_det_derivative.begin(); it != tmp_exo_det_derivative.end(); it++)
+      for (const auto &it : tmp_exo_det_derivative)
         {
-          int lag = it->first.first;
-          int eq = it->first.second.second;
-          int var = it->first.second.first;
+          int lag, eq, var;
+          tie(lag, var, eq) = it.first;
           int eqr = getBlockInitialEquationID(block, eq);
           int varr = getBlockInitialDetExogenousID(block, var);
           if (prev_var != var || prev_lag != lag)
@@ -1654,7 +1637,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
               prev_lag = lag;
               count_col_exo_det++;
             }
-          expr_t id = it->second;
+          expr_t id = it.second;
 
           FNUMEXPR_ fnumexpr(FirstExodetDerivative, eqr, varr, lag);
           fnumexpr.write(code_file, instruction_number);
@@ -1665,11 +1648,10 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
       prev_var = -1;
       prev_lag = -999999999;
       count_col_other_endo = 0;
-      for (map<pair<int, pair<int, int>>, expr_t>::const_iterator it = tmp_other_endo_derivative.begin(); it != tmp_other_endo_derivative.end(); it++)
+      for (const auto &it : tmp_other_endo_derivative)
         {
-          int lag = it->first.first;
-          int eq = it->first.second.second;
-          int var = it->first.second.first;
+          int lag, eq, var;
+          tie(lag, var, eq) = it.first;
           int eqr = getBlockInitialEquationID(block, eq);
           int varr = getBlockInitialOtherEndogenousID(block, var);;
           if (prev_var != var || prev_lag != lag)
@@ -1678,7 +1660,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id
               prev_lag = lag;
               count_col_other_endo++;
             }
-          expr_t id = it->second;
+          expr_t id = it.second;
 
           FNUMEXPR_ fnumexpr(FirstOtherEndoDerivative, eqr, varr, lag);
           fnumexpr.write(code_file, instruction_number);
@@ -1936,11 +1918,11 @@ DynamicModel::Write_Inf_To_Bin_File_Block(const string &basename, const int &num
   unsigned int block_size = getBlockSize(num);
   unsigned int block_mfs = getBlockMfs(num);
   unsigned int block_recursive = block_size - block_mfs;
-  for (auto it = blocks_derivatives[num].begin(); it != (blocks_derivatives[num]).end(); it++)
+  for (const auto &it : blocks_derivatives[num])
     {
-      unsigned int eq = it->first.first;
-      unsigned int var = it->first.second;
-      int lag = it->second.first;
+      unsigned int eq, var;
+      int lag;
+      tie(eq, var, lag, ignore) = it;
       if (lag != 0 && !is_two_boundaries)
         continue;
       if (eq >= block_recursive && var >= block_recursive)
@@ -1985,18 +1967,18 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
       cerr << "Error: Can't open file " << filename << " for writing" << endl;
       exit(EXIT_FAILURE);
     }
-  mDynamicModelFile << "%\n";
-  mDynamicModelFile << "% " << filename << " : Computes dynamic model for Dynare\n";
-  mDynamicModelFile << "%\n";
-  mDynamicModelFile << "% Warning : this file is generated automatically by Dynare\n";
-  mDynamicModelFile << "%           from model file (.mod)\n\n";
-  mDynamicModelFile << "%/\n";
+  mDynamicModelFile << "%" << endl
+                    << "% " << filename << " : Computes dynamic model for Dynare" << endl
+                    << "%" << endl
+                    << "% Warning : this file is generated automatically by Dynare" << endl
+                    << "%           from model file (.mod)" << endl << endl
+                    << "%/" << endl;
 
   int Nb_SGE = 0;
   bool open_par = false;
 
-  mDynamicModelFile << "function [varargout] = dynamic(options_, M_, oo_, varargin)\n";
-  mDynamicModelFile << "  g2=[];g3=[];\n";
+  mDynamicModelFile << "function [varargout] = dynamic(options_, M_, oo_, varargin)" << endl
+                    << "  g2=[];g3=[];" << endl;
   //Temporary variables declaration
   OK = true;
   ostringstream tmp_output;
@@ -2010,16 +1992,16 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
       temporary_term->writeOutput(tmp_output, ExprNodeOutputType::matlabStaticModelSparse, temporary_terms, {});
     }
   if (tmp_output.str().length() > 0)
-    mDynamicModelFile << "  global " << tmp_output.str() << ";\n";
+    mDynamicModelFile << "  global " << tmp_output.str() << ";" << endl;
 
-  mDynamicModelFile << "  T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);\n";
+  mDynamicModelFile << "  T_init=zeros(1,options_.periods+M_.maximum_lag+M_.maximum_lead);" << endl;
   tmp_output.str("");
   for (auto temporary_term : temporary_terms)
     {
       tmp_output << "  ";
       // In the following, "Static" is used to avoid getting the "(it_)" subscripting
       temporary_term->writeOutput(tmp_output, ExprNodeOutputType::matlabStaticModelSparse, temporary_terms, {});
-      tmp_output << "=T_init;\n";
+      tmp_output << "=T_init;" << endl;
     }
   if (tmp_output.str().length() > 0)
     mDynamicModelFile << tmp_output.str();
@@ -2065,30 +2047,30 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
               tmp_eq << " " << getBlockEquationID(block, ik)+1;
             }
         }
-      mDynamicModelFile << "    y_index_eq=[" << tmp_eq.str() << "];\n";
-      mDynamicModelFile << "    y_index=[" << tmp.str() << "];\n";
+      mDynamicModelFile << "    y_index_eq=[" << tmp_eq.str() << "];" << endl
+                        << "    y_index=[" << tmp.str() << "];" << endl;
 
       switch (simulation_type)
         {
         case EVALUATE_FORWARD:
         case EVALUATE_BACKWARD:
-          mDynamicModelFile << "    [y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, 1, it_-1, 1);\n";
-          mDynamicModelFile << "    residual(y_index_eq)=ys(y_index)-y(it_, y_index);\n";
+          mDynamicModelFile << "    [y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, 1, it_-1, 1);" << endl
+                            << "    residual(y_index_eq)=ys(y_index)-y(it_, y_index);" << endl;
           break;
         case SOLVE_FORWARD_SIMPLE:
         case SOLVE_BACKWARD_SIMPLE:
-          mDynamicModelFile << "    [r, y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, it_, 1);\n";
-          mDynamicModelFile << "    residual(y_index_eq)=r;\n";
+          mDynamicModelFile << "    [r, y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, it_, 1);" << endl
+                            << "    residual(y_index_eq)=r;" << endl;
           break;
         case SOLVE_FORWARD_COMPLETE:
         case SOLVE_BACKWARD_COMPLETE:
-          mDynamicModelFile << "    [r, y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, it_, 1);\n";
-          mDynamicModelFile << "    residual(y_index_eq)=r;\n";
+          mDynamicModelFile << "    [r, y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, it_, 1);" << endl
+                            << "    residual(y_index_eq)=r;" << endl;
           break;
         case SOLVE_TWO_BOUNDARIES_COMPLETE:
         case SOLVE_TWO_BOUNDARIES_SIMPLE:
-          mDynamicModelFile << "    [r, y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" <<  block + 1 << "(y, x, params, steady_state, it_-" << max_lag << ", 1, " << max_lag << ", " << block_recursive << "," << "options_.periods" << ");\n";
-          mDynamicModelFile << "    residual(y_index_eq)=r(:,M_.maximum_lag+1);\n";
+          mDynamicModelFile << "    [r, y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_xd, dr(" << count_call << ").g1_o]=" << basename << ".block.dynamic_" <<  block + 1 << "(y, x, params, steady_state, it_-" << max_lag << ", 1, " << max_lag << ", " << block_recursive << "," << "options_.periods" << ");" << endl
+                            << "    residual(y_index_eq)=r(:,M_.maximum_lag+1);" << endl;
           break;
         default:
           break;
@@ -2128,11 +2110,10 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
                     << "  maxit_=options_.simul.maxit;" << endl
                     << "  solve_tolf=options_.solve_tolf;" << endl
                     << "  y=oo_.endo_simul';" << endl
-                    << "  x=oo_.exo_simul;" << endl;
-
-  mDynamicModelFile << "  params=M_.params;\n";
-  mDynamicModelFile << "  steady_state=oo_.steady_state;\n";
-  mDynamicModelFile << "  oo_.deterministic_simulation.status = 0;\n";
+                    << "  x=oo_.exo_simul;" << endl
+                    << "  params=M_.params;" << endl
+                    << "  steady_state=oo_.steady_state;" << endl
+                    << "  oo_.deterministic_simulation.status = 0;" << endl;
   for (block = 0; block < nb_blocks; block++)
     {
       unsigned int block_size = getBlockSize(block);
@@ -2143,165 +2124,155 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
       if ((simulation_type == EVALUATE_FORWARD) && (block_size))
         {
           if (open_par)
-            {
-              mDynamicModelFile << "  end\n";
-            }
-          mDynamicModelFile << "  oo_.deterministic_simulation.status = 1;\n";
-          mDynamicModelFile << "  oo_.deterministic_simulation.error = 0;\n";
-          mDynamicModelFile << "  oo_.deterministic_simulation.iterations = 0;\n";
-          mDynamicModelFile << "  if(isfield(oo_.deterministic_simulation,'block'))\n";
-          mDynamicModelFile << "    blck_num = length(oo_.deterministic_simulation.block)+1;\n";
-          mDynamicModelFile << "  else\n";
-          mDynamicModelFile << "    blck_num = 1;\n";
-          mDynamicModelFile << "  end;\n";
-          mDynamicModelFile << "  oo_.deterministic_simulation.block(blck_num).status = 1;\n";
-          mDynamicModelFile << "  oo_.deterministic_simulation.block(blck_num).error = 0;\n";
-          mDynamicModelFile << "  oo_.deterministic_simulation.block(blck_num).iterations = 0;\n";
-          mDynamicModelFile << "  g1=[];g2=[];g3=[];\n";
-          mDynamicModelFile << "  y=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, 0, y_kmin, periods);\n";
-          mDynamicModelFile << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);\n";
-          mDynamicModelFile << "  if any(isnan(tmp) | isinf(tmp))\n";
-          mDynamicModelFile << "    disp(['Inf or Nan value during the evaluation of block " << block <<"']);\n";
-          mDynamicModelFile << "    oo_.deterministic_simulation.status = 0;\n";
-          mDynamicModelFile << "    oo_.deterministic_simulation.error = 100;\n";
-          mDynamicModelFile << "    varargout{1} = oo_;\n";
-          mDynamicModelFile << "    return;\n";
-          mDynamicModelFile << "  end;\n";
+            mDynamicModelFile << "  end" << endl;
+          mDynamicModelFile << "  oo_.deterministic_simulation.status = 1;" << endl
+                            << "  oo_.deterministic_simulation.error = 0;" << endl
+                            << "  oo_.deterministic_simulation.iterations = 0;" << endl
+                            << "  if(isfield(oo_.deterministic_simulation,'block'))" << endl
+                            << "    blck_num = length(oo_.deterministic_simulation.block)+1;" << endl
+                            << "  else" << endl
+                            << "    blck_num = 1;" << endl
+                            << "  end;" << endl
+                            << "  oo_.deterministic_simulation.block(blck_num).status = 1;" << endl
+                            << "  oo_.deterministic_simulation.block(blck_num).error = 0;" << endl
+                            << "  oo_.deterministic_simulation.block(blck_num).iterations = 0;" << endl
+                            << "  g1=[];g2=[];g3=[];" << endl
+                            << "  y=" << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, 0, y_kmin, periods);" << endl
+                            << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
+                            << "  if any(isnan(tmp) | isinf(tmp))" << endl
+                            << "    disp(['Inf or Nan value during the evaluation of block " << block <<"']);" << endl
+                            << "    oo_.deterministic_simulation.status = 0;" << endl
+                            << "    oo_.deterministic_simulation.error = 100;" << endl
+                            << "    varargout{1} = oo_;" << endl
+                            << "    return;" << endl
+                            << "  end;" << endl;
         }
       else if ((simulation_type == EVALUATE_BACKWARD) && (block_size))
         {
           if (open_par)
-            {
-              mDynamicModelFile << "  end\n";
-            }
-          mDynamicModelFile << "  oo_.deterministic_simulation.status = 1;\n";
-          mDynamicModelFile << "  oo_.deterministic_simulation.error = 0;\n";
-          mDynamicModelFile << "  oo_.deterministic_simulation.iterations = 0;\n";
-          mDynamicModelFile << "  if(isfield(oo_.deterministic_simulation,'block'))\n";
-          mDynamicModelFile << "    blck_num = length(oo_.deterministic_simulation.block)+1;\n";
-          mDynamicModelFile << "  else\n";
-          mDynamicModelFile << "    blck_num = 1;\n";
-          mDynamicModelFile << "  end;\n";
-          mDynamicModelFile << "  oo_.deterministic_simulation.block(blck_num).status = 1;\n";
-          mDynamicModelFile << "  oo_.deterministic_simulation.block(blck_num).error = 0;\n";
-          mDynamicModelFile << "  oo_.deterministic_simulation.block(blck_num).iterations = 0;\n";
-          mDynamicModelFile << "  g1=[];g2=[];g3=[];\n";
-          mDynamicModelFile << "  " << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, 0, y_kmin, periods);\n";
-          mDynamicModelFile << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);\n";
-          mDynamicModelFile << "  if any(isnan(tmp) | isinf(tmp))\n";
-          mDynamicModelFile << "    disp(['Inf or Nan value during the evaluation of block " << block <<"']);\n";
-          mDynamicModelFile << "    oo_.deterministic_simulation.status = 0;\n";
-          mDynamicModelFile << "    oo_.deterministic_simulation.error = 100;\n";
-          mDynamicModelFile << "    varargout{1} = oo_;\n";
-          mDynamicModelFile << "    return;\n";
-          mDynamicModelFile << "  end;\n";
+            mDynamicModelFile << "  end" << endl;
+          mDynamicModelFile << "  oo_.deterministic_simulation.status = 1;" << endl
+                            << "  oo_.deterministic_simulation.error = 0;" << endl
+                            << "  oo_.deterministic_simulation.iterations = 0;" << endl
+                            << "  if(isfield(oo_.deterministic_simulation,'block'))" << endl
+                            << "    blck_num = length(oo_.deterministic_simulation.block)+1;" << endl
+                            << "  else" << endl
+                            << "    blck_num = 1;" << endl
+                            << "  end;" << endl
+                            << "  oo_.deterministic_simulation.block(blck_num).status = 1;" << endl
+                            << "  oo_.deterministic_simulation.block(blck_num).error = 0;" << endl
+                            << "  oo_.deterministic_simulation.block(blck_num).iterations = 0;" << endl
+                            << "  g1=[];g2=[];g3=[];" << endl
+                            << "  " << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, 0, y_kmin, periods);" << endl
+                            << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
+                            << "  if any(isnan(tmp) | isinf(tmp))" << endl
+                            << "    disp(['Inf or Nan value during the evaluation of block " << block <<"']);" << endl
+                            << "    oo_.deterministic_simulation.status = 0;" << endl
+                            << "    oo_.deterministic_simulation.error = 100;" << endl
+                            << "    varargout{1} = oo_;" << endl
+                            << "    return;" << endl
+                            << "  end;" << endl;
         }
       else if ((simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_FORWARD_SIMPLE) && (block_size))
         {
           if (open_par)
-            mDynamicModelFile << "  end\n";
+            mDynamicModelFile << "  end" << endl;
           open_par = false;
-          mDynamicModelFile << "  g1=0;\n";
-          mDynamicModelFile << "  r=0;\n";
+          mDynamicModelFile << "  g1=0;" << endl
+                            << "  r=0;" << endl;
           tmp.str("");
           for (unsigned int ik = block_recursive; ik < block_size; ik++)
-            {
-              tmp << " " << getBlockVariableID(block, ik)+1;
-            }
-          mDynamicModelFile << "  y_index = [" << tmp.str() << "];\n";
+            tmp << " " << getBlockVariableID(block, ik)+1;
+          mDynamicModelFile << "  y_index = [" << tmp.str() << "];" << endl;
           int nze = blocks_derivatives[block].size();
-          mDynamicModelFile << "  if(isfield(oo_.deterministic_simulation,'block'))\n";
-          mDynamicModelFile << "    blck_num = length(oo_.deterministic_simulation.block)+1;\n";
-          mDynamicModelFile << "  else\n";
-          mDynamicModelFile << "    blck_num = 1;\n";
-          mDynamicModelFile << "  end;\n";
-          mDynamicModelFile << "  y = solve_one_boundary('" << basename << ".block.dynamic_" <<  block + 1 << "'"
-                            <<", y, x, params, steady_state, y_index, " << nze
-                            <<", options_.periods, " << blocks_linear[block]
-                            <<", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0);\n";
-          mDynamicModelFile << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);\n";
-          mDynamicModelFile << "  if any(isnan(tmp) | isinf(tmp))\n";
-          mDynamicModelFile << "    disp(['Inf or Nan value during the resolution of block " << block <<"']);\n";
-          mDynamicModelFile << "    oo_.deterministic_simulation.status = 0;\n";
-          mDynamicModelFile << "    oo_.deterministic_simulation.error = 100;\n";
-          mDynamicModelFile << "    varargout{1} = oo_;\n";
-          mDynamicModelFile << "    return;\n";
-          mDynamicModelFile << "  end;\n";
+          mDynamicModelFile << "  if(isfield(oo_.deterministic_simulation,'block'))" << endl
+                            << "    blck_num = length(oo_.deterministic_simulation.block)+1;" << endl
+                            << "  else" << endl
+                            << "    blck_num = 1;" << endl
+                            << "  end;" << endl
+                            << "  y = solve_one_boundary('" << basename << ".block.dynamic_" <<  block + 1 << "'"
+                            << ", y, x, params, steady_state, y_index, " << nze
+                            << ", options_.periods, " << blocks_linear[block]
+                            << ", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0);" << endl
+                            << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
+                            << "  if any(isnan(tmp) | isinf(tmp))" << endl
+                            << "    disp(['Inf or Nan value during the resolution of block " << block <<"']);" << endl
+                            << "    oo_.deterministic_simulation.status = 0;" << endl
+                            << "    oo_.deterministic_simulation.error = 100;" << endl
+                            << "    varargout{1} = oo_;" << endl
+                            << "    return;" << endl
+                            << "  end;" << endl;
         }
       else if ((simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_SIMPLE) && (block_size))
         {
           if (open_par)
-            mDynamicModelFile << "  end\n";
+            mDynamicModelFile << "  end" << endl;
           open_par = false;
-          mDynamicModelFile << "  g1=0;\n";
-          mDynamicModelFile << "  r=0;\n";
+          mDynamicModelFile << "  g1=0;" << endl
+                            << "  r=0;" << endl;
           tmp.str("");
           for (unsigned int ik = block_recursive; ik < block_size; ik++)
-            {
-              tmp << " " << getBlockVariableID(block, ik)+1;
-            }
-          mDynamicModelFile << "  y_index = [" << tmp.str() << "];\n";
+            tmp << " " << getBlockVariableID(block, ik)+1;
+          mDynamicModelFile << "  y_index = [" << tmp.str() << "];" << endl;
           int nze = blocks_derivatives[block].size();
 
-          mDynamicModelFile << "  if(isfield(oo_.deterministic_simulation,'block'))\n";
-          mDynamicModelFile << "    blck_num = length(oo_.deterministic_simulation.block)+1;\n";
-          mDynamicModelFile << "  else\n";
-          mDynamicModelFile << "    blck_num = 1;\n";
-          mDynamicModelFile << "  end;\n";
-          mDynamicModelFile << "  y = solve_one_boundary('" << basename << ".block.dynamic_" <<  block + 1 << "'"
+          mDynamicModelFile << "  if(isfield(oo_.deterministic_simulation,'block'))" << endl
+                            << "    blck_num = length(oo_.deterministic_simulation.block)+1;" << endl
+                            << "  else" << endl
+                            << "    blck_num = 1;" << endl
+                            << "  end;" << endl
+                            << "  y = solve_one_boundary('" << basename << ".block.dynamic_" <<  block + 1 << "'"
                             <<", y, x, params, steady_state, y_index, " << nze
                             <<", options_.periods, " << blocks_linear[block]
-                            <<", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0);\n";
-          mDynamicModelFile << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);\n";
-          mDynamicModelFile << "  if any(isnan(tmp) | isinf(tmp))\n";
-          mDynamicModelFile << "    disp(['Inf or Nan value during the resolution of block " << block <<"']);\n";
-          mDynamicModelFile << "    oo_.deterministic_simulation.status = 0;\n";
-          mDynamicModelFile << "    oo_.deterministic_simulation.error = 100;\n";
-          mDynamicModelFile << "    varargout{1} = oo_;\n";
-          mDynamicModelFile << "    return;\n";
-          mDynamicModelFile << "  end;\n";
+                            <<", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0);" << endl
+                            << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
+                            << "  if any(isnan(tmp) | isinf(tmp))" << endl
+                            << "    disp(['Inf or Nan value during the resolution of block " << block <<"']);" << endl
+                            << "    oo_.deterministic_simulation.status = 0;" << endl
+                            << "    oo_.deterministic_simulation.error = 100;" << endl
+                            << "    varargout{1} = oo_;" << endl
+                            << "    return;" << endl
+                            << "  end;" << endl;
         }
       else if ((simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE) && (block_size))
         {
           if (open_par)
-            mDynamicModelFile << "  end\n";
+            mDynamicModelFile << "  end" << endl;
           open_par = false;
           Nb_SGE++;
           int nze = blocks_derivatives[block].size();
           mDynamicModelFile << "  y_index=[";
           for (unsigned int ik = block_recursive; ik < block_size; ik++)
-            {
-              mDynamicModelFile << " " << getBlockVariableID(block, ik)+1;
-            }
-          mDynamicModelFile << "  ];\n";
-          mDynamicModelFile << "  if(isfield(oo_.deterministic_simulation,'block'))\n";
-          mDynamicModelFile << "    blck_num = length(oo_.deterministic_simulation.block)+1;\n";
-          mDynamicModelFile << "  else\n";
-          mDynamicModelFile << "    blck_num = 1;\n";
-          mDynamicModelFile << "  end;\n";
-          mDynamicModelFile << "  [y oo_] = solve_two_boundaries('" << basename << ".block.dynamic_" <<  block + 1 << "'"
+            mDynamicModelFile << " " << getBlockVariableID(block, ik)+1;
+          mDynamicModelFile << "  ];" << endl
+                            << "  if(isfield(oo_.deterministic_simulation,'block'))" << endl
+                            << "    blck_num = length(oo_.deterministic_simulation.block)+1;" << endl
+                            << "  else" << endl
+                            << "    blck_num = 1;" << endl
+                            << "  end;" << endl
+                            << "  [y oo_] = solve_two_boundaries('" << basename << ".block.dynamic_" <<  block + 1 << "'"
                             <<", y, x, params, steady_state, y_index, " << nze
                             <<", options_.periods, " << max_leadlag_block[block].first
                             <<", " << max_leadlag_block[block].second
                             <<", " << blocks_linear[block]
-                            <<", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, options_, M_, oo_);\n";
-          mDynamicModelFile << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);\n";
-          mDynamicModelFile << "  if any(isnan(tmp) | isinf(tmp))\n";
-          mDynamicModelFile << "    disp(['Inf or Nan value during the resolution of block " << block <<"']);\n";
-          mDynamicModelFile << "    oo_.deterministic_simulation.status = 0;\n";
-          mDynamicModelFile << "    oo_.deterministic_simulation.error = 100;\n";
-          mDynamicModelFile << "    varargout{1} = oo_;\n";
-          mDynamicModelFile << "    return;\n";
-          mDynamicModelFile << "  end;\n";
+                            <<", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, options_, M_, oo_);" << endl
+                            << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
+                            << "  if any(isnan(tmp) | isinf(tmp))" << endl
+                            << "    disp(['Inf or Nan value during the resolution of block " << block <<"']);" << endl
+                            << "    oo_.deterministic_simulation.status = 0;" << endl
+                            << "    oo_.deterministic_simulation.error = 100;" << endl
+                            << "    varargout{1} = oo_;" << endl
+                            << "    return;" << endl
+                            << "  end;" << endl;
         }
     }
   if (open_par)
-    mDynamicModelFile << "  end;\n";
+    mDynamicModelFile << "  end;" << endl;
   open_par = false;
-  mDynamicModelFile << "  oo_.endo_simul = y';\n";
-  mDynamicModelFile << "  varargout{1} = oo_;\n";
-  mDynamicModelFile << "return;\n";
-  mDynamicModelFile << "end" << endl;
+  mDynamicModelFile << "  oo_.endo_simul = y';" << endl
+                    << "  varargout{1} = oo_;" << endl
+                    << "return;" << endl
+                    << "end" << endl;
 
   mDynamicModelFile.close();
 
@@ -3204,35 +3175,32 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
               tmp_s_eq << " " << getBlockEquationID(block, i)+1;
             }
           set<int> exogenous;
-          exogenous.clear();
           for (const auto & it : exo_block[block])
-            for (auto it1 = it.second.begin(); it1 != it.second.end(); it1++)
-              exogenous.insert(*it1);
+            for (int it1 : it.second)
+              exogenous.insert(it1);
           set<int> exogenous_det;
-          exogenous_det.clear();
           for (const auto & it : exo_det_block[block])
-            for (auto it1 = it.second.begin(); it1 != it.second.end(); it1++)
-              exogenous_det.insert(*it1);
+            for (int it1 : it.second)
+              exogenous_det.insert(it1);
           set<int> other_endogenous;
-          other_endogenous.clear();
           for (const auto & it : other_endo_block[block])
-            for (auto it1 = it.second.begin(); it1 != it.second.end(); it1++)
-              other_endogenous.insert(*it1);
-          output << "block_structure.block(" << block+1 << ").Simulation_Type = " << simulation_type << ";\n";
-          output << "block_structure.block(" << block+1 << ").maximum_lag = " << max_lag << ";\n";
-          output << "block_structure.block(" << block+1 << ").maximum_lead = " << max_lead << ";\n";
-          output << "block_structure.block(" << block+1 << ").maximum_endo_lag = " << max_lag_endo << ";\n";
-          output << "block_structure.block(" << block+1 << ").maximum_endo_lead = " << max_lead_endo << ";\n";
-          output << "block_structure.block(" << block+1 << ").maximum_exo_lag = " << max_lag_exo << ";\n";
-          output << "block_structure.block(" << block+1 << ").maximum_exo_lead = " << max_lead_exo << ";\n";
-          output << "block_structure.block(" << block+1 << ").maximum_exo_det_lag = " << max_lag_exo_det << ";\n";
-          output << "block_structure.block(" << block+1 << ").maximum_exo_det_lead = " << max_lead_exo_det << ";\n";
-          output << "block_structure.block(" << block+1 << ").endo_nbr = " << block_size << ";\n";
-          output << "block_structure.block(" << block+1 << ").mfs = " << getBlockMfs(block) << ";\n";
-          output << "block_structure.block(" << block+1 << ").equation = [" << tmp_s_eq.str() << "];\n";
-          output << "block_structure.block(" << block+1 << ").variable = [" << tmp_s.str() << "];\n";
-          output << "block_structure.block(" << block+1 << ").exo_nbr = " << getBlockExoSize(block) << ";\n";
-          output << "block_structure.block(" << block+1 << ").exogenous = [";
+            for (int it1 : it.second)
+              other_endogenous.insert(it1);
+          output << "block_structure.block(" << block+1 << ").Simulation_Type = " << simulation_type << ";" << endl
+                 << "block_structure.block(" << block+1 << ").maximum_lag = " << max_lag << ";" << endl
+                 << "block_structure.block(" << block+1 << ").maximum_lead = " << max_lead << ";" << endl
+                 << "block_structure.block(" << block+1 << ").maximum_endo_lag = " << max_lag_endo << ";" << endl
+                 << "block_structure.block(" << block+1 << ").maximum_endo_lead = " << max_lead_endo << ";" << endl
+                 << "block_structure.block(" << block+1 << ").maximum_exo_lag = " << max_lag_exo << ";" << endl
+                 << "block_structure.block(" << block+1 << ").maximum_exo_lead = " << max_lead_exo << ";" << endl
+                 << "block_structure.block(" << block+1 << ").maximum_exo_det_lag = " << max_lag_exo_det << ";" << endl
+                 << "block_structure.block(" << block+1 << ").maximum_exo_det_lead = " << max_lead_exo_det << ";" << endl
+                 << "block_structure.block(" << block+1 << ").endo_nbr = " << block_size << ";" << endl
+                 << "block_structure.block(" << block+1 << ").mfs = " << getBlockMfs(block) << ";" << endl
+                 << "block_structure.block(" << block+1 << ").equation = [" << tmp_s_eq.str() << "];" << endl
+                 << "block_structure.block(" << block+1 << ").variable = [" << tmp_s.str() << "];" << endl
+                 << "block_structure.block(" << block+1 << ").exo_nbr = " << getBlockExoSize(block) << ";" << endl
+                 << "block_structure.block(" << block+1 << ").exogenous = [";
           int i = 0;
           for (int exogenou : exogenous)
             if (exogenou >= 0)
@@ -3240,9 +3208,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
                 output << " " << exogenou+1;
                 i++;
               }
-          output << "];\n";
-
-          output << "block_structure.block(" << block+1 << ").exogenous_det = [";
+          output << "];" << endl
+                 << "block_structure.block(" << block+1 << ").exogenous_det = [";
           i = 0;
           for (int it_exogenous_det : exogenous_det)
             if (it_exogenous_det >= 0)
@@ -3250,10 +3217,9 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
                 output << " " << it_exogenous_det+1;
                 i++;
               }
-          output << "];\n";
-          output << "block_structure.block(" << block+1 << ").exo_det_nbr = " << i << ";\n";
-
-          output << "block_structure.block(" << block+1 << ").other_endogenous = [";
+          output << "];" << endl
+                 << "block_structure.block(" << block+1 << ").exo_det_nbr = " << i << ";" << endl
+                 << "block_structure.block(" << block+1 << ").other_endogenous = [";
           i = 0;
           for (int other_endogenou : other_endogenous)
             if (other_endogenou >= 0)
@@ -3261,8 +3227,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
                 output << " " << other_endogenou+1;
                 i++;
               }
-          output << "];\n";
-          output << "block_structure.block(" << block+1 << ").other_endogenous_block = [";
+          output << "];" << endl
+                 << "block_structure.block(" << block+1 << ").other_endogenous_block = [";
           i = 0;
           for (int other_endogenou : other_endogenous)
             if (other_endogenou >= 0)
@@ -3279,21 +3245,21 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
                   output << " " << j;
                 i++;
               }
-          output << "];\n";
+          output << "];" << endl;
 
           //vector<int> inter_state_var;
-          output << "block_structure.block(" << block+1 << ").tm1 = zeros(" << i << ", " << state_var.size() << ");\n";
+          output << "block_structure.block(" << block+1 << ").tm1 = zeros(" << i << ", " << state_var.size() << ");" << endl;
           int count_other_endogenous = 1;
           for (int other_endogenou : other_endogenous)
             {
-              for (vector<int>::const_iterator it = state_var.begin(); it != state_var.end(); it++)
+              for (auto it = state_var.begin(); it != state_var.end(); ++it)
                 {
                   //cout << "block = " << block+1 << " state_var = " << *it << " it_other_endogenous=" << *it_other_endogenous + 1 << "\n";
                   if (*it == other_endogenou + 1)
                     {
                       output << "block_structure.block(" << block+1 << ").tm1("
                              << count_other_endogenous << ", "
-                             << it - state_var.begin()+1 << ") = 1;\n";
+                             << it - state_var.begin()+1 << ") = 1;" << endl;
                       /*output << "block_structure.block(" << block+1 << ").tm1("
                         << it - state_var.begin()+1 << ", "
                         << count_other_endogenous << ") = 1;\n";*/
@@ -3303,14 +3269,14 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
               count_other_endogenous++;
             }
 
-          output << "block_structure.block(" << block+1 << ").other_endo_nbr = " << i << ";\n";
+          output << "block_structure.block(" << block+1 << ").other_endo_nbr = " << i << ";" << endl;
 
           tmp_s.str("");
           count_lead_lag_incidence = 0;
           dynamic_jacob_map_t reordered_dynamic_jacobian;
           for (const auto & it : blocks_derivatives[block])
-            reordered_dynamic_jacobian[{ it.second.first, { it.first.second, it.first.first } }] = it.second.second;
-          output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [];\n";
+            reordered_dynamic_jacobian[{ get<2>(it), get<1>(it), get<0>(it) }] = get<3>(it);
+          output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [];" << endl;
           int last_var = -1;
           vector<int> local_state_var;
           vector<int> local_stat_var;
@@ -3318,62 +3284,62 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
           for (int lag = -1; lag < 1+1; lag++)
             {
               last_var = -1;
-              for (dynamic_jacob_map_t::const_iterator it = reordered_dynamic_jacobian.begin(); it != reordered_dynamic_jacobian.end(); it++)
+              for (const auto &it : reordered_dynamic_jacobian)
                 {
-                  if (lag == it->first.first && last_var != it->first.second.first)
+                  if (lag == get<0>(it.first) && last_var != get<1>(it.first))
                     {
                       if (lag == -1)
                         {
-                          local_state_var.push_back(getBlockVariableID(block, it->first.second.first)+1);
+                          local_state_var.push_back(getBlockVariableID(block, get<1>(it.first))+1);
                           n_backward++;
                         }
                       else if (lag == 0)
                         {
-                          if (find(local_state_var.begin(), local_state_var.end(), getBlockVariableID(block, it->first.second.first)+1) == local_state_var.end())
+                          if (find(local_state_var.begin(), local_state_var.end(), getBlockVariableID(block, get<1>(it.first))+1) == local_state_var.end())
                             {
-                              local_stat_var.push_back(getBlockVariableID(block, it->first.second.first)+1);
+                              local_stat_var.push_back(getBlockVariableID(block, get<1>(it.first))+1);
                               n_static++;
                             }
                         }
                       else
                         {
-                          if (find(local_state_var.begin(), local_state_var.end(), getBlockVariableID(block, it->first.second.first)+1) != local_state_var.end())
+                          if (find(local_state_var.begin(), local_state_var.end(), getBlockVariableID(block, get<1>(it.first))+1) != local_state_var.end())
                             {
                               n_backward--;
                               n_mixed++;
                             }
                           else
                             {
-                              if (find(local_stat_var.begin(), local_stat_var.end(), getBlockVariableID(block, it->first.second.first)+1) != local_stat_var.end())
+                              if (find(local_stat_var.begin(), local_stat_var.end(), getBlockVariableID(block, get<1>(it.first))+1) != local_stat_var.end())
                                 n_static--;
                               n_forward++;
                             }
                         }
                       count_lead_lag_incidence++;
-                      for (int i = last_var; i < it->first.second.first-1; i++)
+                      for (int i = last_var; i < get<1>(it.first)-1; i++)
                         tmp_s << " 0";
                       if (tmp_s.str().length())
                         tmp_s << " ";
                       tmp_s << count_lead_lag_incidence;
-                      last_var = it->first.second.first;
+                      last_var = get<1>(it.first);
                     }
                 }
               for (int i = last_var + 1; i < block_size; i++)
                 tmp_s << " 0";
-              output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [ block_structure.block(" << block+1 << ").lead_lag_incidence; " << tmp_s.str() << "]; %lag = " << lag << "\n";
+              output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [ block_structure.block(" << block+1 << ").lead_lag_incidence; " << tmp_s.str() << "]; %lag = " << lag << endl;
               tmp_s.str("");
             }
           vector<int> inter_state_var;
-          for (vector<int>::const_iterator it_l = local_state_var.begin(); it_l != local_state_var.end(); it_l++)
-            for (vector<int>::const_iterator it = state_var.begin(); it != state_var.end(); it++)
+          for (auto it_l = local_state_var.begin(); it_l != local_state_var.end(); ++it_l)
+            for (auto it = state_var.begin(); it != state_var.end(); ++it)
               if (*it == *it_l)
                 inter_state_var.push_back(it - state_var.begin()+1);
           output << "block_structure.block(" << block+1 << ").sorted_col_dr_ghx = [";
-          for (vector<int>::const_iterator it = inter_state_var.begin(); it != inter_state_var.end(); it++)
-            output << *it << " ";
-          output << "];\n";
+          for (int it : inter_state_var)
+            output << it << " ";
+          output << "];" << endl;
           count_lead_lag_incidence = 0;
-          output << "block_structure.block(" << block+1 << ").lead_lag_incidence_other = [];\n";
+          output << "block_structure.block(" << block+1 << ").lead_lag_incidence_other = [];" << endl;
           for (int lag = -1; lag <= 1; lag++)
             {
               tmp_s.str("");
@@ -3383,7 +3349,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
                   for (int i = 0; i < block_size; i++)
                     {
                       unsigned int eq = getBlockEquationID(block, i);
-                      auto it = derivative_other_endo[block].find({ lag, { eq, other_endogenou } });
+                      auto it = derivative_other_endo[block].find({ lag, eq, other_endogenou });
                       if (it != derivative_other_endo[block].end())
                         {
                           count_lead_lag_incidence++;
@@ -3395,33 +3361,33 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
                   if (!done)
                     tmp_s << " 0";
                 }
-              output << "block_structure.block(" << block+1 << ").lead_lag_incidence_other = [ block_structure.block(" << block+1 << ").lead_lag_incidence_other; " << tmp_s.str() << "]; %lag = " << lag << "\n";
+              output << "block_structure.block(" << block+1 << ").lead_lag_incidence_other = [ block_structure.block(" << block+1 << ").lead_lag_incidence_other; " << tmp_s.str() << "]; %lag = " << lag << endl;
             }
-          output << "block_structure.block(" << block+1 << ").n_static = " << n_static << ";\n";
-          output << "block_structure.block(" << block+1 << ").n_forward = " << n_forward << ";\n";
-          output << "block_structure.block(" << block+1 << ").n_backward = " << n_backward << ";\n";
-          output << "block_structure.block(" << block+1 << ").n_mixed = " << n_mixed << ";\n";
+          output << "block_structure.block(" << block+1 << ").n_static = " << n_static << ";" << endl
+                 << "block_structure.block(" << block+1 << ").n_forward = " << n_forward << ";" << endl
+                 << "block_structure.block(" << block+1 << ").n_backward = " << n_backward << ";" << endl
+                 << "block_structure.block(" << block+1 << ").n_mixed = " << n_mixed << ";" << endl;
         }
-      output << modstruct << "block_structure.block = block_structure.block;\n";
+      output << modstruct << "block_structure.block = block_structure.block;" << endl;
       string cst_s;
       int nb_endo = symbol_table.endo_nbr();
       output << modstruct << "block_structure.variable_reordered = [";
       for (int i = 0; i < nb_endo; i++)
         output << " " << variable_reordered[i]+1;
-      output << "];\n";
+      output << "];" << endl;
       output << modstruct << "block_structure.equation_reordered = [";
       for (int i = 0; i < nb_endo; i++)
         output << " " << equation_reordered[i]+1;
-      output << "];\n";
+      output << "];" << endl;
       vector<int> variable_inv_reordered(nb_endo);
 
       for (int i = 0; i < nb_endo; i++)
         variable_inv_reordered[variable_reordered[i]] = i;
 
-      for (vector<int>::const_iterator it = state_var.begin(); it != state_var.end(); it++)
-        state_equ.push_back(equation_reordered[variable_inv_reordered[*it - 1]]+1);
+      for (int it : state_var)
+        state_equ.push_back(equation_reordered[variable_inv_reordered[it - 1]]+1);
 
-      map<pair< int, pair<int, int>>,  int>  lag_row_incidence;
+      map<tuple<int, int, int>,  int> lag_row_incidence;
       for (const auto & first_derivative : derivatives[1])
         {
           int deriv_id = first_derivative.first[1];
@@ -3431,23 +3397,23 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
               int symb = getSymbIDByDerivID(deriv_id);
               int var = symbol_table.getTypeSpecificID(symb);
               int lag = getLagByDerivID(deriv_id);
-              lag_row_incidence[{ lag, { eq, var } }] = 1;
+              lag_row_incidence[{ lag, eq, var }] = 1;
             }
         }
       int prev_lag = -1000000;
-      for (map<pair< int, pair<int, int>>,  int>::const_iterator it = lag_row_incidence.begin(); it != lag_row_incidence.end(); it++)
+      for (const auto &it : lag_row_incidence)
         {
-          if (prev_lag != it->first.first)
+          if (prev_lag != get<0>(it.first))
             {
               if (prev_lag != -1000000)
-                output << "];\n";
-              prev_lag = it->first.first;
-              output << modstruct << "block_structure.incidence(" << max_endo_lag+it->first.first+1 << ").lead_lag = " << prev_lag << ";\n";
-              output << modstruct << "block_structure.incidence(" << max_endo_lag+it->first.first+1 << ").sparse_IM = [";
+                output << "];" << endl;
+              prev_lag = get<0>(it.first);
+              output << modstruct << "block_structure.incidence(" << max_endo_lag+get<0>(it.first)+1 << ").lead_lag = " << prev_lag << ";" << endl
+                     << modstruct << "block_structure.incidence(" << max_endo_lag+get<0>(it.first)+1 << ").sparse_IM = [";
             }
-          output << it->first.second.first+1 << " " << it->first.second.second+1 << ";\n";
+          output << get<1>(it.first)+1 << " " << get<2>(it.first)+1 << ";" << endl;
         }
-      output << "];\n";
+      output << "];" << endl;
       if (estimation_present)
         {
           ofstream KF_index_file;
@@ -3456,8 +3422,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
           KF_index_file.open(main_name, ios::out | ios::binary | ios::ate);
           int n_obs = symbol_table.observedVariablesNbr();
           int n_state = state_var.size();
-          for (vector<int>::const_iterator it = state_var.begin(); it != state_var.end(); it++)
-            if (symbol_table.isObservedVariable(symbol_table.getID(SymbolType::endogenous, *it-1)))
+          for (int it : state_var)
+            if (symbol_table.isObservedVariable(symbol_table.getID(SymbolType::endogenous, it-1)))
               n_obs--;
 
           int n = n_obs + n_state;
@@ -3478,19 +3444,19 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
               for (int i = 0; i < block_size; i++)
                 {
                   int var = getBlockVariableID(block, i);
-                  vector<int>::const_iterator it_state_var = find(state_var.begin(), state_var.end(), var+1);
+                  auto it_state_var = find(state_var.begin(), state_var.end(), var+1);
                   if (it_state_var != state_var.end())
                     nze++;
                 }
               if (block == 0)
                 {
                   set<pair<int, int>> row_state_var_incidence;
-                  for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+                  for (const auto &it : blocks_derivatives[block])
                     {
-                      vector<int>::const_iterator it_state_var = find(state_var.begin(), state_var.end(), getBlockVariableID(block, it->first.second)+1);
+                      auto it_state_var = find(state_var.begin(), state_var.end(), getBlockVariableID(block, get<1>(it))+1);
                       if (it_state_var != state_var.end())
                         {
-                          vector<int>::const_iterator it_state_equ = find(state_equ.begin(), state_equ.end(), getBlockEquationID(block, it->first.first)+1);
+                          auto it_state_equ = find(state_equ.begin(), state_equ.end(), getBlockEquationID(block, get<0>(it))+1);
                           if (it_state_equ != state_equ.end())
                             row_state_var_incidence.emplace(it_state_equ - state_equ.begin(), it_state_var - state_var.begin());
                         }
@@ -3597,8 +3563,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
     }
 
   output << modstruct << "state_var = [";
-  for (vector<int>::const_iterator it=state_var.begin(); it != state_var.end(); it++)
-    output << *it << (julia ? "," : " ");
+  for (int it : state_var)
+    output << it << (julia ? "," : " ");
   output << "];" << endl;
 
   // Writing initialization for some other variables
@@ -3651,10 +3617,10 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
     it->ExprNode::writeOutput(output, ExprNodeOutputType::matlabDynamicModel);
 }
 
-map<pair<int, pair<int, int >>, expr_t>
+map<tuple<int, int, int>, expr_t>
 DynamicModel::collect_first_order_derivatives_endogenous()
 {
-  map<pair<int, pair<int, int >>, expr_t> endo_derivatives;
+  map<tuple<int, int, int>, expr_t> endo_derivatives;
   for (auto & first_derivative : derivatives[1])
     {
       if (getTypeByDerivID(first_derivative.first[1]) == SymbolType::endogenous)
@@ -3662,7 +3628,7 @@ DynamicModel::collect_first_order_derivatives_endogenous()
           int eq = first_derivative.first[0];
           int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first[1]));
           int lag = getLagByDerivID(first_derivative.first[1]);
-          endo_derivatives[{ eq, { var, lag } }] = first_derivative.second;
+          endo_derivatives[{ eq, var, lag }] = first_derivative.second;
         }
     }
   return endo_derivatives;
@@ -4406,14 +4372,14 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
     }
 
   jacob_map_t contemporaneous_jacobian, static_jacobian;
-  map<pair<int, pair<int, int> >, expr_t> first_order_endo_derivatives;
+  map<tuple<int, int, int>, expr_t> first_order_endo_derivatives;
   // for each block contains pair<Size, Feddback_variable>
   vector<pair<int, int> > blocks;
   vector<unsigned int> n_static, n_forward, n_backward, n_mixed;
 
   if (linear_decomposition)
     {
-      map<pair<int, pair<int, int>>, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
+      first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
       is_equation_linear = equationLinear(first_order_endo_derivatives);
 
       evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false);
@@ -4462,7 +4428,7 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
       equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
 
       if (!nopreprocessoroutput)
-        cout << "Finding the optimal block decomposition of the model ...\n";
+        cout << "Finding the optimal block decomposition of the model ..." << endl;
 
       lag_lead_vector_t equation_lag_lead, variable_lag_lead;
 
@@ -4484,15 +4450,15 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
       if (!no_tmp_terms)
         computeTemporaryTermsOrdered();
       int k = 0;
-      equation_block = vector<int>(equations.size());
-      variable_block_lead_lag = vector< pair< int, pair< int, int>>>(equations.size());
+      equation_block.resize(equations.size());
+      variable_block_lead_lag = vector<tuple<int, int, int>>(equations.size());
       for (unsigned int i = 0; i < getNbBlocks(); i++)
         {
           for (unsigned int j = 0; j < getBlockSize(i); j++)
             {
               equation_block[equation_reordered[k]] = i;
               int l = variable_reordered[k];
-              variable_block_lead_lag[l] = { i, { variable_lag_lead[l].first, variable_lag_lead[l].second } };
+              variable_block_lead_lag[l] = { i, variable_lag_lead[l].first, variable_lag_lead[l].second };
               k++;
             }
         }
@@ -4610,12 +4576,11 @@ DynamicModel::writeRevXrefs(ostream &output, const map<pair<int, int>, set<int>>
     }
 }
 
-map<pair<pair<int, pair<int, int>>, pair<int, int>>, int>
+map<tuple<int, int, int, int, int>, int>
 DynamicModel::get_Derivatives(int block)
 {
   int max_lag, max_lead;
-  map<pair<pair<int, pair<int, int>>, pair<int, int>>, int> Derivatives;
-  Derivatives.clear();
+  map<tuple<int, int, int, int, int>, int> Derivatives;
   BlockSimulationType simulation_type = getBlockSimulationType(block);
   if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
     {
@@ -4638,10 +4603,10 @@ DynamicModel::get_Derivatives(int block)
           for (int var = 0; var < block_size; var++)
             {
               int varr = getBlockVariableID(block, var);
-              if (dynamic_jacobian.find({ lag, { eqr, varr } }) != dynamic_jacobian.end())
+              if (dynamic_jacobian.find({ lag, eqr, varr }) != dynamic_jacobian.end())
                 {
                   bool OK = true;
-                  map<pair<pair<int, pair<int, int>>, pair<int, int>>, int>::const_iterator its = Derivatives.find({ { lag, { eq, var } }, { eqr, varr } });
+                  auto its = Derivatives.find({ lag, eq, var, eqr, varr });
                   if (its != Derivatives.end())
                     {
                       if (its->second == 2)
@@ -4652,10 +4617,10 @@ DynamicModel::get_Derivatives(int block)
                     {
                       if (getBlockEquationType(block, eq) == E_EVALUATE_S && eq < block_nb_recursive)
                         //It's a normalized equation, we have to recompute the derivative using chain rule derivative function
-                        Derivatives[{ { lag, { eq, var } }, { eqr, varr } }] = 1;
+                        Derivatives[{ lag, eq, var, eqr, varr }] = 1;
                       else
                         //It's a feedback equation we can use the derivatives
-                        Derivatives[{ { lag, { eq, var } }, { eqr, varr } }] = 0;
+                        Derivatives[{ lag, eq, var, eqr, varr }] = 0;
                     }
                   if (var < block_nb_recursive)
                     {
@@ -4664,15 +4629,15 @@ DynamicModel::get_Derivatives(int block)
                         {
                           int varrs = getBlockVariableID(block, vars);
                           //A new derivative needs to be computed using the chain rule derivative function (a feedback variable appears in a recursive equation)
-                          if (Derivatives.find({ { lag, { var, vars } }, { eqs, varrs } }) != Derivatives.end())
-                            Derivatives[{ { lag, { eq, vars } }, { eqr, varrs } }] = 2;
+                          if (Derivatives.find({ lag, var, vars, eqs, varrs }) != Derivatives.end())
+                            Derivatives[{ lag, eq, vars, eqr, varrs }] = 2;
                         }
                     }
                 }
             }
         }
     }
-  return (Derivatives);
+  return Derivatives;
 }
 
 void
@@ -4696,30 +4661,24 @@ DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_endo_derivat
           else
             recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
         }
-      map<pair<pair<int, pair<int, int>>, pair<int, int>>, int> Derivatives = get_Derivatives(block);
-      map<pair<pair<int, pair<int, int>>, pair<int, int>>, int>::const_iterator it = Derivatives.begin();
-      for (int i = 0; i < (int) Derivatives.size(); i++)
+      auto Derivatives = get_Derivatives(block);
+      for (const auto &it : Derivatives)
         {
-          int Deriv_type = it->second;
-          pair<pair<int, pair<int, int>>, pair<int, int>> it_l(it->first);
-          it++;
-          int lag = it_l.first.first;
-          int eq = it_l.first.second.first;
-          int var = it_l.first.second.second;
-          int eqr = it_l.second.first;
-          int varr = it_l.second.second;
+          int Deriv_type = it.second;
+          int lag, eq, var, eqr, varr;
+          tie(lag, eq, var, eqr, varr) = it.first;
           if (Deriv_type == 0)
-            first_chain_rule_derivatives[{ eqr, { varr, lag } }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
+            first_chain_rule_derivatives[{ eqr, varr, lag }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
           else if (Deriv_type == 1)
-            first_chain_rule_derivatives[{ eqr, { varr, lag } }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
+            first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
           else if (Deriv_type == 2)
             {
               if (getBlockEquationType(block, eq) == E_EVALUATE_S && eq < block_nb_recursives)
-                first_chain_rule_derivatives[{ eqr, { varr, lag } }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
+                first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
               else
-                first_chain_rule_derivatives[{ eqr, { varr, lag } }] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
+                first_chain_rule_derivatives[{ eqr, varr, lag }] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
             }
-          tmp_derivatives.emplace_back(make_pair(eq, var), make_pair(lag, first_chain_rule_derivatives[{ eqr, { varr, lag } }]));
+          tmp_derivatives.emplace_back(eq, var, lag, first_chain_rule_derivatives[{ eqr, varr, lag }]);
         }
       blocks_endo_derivatives[block] = tmp_derivatives;
     }
@@ -4729,10 +4688,8 @@ void
 DynamicModel::collect_block_first_order_derivatives()
 {
   //! vector for an equation or a variable indicates the block number
-  vector<int> equation_2_block, variable_2_block;
+  vector<int> equation_2_block(equation_reordered.size()), variable_2_block(variable_reordered.size());
   unsigned int nb_blocks = getNbBlocks();
-  equation_2_block = vector<int>(equation_reordered.size());
-  variable_2_block = vector<int>(variable_reordered.size());
   for (unsigned int block = 0; block < nb_blocks; block++)
     {
       unsigned int block_size = getBlockSize(block);
@@ -4774,7 +4731,7 @@ DynamicModel::collect_block_first_order_derivatives()
               if (lag > 0 && lag > endo_max_leadlag_block[block_eq].second)
                 endo_max_leadlag_block[block_eq] = { endo_max_leadlag_block[block_eq].first, lag };
               tmp_derivative = derivative_endo[block_eq];
-              tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
+              tmp_derivative[{ lag, eq, var }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
               derivative_endo[block_eq] = tmp_derivative;
             }
           else
@@ -4785,7 +4742,7 @@ DynamicModel::collect_block_first_order_derivatives()
                 other_endo_max_leadlag_block[block_eq] = { other_endo_max_leadlag_block[block_eq].first, lag };
               tmp_derivative = derivative_other_endo[block_eq];
               {
-                map< int, map<int, int>>::const_iterator it = block_other_endo_index.find(block_eq);
+                auto it = block_other_endo_index.find(block_eq);
                 if (it == block_other_endo_index.end())
                   block_other_endo_index[block_eq][var] = 0;
                 else
@@ -4798,7 +4755,7 @@ DynamicModel::collect_block_first_order_derivatives()
                       }
                   }
               }
-              tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
+              tmp_derivative[{ lag, eq, var }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
               derivative_other_endo[block_eq] = tmp_derivative;
               lag_var = other_endo_block[block_eq];
               if (lag_var.find(lag) == lag_var.end())
@@ -4814,7 +4771,7 @@ DynamicModel::collect_block_first_order_derivatives()
             exo_max_leadlag_block[block_eq] = { exo_max_leadlag_block[block_eq].first, lag };
           tmp_derivative = derivative_exo[block_eq];
           {
-            map< int, map<int, int>>::const_iterator it = block_exo_index.find(block_eq);
+            auto it = block_exo_index.find(block_eq);
             if (it == block_exo_index.end())
               block_exo_index[block_eq][var] = 0;
             else
@@ -4827,7 +4784,7 @@ DynamicModel::collect_block_first_order_derivatives()
                   }
               }
           }
-          tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }];
+          tmp_derivative[{ lag, eq, var }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }];
           derivative_exo[block_eq] = tmp_derivative;
           lag_var = exo_block[block_eq];
           if (lag_var.find(lag) == lag_var.end())
@@ -4842,7 +4799,7 @@ DynamicModel::collect_block_first_order_derivatives()
             exo_det_max_leadlag_block[block_eq] = { exo_det_max_leadlag_block[block_eq].first, lag };
           tmp_derivative = derivative_exo_det[block_eq];
           {
-            map< int, map<int, int>>::const_iterator it = block_det_exo_index.find(block_eq);
+            auto it = block_det_exo_index.find(block_eq);
             if (it == block_det_exo_index.end())
               block_det_exo_index[block_eq][var] = 0;
             else
@@ -4855,7 +4812,7 @@ DynamicModel::collect_block_first_order_derivatives()
                   }
               }
           }
-          tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }];
+          tmp_derivative[{ lag, eq, var }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::exogenous, var), lag) }];
           derivative_exo_det[block_eq] = tmp_derivative;
           lag_var = exo_det_block[block_eq];
           if (lag_var.find(lag) == lag_var.end())
@@ -4883,10 +4840,10 @@ DynamicModel::collectBlockVariables()
       int prev_lag = -999999999;
       int count_col_exo = 0;
       var_t tmp_var_exo;
-      for (lag_var_t::const_iterator it = exo_block[block].begin(); it != exo_block[block].end(); it++)
+      for (const auto &it : exo_block[block])
         {
-          int lag = it->first;
-          for (int var : it->second)
+          int lag = it.first;
+          for (int var : it.second)
             {
               tmp_var_exo.insert(var);
               if (prev_var != var || prev_lag != lag)
@@ -5070,9 +5027,8 @@ void
 DynamicModel::toNonlinearPart(DynamicModel &non_linear_equations_dynamic_model) const
 {
   // Convert model local variables (need to be done first)
-  for (map<int, expr_t>::const_iterator it = local_variables_table.begin();
-       it != local_variables_table.end(); it++)
-    non_linear_equations_dynamic_model.AddLocalVariable(it->first, it->second);
+  for (const auto & it : local_variables_table)
+    non_linear_equations_dynamic_model.AddLocalVariable(it.first, it.second);
 }
 
 bool
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index d84914b51527b53ff6f0e408f1390efaf0e7dc9b..dedafa40cc94934bbdf249f917917ba752d9ecba 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -99,8 +99,8 @@ private:
 
   vector<temporary_terms_inuse_t> v_temporary_terms_inuse;
 
-  //! Store the derivatives or the chainrule derivatives:map<pair< equation, pair< variable, lead_lag >, expr_t>
-  using first_chain_rule_derivatives_t = map< pair< int, pair< int, int>>, expr_t>;
+  //! Store the derivatives or the chainrule derivatives:map<tuple<equation, variable, lead_lag>, expr_t>
+  using first_chain_rule_derivatives_t = map<tuple<int, int, int>, expr_t>;
   first_chain_rule_derivatives_t first_chain_rule_derivatives;
 
   //! Writes dynamic model file (Matlab version)
@@ -135,7 +135,7 @@ private:
   //void evaluateJacobian(const eval_context_t &eval_context, jacob_map *j_m, bool dynamic);
 
   //! return a map on the block jacobian
-  map<pair<pair<int, pair<int, int>>, pair<int, int>>, int> get_Derivatives(int block);
+  map<tuple<int, int, int, int, int>, int> get_Derivatives(int block);
   //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
   void computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives);
 
@@ -163,7 +163,7 @@ private:
   //! Computes derivatives of the Jacobian w.r. to trend vars and tests that they are equal to zero
   void testTrendDerivativesEqualToZero(const eval_context_t &eval_context);
   //! Collect only the first derivatives
-  map<pair<int, pair<int, int>>, expr_t> collect_first_order_derivatives_endogenous();
+  map<tuple<int, int, int>, expr_t> collect_first_order_derivatives_endogenous();
 
   //! Allocates the derivation IDs for all dynamic variables of the model
   /*! Also computes max_{endo,exo}_{lead_lag}, and initializes dynJacobianColsNbr to the number of dynamic endos */
@@ -200,8 +200,8 @@ private:
   //! Vector indicating if the block is linear in endogenous variable (true) or not (false)
   vector<bool> blocks_linear;
 
-  //! Map the derivatives for a block pair<lag, make_pair(make_pair(eq, var)), expr_t>
-  using derivative_t = map<pair< int, pair<int, int>>, expr_t>;
+  //! Map the derivatives for a block tuple<lag, eq, var>
+  using derivative_t = map<tuple<int, int, int>, expr_t>;
   //! Vector of derivative for each blocks
   vector<derivative_t> derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
 
@@ -216,8 +216,8 @@ private:
   map< int, map<int, int>> block_exo_index, block_det_exo_index, block_other_endo_index;
 
   //! for each block described the number of static, forward, backward and mixed variables in the block
-  /*! pair< pair<static, forward>, pair<backward,mixed>> */
-  vector<pair< pair<int, int>, pair<int, int>>> block_col_type;
+  /*! tuple<static, forward, backward, mixed> */
+  vector<tuple<int, int, int, int>> block_col_type;
 
   //! Help computeXrefs to compute the reverse references (i.e. param->eqs, endo->eqs, etc)
   void computeRevXref(map<pair<int, int>, set<int>> &xrefset, const set<pair<int, int>> &eiref, int eqn);
@@ -226,7 +226,7 @@ private:
   void writeRevXrefs(ostream &output, const map<pair<int, int>, set<int>> &xrefmap, const string &type) const;
 
   //! List for each variable its block number and its maximum lag and lead inside the block
-  vector<pair<int, pair<int, int>>> variable_block_lead_lag;
+  vector<tuple<int, int, int>> variable_block_lead_lag;
   //! List for each equation its block number
   vector<int> equation_block;
 
@@ -481,19 +481,19 @@ public:
   BlockSimulationType
   getBlockSimulationType(int block_number) const override
   {
-    return (block_type_firstequation_size_mfs[block_number].first.first);
+    return (get<0>(block_type_firstequation_size_mfs[block_number]));
   };
   //! Return the first equation number of a block
   unsigned int
   getBlockFirstEquation(int block_number) const override
   {
-    return (block_type_firstequation_size_mfs[block_number].first.second);
+    return (get<1>(block_type_firstequation_size_mfs[block_number]));
   };
   //! Return the size of the block block_number
   unsigned int
   getBlockSize(int block_number) const override
   {
-    return (block_type_firstequation_size_mfs[block_number].second.first);
+    return (get<2>(block_type_firstequation_size_mfs[block_number]));
   };
   //! Return the number of exogenous variable in the block block_number
   unsigned int
@@ -511,7 +511,7 @@ public:
   unsigned int
   getBlockMfs(int block_number) const override
   {
-    return (block_type_firstequation_size_mfs[block_number].second.second);
+    return (get<3>(block_type_firstequation_size_mfs[block_number]));
   };
   //! Return the maximum lag in a block
   unsigned int
@@ -529,37 +529,37 @@ public:
   EquationType
   getBlockEquationType(int block_number, int equation_number) const override
   {
-    return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first);
+    return (equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first);
   };
   //! Return true if the equation has been normalized
   bool
   isBlockEquationRenormalized(int block_number, int equation_number) const override
   {
-    return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first == E_EVALUATE_S);
+    return (equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first == E_EVALUATE_S);
   };
   //! Return the expr_t of the equation equation_number belonging to the block block_number
   expr_t
   getBlockEquationExpr(int block_number, int equation_number) const override
   {
-    return (equations[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]]);
+    return (equations[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]]);
   };
   //! Return the expr_t of the renormalized equation equation_number belonging to the block block_number
   expr_t
   getBlockEquationRenormalizedExpr(int block_number, int equation_number) const override
   {
-    return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].second);
+    return (equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].second);
   };
   //! Return the original number of equation equation_number belonging to the block block_number
   int
   getBlockEquationID(int block_number, int equation_number) const override
   {
-    return (equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]);
+    return (equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]);
   };
   //! Return the original number of variable variable_number belonging to the block block_number
   int
   getBlockVariableID(int block_number, int variable_number) const override
   {
-    return (variable_reordered[block_type_firstequation_size_mfs[block_number].first.second+variable_number]);
+    return (variable_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+variable_number]);
   };
   //! Return the original number of the exogenous variable varexo_number belonging to the block block_number
   int
@@ -572,19 +572,19 @@ public:
   int
   getBlockInitialEquationID(int block_number, int equation_number) const override
   {
-    return ((int) inv_equation_reordered[equation_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
+    return ((int) inv_equation_reordered[equation_number] - (int) get<1>(block_type_firstequation_size_mfs[block_number]));
   };
   //! Return the position of variable_number in the block number belonging to the block block_number
   int
   getBlockInitialVariableID(int block_number, int variable_number) const override
   {
-    return ((int) inv_variable_reordered[variable_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
+    return ((int) inv_variable_reordered[variable_number] - (int) get<1>(block_type_firstequation_size_mfs[block_number]));
   };
   //! Return the block number containing the endogenous variable variable_number
   int
   getBlockVariableID(int variable_number) const
   {
-    return (variable_block_lead_lag[variable_number].first);
+    return (get<0>(variable_block_lead_lag[variable_number]));
   };
   //! Return the position of the exogenous variable_number in the block number belonging to the block block_number
   int
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index dcdb505966a35cabf2958cb74c007604cb98c7f0..2bd75d68a489ddb21ff754f6d896f864e1ebaf3f 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -185,7 +185,7 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
 #if 1
   bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]);
 #else // Alternative way to compute normalization, by giving an initial matching using natural normalizations
-  fill(mate_map.begin(), mate_map.end(), graph_traits<BipartiteGraph>::null_vertex());
+  fill(mate_map.begin(), mate_map.end(), boost::graph_traits<BipartiteGraph>::null_vertex());
 
   multimap<int, int> natural_endo2eqs;
   computeNormalizedEquations(natural_endo2eqs);
@@ -201,15 +201,14 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
       put(&mate_map[0], n+j, i);
     }
 
-  edmonds_augmenting_path_finder<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type> augmentor(g, &mate_map[0], get(vertex_index, g));
-  bool not_maximum_yet = true;
-  while (not_maximum_yet)
+  boost::edmonds_augmenting_path_finder<BipartiteGraph, int *, boost::property_map<BipartiteGraph, boost::vertex_index_t>::type> augmentor(g, &mate_map[0], get(boost::vertex_index, g));
+  while (augmentor.augment_matching())
     {
-      not_maximum_yet = augmentor.augment_matching();
-    }
+    };
+
   augmentor.get_current_matching(&mate_map[0]);
 
-  bool check = maximum_cardinality_matching_verifier<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type>::verify_matching(g, &mate_map[0], get(vertex_index, g));
+  bool check = boost::maximum_cardinality_matching_verifier<BipartiteGraph, int *, boost::property_map<BipartiteGraph, boost::vertex_index_t>::type>::verify_matching(g, &mate_map[0], get(boost::vertex_index, g));
 #endif
 
   assert(check);
@@ -222,7 +221,7 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
 
   // Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
   endo2eq.resize(equations.size());
-  transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus<int>(), n));
+  transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), [=](int i) { return i-n; });
 
 #ifdef DEBUG
   multimap<int, int> natural_endo2eqs;
@@ -237,9 +236,9 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
 
       n1++;
 
-      pair<multimap<int, int>::const_iterator, multimap<int, int>::const_iterator> x = natural_endo2eqs.equal_range(i);
-      if (find_if(x.first, x.second, compose1(bind2nd(equal_to<int>(), endo2eq[i]), select2nd<multimap<int, int>::value_type>())) == x.second)
-        cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
+      auto x = natural_endo2eqs.equal_range(i);
+      if (find_if(x.first, x.second, [=](auto y) { return y.second == endo2eq[i]; }) == x.second)
+        cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, i))
              << " not used." << endl;
       else
         n2++;
@@ -274,9 +273,9 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
   //jacob_map normalized_contemporaneous_jacobian;
   jacob_map_t normalized_contemporaneous_jacobian(contemporaneous_jacobian);
   vector<double> max_val(n, 0.0);
-  for (jacob_map_t::const_iterator iter = contemporaneous_jacobian.begin(); iter != contemporaneous_jacobian.end(); iter++)
-    if (fabs(iter->second) > max_val[iter->first.first])
-      max_val[iter->first.first] = fabs(iter->second);
+  for (const auto &it : contemporaneous_jacobian)
+    if (fabs(it.second) > max_val[it.first.first])
+      max_val[it.first.first] = fabs(it.second);
 
   for (auto & iter : normalized_contemporaneous_jacobian)
     iter.second /= max_val[iter.first.first];
@@ -324,22 +323,22 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
       if (check)
         {
           // Update the jacobian matrix
-          for (jacob_map_t::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
+          for (const auto &it : tmp_normalized_contemporaneous_jacobian)
             {
-              if (static_jacobian.find({ it->first.first, it->first.second }) == static_jacobian.end())
-                static_jacobian[{ it->first.first, it->first.second }] = 0;
-              if (dynamic_jacobian.find({ 0, { it->first.first, it->first.second } }) == dynamic_jacobian.end())
-                dynamic_jacobian[{ 0, { it->first.first, it->first.second } }] = nullptr;
-              if (contemporaneous_jacobian.find({ it->first.first, it->first.second }) == contemporaneous_jacobian.end())
-                contemporaneous_jacobian[{ it->first.first, it->first.second }] = 0;
+              if (static_jacobian.find({ it.first.first, it.first.second }) == static_jacobian.end())
+                static_jacobian[{ it.first.first, it.first.second }] = 0;
+              if (dynamic_jacobian.find({ 0, it.first.first, it.first.second }) == dynamic_jacobian.end())
+                dynamic_jacobian[{ 0, it.first.first, it.first.second }] = nullptr;
+              if (contemporaneous_jacobian.find({ it.first.first, it.first.second }) == contemporaneous_jacobian.end())
+                contemporaneous_jacobian[{ it.first.first, it.first.second }] = 0;
               try
                 {
-                  if (derivatives[1].find({ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }) == derivatives[1].end())
-                    derivatives[1][{ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }] = Zero;
+                  if (derivatives[1].find({ it.first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it.first.second), 0) }) == derivatives[1].end())
+                    derivatives[1][{ it.first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it.first.second), 0) }] = Zero;
                 }
               catch (DataTree::UnknownDerivIDException &e)
                 {
-                  cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it->first.second))
+                  cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it.first.second))
                        << " does not appear at the current period (i.e. with no lead and no lag); this case is not handled by the 'block' option of the 'model' block." << endl;
                   exit(EXIT_FAILURE);
                 }
@@ -382,14 +381,13 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m
 {
   int nb_elements_contemparenous_Jacobian = 0;
   set<vector<int>> jacobian_elements_to_delete;
-  for (auto it = derivatives[1].begin();
-       it != derivatives[1].end(); it++)
+  for (const auto &it : derivatives[1])
     {
-      int deriv_id = it->first[1];
+      int deriv_id = it.first[1];
       if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
         {
-          expr_t Id = it->second;
-          int eq = it->first[0];
+          expr_t Id = it.second;
+          int eq = it.first[0];
           int symb = getSymbIDByDerivID(deriv_id);
           int var = symbol_table.getTypeSpecificID(symb);
           int lag = getLagByDerivID(deriv_id);
@@ -426,7 +424,7 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m
                 static_jacobian[{ eq, var }] += val;
               else
                 static_jacobian[{ eq, var }] = val;
-              dynamic_jacobian[{ lag, { eq, var } }] = Id;
+              dynamic_jacobian[{ lag, eq, var }] = Id;
             }
         }
     }
@@ -452,17 +450,17 @@ ModelTree::select_non_linear_equations_and_variables(vector<bool> is_equation_li
   /*equation_reordered.resize(equations.size());
   variable_reordered.resize(equations.size());*/
   unsigned int num = 0;
-  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++)
-    if (!is_equation_linear[*it])
+  for (auto it : endo2eq)
+    if (!is_equation_linear[it])
       num++;
   vector<int> endo2block = vector<int>(endo2eq.size(), 1);
   vector<pair<set<int>, pair<set<int>, vector<int> > > > components_set(num);
   int i = 0, j = 0;
-  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, j++)
+  for (auto it : endo2eq)
     {
-      if (!is_equation_linear[*it])
+      if (!is_equation_linear[it])
         {
-          equation_reordered[i] = *it;
+          equation_reordered[i] = it;
           variable_reordered[i] = j;
           endo2block[j] = 0;
           components_set[endo2block[j]].first.insert(i);
@@ -477,6 +475,7 @@ ModelTree::select_non_linear_equations_and_variables(vector<bool> is_equation_li
           cout.flush();
           */
           i++;
+          j++;
         }
     }
 /*  for (unsigned int j = 0; j < is_equation_linear.size() ; j++)
@@ -501,9 +500,9 @@ ModelTree::select_non_linear_equations_and_variables(vector<bool> is_equation_li
     }
   cout.flush();
   int nb_endo = is_equation_linear.size();
-  vector<pair<int, int> > blocks = vector<pair<int, int> >(1, make_pair(i, i));
-  inv_equation_reordered = vector<int>(nb_endo);
-  inv_variable_reordered = vector<int>(nb_endo);
+  vector<pair<int, int>> blocks(1, make_pair(i, i));
+  inv_equation_reordered.resize(nb_endo);
+  inv_variable_reordered.resize(nb_endo);
   for (int i = 0; i < nb_endo; i++)
     {
       inv_variable_reordered[variable_reordered[i]] = i;
@@ -528,7 +527,7 @@ ModelTree::computeNaturalNormalization()
         if (result.size() == 1 && result.begin()->second == 0)
           {
             //check if the endogenous variable has not been already used in an other match !
-             vector<int>::iterator it = find(endo2eq.begin(), endo2eq.end(), result.begin()->first);
+             auto it = find(endo2eq.begin(), endo2eq.end(), result.begin()->first);
              if (it == endo2eq.end())
                endo2eq[result.begin()->first] = eq;
              else
@@ -547,15 +546,15 @@ ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, ve
   vector<int> eq2endo(equations.size(), 0);
   equation_reordered.resize(equations.size());
   variable_reordered.resize(equations.size());
-  bool *IM;
   int n = equations.size();
-  IM = (bool *) calloc(n*n, sizeof(bool));
+  vector<bool> IM(n*n);
   int i = 0;
-  for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, i++)
+  for (auto it : endo2eq)
     {
-      eq2endo[*it] = i;
+      eq2endo[it] = i;
       equation_reordered[i] = i;
-      variable_reordered[*it] = i;
+      variable_reordered[it] = i;
+      i++;
     }
   if (cutoff == 0)
     {
@@ -657,11 +656,10 @@ ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, ve
         }
       epilogue = tmp_epilogue;
     }
-  free(IM);
 }
 
 equation_type_and_normalized_equation_t
-ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int>>, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const
+ModelTree::equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const
 {
   expr_t lhs;
   BinaryOpNode *eq_node;
@@ -674,7 +672,7 @@ ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int>>, expr_t
       eq_node = equations[eq];
       lhs = eq_node->get_arg1();
       Equation_Simulation_Type = E_SOLVE;
-      auto derivative = first_order_endo_derivatives.find({ eq, { var, 0 } });
+      auto derivative = first_order_endo_derivatives.find({ eq, var, 0 });
       pair<bool, expr_t> res;
       if (derivative != first_order_endo_derivatives.end())
         {
@@ -683,9 +681,7 @@ ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int>>, expr_t
           auto d_endo_variable = result.find({ var, 0 });
           //Determine whether the equation could be evaluated rather than to be solved
           if (lhs->isVariableNodeEqualTo(SymbolType::endogenous, Index_Var_IM[i], 0) && derivative->second->isNumConstNodeEqualTo(1))
-            {
-              Equation_Simulation_Type = E_EVALUATE;
-            }
+            Equation_Simulation_Type = E_EVALUATE;
           else
             {
               vector<pair<int, pair<expr_t, expr_t>>> List_of_Op_RHS;
@@ -704,7 +700,7 @@ ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int>>, expr_t
         }
       V_Equation_Simulation_Type[eq] = { Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second) };
     }
-  return (V_Equation_Simulation_Type);
+  return V_Equation_Simulation_Type;
 }
 
 void
@@ -734,9 +730,8 @@ ModelTree::getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian
     }
   for (const auto & it : dynamic_jacobian)
     {
-      int lag = it.first.first;
-      int j_1 = it.first.second.first;
-      int i_1 = it.first.second.second;
+      int lag, j_1, i_1;
+      tie(lag, j_1, i_1) = it.first;
       if (variable_blck[i_1] == equation_blck[j_1])
         {
           if (lag > variable_lead_lag[i_1].second)
@@ -786,12 +781,12 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
     }
   else
     tmp_normalized_contemporaneous_jacobian = static_jacobian;
-  for (jacob_map_t::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
-    if (reverse_equation_reordered[it->first.first] >= (int) prologue && reverse_equation_reordered[it->first.first] < (int) (nb_var - epilogue)
-        && reverse_variable_reordered[it->first.second] >= (int) prologue && reverse_variable_reordered[it->first.second] < (int) (nb_var - epilogue)
-        && it->first.first != endo2eq[it->first.second])
-      add_edge(vertex(reverse_equation_reordered[endo2eq[it->first.second]]-prologue, G2),
-               vertex(reverse_equation_reordered[it->first.first]-prologue, G2),
+  for (const auto &it : tmp_normalized_contemporaneous_jacobian)
+    if (reverse_equation_reordered[it.first.first] >= (int) prologue && reverse_equation_reordered[it.first.first] < (int) (nb_var - epilogue)
+        && reverse_variable_reordered[it.first.second] >= (int) prologue && reverse_variable_reordered[it.first.second] < (int) (nb_var - epilogue)
+        && it.first.first != endo2eq[it.first.second])
+      add_edge(vertex(reverse_equation_reordered[endo2eq[it.first.second]]-prologue, G2),
+               vertex(reverse_equation_reordered[it.first.first]-prologue, G2),
                G2);
 
   vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
@@ -832,12 +827,12 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
   //   - first set = equations belonging to the block,
   //   - second set = the feeback variables,
   //   - third vector = the reordered non-feedback variables.
-  vector<pair<set<int>, pair<set<int>, vector<int>>>> components_set(num);
+  vector<tuple<set<int>, set<int>, vector<int>>> components_set(num);
   for (unsigned int i = 0; i < endo2block.size(); i++)
     {
       endo2block[i] = unordered2ordered[endo2block[i]];
       blocks[endo2block[i]].first++;
-      components_set[endo2block[i]].first.insert(i);
+      get<0>(components_set[endo2block[i]]).insert(i);
     }
 
   getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
@@ -885,12 +880,12 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
 
   for (int i = 0; i < num; i++)
     {
-      AdjacencyList_t G = extract_subgraph(G2, components_set[i].first);
+      AdjacencyList_t G = extract_subgraph(G2, get<0>(components_set[i]));
       set<int> feed_back_vertices;
       //Print(G);
       AdjacencyList_t G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
       auto v_index = get(boost::vertex_index, G);
-      components_set[i].second.first = feed_back_vertices;
+      get<1>(components_set[i]) = feed_back_vertices;
       blocks[i].second = feed_back_vertices.size();
       vector<int> Reordered_Vertice;
       Reorder_the_recursive_variables(G, feed_back_vertices, Reordered_Vertice);
@@ -898,7 +893,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
       //First we have the recursive equations conditional on feedback variables
       for (int j = 0; j < 4; j++)
         {
-          for (int & its : Reordered_Vertice)
+          for (int its : Reordered_Vertice)
             {
               bool something_done = false;
               if      (j == 2 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second != 0)
@@ -929,7 +924,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
                 }
             }
         }
-      components_set[i].second.second = Reordered_Vertice;
+      get<2>(components_set[i]) = Reordered_Vertice;
       //Second we have the equations related to the feedback variables
       for (int j = 0; j < 4; j++)
         {
@@ -978,8 +973,8 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
         n_static[prologue+num+i]++;
     }
 
-  inv_equation_reordered = vector<int>(nb_var);
-  inv_variable_reordered = vector<int>(nb_var);
+  inv_equation_reordered.resize(nb_var);
+  inv_variable_reordered.resize(nb_var);
   for (int i = 0; i < nb_var; i++)
     {
       inv_variable_reordered[variable_reordered[i]] = i;
@@ -1017,7 +1012,7 @@ ModelTree::printBlockDecomposition(const vector<pair<int, int>> &blocks) const
 }
 
 block_type_firstequation_size_mfs_t
-ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<pair< pair<int, int>, pair<int, int>>> &block_col_type, bool linear_decomposition)
+ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<tuple<int, int, int, int>> &block_col_type, bool linear_decomposition)
 {
   int i = 0;
   int count_equ = 0, blck_count_simult = 0;
@@ -1063,7 +1058,7 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
               auto it1 = find(variable_reordered.begin()+first_count_equ, variable_reordered.begin()+(first_count_equ+Blck_Size), curr_variable);
               if (linear_decomposition)
                 {
-                  if (dynamic_jacobian.find(make_pair(curr_lag, make_pair(equation_reordered[count_equ], curr_variable))) != dynamic_jacobian.end())
+                  if (dynamic_jacobian.find({ curr_lag, equation_reordered[count_equ], curr_variable }) != dynamic_jacobian.end())
                     {
                       if (curr_lag > Lead)
                         Lead = curr_lag;
@@ -1074,7 +1069,7 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
               else
                 {
                   if (it1 != variable_reordered.begin()+(first_count_equ+Blck_Size))
-                    if (dynamic_jacobian.find({ curr_lag, { equation_reordered[count_equ], curr_variable } }) != dynamic_jacobian.end())
+                    if (dynamic_jacobian.find({ curr_lag, equation_reordered[count_equ], curr_variable }) != dynamic_jacobian.end())
                       {
                         if (curr_lag > Lead)
                           Lead = curr_lag;
@@ -1121,17 +1116,17 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
           if (i > 0)
             {
               bool is_lead = false, is_lag = false;
-              int c_Size = (block_type_size_mfs[block_type_size_mfs.size()-1]).second.first;
-              int first_equation = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.second;
+              int c_Size = get<2>(block_type_size_mfs[block_type_size_mfs.size()-1]);
+              int first_equation = get<1>(block_type_size_mfs[block_type_size_mfs.size()-1]);
               if (c_Size > 0 && ((prev_Type ==  EVALUATE_FORWARD && Simulation_Type == EVALUATE_FORWARD && !is_lead)
                                  || (prev_Type ==  EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD && !is_lag)))
                 {
                   for (int j = first_equation; j < first_equation+c_Size; j++)
                     {
-                      auto it = dynamic_jacobian.find({ -1, { equation_reordered[eq], variable_reordered[j] } });
+                      auto it = dynamic_jacobian.find({ -1, equation_reordered[eq], variable_reordered[j] });
                       if (it != dynamic_jacobian.end())
                         is_lag = true;
-                      it = dynamic_jacobian.find({ +1, { equation_reordered[eq], variable_reordered[j] } });
+                      it = dynamic_jacobian.find({ +1, equation_reordered[eq], variable_reordered[j] });
                       if (it != dynamic_jacobian.end())
                         is_lead = true;
                     }
@@ -1140,55 +1135,55 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
                   || (prev_Type ==  EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD && !is_lag))
                 {
                   //merge the current block with the previous one
-                  BlockSimulationType c_Type = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.first;
+                  BlockSimulationType c_Type = get<0>(block_type_size_mfs[block_type_size_mfs.size()-1]);
                   c_Size++;
-                  block_type_size_mfs[block_type_size_mfs.size()-1] = { { c_Type, first_equation }, { c_Size, c_Size } };
+                  block_type_size_mfs[block_type_size_mfs.size()-1] = { c_Type, first_equation, c_Size, c_Size };
                   if (block_lag_lead[block_type_size_mfs.size()-1].first > Lag)
                     Lag = block_lag_lead[block_type_size_mfs.size()-1].first;
                   if (block_lag_lead[block_type_size_mfs.size()-1].second > Lead)
                     Lead = block_lag_lead[block_type_size_mfs.size()-1].second;
                   block_lag_lead[block_type_size_mfs.size()-1] = { Lag, Lead };
-                  pair< pair< unsigned int, unsigned int>, pair<unsigned int, unsigned int>> tmp = block_col_type[block_col_type.size()-1];
-                  block_col_type[block_col_type.size()-1] = { { tmp.first.first+l_n_static, tmp.first.second+l_n_forward }, { tmp.second.first+l_n_backward, tmp.second.second+l_n_mixed } };
+                  auto tmp = block_col_type[block_col_type.size()-1];
+                  block_col_type[block_col_type.size()-1] = { get<0>(tmp)+l_n_static, get<1>(tmp)+l_n_forward, get<2>(tmp)+l_n_backward, get<3>(tmp)+l_n_mixed };
                 }
               else
                 {
-                  block_type_size_mfs.emplace_back(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size));
+                  block_type_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
                   block_lag_lead.emplace_back(Lag, Lead);
-                  block_col_type.emplace_back(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed));
+                  block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
                 }
             }
           else
             {
-              block_type_size_mfs.emplace_back(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size));
+              block_type_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
               block_lag_lead.emplace_back(Lag, Lead);
-              block_col_type.emplace_back(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed));
+              block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
             }
         }
       else
         {
-          block_type_size_mfs.emplace_back(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size));
+          block_type_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
           block_lag_lead.emplace_back(Lag, Lead);
-          block_col_type.emplace_back(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed));
+          block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
         }
       prev_Type = Simulation_Type;
       eq += Blck_Size;
     }
-  return (block_type_size_mfs);
+  return block_type_size_mfs;
 }
 
 vector<bool>
-ModelTree::equationLinear(map<pair<int, pair<int, int> >, expr_t> first_order_endo_derivatives) const
+ModelTree::equationLinear(map<tuple<int, int, int>, expr_t> first_order_endo_derivatives) const
 {
   vector<bool> is_linear(symbol_table.endo_nbr(), true);
-  for (map<pair<int, pair<int, int> >, expr_t>::const_iterator it = first_order_endo_derivatives.begin(); it != first_order_endo_derivatives.end(); it++)
+  for (const auto &it : first_order_endo_derivatives)
     {
-       expr_t Id = it->second;
-       set<pair<int, int> > endogenous;
+       expr_t Id = it.second;
+       set<pair<int, int>> endogenous;
        Id->collectEndogenous(endogenous);
        if (endogenous.size() > 0)
          {
-           int eq = it->first.first;
+           int eq = get<0>(it.first);
            is_linear[eq] = false;
          }
     }
@@ -1208,12 +1203,12 @@ ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vec
       int first_variable_position = getBlockFirstEquation(block);
       if (simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
         {
-          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = derivatives_block.begin(); it != derivatives_block.end(); it++)
+          for (const auto &it : derivatives_block)
             {
-              int lag = it->second.first;
+              int lag = get<2>(it);
               if (lag == 0)
                 {
-                  expr_t Id = it->second.second;
+                  expr_t Id = get<3>(it);
                   set<pair<int, int>> endogenous;
                   Id->collectEndogenous(endogenous);
                   if (endogenous.size() > 0)
@@ -1232,10 +1227,10 @@ ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vec
         }
       else if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
         {
-          for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = derivatives_block.begin(); it != derivatives_block.end(); it++)
+          for (const auto &it : derivatives_block)
             {
-              int lag = it->second.first;
-              expr_t Id = it->second.second; //
+              int lag = get<2>(it);
+              expr_t Id = get<3>(it); //
               set<pair<int, int>> endogenous;
               Id->collectEndogenous(endogenous);
               if (endogenous.size() > 0)
@@ -1254,7 +1249,7 @@ ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vec
     the_end:
       ;
     }
-  return (blocks_linear);
+  return blocks_linear;
 }
 
 int
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index aec7fe645768450d852bf16037733bf425b81e37..882a06d67416c1d19c19a0b17e020cf742a9aac5 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -46,16 +46,16 @@ auto vectorToTuple(const vector<T>& v) {
 }
 
 //! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a expr_t on the new normalized equation
-using equation_type_and_normalized_equation_t = vector<pair<EquationType, expr_t >>;
+using equation_type_and_normalized_equation_t = vector<pair<EquationType, expr_t>>;
 
 //! Vector describing variables: max_lag in the block, max_lead in the block
 using lag_lead_vector_t = vector<pair< int, int>>;
 
-//! for each block contains pair< pair<Simulation_Type, first_equation>, pair < Block_Size, Recursive_part_Size >>
-using block_type_firstequation_size_mfs_t = vector<pair< pair< BlockSimulationType, int>, pair<int, int>>>;
+//! for each block contains tuple<Simulation_Type, first_equation, Block_Size, Recursive_part_Size>
+using block_type_firstequation_size_mfs_t = vector<tuple<BlockSimulationType, int, int, int>>;
 
-//! for a block contains derivatives pair< pair<block_equation_number, block_variable_number> , pair<lead_lag, expr_t>>
-using block_derivatives_equation_variable_laglead_nodeid_t = vector< pair<pair<int, int>, pair< int, expr_t >>>;
+//! for a block contains derivatives tuple<block_equation_number, block_variable_number, lead_lag, expr_t>
+using block_derivatives_equation_variable_laglead_nodeid_t = vector<tuple<int, int, int, expr_t>>;
 
 //! for all blocks derivatives description
 using blocks_derivatives_t = vector<block_derivatives_equation_variable_laglead_nodeid_t>;
@@ -194,7 +194,7 @@ protected:
 
   //! Sparse matrix of double to store the values of the Jacobian
   /*! First index is lag, second index is equation number, third index is endogenous type specific ID */
-  using dynamic_jacob_map_t = map<pair<int, pair<int, int>>, expr_t>;
+  using dynamic_jacob_map_t = map<tuple<int, int, int>, expr_t>;
 
   //! Normalization of equations
   /*! Maps endogenous type specific IDs to equation numbers */
@@ -234,15 +234,15 @@ protected:
   //! Search the equations and variables belonging to the prologue and the epilogue of the model
   void computePrologueAndEpilogue(const jacob_map_t &static_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered);
   //! Determine the type of each equation of model and try to normalized the unnormalized equation using computeNormalizedEquations
-  equation_type_and_normalized_equation_t equationTypeDetermination(const map<pair<int, pair<int, int>>, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const;
+  equation_type_and_normalized_equation_t equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const;
   //! Compute the block decomposition and for a non-recusive block find the minimum feedback set
   void computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered, lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead_t, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed) const;
   //! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block
-  block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<pair< pair<int, int>, pair<int, int>>> &block_col_type, bool linear_decomposition);
+  block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<tuple<int, int, int, int>> &block_col_type, bool linear_decomposition);
   //! Determine the maximum number of lead and lag for the endogenous variable in a bloc
   void getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian, const vector<int> &components_set, int nb_blck_sim, lag_lead_vector_t &equation_lead_lag, lag_lead_vector_t &variable_lead_lag, const vector<int> &equation_reordered, const vector<int> &variable_reordered) const;
   //! For each equation determine if it is linear or not
-  vector<bool> equationLinear(map<pair<int, pair<int, int> >, expr_t> first_order_endo_derivatives) const;
+  vector<bool> equationLinear(map<tuple<int, int, int>, expr_t> first_order_endo_derivatives) const;
   //! Print an abstract of the block structure of the model
   void printBlockDecomposition(const vector<pair<int, int>> &blocks) const;
   //! Determine for each block if it is linear or not
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 584aad0547a6bbe7f111654601fb3d69e3bf8224..e9afe202795ae42b07a4f1e6797be834149fcfab 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -61,7 +61,7 @@ StaticModel::copyHelper(const StaticModel &m)
     {
       block_derivatives_equation_variable_laglead_nodeid_t v;
       for (const auto &it2 : it)
-        v.push_back(make_pair(it2.first, make_pair(it2.second.first, f(it2.second.second))));
+        v.emplace_back(get<0>(it2), get<1>(it2), get<2>(it2), f(get<3>(it2)));
       blocks_derivatives.push_back(v);
     }
 
@@ -100,12 +100,7 @@ StaticModel::StaticModel(const StaticModel &m) :
   global_temporary_terms {m.global_temporary_terms},
   block_type_firstequation_size_mfs {m.block_type_firstequation_size_mfs},
   blocks_linear {m.blocks_linear},
-  other_endo_block {m.other_endo_block},
-  exo_block {m.exo_block},
-  exo_det_block {m.exo_det_block},
   block_col_type {m.block_col_type},
-  variable_block_lead_lag {m.variable_block_lead_lag},
-  equation_block {m.equation_block},
   endo_max_leadlag_block {m.endo_max_leadlag_block},
   other_endo_max_leadlag_block {m.other_endo_max_leadlag_block},
   exo_max_leadlag_block {m.exo_max_leadlag_block},
@@ -145,12 +140,7 @@ StaticModel::operator=(const StaticModel &m)
   derivative_exo.clear();
   derivative_exo_det.clear();
 
-  other_endo_block = m.other_endo_block;
-  exo_block = m.exo_block;
-  exo_det_block = m.exo_det_block;
   block_col_type = m.block_col_type;
-  variable_block_lead_lag = m.variable_block_lead_lag;
-  equation_block = m.equation_block;
   endo_max_leadlag_block = m.endo_max_leadlag_block;
   other_endo_max_leadlag_block = m.other_endo_max_leadlag_block;
   exo_max_leadlag_block = m.exo_max_leadlag_block;
@@ -231,7 +221,7 @@ StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_nu
 void
 StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
 {
-  auto it = first_chain_rule_derivatives.find({ eqr, { varr, lag } });
+  auto it = first_chain_rule_derivatives.find({ eqr, varr, lag });
   if (it != first_chain_rule_derivatives.end())
     (it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
   else
@@ -249,7 +239,6 @@ StaticModel::computeTemporaryTermsOrdered()
   BinaryOpNode *eq_node;
   first_chain_rule_derivatives_t::const_iterator it_chr;
   ostringstream tmp_s;
-  v_temporary_terms.clear();
   map_idx.clear();
 
   unsigned int nb_blocks = getNbBlocks();
@@ -266,11 +255,8 @@ StaticModel::computeTemporaryTermsOrdered()
   for (unsigned int block = 0; block < nb_blocks; block++)
     {
       map<expr_t, int> reference_count_local;
-      reference_count_local.clear();
       map<expr_t, pair<int, int>> first_occurence_local;
-      first_occurence_local.clear();
       temporary_terms_t temporary_terms_l;
-      temporary_terms_l.clear();
 
       unsigned int block_size = getBlockSize(block);
       unsigned int block_nb_mfs = getBlockMfs(block);
@@ -287,14 +273,12 @@ StaticModel::computeTemporaryTermsOrdered()
               eq_node->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local,  i);
             }
         }
-      for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+      for (const auto &it : blocks_derivatives[block])
         {
-          expr_t id = it->second.second;
+          expr_t id = get<3>(it);
           id->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local,  block_size-1);
         }
-      set<int> temporary_terms_in_use;
-      temporary_terms_in_use.clear();
-      v_temporary_terms_inuse[block] = temporary_terms_in_use;
+      v_temporary_terms_inuse[block] = {};
       computeTemporaryTermsMapping(temporary_terms_l, map_idx2[block]);
     }
 
@@ -316,9 +300,9 @@ StaticModel::computeTemporaryTermsOrdered()
               eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
             }
         }
-      for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+      for (const auto &it : blocks_derivatives[block])
         {
-          expr_t id = it->second.second;
+          expr_t id = get<3>(it);
           id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
         }
     }
@@ -340,15 +324,14 @@ StaticModel::computeTemporaryTermsOrdered()
               eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
             }
         }
-      for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+      for (const auto &it : blocks_derivatives[block])
         {
-          expr_t id = it->second.second;
+          expr_t id = get<3>(it);
           id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
         }
       for (int i = 0; i < (int) getBlockSize(block); i++)
-        for (auto it = v_temporary_terms[block][i].begin();
-             it != v_temporary_terms[block][i].end(); it++)
-          (*it)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+        for (const auto &it : v_temporary_terms[block][i])
+          it->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
       v_temporary_terms_inuse[block] = temporary_terms_in_use;
     }
   computeTemporaryTermsMapping(temporary_terms, map_idx);
@@ -396,16 +379,16 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
       tmp1_output.str("");
       tmp1_output << packageDir(basename + ".block") << "/static_" << block+1 << ".m";
       output.open(tmp1_output.str(), ios::out | ios::binary);
-      output << "%\n";
-      output << "% " << tmp1_output.str() << " : Computes static model for Dynare\n";
-      output << "%\n";
-      output << "% Warning : this file is generated automatically by Dynare\n";
-      output << "%           from model file (.mod)\n\n";
-      output << "%/\n";
+      output << "%" << endl
+             << "% " << tmp1_output.str() << " : Computes static model for Dynare" << endl
+             << "%" << endl
+             << "% Warning : this file is generated automatically by Dynare" << endl
+             << "%           from model file (.mod)" << endl << endl
+             << "%/" << endl;
       if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
-        output << "function y = static_" << block+1 << "(y, x, params)\n";
+        output << "function y = static_" << block+1 << "(y, x, params)" << endl;
       else
-        output << "function [residual, y, g1] = static_" << block+1 << "(y, x, params)\n";
+        output << "function [residual, y, g1] = static_" << block+1 << "(y, x, params)" << endl;
 
       BlockType block_type;
       if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
@@ -436,11 +419,11 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
           tmp_output.str("");
           for (int it : v_temporary_terms_inuse[block])
             tmp_output << " T" << it;
-          output << "  global" << tmp_output.str() << ";\n";
+          output << "  global" << tmp_output.str() << ";" << endl;
         }
 
       if (simulation_type != EVALUATE_BACKWARD && simulation_type != EVALUATE_FORWARD)
-        output << "  residual=zeros(" << block_mfs << ",1);\n";
+        output << "  residual=zeros(" << block_mfs << ",1);" << endl;
 
       // The equations
       temporary_terms_idxs_t temporary_terms_idxs;
@@ -449,7 +432,6 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
           if (!global_temporary_terms)
             local_temporary_terms = v_temporary_terms[block][i];
           temporary_terms_t tt2;
-          tt2.clear();
           if (v_temporary_terms[block].size())
             {
               output << "  " << "% //Temporary variables" << endl;
@@ -498,7 +480,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
                   if (isBlockEquationRenormalized(block, i))
                     {
                       rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
-                      output << "\n  ";
+                      output << endl << "  ";
                       tmp_output.str("");
                       eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
                       lhs = eq_node->get_arg1();
@@ -510,10 +492,10 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
                 }
               else
                 {
-                  cerr << "Type mismatch for equation " << equation_ID+1  << "\n";
+                  cerr << "Type mismatch for equation " << equation_ID+1  << endl;
                   exit(EXIT_FAILURE);
                 }
-              output << ";\n";
+              output << ";" << endl;
               break;
             case SOLVE_BACKWARD_SIMPLE:
             case SOLVE_FORWARD_SIMPLE:
@@ -531,7 +513,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
               output << tmp_output.str();
               output << ") - (";
               rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
-              output << ");\n";
+              output << ");" << endl;
             }
         }
       // The Jacobian if we have to solve the block
@@ -544,13 +526,13 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
         case SOLVE_FORWARD_SIMPLE:
         case SOLVE_BACKWARD_COMPLETE:
         case SOLVE_FORWARD_COMPLETE:
-          for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+          for (const auto &it : blocks_derivatives[block])
             {
-              unsigned int eq = it->first.first;
-              unsigned int var = it->first.second;
+              unsigned int eq, var;
+              expr_t id;
+              tie(eq, var, ignore, id) = it;
               unsigned int eqr = getBlockEquationID(block, eq);
               unsigned int varr = getBlockVariableID(block, var);
-              expr_t id = it->second.second;
               output << "    g1(" << eq+1-block_recursive << ", " << var+1-block_recursive << ") = ";
               id->writeOutput(output, local_output_type, local_temporary_terms, {});
               output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
@@ -626,8 +608,7 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
   fjmp_if_eval.write(code_file, instruction_number);
   int prev_instruction_number = instruction_number;
 
-  vector<vector<pair<int, int>>> my_derivatives;
-  my_derivatives.resize(symbol_table.endo_nbr());
+  vector<vector<pair<int, int>>> my_derivatives(symbol_table.endo_nbr());
   count_u = symbol_table.endo_nbr();
   for (const auto & first_derivative : derivatives[1])
     {
@@ -657,8 +638,7 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
       fldr.write(code_file, instruction_number);
       if (my_derivatives[i].size())
         {
-          for (vector<pair<int, int>>::const_iterator it = my_derivatives[i].begin();
-               it != my_derivatives[i].end(); it++)
+          for (auto it = my_derivatives[i].begin(); it != my_derivatives[i].end(); it++)
             {
               FLDSU_ fldsu(it->second);
               fldsu.write(code_file, instruction_number);
@@ -690,10 +670,7 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
   code_file.seekp(pos3);
   prev_instruction_number = instruction_number;
 
-  temporary_terms_t tt2;
-  tt2.clear();
-  temporary_terms_t tt3;
-  tt3.clear();
+  temporary_terms_t tt2, tt3;
 
   // The Jacobian if we have to solve the block determinsitic bloc
   for (const auto & first_derivative : derivatives[1])
@@ -820,7 +797,6 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
         {
           //The Temporary terms
           temporary_terms_t tt2;
-          tt2.clear();
           if (v_temporary_terms[block].size())
             {
               for (auto it : v_temporary_terms[block][i])
@@ -919,10 +895,10 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
             case SOLVE_BACKWARD_COMPLETE:
             case SOLVE_FORWARD_COMPLETE:
               count_u = feedback_variables.size();
-              for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+              for (const auto &it : blocks_derivatives[block])
                 {
-                  unsigned int eq = it->first.first;
-                  unsigned int var = it->first.second;
+                  unsigned int eq, var;
+                  tie(eq, var, ignore, ignore) = it;
                   unsigned int eqr = getBlockEquationID(block, eq);
                   unsigned int varr = getBlockVariableID(block, var);
                   if (eq >= block_recursive && var >= block_recursive)
@@ -1005,32 +981,28 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
       code_file.seekp(pos3);
       prev_instruction_number = instruction_number;
 
-      temporary_terms_t tt2;
-      tt2.clear();
-      temporary_terms_t tt3;
-      tt3.clear();
+      temporary_terms_t tt2, tt3;
       deriv_node_temp_terms_t tef_terms2;
 
       for (i = 0; i < (int) block_size; i++)
         {
           if (v_temporary_terms_local[block].size())
             {
-              for (auto it = v_temporary_terms_local[block][i].begin();
-                   it != v_temporary_terms_local[block][i].end(); it++)
+              for (const auto &it : v_temporary_terms_local[block][i])
                 {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != nullptr)
-                    (*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms2);
+                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
+                    it->compileExternalFunctionOutput(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms2);
 
-                  FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx2[block].find((*it)->idx)->second));
+                  FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx2[block].find(it->idx)->second));
                   fnumexpr.write(code_file, instruction_number);
 
-                  (*it)->compile(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms);
+                  it->compile(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms);
 
-                  FSTPST_ fstpst((int)(map_idx2[block].find((*it)->idx)->second));
+                  FSTPST_ fstpst((int)(map_idx2[block].find(it->idx)->second));
                   fstpst.write(code_file, instruction_number);
                   // Insert current node into tt2
-                  tt3.insert(*it);
-                  tt2.insert(*it);
+                  tt3.insert(it);
+                  tt2.insert(it);
                 }
             }
 
@@ -1113,10 +1085,10 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
         case SOLVE_BACKWARD_COMPLETE:
         case SOLVE_FORWARD_COMPLETE:
           count_u = feedback_variables.size();
-          for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+          for (const auto &it : blocks_derivatives[block])
             {
-              unsigned int eq = it->first.first;
-              unsigned int var = it->first.second;
+              unsigned int eq, var;
+              tie(eq, var, ignore, ignore) = it;
               unsigned int eqr = getBlockEquationID(block, eq);
               unsigned int varr = getBlockVariableID(block, var);
               FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, 0);
@@ -1165,10 +1137,10 @@ StaticModel::Write_Inf_To_Bin_File_Block(const string &basename, const int &num,
   unsigned int block_size = getBlockSize(num);
   unsigned int block_mfs = getBlockMfs(num);
   unsigned int block_recursive = block_size - block_mfs;
-  for (auto it = blocks_derivatives[num].begin(); it != (blocks_derivatives[num]).end(); it++)
+  for (const auto &it : blocks_derivatives[num])
     {
-      unsigned int eq = it->first.first;
-      unsigned int var = it->first.second;
+      unsigned int eq, var;
+      tie(eq, var, ignore, ignore) = it;
       int lag = 0;
       if (eq >= block_recursive && var >= block_recursive)
         {
@@ -1196,10 +1168,10 @@ StaticModel::Write_Inf_To_Bin_File_Block(const string &basename, const int &num,
   SaveCode.close();
 }
 
-map<pair<int, pair<int, int >>, expr_t>
+map<tuple<int, int, int>, expr_t>
 StaticModel::collect_first_order_derivatives_endogenous()
 {
-  map<pair<int, pair<int, int >>, expr_t> endo_derivatives;
+  map<tuple<int, int, int>, expr_t> endo_derivatives;
   for (auto & first_derivative : derivatives[1])
     {
       if (getTypeByDerivID(first_derivative.first[1]) == SymbolType::endogenous)
@@ -1207,7 +1179,7 @@ StaticModel::collect_first_order_derivatives_endogenous()
           int eq = first_derivative.first[0];
           int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first[1]));
           int lag = 0;
-          endo_derivatives[{ eq, { var, lag } }] = first_derivative.second;
+          endo_derivatives[{ eq, var, lag }] = first_derivative.second;
         }
     }
   return endo_derivatives;
@@ -1270,12 +1242,12 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
 
       computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered);
 
-      map<pair<int, pair<int, int>>, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
+      map<tuple<int, int, int>, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
 
       equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
 
       if (!nopreprocessoroutput)
-        cout << "Finding the optimal block decomposition of the model ...\n";
+        cout << "Finding the optimal block decomposition of the model ..." << endl;
 
       lag_lead_vector_t equation_lag_lead, variable_lag_lead;
 
@@ -2209,7 +2181,7 @@ StaticModel::writeStaticBlockMFSFile(const string &basename) const
   output << "function [residual, g1, y, var_index] = static(nblock, y, x, params)" << endl
          << "  residual = [];" << endl
          << "  g1 = [];" << endl
-         << "  var_index = [];\n" << endl
+         << "  var_index = [];" << endl << endl
          << "  switch nblock" << endl;
 
   unsigned int nb_blocks = getNbBlocks();
@@ -2225,16 +2197,16 @@ StaticModel::writeStaticBlockMFSFile(const string &basename) const
 
       if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
         {
-          output << "      y_tmp = " << basename << ".block.static_" << b+1 << "(y, x, params);\n";
+          output << "      y_tmp = " << basename << ".block.static_" << b+1 << "(y, x, params);" << endl;
           ostringstream tmp;
           for (int i = 0; i < (int) getBlockSize(b); i++)
             tmp << " " << getBlockVariableID(b, i)+1;
-          output << "      var_index = [" << tmp.str() << "];\n";
-          output << "      residual  = y(var_index) - y_tmp(var_index);\n";
-          output << "      y = y_tmp;\n";
+          output << "      var_index = [" << tmp.str() << "];" << endl
+                 << "      residual  = y(var_index) - y_tmp(var_index);" << endl
+                 << "      y = y_tmp;" << endl;
         }
       else
-        output << "      [residual, y, g1] = " << basename << ".block.static_" << b+1 << "(y, x, params);\n";
+        output << "      [residual, y, g1] = " << basename << ".block.static_" << b+1 << "(y, x, params);" << endl;
 
     }
   output << "  end" << endl
@@ -2267,23 +2239,23 @@ StaticModel::writeOutput(ostream &output, bool block) const
           tmp_s << " " << getBlockVariableID(b, i)+1;
           tmp_s_eq << " " << getBlockEquationID(b, i)+1;
         }
-      output << "block_structure_stat.block(" << b+1 << ").Simulation_Type = " << simulation_type << ";\n";
-      output << "block_structure_stat.block(" << b+1 << ").endo_nbr = " << block_size << ";\n";
-      output << "block_structure_stat.block(" << b+1 << ").mfs = " << getBlockMfs(block) << ";\n";
-      output << "block_structure_stat.block(" << b+1 << ").equation = [" << tmp_s_eq.str() << "];\n";
-      output << "block_structure_stat.block(" << b+1 << ").variable = [" << tmp_s.str() << "];\n";
+      output << "block_structure_stat.block(" << b+1 << ").Simulation_Type = " << simulation_type << ";" << endl
+             << "block_structure_stat.block(" << b+1 << ").endo_nbr = " << block_size << ";" << endl
+             << "block_structure_stat.block(" << b+1 << ").mfs = " << getBlockMfs(block) << ";" << endl
+             << "block_structure_stat.block(" << b+1 << ").equation = [" << tmp_s_eq.str() << "];" << endl
+             << "block_structure_stat.block(" << b+1 << ").variable = [" << tmp_s.str() << "];" << endl;
     }
-  output << "M_.block_structure_stat.block = block_structure_stat.block;\n";
+  output << "M_.block_structure_stat.block = block_structure_stat.block;" << endl;
   string cst_s;
   int nb_endo = symbol_table.endo_nbr();
   output << "M_.block_structure_stat.variable_reordered = [";
   for (int i = 0; i < nb_endo; i++)
     output << " " << variable_reordered[i]+1;
-  output << "];\n";
-  output << "M_.block_structure_stat.equation_reordered = [";
+  output << "];" << endl
+         << "M_.block_structure_stat.equation_reordered = [";
   for (int i = 0; i < nb_endo; i++)
     output << " " << equation_reordered[i]+1;
-  output << "];\n";
+  output << "];" << endl;
 
   map<pair<int, int>,  int>  row_incidence;
   for (const auto & first_derivative : derivatives[1])
@@ -2299,11 +2271,9 @@ StaticModel::writeOutput(ostream &output, bool block) const
         }
     }
   output << "M_.block_structure_stat.incidence.sparse_IM = [";
-  for (map<pair< int, int >,  int>::const_iterator it = row_incidence.begin(); it != row_incidence.end(); it++)
-    {
-      output << it->first.first+1 << " " << it->first.second+1 << ";\n";
-    }
-  output << "];\n";
+  for (const auto &it : row_incidence)
+    output << it.first.first+1 << " " << it.first.second+1 << ";" << endl;
+  output << "];" << endl;
 }
 
 SymbolType
@@ -2352,11 +2322,10 @@ StaticModel::addAllParamDerivId(set<int> &deriv_id_set)
     deriv_id_set.insert(i + symbol_table.endo_nbr());
 }
 
-map<pair<pair<int, pair<int, int>>, pair<int, int>>, int>
+map<tuple<int, int, int, int, int>, int>
 StaticModel::get_Derivatives(int block)
 {
-  map<pair<pair<int, pair<int, int>>, pair<int, int>>, int> Derivatives;
-  Derivatives.clear();
+  map<tuple<int, int, int, int, int>, int> Derivatives;
   int block_size = getBlockSize(block);
   int block_nb_recursive = block_size - getBlockMfs(block);
   int lag = 0;
@@ -2366,10 +2335,10 @@ StaticModel::get_Derivatives(int block)
       for (int var = 0; var < block_size; var++)
         {
           int varr = getBlockVariableID(block, var);
-          if (dynamic_jacobian.find({ lag, { eqr, varr } }) != dynamic_jacobian.end())
+          if (dynamic_jacobian.find({ lag, eqr, varr }) != dynamic_jacobian.end())
             {
               bool OK = true;
-              map<pair<pair<int, pair<int, int>>, pair<int, int>>, int>::const_iterator its = Derivatives.find({ { lag, { eq, var } }, { eqr, varr } });
+              auto its = Derivatives.find({ lag, eq, var, eqr, varr });
               if (its != Derivatives.end())
                 {
                   if (its->second == 2)
@@ -2380,10 +2349,10 @@ StaticModel::get_Derivatives(int block)
                 {
                   if (getBlockEquationType(block, eq) == E_EVALUATE_S && eq < block_nb_recursive)
                     //It's a normalized equation, we have to recompute the derivative using chain rule derivative function
-                    Derivatives[{ { lag, { eq, var } }, { eqr, varr } }] = 1;
+                    Derivatives[{ lag, eq, var, eqr, varr }] = 1;
                   else
                     //It's a feedback equation we can use the derivatives
-                    Derivatives[{ { lag, { eq, var } }, { eqr, varr } }] = 0;
+                    Derivatives[{ lag, eq, var, eqr, varr }] = 0;
                 }
               if (var < block_nb_recursive)
                 {
@@ -2392,15 +2361,15 @@ StaticModel::get_Derivatives(int block)
                     {
                       int varrs = getBlockVariableID(block, vars);
                       //A new derivative needs to be computed using the chain rule derivative function (a feedback variable appears in a recursive equation)
-                      if (Derivatives.find({ { lag, { var, vars } }, { eqs, varrs } }) != Derivatives.end())
-                        Derivatives[{ { lag, { eq, vars } }, { eqr, varrs } }] = 2;
+                      if (Derivatives.find({ lag, var, vars, eqs, varrs }) != Derivatives.end())
+                        Derivatives[{ lag, eq, vars, eqr, varrs }] = 2;
                     }
                 }
             }
         }
     }
 
-  return (Derivatives);
+  return Derivatives;
 }
 
 void
@@ -2427,30 +2396,24 @@ StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
               else
                 recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
             }
-          map<pair<pair<int, pair<int, int>>, pair<int, int>>, int> Derivatives = get_Derivatives(block);
-          map<pair<pair<int, pair<int, int>>, pair<int, int>>, int>::const_iterator it = Derivatives.begin();
-          for (int i = 0; i < (int) Derivatives.size(); i++)
+          auto Derivatives = get_Derivatives(block);
+          for (const auto &it : Derivatives)
             {
-              int Deriv_type = it->second;
-              pair<pair<int, pair<int, int>>, pair<int, int>> it_l(it->first);
-              it++;
-              int lag = it_l.first.first;
-              int eq = it_l.first.second.first;
-              int var = it_l.first.second.second;
-              int eqr = it_l.second.first;
-              int varr = it_l.second.second;
+              int Deriv_type = it.second;
+              int lag, eq, var, eqr, varr;
+              tie(lag, eq, var, eqr, varr) = it.first;
               if (Deriv_type == 0)
-                first_chain_rule_derivatives[{ eqr, { varr, lag } }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
+                first_chain_rule_derivatives[{ eqr, varr, lag }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
               else if (Deriv_type == 1)
-                first_chain_rule_derivatives[{ eqr, { varr, lag } }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
+                first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
               else if (Deriv_type == 2)
                 {
                   if (getBlockEquationType(block, eq) == E_EVALUATE_S && eq < block_nb_recursives)
-                    first_chain_rule_derivatives[{ eqr, { varr, lag } }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
+                    first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
                   else
-                    first_chain_rule_derivatives[{ eqr, { varr, lag } }] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
+                    first_chain_rule_derivatives[{ eqr, varr, lag }] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
                 }
-              tmp_derivatives.emplace_back(make_pair(eq, var), make_pair(lag, first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))]));
+              tmp_derivatives.emplace_back(eq, var, lag, first_chain_rule_derivatives[{ eqr, varr, lag }]);
             }
         }
       else
@@ -2472,8 +2435,8 @@ StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
                   expr_t d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), 0), recursive_variables);
                   if (d1 == Zero)
                     continue;
-                  first_chain_rule_derivatives[{ eqr, { varr, 0 } }] = d1;
-                  tmp_derivatives.emplace_back(make_pair(eq, var), make_pair(0, first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, 0))]));
+                  first_chain_rule_derivatives[{ eqr, varr, 0 }] = d1;
+                  tmp_derivatives.emplace_back(eq, var, 0, first_chain_rule_derivatives[{ eqr, varr, 0 }]);
                 }
             }
         }
@@ -2485,10 +2448,8 @@ void
 StaticModel::collect_block_first_order_derivatives()
 {
   //! vector for an equation or a variable indicates the block number
-  vector<int> equation_2_block, variable_2_block;
+  vector<int> equation_2_block(equation_reordered.size()), variable_2_block(variable_reordered.size());
   unsigned int nb_blocks = getNbBlocks();
-  equation_2_block = vector<int>(equation_reordered.size());
-  variable_2_block = vector<int>(variable_reordered.size());
   for (unsigned int block = 0; block < nb_blocks; block++)
     {
       unsigned int block_size = getBlockSize(block);
@@ -2517,7 +2478,7 @@ StaticModel::collect_block_first_order_derivatives()
       if (getTypeByDerivID(first_derivative.first[1]) == SymbolType::endogenous && block_eq == block_var)
         {
           tmp_derivative = derivative_endo[block_eq];
-          tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
+          tmp_derivative[{ lag, eq, var }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
           derivative_endo[block_eq] = tmp_derivative;
         }
     }
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index de40a312c6de5b062f32405d85291ea31225c78e..59bfb319145f621930c35bc805fb0a2a09881cc7 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -42,7 +42,7 @@ private:
 
   vector<temporary_terms_inuse_t> v_temporary_terms_inuse;
 
-  using first_chain_rule_derivatives_t = map< pair< int, pair< int, int>>, expr_t>;
+  using first_chain_rule_derivatives_t = map<tuple<int, int, int>, expr_t>;
   first_chain_rule_derivatives_t first_chain_rule_derivatives;
 
   //! Writes static model file (standard Matlab version)
@@ -99,11 +99,11 @@ private:
   //! Compute the column indices of the static Jacobian
   void computeStatJacobianCols();
   //! return a map on the block jacobian
-  map<pair<pair<int, pair<int, int>>, pair<int, int>>, int> get_Derivatives(int block);
+  map<tuple<int, int, int, int, int>, int> get_Derivatives(int block);
   //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
   void computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives);
   //! Collect only the first derivatives
-  map<pair<int, pair<int, int>>, expr_t> collect_first_order_derivatives_endogenous();
+  map<tuple<int, int, int>, expr_t> collect_first_order_derivatives_endogenous();
 
   //! Collecte the derivatives w.r. to endogenous of the block, to endogenous of previouys blocks and to exogenous
   void collect_block_first_order_derivatives();
@@ -126,24 +126,18 @@ private:
   //! Vector indicating if the block is linear in endogenous variable (true) or not (false)
   vector<bool> blocks_linear;
 
-  //! Map the derivatives for a block pair<lag, make_pair(make_pair(eq, var)), expr_t>
-  using derivative_t = map<pair< int, pair<int, int>>, expr_t>;
+  //! Map the derivatives for a block tuple<lag, eq, var>
+  using derivative_t = map<tuple<int, int, int>, expr_t>;
   //! Vector of derivative for each blocks
   vector<derivative_t> derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
 
   //!List for each block and for each lag-leag all the other endogenous variables and exogenous variables
   using var_t = set<int>;
   using lag_var_t = map<int, var_t>;
-  vector<lag_var_t> other_endo_block, exo_block, exo_det_block;
 
   //! for each block described the number of static, forward, backward and mixed variables in the block
-  /*! pair< pair<static, forward>, pair<backward,mixed>> */
-  vector<pair< pair<int, int>, pair<int, int>>> block_col_type;
-
-  //! List for each variable its block number and its maximum lag and lead inside the block
-  vector<pair<int, pair<int, int>>> variable_block_lead_lag;
-  //! List for each equation its block number
-  vector<int> equation_block;
+  /*! tuple<static, forward, backward, mixed> */
+  vector<tuple<int, int, int, int>> block_col_type;
 
   //!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
   vector<pair<int, int>> endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
@@ -246,19 +240,19 @@ public:
   BlockSimulationType
   getBlockSimulationType(int block_number) const override
   {
-    return (block_type_firstequation_size_mfs[block_number].first.first);
+    return (get<0>(block_type_firstequation_size_mfs[block_number]));
   };
   //! Return the first equation number of a block
   unsigned int
   getBlockFirstEquation(int block_number) const override
   {
-    return (block_type_firstequation_size_mfs[block_number].first.second);
+    return (get<1>(block_type_firstequation_size_mfs[block_number]));
   };
   //! Return the size of the block block_number
   unsigned int
   getBlockSize(int block_number) const override
   {
-    return (block_type_firstequation_size_mfs[block_number].second.first);
+    return (get<2>(block_type_firstequation_size_mfs[block_number]));
   };
   //! Return the number of exogenous variable in the block block_number
   unsigned int
@@ -276,7 +270,7 @@ public:
   unsigned int
   getBlockMfs(int block_number) const override
   {
-    return (block_type_firstequation_size_mfs[block_number].second.second);
+    return (get<3>(block_type_firstequation_size_mfs[block_number]));
   };
   //! Return the maximum lag in a block
   unsigned int
@@ -294,37 +288,37 @@ public:
   EquationType
   getBlockEquationType(int block_number, int equation_number) const override
   {
-    return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first);
+    return (equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first);
   };
   //! Return true if the equation has been normalized
   bool
   isBlockEquationRenormalized(int block_number, int equation_number) const override
   {
-    return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first == E_EVALUATE_S);
+    return (equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first == E_EVALUATE_S);
   };
   //! Return the expr_t of the equation equation_number belonging to the block block_number
   expr_t
   getBlockEquationExpr(int block_number, int equation_number) const override
   {
-    return (equations[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]]);
+    return (equations[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]]);
   };
   //! Return the expr_t of the renormalized equation equation_number belonging to the block block_number
   expr_t
   getBlockEquationRenormalizedExpr(int block_number, int equation_number) const override
   {
-    return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].second);
+    return (equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].second);
   };
   //! Return the original number of equation equation_number belonging to the block block_number
   int
   getBlockEquationID(int block_number, int equation_number) const override
   {
-    return (equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]);
+    return (equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]);
   };
   //! Return the original number of variable variable_number belonging to the block block_number
   int
   getBlockVariableID(int block_number, int variable_number) const override
   {
-    return (variable_reordered[block_type_firstequation_size_mfs[block_number].first.second+variable_number]);
+    return (variable_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+variable_number]);
   };
   //! Return the original number of the exogenous variable varexo_number belonging to the block block_number
   int
@@ -336,13 +330,13 @@ public:
   int
   getBlockInitialEquationID(int block_number, int equation_number) const override
   {
-    return ((int) inv_equation_reordered[equation_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
+    return ((int) inv_equation_reordered[equation_number] - (int) get<1>(block_type_firstequation_size_mfs[block_number]));
   };
   //! Return the position of variable_number in the block number belonging to the block block_number
   int
   getBlockInitialVariableID(int block_number, int variable_number) const override
   {
-    return ((int) inv_variable_reordered[variable_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
+    return ((int) inv_variable_reordered[variable_number] - (int) get<1>(block_type_firstequation_size_mfs[block_number]));
   };
   //! Return the position of variable_number in the block number belonging to the block block_number
   int