v4 matlab+preprocessor: allow syntax of homotopy_setup without initial value

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1779 ac1d8469-bf42-47a9-8791-bf33cf982152

matlab/homotopy1.m

@@ -13,6 +13,8 @@ function homotopy1(values, step_nbr)
 %                   exogenous deterministic, 4 for parameters)
 %                   Column 2 is symbol integer identifier.
 %                   Column 3 is initial value, and column 4 is final value.
+%                   Column 3 can contain NaNs, in which case previous
+%                   initialization of variable will be used as initial value.
 %    step_nbr:      number of steps for homotopy
 %
 % OUTPUTS
@@ -36,13 +38,23 @@ function homotopy1(values, step_nbr)
     error('HOMOTOPY: incorrect variable types specified')
   end
 
-  if any(values(:,3) == values(:,4))
+  % Construct vector of starting values, using previously initialized values
+  % when initial value has not been given in homotopy_setup block
+  oldvalues = values(:,3);
+  ipn = find(values(:,1) == 4 & isnan(oldvalues));
+  oldvalues(ipn) = M_.params(values(ipn, 2));
+  ixn = find(values(:,1) == 1 & isnan(oldvalues));
+  oldvalues(ixn) = oo_.exo_steady_state(values(ixn, 2));
+  ixdn = find(values(:,1) == 2 & isnan(oldvalues));
+  oldvalues(ixdn) = oo_.exo_det_steady_state(values(ixdn, 2));
+  
+  if any(oldvalues == values(:,4))
     error('HOMOTOPY: initial and final values should be different')
   end
   
   points = zeros(nv, step_nbr+1);
   for i = 1:nv
-    points(i,:) = values(i,3):(values(i,4)-values(i,3))/step_nbr:values(i,4);
+    points(i,:) = oldvalues(i):(values(i,4)-oldvalues(i))/step_nbr:values(i,4);
   end
   
   for i=1:step_nbr+1

matlab/homotopy2.m

@@ -15,6 +15,8 @@ function homotopy2(values, step_nbr)
 %                   exogenous deterministic, 4 for parameters)
 %                   Column 2 is symbol integer identifier.
 %                   Column 3 is initial value, and column 4 is final value.
+%                   Column 3 can contain NaNs, in which case previous
+%                   initialization of variable will be used as initial value.
 %    step_nbr:      number of steps for homotopy
 %
 % OUTPUTS
@@ -30,27 +32,41 @@ function homotopy2(values, step_nbr)
 
   nv = size(values, 1);
   
-  % Initialize all variables with initial value
+  oldvalues = values(:,3);
+  
+  % Initialize all variables with initial value, or the reverse...
   for i = 1:nv
     switch values(i,1)
      case 1
-      oo_.exo_steady_state(values(i,2)) = values(i,3);
+      if isnan(oldvalues(i))
+        oldvalues(i) = oo_.exo_steady_state(values(i,2));
+      else
+        oo_.exo_steady_state(values(i,2)) = oldvalues(i);
+      end
      case 2
-      oo_.exo_det_steady_state(values(i,2)) = values(i,3);
+      if isnan(oldvalues(i))
+        oldvalues(i) = oo_.exo_det_steady_state(values(i,2));
+      else
+        oo_.exo_det_steady_state(values(i,2)) = oldvalues(i);
+      end
      case 4
-      M_.params(values(i,2)) = values(i,3);
+      if isnan(oldvalues(i))
+        oldvalues(i) = M_.params(values(i,2));
+      else
+        M_.params(values(i,2)) = oldvalues(i);
+      end
      otherwise
       error('HOMOTOPY: incorrect variable types specified')
     end
   end
 
-  if any(values(:,3) == values(:,4))
+  if any(oldvalues == values(:,4))
     error('HOMOTOPY: initial and final values should be different')
   end
   
   % Actually do the homotopy
   for i = 1:nv
-    for v = values(i,3):(values(i,4)-values(i,3))/step_nbr:values(i,4)
+    for v = oldvalues(i):(values(i,4)-oldvalues(i))/step_nbr:values(i,4)
       switch values(i,1)
        case 1
         oo_.exo_steady_state(values(i,2)) = v;

matlab/homotopy3.m

@@ -15,6 +15,8 @@ function homotopy3(values, step_nbr)
 %                   exogenous deterministic, 4 for parameters)
 %                   Column 2 is symbol integer identifier.
 %                   Column 3 is initial value, and column 4 is final value.
