From 153c7b2542a0412f0be7f71657574fbf04c0df10 Mon Sep 17 00:00:00 2001
From: ferhat <ferhat@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Fri, 11 Sep 2009 17:06:54 +0000
Subject: [PATCH] New option: "evaluate' to evaluate the static or the dynamic
 model

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2919 ac1d8469-bf42-47a9-8791-bf33cf982152
---
 mex/sources/bytecode/Interpreter.cc | 297 ++++++++++++++++++++++------
 mex/sources/bytecode/Interpreter.hh |   7 +-
 mex/sources/bytecode/bytecode.cc    | 103 +++++++---
 3 files changed, 324 insertions(+), 83 deletions(-)

diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index 6441a339d..dc145cd66 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -128,7 +128,7 @@ Interpreter::log1(double a)
 
 
 void
-Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
+Interpreter::compute_block_time(int Per_u_, bool evaluate) /*throw(EvalException)*/
 {
   int var, lag, op;
   ostringstream tmp_out;
@@ -154,7 +154,10 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
                 case eEndogenous :
                   var=get_code_int;
                   lag=get_code_int;
-                  Stack.push(y[(it_+lag)*y_size+var]);
+                  if(evaluate)
+                    Stack.push(ya[(it_+lag)*y_size+var]);
+                  else
+                    Stack.push(y[(it_+lag)*y_size+var]);
 #ifdef DEBUG
                   tmp_out << " y[" << it_+lag << ", " << var << "](" << y[(it_+lag)*y_size+var] << ")";
 #endif
@@ -189,7 +192,10 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
                   break;
                 case eEndogenous :
                   var=get_code_int;
-                  Stack.push(y[var]);
+                  if(evaluate)
+                    Stack.push(ya[var]);
+                  else
+                    Stack.push(y[var]);
 #ifdef DEBUG
                   tmp_out << " y[" << var << "](" << y[var] << ")";
 #endif
@@ -640,6 +646,178 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
     }
 }
 
