Bytecode: rework message verbosity

Now uses options_.verbosity to decide what to print:
– if options_.verbosity == 0, prints nothing
– if options_.verbosity >= 1, prints iteration counter and duration, and fatal errors
– if options_.verbosity >= 2, additionally print floating point exceptions and
  details about algorithmic decisions

matlab/perfect-foresight-models/perfect_foresight_solver_core.m

@@ -65,7 +65,9 @@ if options_.block
             y = bytecode('dynamic', 'block_decomposed', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state,1, periods+2), periods);
             success = true;
         catch ME
-            disp(ME.message)
+            if options_.verbosity >= 1
+                disp(ME.message)
+            end
             if options_.no_homotopy
                 error('Error in bytecode')
             end
@@ -81,7 +83,9 @@ else
             y = bytecode('dynamic', oo_.endo_simul, oo_.exo_simul, M_.params, repmat(oo_.steady_state, 1, periods+2), periods);
             success = true;
         catch ME
-            disp(ME.message)
+            if options_.verbosity >= 1
+                disp(ME.message)
+            end
             if options_.no_homotopy
                 error('Error in bytecode')
             end

mex/sources/bytecode/Interpreter.cc

@@ -32,9 +32,9 @@ Interpreter::Interpreter(Evaluate &evaluator_arg, double *params_arg, double *y_
                          size_t nb_row_x_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg,
                          int maxit_arg_, double solve_tolf_arg, int y_decal_arg, double markowitz_c_arg,
                          string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg,
-                         bool global_temporary_terms_arg, bool print_arg, bool print_error_arg, mxArray *GlobalTemporaryTerms_arg,
-                         bool steady_state_arg, bool block_decomposed_arg, bool print_it_arg, int col_x_arg, int col_y_arg, const BasicSymbolTable &symbol_table_arg)
-: dynSparseMatrix {evaluator_arg, y_size_arg, y_kmin_arg, y_kmax_arg, print_it_arg, steady_state_arg, block_decomposed_arg, periods_arg, minimal_solving_periods_arg, symbol_table_arg, print_error_arg}
+                         bool global_temporary_terms_arg, bool print_arg, mxArray *GlobalTemporaryTerms_arg,
+                         bool steady_state_arg, bool block_decomposed_arg, int col_x_arg, int col_y_arg, const BasicSymbolTable &symbol_table_arg, int verbosity_arg)
+: dynSparseMatrix {evaluator_arg, y_size_arg, y_kmin_arg, y_kmax_arg, steady_state_arg, block_decomposed_arg, periods_arg, minimal_solving_periods_arg, symbol_table_arg, verbosity_arg}
 {
   params = params_arg;
   y = y_arg;
@@ -60,7 +60,6 @@ Interpreter::Interpreter(Evaluate &evaluator_arg, double *params_arg, double *y_
   col_x = col_x_arg;
   col_y = col_y_arg;
   GlobalTemporaryTerms = GlobalTemporaryTerms_arg;
-  print_it = print_it_arg;
 }
 
 void
@@ -107,13 +106,14 @@ Interpreter::solve_simple_one_periods()
               if (!isfinite(res1))
                 {
                   slowc /= 1.5;
-                  mexPrintf("Reducing the path length in Newton step slowc=%f\n", slowc);
+                  if (verbosity >= 2)
+                    mexPrintf("Reducing the path length in Newton step slowc=%f\n", slowc);
                   feclearexcept(FE_ALL_EXCEPT);
                   y[Block_Contain[0].Variable + Per_y_] = ya - slowc * (r[0] / g1[0]);
                   if (fetestexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW))
                     {
                       res1 = numeric_limits<double>::quiet_NaN();
-                      if (print_error)
+                      if (verbosity >= 1)
                         mexPrintf("      Singularity in block %d", block_num+1);
                     }
                 }
@@ -129,7 +129,7 @@ Interpreter::solve_simple_one_periods()
       if (fetestexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW))
         {
           res1 = numeric_limits<double>::quiet_NaN();
-          if (print_error)
+          if (verbosity >= 1)
             mexPrintf("      Singularity in block %d", block_num+1);
         }
       iter++;