+%                   Column 3 can contain NaNs, in which case previous
+%                   initialization of variable will be used as initial value.
 %    step_nbr:      maximum number of steps to try before aborting
 %
 % OUTPUTS
@@ -39,7 +41,16 @@ function homotopy3(values, step_nbr)
     error('HOMOTOPY: incorrect variable types specified')
   end
 
+  % Construct vector of starting values, using previously initialized values
+  % when initial value has not been given in homotopy_setup block
   oldvalues = values(:,3);
+  ipn = find(values(:,1) == 4 & isnan(oldvalues));
+  oldvalues(ipn) = M_.params(values(ipn, 2));
+  ixn = find(values(:,1) == 1 & isnan(oldvalues));
+  oldvalues(ixn) = oo_.exo_steady_state(values(ixn, 2));
+  ixdn = find(values(:,1) == 2 & isnan(oldvalues));
+  oldvalues(ixdn) = oo_.exo_det_steady_state(values(ixdn, 2));
+
   targetvalues = values(:,4);
 
   if min(abs(targetvalues - oldvalues)) < tol

preprocessor/DynareBison.yy

@@ -1250,7 +1250,10 @@ homotopy_list : homotopy_item
               ;
 
 homotopy_item : NAME COMMA expression COMMA expression ';'
-              { driver.homotopy_val($1,$3,$5);};
+                { driver.homotopy_val($1, $3, $5);}
+              | NAME COMMA expression ';'
+                { driver.homotopy_val($1, NULL, $3);}
+              ;
 
 number : INT_NUMBER
        | FLOAT_NUMBER

preprocessor/NumericalInitialization.cc

@@ -179,7 +179,10 @@ HomotopyStatement::writeOutput(ostream &output, const string &basename) const
       const int id = symbol_table.getID(name) + 1;
 
       output << "options_.homotopy_values = vertcat(options_.homotopy_values, [ " << type << ", " << id << ", ";
+      if (expression1 != NULL)
         expression1->writeOutput(output);
+      else
+        output << "NaN";
       output << ", ";
       expression2->writeOutput(output);
       output << "]);" << endl;

preprocessor/ParsingDriver.cc

@@ -347,10 +347,7 @@ ParsingDriver::homotopy_val(string *name, NodeID val1, NodeID val2)
       && type != eExogenousDet)
     error("homotopy_val: " + *name + " should be a parameter or exogenous variable");
 
-  if (homotopy_values.find(*name) != homotopy_values.end())
-    error("homotopy_val: " + *name +" declared twice");
-
-  homotopy_values[*name] = make_pair(val1, val2);
+  homotopy_values.push_back(make_pair(*name, make_pair(val1, val2)));
 
   delete name;
 }

preprocessor/include/NumericalInitialization.hh

@@ -95,11 +95,9 @@ public:
 class HomotopyStatement : public Statement
 {
 public:
-  /*!
-    Contrary to Initval and Endval, we use a map, since it is impossible to reuse
-    a given initialization value in a second initialization inside the block.
-  */
-  typedef map<string,pair<NodeID,NodeID> > homotopy_values_type;
+  //! Stores the declarations of homotopy_setup
+  /*! Order matter so we use a vector. First NodeID can be NULL if no initial value given. */
+  typedef vector<pair<string, pair<NodeID, NodeID> > > homotopy_values_type;
 private:
   const homotopy_values_type homotopy_values;
   const SymbolTable &symbol_table;

preprocessor/include/ParsingDriver.hh

@@ -209,7 +209,8 @@ public:
   void init_val(string *name, NodeID rhs);
   //! Writes an histval block
   void hist_val(string *name, string *lag, NodeID rhs);
-  //! Writes an homotopy_setup block
+  //! Adds an entry in a homotopy_setup block
+  /*! Second argument "val1" can be NULL if no initial value provided */
   void homotopy_val(string *name, NodeID val1, NodeID val2);
   //! Writes end of an initval block
   void end_initval();

tests/homotopy/ramst_homotopy.mod

@@ -7,13 +7,14 @@
 initval;
 k = 12.75;
 c = 1.5;
+x = 1;
 end;
 
 homotopy_setup;
 bet, 0.05, 0.1;
-x, 1, 2;
+x, 2;
 end;
 
-steady(homotopy_mode = 1, homotopy_steps = 50);
+steady(homotopy_mode = 2, homotopy_steps = 50);
 //steady(homotopy_mode = 2, homotopy_steps = 50);
 //steady(homotopy_mode = 3, homotopy_steps = 50);