+void
+Interpreter::evaluate_a_block(int size, int type, string bin_basename, bool steady_state, int block_num)
+{
+  char *begining;
+  int max_lag_plus_max_lead_plus_1;
+  int symbol_table_endo_nbr;
+  int Block_List_Max_Lag;
+  int Block_List_Max_Lead;
+  int u_count_int;
+  bool is_linear;
+  //mexPrintf("Block = %d type = %d, size = %d, steady_state = %d\n",block_num, type, size, steady_state);
+  //mexPrintf("y_kmin=%d periods=%d\n",y_kmin, periods);
+  //mexEvalString("drawnow;");
+  if(type == SOLVE_FORWARD_COMPLETE || type == SOLVE_BACKWARD_COMPLETE || type == SOLVE_TWO_BOUNDARIES_SIMPLE || type == SOLVE_TWO_BOUNDARIES_COMPLETE)
+    {
+      is_linear=get_code_bool;
+      max_lag_plus_max_lead_plus_1=get_code_int;
+      symbol_table_endo_nbr=get_code_int;
+      Block_List_Max_Lag=get_code_int;
+      Block_List_Max_Lead=get_code_int;
+      u_count_int=get_code_int;
+    }
+  switch (type)
+    {
+      case EVALUATE_FORWARD :
+        if(steady_state)
+					compute_block_time(0, true);
+				else
+				  {
+            begining=get_code_pointer;
+            for (it_=y_kmin;it_<periods+y_kmin;it_++)
+              {
+                //mexPrintf("it_=%d, y_size=%d\n",it_, y_size);
+                set_code_pointer(begining);
+                Per_y_=it_*y_size;
+                compute_block_time(0, true);
+              }
+				  }
+        break;
+      case SOLVE_FORWARD_SIMPLE :
+        g1=(double*)mxMalloc(size*size*sizeof(double));
+				r=(double*)mxMalloc(size*sizeof(double));
+				if(steady_state)
+				  {
+					  compute_block_time(0, true);
+					  for(int j=0; j<size; j++)
+              y[Block_Contain[j].Variable] += r[j];
+				  }
+				else
+				  {
+            begining=get_code_pointer;
+            for (it_=y_kmin;it_<periods+y_kmin;it_++)
+              {
+                //mexPrintf("it_=%d, y_size=%d\n",it_, y_size);
+                set_code_pointer(begining);
+                Per_y_=it_*y_size;
+                compute_block_time(0, true);
+                for(int j=0; j<size; j++)
+                  y[it_*y_size+Block_Contain[j].Variable] += r[j];
+              }
+				  }
+        mxFree(g1);
+        mxFree(r);
+        break;
+      case SOLVE_FORWARD_COMPLETE :
+        fixe_u(&u, u_count_int, u_count_int);
+        Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false);
+        r=(double*)mxMalloc(size*sizeof(double));
+				if(steady_state)
+				  {
+					  compute_block_time(0, true);
+					  for(int j=0; j<size; j++)
+              y[Block_Contain[j].Variable] += r[j];
+				  }
+				else
+				  {
+            begining=get_code_pointer;
+            for (it_=y_kmin;it_<periods+y_kmin;it_++)
+              {
+                //mexPrintf("it_=%d, y_size=%d\n",it_, y_size);
+                set_code_pointer(begining);
+                Per_y_=it_*y_size;
+                compute_block_time(0, true);
+                for(int j=0; j<size; j++)
+                  y[it_*y_size+Block_Contain[j].Variable] += r[j];
+              }
+				  }
+        mxFree(r);
+        break;
+      case EVALUATE_BACKWARD :
+        if(steady_state)
+          compute_block_time(0, true);
+				else
+				  {
+            begining=get_code_pointer;
+            for (it_=periods+y_kmin-1;it_>=y_kmin;it_--)
+              {
+                set_code_pointer(begining);
+                Per_y_=it_*y_size;
+                compute_block_time(0, true);
+              }
+				  }
+        break;
+      case SOLVE_BACKWARD_SIMPLE :
+        g1=(double*)mxMalloc(size*size*sizeof(double));
+				r=(double*)mxMalloc(size*sizeof(double));
+        if(steady_state)
+          {
+            compute_block_time(0, true);
+            for(int j=0; j<size; j++)
+              y[Block_Contain[j].Variable] += r[j];
+          }
+				else
+				  {
+            begining=get_code_pointer;
+            for (it_=periods+y_kmin-1;it_>=y_kmin;it_--)
+              {
+                set_code_pointer(begining);
+                Per_y_=it_*y_size;
+                compute_block_time(0,true);
+                for(int j=0; j<size; j++)
+                  y[it_*y_size+Block_Contain[j].Variable] += r[j];
+              }
+				  }
+        mxFree(g1);
+        mxFree(r);
+        break;
+      case SOLVE_BACKWARD_COMPLETE :
+        fixe_u(&u, u_count_int, u_count_int);
+        Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false);
+        r=(double*)mxMalloc(size*sizeof(double));
+        if(steady_state)
+          {
+            compute_block_time(0, true);
+            for(int j=0; j<size; j++)
+              y[Block_Contain[j].Variable] += r[j];
+          }
+				else
+				  {
+            begining=get_code_pointer;
+            for (it_=periods+y_kmin-1;it_>=y_kmin;it_--)
+              {
+                set_code_pointer(begining);
+                Per_y_=it_*y_size;
+                compute_block_time(0,true);
+                for(int j=0; j<size; j++)
+                  y[it_*y_size+Block_Contain[j].Variable] += r[j];
+              }
+				  }
+        mxFree(r);
+        break;
+      case SOLVE_TWO_BOUNDARIES_SIMPLE :
+      case SOLVE_TWO_BOUNDARIES_COMPLETE:
+        fixe_u(&u, u_count_int, u_count_int);
+        Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state, true);
+        u_count=u_count_int*(periods+y_kmax+y_kmin);
+        r=(double*)mxMalloc(size*sizeof(double));
+        begining=get_code_pointer;
+        for (it_=y_kmin;it_<periods+y_kmin;it_++)
+          {
+            Per_u_=(it_-y_kmin)*u_count_int;
+            Per_y_=it_*y_size;
+            set_code_pointer(begining);
+            compute_block_time(Per_u_, true);
+            for(int j=0; j<size; j++)
+              y[it_*y_size+Block_Contain[j].Variable] += r[j];
+          }
+        mxFree(r);
+        break;
+    }
+}
+
 bool
 Interpreter::simulate_a_block(int size,int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num)
 {
@@ -659,7 +837,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
     {
       case EVALUATE_FORWARD :
         if(steady_state)
-					compute_block_time(0);
+					compute_block_time(0, false);
 				else
 				  {
             begining=get_code_pointer;
@@ -667,13 +845,13 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
               {
                 set_code_pointer(begining);
                 Per_y_=it_*y_size;
-                compute_block_time(0);
+                compute_block_time(0, false);
               }
 				  }
         break;
       case EVALUATE_BACKWARD :
         if(steady_state)
-          compute_block_time(0);
+          compute_block_time(0, false);
 				else
 				  {
             begining=get_code_pointer;
@@ -681,7 +859,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
               {
                 set_code_pointer(begining);
                 Per_y_=it_*y_size;
-                compute_block_time(0);
+                compute_block_time(0, false);
               }
 				  }
         break;
@@ -697,50 +875,50 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
               {
                 set_code_pointer(begining);
                 Per_y_=it_*y_size;
-                    compute_block_time(0);
-                    y[Block_Contain[0].Variable] += -r[0]/g1[0];
+                compute_block_time(0, false);
+                y[Block_Contain[0].Variable] += -r[0]/g1[0];
+                double rr;
+							  rr=r[0];
+                cvg=((rr*rr)<solve_tolf);
+                iter++;
+              }
+            if (!cvg)
+              {
+                mexPrintf("Convergence not achieved in block %d, after %d iterations\n",Block_Count,iter);
+                /*mexEvalString("st=fclose('all');clear all;");
+                mexErrMsgTxt("End of simulate");*/
+                return false;
+              }
+          }
+			  else
+				  {
+            for (it_=y_kmin;it_<periods+y_kmin;it_++)
+              {
+                cvg=false;
+                iter=0;
+                Per_y_=it_*y_size;
+                while (!(cvg||(iter>maxit_)))
+                  {
+                    set_code_pointer(begining);
+                    Per_y_=it_*y_size;
+                    compute_block_time(0, false);
+                    y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0];
                     double rr;
-							    	rr=r[0];
+                    if(fabs(1+y[Per_y_+Block_Contain[0].Variable])>eps)
+				    				  rr=r[0]/(1+y[Per_y_+Block_Contain[0].Variable]);
+						    		else
+								    	rr=r[0];
                     cvg=((rr*rr)<solve_tolf);
                     iter++;
                   }
                 if (!cvg)
                   {
-                    mexPrintf("Convergence not achieved in block %d, after %d iterations\n",Block_Count,iter);
-                    /*mexEvalString("st=fclose('all');clear all;");
-                    mexErrMsgTxt("End of simulate");*/
-                    return false;
+                    mexPrintf("Convergence not achieved in block %d, at time %d after %d iterations\n",Block_Count,it_,iter);
+                    mexEvalString("st=fclose('all');clear all;");
+                    mexErrMsgTxt("End of simulate");
                   }
-              }
-				   else
-				     {
-                for (it_=y_kmin;it_<periods+y_kmin;it_++)
-                  {
-                    cvg=false;
-                    iter=0;
-                    Per_y_=it_*y_size;
-                    while (!(cvg||(iter>maxit_)))
-                      {
-                        set_code_pointer(begining);
-                        Per_y_=it_*y_size;
-                        compute_block_time(0);
-                        y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0];
-                        double rr;
-                        if(fabs(1+y[Per_y_+Block_Contain[0].Variable])>eps)
-				    					    rr=r[0]/(1+y[Per_y_+Block_Contain[0].Variable]);
-						    		    else
-								    	    rr=r[0];
-                        cvg=((rr*rr)<solve_tolf);
-                        iter++;
-                      }
-                    if (!cvg)
-                      {
-                        mexPrintf("Convergence not achieved in block %d, at time %d after %d iterations\n",Block_Count,it_,iter);
-                        mexEvalString("st=fclose('all');clear all;");
-                        mexErrMsgTxt("End of simulate");
-                     }
-								}
-					 }
+							}
+          }
         mxFree(g1);
         mxFree(r);
         break;