@@ -518,7 +518,8 @@ Interpreter::simulate_a_block(const vector_table_conditional_local_type &vector_
 #endif
       if (steady_state)
         {
-          mexPrintf("SOLVE TWO BOUNDARIES in a steady state model: impossible case\n");
+          if (verbosity >= 1)
+            mexPrintf("SOLVE TWO BOUNDARIES in a steady state model: impossible case\n");
           return ERROR_ON_EXIT;
         }
       if (vector_table_conditional_local.size())
@@ -728,8 +729,8 @@ Interpreter::MainLoop(const string &bin_basename, bool evaluate, int block, bool
       else
         {
 #ifdef DEBUG
-          mexPrintf("endo in block %d, size=%d, type=%d, steady_state=%d, print_it=%d, is_linear=%d, endo_nbr=%d, u_count_int=%d\n",
-                    current_block+1, size, type, steady_state, print_it, is_linear, symbol_table_endo_nbr, u_count_int);
+          mexPrintf("endo in block %d, size=%d, type=%d, steady_state=%d, is_linear=%d, endo_nbr=%d, u_count_int=%d\n",
+                    current_block+1, size, type, steady_state, is_linear, symbol_table_endo_nbr, u_count_int);
 #endif
           bool result;
           if (sconstrained_extended_path.size())
@@ -822,8 +823,8 @@ Interpreter::extended_path(const string &file_name, bool evaluate, int block, in
     y_save[i] = y[i];
   if (endo_name_length_l < 8)
     endo_name_length_l = 8;
-  bool old_print_it = print_it;
-  print_it = false;
+  int old_verbosity {verbosity};
+  verbosity = 0;
   ostringstream res1;
   res1 << std::scientific << 2.54656875434865131;
   int real_max_length = res1.str().length();
@@ -832,7 +833,7 @@ Interpreter::extended_path(const string &file_name, bool evaluate, int block, in
   string line;
   line.insert(line.begin(), table_length, '-');
   line.insert(line.length(), "\n");
-  if (old_print_it)
+  if (old_verbosity >= 1)
     {
       mexPrintf("\nExtended Path simulation:\n");
       mexPrintf("-------------------------\n");
@@ -846,7 +847,7 @@ Interpreter::extended_path(const string &file_name, bool evaluate, int block, in
   for (int t = 0; t < nb_periods; t++)
     {
       previous_block_exogenous.clear();
-      if (old_print_it)
+      if (old_verbosity >= 1)
         {
           mexPrintf("|%s|", elastic(dates[t], date_length+2, false).c_str());
           mexEvalString("drawnow;");
@@ -871,7 +872,7 @@ Interpreter::extended_path(const string &file_name, bool evaluate, int block, in
             x[y_kmin + j * nb_row_x] = x_save[t + 1 + y_kmin + j * nb_row_x];
         }
 
-      if (old_print_it)
+      if (old_verbosity >= 1)
         {
           ostringstream res1;
           res1 << std::scientific << max_res;
@@ -880,7 +881,7 @@ Interpreter::extended_path(const string &file_name, bool evaluate, int block, in
           mexEvalString("drawnow;");
         }
     }
-  print_it = old_print_it;
+  verbosity = old_verbosity;
   for (int i = 0; i < y_size * col_y; i++)
     y[i] = y_save[i];
   for (int j = 0; j < col_x * nb_row_x; j++)

mex/sources/bytecode/Interpreter.hh

@@ -38,7 +38,7 @@ class Interpreter : public dynSparseMatrix
 private:
   vector<int> previous_block_exogenous;
   bool global_temporary_terms;
-  bool print;
+  bool print; // Whether the “print” command is requested
   int col_x, col_y;
   vector<double> residual;
   void evaluate_over_periods(bool forward);
@@ -56,8 +56,8 @@ public:
               size_t nb_row_x_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg,
               int maxit_arg_, double solve_tolf_arg, int y_decal_arg, double markowitz_c_arg,
               string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg,
-              bool global_temporary_terms_arg, bool print_arg, bool print_error_arg, mxArray *GlobalTemporaryTerms_arg,
-              bool steady_state_arg, bool block_decomposed_arg, bool print_it_arg, int col_x_arg, int col_y_arg, const BasicSymbolTable &symbol_table_arg);
+              bool global_temporary_terms_arg, bool print_arg, mxArray *GlobalTemporaryTerms_arg,
+              bool steady_state_arg, bool block_decomposed_arg, int col_x_arg, int col_y_arg, const BasicSymbolTable &symbol_table_arg, int verbosity_arg);
   pair<bool, vector<int>> extended_path(const string &file_name, bool evaluate, int block, int nb_periods, const vector<s_plan> &sextended_path, const vector<s_plan> &sconstrained_extended_path, const vector<string> &dates, const table_conditional_global_type &table_conditional_global);
   pair<bool, vector<int>> compute_blocks(const string &file_name, bool evaluate, int block);
   void check_for_controlled_exo_validity(int current_block, const vector<s_plan> &sconstrained_extended_path);

mex/sources/bytecode/SparseMatrix.cc

@@ -24,20 +24,18 @@
 
 #include "SparseMatrix.hh"
 
-dynSparseMatrix::dynSparseMatrix(Evaluate &evaluator_arg, int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, bool block_decomposed_arg, int periods_arg,
-                                 int minimal_solving_periods_arg, const BasicSymbolTable &symbol_table_arg,
-                                 bool print_error_arg) :
+dynSparseMatrix::dynSparseMatrix(Evaluate &evaluator_arg, int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool steady_state_arg, bool block_decomposed_arg, int periods_arg,
+                                 int minimal_solving_periods_arg, const BasicSymbolTable &symbol_table_arg, int verbosity_arg) :
   symbol_table {symbol_table_arg},
   steady_state {steady_state_arg},
   block_decomposed {block_decomposed_arg},
   evaluator {evaluator_arg},
   minimal_solving_periods {minimal_solving_periods_arg},
-  print_it {print_it_arg},
   y_size {y_size_arg},
   y_kmin {y_kmin_arg},
   y_kmax {y_kmax_arg},
   periods {periods_arg},
-  print_error {print_error_arg}
+  verbosity {verbosity_arg}
 {
   pivotva = nullptr;
   g_save_op = nullptr;
@@ -1812,7 +1810,7 @@ dynSparseMatrix::compute_block_time(int Per_u_, bool evaluate, bool no_derivativ
   catch (FloatingPointException &e)
     {
       res1 = numeric_limits<double>::quiet_NaN();
-      if (print_error)
+      if (verbosity >= 2)
         mexPrintf("%s\n      %s\n", e.message.c_str(), e.location.c_str());
     }
 }
@@ -1894,7 +1892,8 @@ dynSparseMatrix::compute_complete(double lambda, double *crit)
         }
       it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
     }
-  mexPrintf("  lambda=%e, res2=%e\n", lambda, res2_);
+  if (verbosity >= 2)
+    mexPrintf("  lambda=%e, res2=%e\n", lambda, res2_);
   *crit = res2_/2;
   return true;
 }
@@ -1908,7 +1907,8 @@ dynSparseMatrix::mnbrak(double *ax, double *bx, double *cx, double *fa, double *
 
   auto sign = [](double a, double b) { return b >= 0.0 ? fabs(a) : -fabs(a); };
 
-  mexPrintf("bracketing *ax=%f, *bx=%f\n", *ax, *bx);
+  if (verbosity >= 2)
+    mexPrintf("bracketing *ax=%f, *bx=%f\n", *ax, *bx);
   if (!compute_complete(*ax, fa))
     return false;
   if (!compute_complete(*bx, fb))
@@ -1993,7 +1993,8 @@ dynSparseMatrix::golden(double ax, double bx, double cx, double tol, double solv
 {
   const double R = 0.61803399;
   const double C = (1.0-R);
-  mexPrintf("golden\n");
+  if (verbosity >= 2)
+    mexPrintf("golden\n");
   int iter = 0, max_iter = 100;
   double f1, f2, x1, x2;
   double x0 = ax;
@@ -2890,10 +2891,13 @@ dynSparseMatrix::Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, i
           mxFree(bc);
           if (steady_state)
             {
-              if (blck > 1)
-                mexPrintf("Error: singular system in Simulate_NG in block %d\n", blck+1);
-              else
-                mexPrintf("Error: singular system in Simulate_NG");
+              if (verbosity >= 1)
+                {
+                  if (blck > 1)
+                    mexPrintf("Error: singular system in Simulate_NG in block %d\n", blck+1);
+                  else
+                    mexPrintf("Error: singular system in Simulate_NG");
+                }
               return true;
             }
           else
@@ -3217,9 +3221,9 @@ dynSparseMatrix::Solve_ByteCode_Symbolic_Sparse_GaussianElimination(int Size, bo
                         NR_max = NR[j];
                       }
                   }
-              if (fabs(piv) < eps)
+              if (fabs(piv) < eps && verbosity >= 1)
                 mexPrintf("==> Error NR_max=%d, N_max=%d and piv=%f, piv_abs=%f, markovitz_max=%f\n", NR_max, N_max, piv, piv_abs, markovitz_max);
-              if (NR_max == 0)
+              if (NR_max == 0 && verbosity >= 1)
                 mexPrintf("==> Error NR_max=0 and piv=%f, markovitz_max=%f\n", piv, markovitz_max);
               pivot[i] = pivj;
               pivot_save[i] = pivj;
@@ -3676,7 +3680,8 @@ dynSparseMatrix::Check_and_Correct_Previous_Iteration(int y_size, int size)
             }
           compute_complete(true, res1, res2, max_res, max_res_idx);
         }
-      mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save);
+      if (verbosity >= 2)
+        mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n", slowc_save);
       for (int i = 0; i < size; i++)
         {
           int eq = index_vara[i];
@@ -3726,11 +3731,14 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int size)
 #endif
       if (steady_state)
         {
-          if (iter == 0)
-            mexPrintf(" the initial values of endogenous variables are too far from the solution.\nChange them!\n");
-          else
-            mexPrintf(" dynare cannot improve the simulation in block %d at time %d (variable %d)\n", block_num+1, it_+1, index_vara[max_res_idx]+1);
-          mexEvalString("drawnow;");
+          if (verbosity >= 1)
+            {
+              if (iter == 0)
+                mexPrintf(" the initial values of endogenous variables are too far from the solution.\nChange them!\n");
+              else
+                mexPrintf(" dynare cannot improve the simulation in block %d at time %d (variable %d)\n", block_num+1, it_+1, index_vara[max_res_idx]+1);
+              mexEvalString("drawnow;");
+            }
         }
       else
         {
@@ -3743,7 +3751,7 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int size)
         }
     }
 
-  if (print_it)
+  if (verbosity >= 1)
     {
       if (steady_state)
         {
@@ -3848,7 +3856,7 @@ dynSparseMatrix::solve_linear(int block_num, int y_size, int size, int iter)
       if (iter)
         Check_and_Correct_Previous_Iteration(y_size, size);
       bool singular_system = Simulate_One_Boundary(block_num, y_size, size);
-      if (singular_system)
+      if (singular_system && verbosity >= 1)
         Singular_display(block_num, size);
     }
   return cvg;
@@ -3984,20 +3992,12 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
   double *Ax = nullptr, *b;
   SuiteSparse_long *Ap = nullptr, *Ai = nullptr;
 
-  if (iter > 0)
-    {
-      if (print_it)
-        {
-          mexPrintf("Sim : %f ms\n", (1000.0*(static_cast<double>(clock())-static_cast<double>(time00)))/static_cast<double>(CLOCKS_PER_SEC));
-          mexEvalString("drawnow;");
-        }
-      time00 = clock();
-    }
   if (isnan(res1) || isinf(res1) || (res2 > 12*g0 && iter > 0))
     {
       if (iter == 0 || fabs(slowc_save) < 1e-8)
         {
-          mexPrintf("res1 = %f, res2 = %f g0 = %f iter = %d\n", res1, res2, g0, iter);
+          if (verbosity >= 2)
+            mexPrintf("res1 = %f, res2 = %f g0 = %f iter = %d\n", res1, res2, g0, iter);
           for (int j = 0; j < y_size; j++)
             {
               bool select = false;
@@ -4007,10 +4007,13 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
                     select = true;
                     break;
                   }
-              if (select)
-                mexPrintf("-> variable %s (%d) at time %d = %f direction = %f\n", symbol_table.getName(SymbolType::endogenous, j).c_str(), j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
-              else
-                mexPrintf("   variable %s (%d) at time %d = %f direction = %f\n", symbol_table.getName(SymbolType::endogenous, j).c_str(), j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
+              if (verbosity >= 2)
+                {
+                  if (select)
+                    mexPrintf("-> variable %s (%d) at time %d = %f direction = %f\n", symbol_table.getName(SymbolType::endogenous, j).c_str(), j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
+                  else
+                    mexPrintf("   variable %s (%d) at time %d = %f direction = %f\n", symbol_table.getName(SymbolType::endogenous, j).c_str(), j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]);
+                }
             }
           if (iter == 0)
             throw FatalException{"In Simulate_Newton_Two_Boundaries, the initial values of endogenous variables are too far from the solution. Change them!"};
@@ -4060,7 +4063,7 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
           prev_slowc_save = slowc_save;
           slowc_save /= 1.05;
         }
-      if (print_it)
+      if (verbosity >= 2)
         {
           if (isnan(res1) || isinf(res1))
             mexPrintf("The model cannot be evaluated, trying to correct it using slowc=%f\n", slowc_save);
@@ -4077,7 +4080,8 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
     {
       if (alt_symbolic && alt_symbolic_count < alt_symbolic_count_max)
         {
-          mexPrintf("Pivoting method will be applied only to the first periods.\n");
+          if (verbosity >= 2)
+            mexPrintf("Pivoting method will be applied only to the first periods.\n");
           alt_symbolic = false;
           symbolic = true;
           markowitz_c = markowitz_c_s;
@@ -4087,7 +4091,8 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
         {
           if (restart > 2)
             {
-              mexPrintf("Divergence or slowdown occurred during simulation.\nIn the next iteration, pivoting method will be applied to all periods.\n");
+              if (verbosity >= 2)
+                mexPrintf("Divergence or slowdown occurred during simulation.\nIn the next iteration, pivoting method will be applied to all periods.\n");
               symbolic = false;
               alt_symbolic = true;
               markowitz_c_s = markowitz_c;
@@ -4095,7 +4100,8 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
             }
           else
             {
-              mexPrintf("Divergence or slowdown occurred during simulation.\nIn the next iteration, pivoting method will be applied for a longer period.\n");
+              if (verbosity >= 2)
+                mexPrintf("Divergence or slowdown occurred during simulation.\nIn the next iteration, pivoting method will be applied for a longer period.\n");
               start_compare = min(tbreak_g, periods);
               restart++;
             }
@@ -4107,7 +4113,7 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
         }
     }
   res1a = res1;
-  if (print_it)
+  if (verbosity >= 1)
     {
       if (iter == 0)
         {
@@ -4186,7 +4192,7 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
       else if (stack_solve_algo == 5)
         Solve_ByteCode_Symbolic_Sparse_GaussianElimination(Size, symbolic, blck);
     }
-  if (print_it)
+  if (verbosity >= 1)
     {
       clock_t t2 = clock();
       mexPrintf("(** %f milliseconds **)\n", 1000.0*(static_cast<double>(t2) - static_cast<double>(t1))/static_cast<double>(CLOCKS_PER_SEC));
@@ -4202,11 +4208,13 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
       if (!golden(ax, bx, cx, 1e-1, solve_tolf, &xmin))
         return;
       slowc = xmin;
-      clock_t t3 = clock();
-      mexPrintf("(** %f milliseconds **)\n", 1000.0*(static_cast<double>(t3) - static_cast<double>(t2))/static_cast<double>(CLOCKS_PER_SEC));
-      mexEvalString("drawnow;");
+      if (verbosity >= 1)
+        {
+          clock_t t3 = clock();
+          mexPrintf("(** %f milliseconds **)\n", 1000.0*(static_cast<double>(t3) - static_cast<double>(t2))/static_cast<double>(CLOCKS_PER_SEC));
+          mexEvalString("drawnow;");
+        }
     }
-  time00 = clock();
   if (tbreak_g == 0)
     tbreak_g = periods;
 }

mex/sources/bytecode/SparseMatrix.hh

@@ -68,7 +68,7 @@ constexpr double mem_increasing_factor = 1.1;
 class dynSparseMatrix
 {
 public:
-  dynSparseMatrix(Evaluate &evaluator_arg, int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, bool block_decomposed_arg, int periods_arg, int minimal_solving_periods_arg, const BasicSymbolTable &symbol_table_arg, bool print_error_arg);
+  dynSparseMatrix(Evaluate &evaluator_arg, int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool steady_state_arg, bool block_decomposed_arg, int periods_arg, int minimal_solving_periods_arg, const BasicSymbolTable &symbol_table_arg, int verbosity_arg);
   void Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin, int y_kmax, int Size, int periods, bool cvg, int minimal_solving_periods, int stack_solve_algo, const vector_table_conditional_local_type &vector_table_conditional_local);
   void Simulate_Newton_One_Boundary(bool forward);
   void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1);
@@ -155,7 +155,6 @@ protected:
   fstream SaveCode;
   string filename;
   int max_u, min_u;
-  clock_t time00;
 
   Mem_Mngr mem_mngr;
   vector<int> u_liste;
@@ -196,7 +195,6 @@ protected:
 
   int stack_solve_algo, solve_algo;
   int minimal_solving_periods;
-  bool print_it;
   int Per_u_, Per_y_;
   int maxit_;
   double *direction;
@@ -230,12 +228,12 @@ protected:
   int u_count_int;
   vector<Block_contain_type> Block_Contain;
 
+  int verbosity; // Corresponds to options_.verbosity
+
   void compute_block_time(int Per_u_, bool evaluate, bool no_derivatives);
   bool compute_complete(bool no_derivatives, double &res1, double &res2, double &max_res, int &max_res_idx);
 
   bool compute_complete(double lambda, double *crit);
-
-  bool print_error; // Whether to stop processing on floating point exceptions
 };
 
 #endif // _SPARSEMATRIX_HH

mex/sources/bytecode/bytecode.cc

@@ -17,7 +17,6 @@
  * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
  */
 
-#include <ctime>
 #include <cmath>
 #include <type_traits>
 #include <algorithm>
@@ -63,7 +62,6 @@ Get_Arguments_and_global_variables(int nrhs,
                                    bool &evaluate, int &block,
                                    mxArray *M_[], mxArray *oo_[], mxArray *options_[], bool &global_temporary_terms,
                                    bool &print,
-                                   bool &print_error,
                                    mxArray *GlobalTemporaryTerms[],
                                    string *plan_struct_name, string *pfplan_struct_name, bool *extended_path, mxArray *ep_struct[])
 {
@@ -126,8 +124,6 @@ Get_Arguments_and_global_variables(int nrhs,
           global_temporary_terms = true;
         else if (Get_Argument(prhs[i]) == "print")
           print = true;
-        else if (Get_Argument(prhs[i]) == "no_print_error")
-          print_error = false;
         else
           {
             pos = 0;
@@ -215,7 +211,8 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   double *yd = nullptr, *xd = nullptr;
   int count_array_argument = 0;
   bool global_temporary_terms = false;
-  bool print = false, print_error = true, print_it = false;
+  bool print = false; // Whether the “print” command is requested
+  int verbosity {1}; // Corresponds to options_.verbosity
   double *steady_yd = nullptr;
   string plan, pfplan;
   bool extended_path;
@@ -244,7 +241,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
                                          &block_structur,
                                          steady_state, block_decomposed, evaluate, block,
                                          &M_, &oo_, &options_, global_temporary_terms,
-                                         print, print_error, &GlobalTemporaryTerms,
+                                         print, &GlobalTemporaryTerms,
                                          &plan, &pfplan, &extended_path, &extended_path_struct);
     }
   catch (GeneralException &feh)
@@ -449,19 +446,20 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
                 splan[i].per_value[j] = { ceil(per_value[j]), per_value[j + num_shocks] };
             }
         }
-      int i = 0;
-      for (auto & it : splan)
-        {
-          mexPrintf("----------------------------------------------------------------------------------------------------\n");
-          mexPrintf("surprise #%d\n", i+1);
-          if (it.exo.length())
-            mexPrintf(" plan fliping var=%s (%d) exo=%s (%d) for the following periods and with the following values:\n", it.var.c_str(), it.var_num, it.exo.c_str(), it.exo_num);
-          else
-            mexPrintf(" plan shocks on var=%s for the following periods and with the following values:\n", it.var.c_str());
-          for (auto &[period, value]: it.per_value)
-            mexPrintf("  %3d %10.5f\n", period, value);
-          i++;
-        }
+      if (verbosity >= 1)
+        for (int i {0};
+             auto & it : splan)
+          {
+            mexPrintf("----------------------------------------------------------------------------------------------------\n");
+            mexPrintf("surprise #%d\n", i+1);
+            if (it.exo.length())
+              mexPrintf(" plan fliping var=%s (%d) exo=%s (%d) for the following periods and with the following values:\n", it.var.c_str(), it.var_num, it.exo.c_str(), it.exo_num);
+            else
+              mexPrintf(" plan shocks on var=%s for the following periods and with the following values:\n", it.var.c_str());
+            for (auto &[period, value]: it.per_value)
+              mexPrintf("  %3d %10.5f\n", period, value);
+            i++;
+          }
     }
 
   if (pfplan.length() > 0)
@@ -509,19 +507,20 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
                 spfplan[i].per_value[j] = { ceil(per_value[j]), per_value[j+ num_shocks] };
             }
         }
-      int i = 0;
-      for (auto & it : spfplan)
-        {
-          mexPrintf("----------------------------------------------------------------------------------------------------\n");
-          mexPrintf("perfect foresight #%d\n", i+1);
-          if (it.exo.length())
-            mexPrintf(" plan flipping var=%s (%d) exo=%s (%d) for the following periods and with the following values:\n", it.var.c_str(), it.var_num, it.exo.c_str(), it.exo_num);
-          else
-            mexPrintf(" plan shocks on var=%s (%d) for the following periods and with the following values:\n", it.var.c_str(), it.var_num);
-          for (auto &[period, value] : it.per_value)
-            mexPrintf("  %3d %10.5f\n", period, value);
-          i++;
-        }
+      if (verbosity >= 1)
+        for (int i {0};
+             auto & it : spfplan)
+          {
+            mexPrintf("----------------------------------------------------------------------------------------------------\n");
+            mexPrintf("perfect foresight #%d\n", i+1);
+            if (it.exo.length())
+              mexPrintf(" plan flipping var=%s (%d) exo=%s (%d) for the following periods and with the following values:\n", it.var.c_str(), it.var_num, it.exo.c_str(), it.exo_num);
+            else
+              mexPrintf(" plan shocks on var=%s (%d) for the following periods and with the following values:\n", it.var.c_str(), it.var_num);
+            for (auto &[period, value] : it.per_value)
+              mexPrintf("  %3d %10.5f\n", period, value);
+            i++;
+          }
     }
 
   int field_steady_state = mxGetFieldNumber(oo_, "steady_state");