@@ -756,7 +934,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
               {
                 set_code_pointer(begining);
                 Per_y_=it_*y_size;
-                    compute_block_time(0);
+                    compute_block_time(0, false);
                     y[Block_Contain[0].Variable] += -r[0]/g1[0];
                     double rr;
     					    	rr=r[0];
@@ -782,7 +960,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
                       {
                         set_code_pointer(begining);
                         Per_y_=it_*y_size;
-                        compute_block_time(0);
+                        compute_block_time(0, false);
                         y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0];
                         double rr;
                         if(fabs(1+y[Per_y_+Block_Contain[0].Variable])>eps)
@@ -833,7 +1011,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
                     res2=0;
 						    		res1=0;
                     max_res=0;
-                    compute_block_time(0);
+                    compute_block_time(0, false);
                     /*if (isnan(res1)||isinf(res1))
                     {
                     memcpy(y, y_save, y_size*sizeof(double)*(periods+y_kmax+y_kmin));
@@ -875,7 +1053,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
                 iter = 0;
                 res1=res2=max_res=0;max_res_idx=0;
                 error_not_printed = true;
-                compute_block_time(0);
+                compute_block_time(0, false);
                 cvg=false;
                 result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true);
               }
@@ -899,7 +1077,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
                         res2=0;
 						    				res1=0;
                         max_res=0;
-                        compute_block_time(0);
+                        compute_block_time(0, false);
                         /*if (isnan(res1)||isinf(res1))
                           {
                             memcpy(y, y_save, y_size*sizeof(double)*(periods+y_kmax+y_kmin));
@@ -945,7 +1123,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
                     iter = 0;
                     res1=res2=max_res=0;max_res_idx=0;
                     error_not_printed = true;
-                    compute_block_time(0);
+                    compute_block_time(0, false);
                     cvg=false;
                     result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false);
                   }
@@ -987,7 +1165,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
                     res2=0;
 						    		res1=0;
                     max_res=0;
-                    compute_block_time(0);
+                    compute_block_time(0, false);
                     /*if (isnan(res1)||isinf(res1))
                     {
                     memcpy(y, y_save, y_size*sizeof(double)*(periods+y_kmax+y_kmin));
@@ -1028,7 +1206,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
                 iter = 0;
                 res1=res2=max_res=0;max_res_idx=0;
                 error_not_printed = true;
-                compute_block_time(0);
+                compute_block_time(0, false);
                 cvg=false;
                 result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true);
               }
@@ -1050,7 +1228,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
                         res2=0;
                         res1=0;
                         max_res=0;
-                        compute_block_time(0);
+                        compute_block_time(0, false);
                         /*if (isnan(res1)||isinf(res1))
                           {
                             memcpy(y, y_save, y_size*sizeof(double)*(periods+y_kmax+y_kmin));
@@ -1094,7 +1272,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
                     set_code_pointer(begining);
                     Per_y_=it_*y_size;
                     error_not_printed = true;
-                    compute_block_time(0);
+                    compute_block_time(0, false);
                     cvg=false;
                     result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false);
                   }
@@ -1147,7 +1325,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
                     Per_u_=(it_-y_kmin)*u_count_int;
                     Per_y_=it_*y_size;
                     set_code_pointer(begining);
-                    compute_block_time(Per_u_);
+                    compute_block_time(Per_u_, false);
                     if (isnan(res1)||isinf(res1))
                       {
                         memcpy(y, y_save, y_size*sizeof(double)*(periods+y_kmax+y_kmin));
@@ -1192,7 +1370,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
                 Per_u_=(it_-y_kmin)*u_count_int;
                 Per_y_=it_*y_size;
                 set_code_pointer(begining);
-                compute_block_time(Per_u_);
+                compute_block_time(Per_u_, false);
                 for (i=0; i< size; i++)
                   {
                     double rr;
@@ -1226,7 +1404,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
 }
 
 bool
-Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_state)
+Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate)
 {
   ifstream CompiledCode;
   bool result = true;
@@ -1283,7 +1461,10 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
                 lBlock_Contain.Own_Derivative=get_code_int;
                 Block_Contain.push_back(lBlock_Contain);
               }
-            result = simulate_a_block(lBlock.size,lBlock.type, file_name, bin_basename,true, steady_state, Block_Count);
+            if(evaluate)
+              evaluate_a_block(lBlock.size,lBlock.type, bin_basename, steady_state, Block_Count);
+            else
+              result = simulate_a_block(lBlock.size,lBlock.type, file_name, bin_basename,true, steady_state, Block_Count);
             if(!result)
               go_on = false;
             break;
diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh
index 10f9256d7..aa4d59bd6 100644
--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
@@ -68,8 +68,9 @@ class Interpreter : SparseMatrix
   protected :
     double pow1(double a, double b);
     double log1(double a);
-    void compute_block_time(int Per_u_);
-    bool simulate_a_block(int size,int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num);
+    void compute_block_time(int Per_u_, bool evaluate);
+    void evaluate_a_block(int size, int type, string bin_basename, bool steady_state, int block_num);
+    bool simulate_a_block(int size, int type, string file_name, string bin_basename, bool Gaussian_Elimination, bool steady_state, int block_num);
     double *T;
     vector<Block_contain_type> Block_Contain;
     vector<Block_type> Block;
@@ -92,7 +93,7 @@ class Interpreter : SparseMatrix
     Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *direction_arg, int y_size_arg, int nb_row_x_arg,
                 int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, int maxit_arg_, double solve_tolf_arg, int size_o_direction_arg,
                 double slowc_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg);
-    bool compute_blocks(string file_name, string bin_basename, bool steady_state);
+    bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate);
 };
 
 
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index 84e211e37..5e83d8d62 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -46,12 +46,21 @@ max(int a, int b)
 using namespace std;
 #include <sstream>
 #include "mex_interface.hh"
+string
+Get_Argument(const char *argv)
+{
+  //mexPrintf("number=%d\n",number);
+  string f(argv);
+  return f;
+}
+
 
 int
 main( int argc, const char* argv[] )
 {
   FILE *fid;
   bool steady_state = false;
+  bool evaluate = false;
   printf("argc=%d\n",argc);
   if(argc<2)
     {
@@ -70,11 +79,23 @@ main( int argc, const char* argv[] )
 
   string file_name(argv[1]);
 
-  if(argc>2)
+  for(i=2;i<argc; i++)
     {
-      string f(argv[1]);
-      if(f == "steady_state")
+      if(Get_Argument(argv[i])=="steady_state")
         steady_state = true;
+      else if(Get_Argument(argv[i])=="dynamic")
+        steady_state = false;
+      else if(Get_Argument(argv[i])=="evaluate")
+        evaluate = true;
+      else
+        {
+          mexPrintf("Unknown argument : ");
+          mexEvalString("st=fclose('all');clear all;");
+          string f;
+          f = Get_Argument(argv[i]);
+          f.append("\n");
+          mexErrMsgTxt(f.c_str());
+        }
     }
   fid = fopen(tmp_out.str().c_str(),"r");
   int periods;
@@ -147,9 +168,10 @@ main( int argc, const char* argv[] )
   clock_t t0= clock();
   Interpreter interprete(params, y, ya, x, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name);
   string f(file_name);
-  interprete.compute_blocks(f, f, steady_state);
+  interprete.compute_blocks(f, f, steady_state, evaluate);
   clock_t t1= clock();
-  mexPrintf("Simulation Time=%f milliseconds\n",1000.0*(double(t1)-double(t0))/double(CLOCKS_PER_SEC));
+  if(!evaluate)
+    mexPrintf("Simulation Time=%f milliseconds\n",1000.0*(double(t1)-double(t0))/double(CLOCKS_PER_SEC));
   if(x)
     mxFree(x);
   if(y)
@@ -162,28 +184,51 @@ main( int argc, const char* argv[] )
 }
 
 #else
+
+string
+Get_Argument(const mxArray *prhs)
+{
+  //mexPrintf("number=%d\n",number);
+  const mxArray *mxa = prhs;
+  int buflen=mxGetM(mxa) * mxGetN(mxa) + 1;
+  char *first_argument;
+  first_argument=(char*)mxCalloc(buflen, sizeof(char));
+  int status = mxGetString(mxa, first_argument, buflen);
+  if (status != 0)
+    mexWarnMsgTxt("Not enough space. The first argument is truncated.");
+  string f(first_argument);
+  mxFree(first_argument);
+  return f;
+}
 /* The gateway routine */
 void
 mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 {
   mxArray *M_, *oo_, *options_;
   int i, row_y, col_y, row_x, col_x, nb_row_xd;
+  int steady_row_y, steady_col_y, steady_row_x, steady_col_x, steady_nb_row_xd;
   int y_kmin=0, y_kmax=0, y_decal=0, periods=1;
   double * pind ;
   double *direction;
   bool steady_state = false;
-  if(nrhs>0)
+  bool evaluate = false;
+  for(i=0;i<nrhs; i++)
     {
-    	const mxArray *mxa = prhs[0];
-      int buflen=mxGetM(mxa) * mxGetN(mxa) + 1;
-      char *first_argument;
-      first_argument=(char*)mxCalloc(buflen, sizeof(char));
-      int status = mxGetString(mxa, first_argument, buflen);
-      if (status != 0)
-        mexWarnMsgTxt("Not enough space. The first argument is truncated.");
-      string f(first_argument);
-      if(f == "steady_state")
+      if(Get_Argument(prhs[i])=="steady_state")
         steady_state = true;
+      else if(Get_Argument(prhs[i])=="dynamic")
+        steady_state = false;
+      else if(Get_Argument(prhs[i])=="evaluate")
+        evaluate = true;
+      else
+        {
+          mexPrintf("Unknown argument : ");
+          mexEvalString("st=fclose('all');clear all;");
+          string f;
+          f = Get_Argument(prhs[i]);
+          f.append("\n");
+          mexErrMsgTxt(f.c_str());
+        }
     }
   M_ = mexGetVariable("global","M_");
   if (M_ == NULL )
@@ -209,7 +254,7 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     }
   //mexPrintf("ok0\n");
   double * params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"params")));
-  double *yd, *xd;
+  double *yd, *xd, *steady_yd, *steady_xd ;
   if(!steady_state)
     {
       yd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"endo_simul")));
@@ -224,6 +269,14 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
       y_kmax=int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"maximum_lead"))))));
       y_decal=max(0,y_kmin-int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"maximum_endo_lag")))))));
       periods=int(floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_,"periods"))))));
+
+      steady_yd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"steady_state")));
+      steady_row_y=mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"steady_state")));
+      steady_col_y=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"steady_state")));;
+      steady_xd= mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_steady_state")));
+      steady_row_x=mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_steady_state")));
+      steady_col_x=mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_,"exo_steady_state")));
+      steady_nb_row_xd=int(floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,"exo_det_nbr"))))));
     }
 	else
 	  {
@@ -262,27 +315,33 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   for (i=0;i<row_x*col_x;i++)
      x[i]=double(xd[i]);
   for (i=0;i<row_y*col_y;i++)
-    y[i]=double(yd[i]);
+    {
+      y[i]  = double(yd[i]);
+      ya[i] = double(yd[i]);
+    }
   int y_size=row_y;
   int nb_row_x=row_x;
 
   /*int it_ = y_kmin;
   for (int j = 0; j < y_size; j++)
 		mexPrintf("   variable %d at time %d and %d = %f\n", j+1, it_, it_+1, y[j+it_*y_size]);*/
-
   clock_t t0= clock();
   Interpreter interprete(params, y, ya, x, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name);
   string f(fname);
-  bool result = interprete.compute_blocks(f, f, steady_state);
+  bool result = interprete.compute_blocks(f, f, steady_state, evaluate);
   clock_t t1= clock();
-  if(!steady_state)
+  if(!steady_state && !evaluate)
     mexPrintf("Simulation Time=%f milliseconds\n",1000.0*(double(t1)-double(t0))/double(CLOCKS_PER_SEC));
   if (nlhs>0)
     {
       plhs[0] = mxCreateDoubleMatrix(row_y, col_y, mxREAL);
       pind = mxGetPr(plhs[0]);
-      for (i=0;i<row_y*col_y;i++)
-        pind[i]=y[i];
+      if(evaluate)
+        for (i=0;i<row_y*col_y;i++)
+          pind[i]=y[i]-ya[i];
+      else
+        for (i=0;i<row_y*col_y;i++)
+          pind[i]=y[i];
 			if(nlhs>1)
 			  {
 			    plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);