@@ -601,13 +600,11 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
         }
     }
   int field = mxGetFieldNumber(options_, "verbosity");
-  int verbose = 0;
   if (field >= 0)
-    verbose = static_cast<int>(*mxGetPr((mxGetFieldByNumber(options_, 0, field))));
+    verbosity = static_cast<int>(mxGetScalar(mxGetFieldByNumber(options_, 0, field)));
   else
     mexErrMsgTxt("verbosity is not a field of options_");
-  if (verbose)
-    print_it = true;
+
   if (!steady_state)
     field = mxGetFieldNumber(options_, "simul");
   else
@@ -710,8 +707,6 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   copy_n(yd, row_y*col_y, y);
   copy_n(yd, row_y*col_y, ya);
 
-  clock_t t0 = clock();
-
   const filesystem::path codfile {file_name + "/model/bytecode/" + (block_decomposed ? "block/" : "")
     + (steady_state ? "static" : "dynamic") + ".cod"};
   Evaluate evaluator {codfile, steady_state, symbol_table};
@@ -719,8 +714,8 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   Interpreter interprete {evaluator, params, y, ya, x, steady_yd, direction, row_y, row_x,
                           periods, y_kmin, y_kmax, maxit_, solve_tolf, y_decal,
                           markowitz_c, file_name, minimal_solving_periods, stack_solve_algo,
-                          solve_algo, global_temporary_terms, print, print_error, GlobalTemporaryTerms,
-                          steady_state, block_decomposed, print_it, col_x, col_y, symbol_table};
+                          solve_algo, global_temporary_terms, print, GlobalTemporaryTerms,
+                          steady_state, block_decomposed, col_x, col_y, symbol_table, verbosity};
   double *pind;
   bool r;
   vector<int> blocks;
@@ -739,10 +734,6 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
       mexErrMsgTxt(feh.message.c_str());
     }
 
-  clock_t t1 = clock();
-  if (!steady_state && !evaluate && print)
-    mexPrintf("Simulation Time=%f milliseconds\n",
-              1000.0*(static_cast<double>(t1)-static_cast<double>(t0))/static_cast<double>(CLOCKS_PER_SEC));
   bool dont_store_a_structure = false;
   if (nlhs > 0)
     {