Commit 75fb140b authored by sebastien's avatar sebastien
Browse files

v4 matlab+preprocessor:

* removed global variables from "static" and "dynamic" files (in standard and USE_DLL modes only)
* added extra arguments when those functions are called from M-files


git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1859 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 21c98f8b
......@@ -79,7 +79,7 @@ if options_.ramsey_policy
k1 = find(M_.orig_model.lead_lag_incidence');
y = repmat(oo_.dr.ys(1:M_.orig_model.endo_nbr),1,M_.orig_model.maximum_lag+M_.orig_model.maximum_lead+1);
[f,fJ,fh] = feval([M_.fname '_dynamic'],y(k1),zeros(1,M_.exo_nbr));
[f,fJ,fh] = feval([M_.fname '_dynamic'],y(k1),zeros(1,M_.exo_nbr), M_.params, it_);
% looking for dynamic variables that are both in the original model
% and in the optimal policy model
......@@ -94,7 +94,7 @@ if options_.ramsey_policy
lead_lag_incidence1 = M_.orig_model.lead_lag_incidence;
y = repmat(oo_.dr.ys(1:M_.orig_model.endo_nbr),1,M_.orig_model.maximum_lag+M_.orig_model.maximum_lead+1);
k1 = find(M_.orig_model.lead_lag_incidence');
[f,fj,fh] = feval([M_.fname '_dynamic'],y(k1),zeros(1,M_.exo_nbr));
[f,fj,fh] = feval([M_.fname '_dynamic'],y(k1),zeros(1,M_.exo_nbr), M_.params, it_);
nrj = size(fj,1);
iy = M_.lead_lag_incidence;
......@@ -174,11 +174,11 @@ else
z = z(iyr0) ;
if options_.order == 1
[junk,jacobia_] = feval([M_.fname '_dynamic'],z,[oo_.exo_simul ...
oo_.exo_det_simul]);
oo_.exo_det_simul], M_.params, it_);
elseif options_.order == 2
[junk,jacobia_,hessian] = feval([M_.fname '_dynamic'],z,...
[oo_.exo_simul ...
oo_.exo_det_simul]);
oo_.exo_det_simul], M_.params, it_);
end
end
......
......@@ -20,7 +20,7 @@ global it_ M_ oo_
n1 = size(x,1) - M_.exo_nbr;
oo_.exo_simul(it_+M_.maximum_lag-M_.maximum_lag,:) = x(n1+1:end)';
fh = str2func([M_.fname '_static']);
y=feval(fh,x(1:n1),oo_.exo_simul);
y=feval(fh,x(1:n1),oo_.exo_simul, M_.params);
......@@ -66,7 +66,7 @@ function homotopy1(values, step_nbr)
oo_.steady_state,...
options_.jacobian_flag, ...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]);
oo_.exo_det_steady_state], M_.params);
if check
error('HOMOTOPY didn''t succeed')
......
......@@ -80,7 +80,7 @@ function homotopy2(values, step_nbr)
oo_.steady_state,...
options_.jacobian_flag, ...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]);
oo_.exo_det_steady_state], M_.params);
if check
error('HOMOTOPY didn''t succeed')
......
......@@ -75,7 +75,7 @@ function homotopy3(values, step_nbr)
oo_.steady_state,...
options_.jacobian_flag, ...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]);
oo_.exo_det_steady_state], M_.params);
if check
inc = inc/2;
......
......@@ -60,7 +60,7 @@ function initial_estimation_checks(xparam1,gend,data)
check1 = max(abs(feval([M_.fname '_static'],...
oo_.steady_state,...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]))) > options_.dynatol ;
oo_.exo_det_steady_state], M_.params))) > options_.dynatol ;
if check1
error(['The seadystate values returned by ' M_.fname ...
'_steadystate.m don''t solve the static model!' ])
......
......@@ -26,8 +26,8 @@ function osr1(i_params,i_var,weights)
% check if ys is steady state
fh = str2func([M_.fname '_static']);
if max(abs(feval(fh,oo_.steady_state,exe))) > options_.dynatol
[oo_.dr.ys, check] = dynare_solve([M_.fname '_static'],oo_.steady_state,1,exe);
if max(abs(feval(fh, oo_.steady_state, exe, M_.params))) > options_.dynatol
[oo_.dr.ys, check] = dynare_solve([M_.fname '_static'], oo_.steady_state, 1, exe, M_.params);
if check
error('OLR: convergence problem in DYNARE_SOLVE')
end
......
......@@ -54,12 +54,12 @@ function J = ramsey_dynamic(ys,lbar)
k = find(i_leadlag');
it_ = 1;
% retrieving derivatives of the objective function
[U,Uy,Uyy] = feval([fname '_objective_static'],ys,zeros(1,exo_nbr));
[U,Uy,Uyy] = feval([fname '_objective_static'],ys,zeros(1,exo_nbr), M_.params);
Uy = Uy';
Uyy = reshape(Uyy,endo_nbr,endo_nbr);
% retrieving derivatives of original model
[f,fJ,fH] = feval([fname '_dynamic'],y(k),zeros(1,exo_nbr));
[f,fJ,fH] = feval([fname '_dynamic'],y(k),zeros(1,exo_nbr), M_.params, it_);
instr_nbr = endo_nbr - size(f,1);
mult_nbr = endo_nbr-instr_nbr;
......
......@@ -43,7 +43,7 @@ function info = ramsey_policy(var_list)
ghs22 = [ghs2(nstatic+(1:npred),:); zeros(nspred-npred,size(ghs2,2))];
fname = [M_.fname '_objective_static'];
[ubar,uj,uh] = feval(fname,dr.ys(1:endo_nbr1),zeros(exo_nbr1,1));
[ubar,uj,uh] = feval(fname,dr.ys(1:endo_nbr1),zeros(exo_nbr1,1), M_.params);
bet = options_.planner_discount;
wbar = ubar/(1-bet);
......
......@@ -38,7 +38,7 @@ function [resids, rJ,mult] = ramsey_static(x)
% value and Jacobian of objective function
ex = zeros(1,M_.exo_nbr);
[U,Uy,Uyy] = feval([fname '_objective_static'],x(i_endo),ex);
[U,Uy,Uyy] = feval([fname '_objective_static'],x(i_endo),ex, M_.params);
Uy = Uy';
Uyy = reshape(Uyy,endo_nbr,endo_nbr);
......@@ -47,7 +47,7 @@ function [resids, rJ,mult] = ramsey_static(x)
k = find(i_lag');
it_ = 1;
% [f,fJ,fH] = feval([fname '_dynamic'],y(k),ex);
[f,fJ] = feval([fname '_dynamic'],y(k),ex);
[f,fJ] = feval([fname '_dynamic'],y(k),ex, M_.params, it_);
% indices of Lagrange multipliers
inst_nbr = endo_nbr - size(f,1);
i_mult = [endo_nbr+1:2*endo_nbr-inst_nbr]';
......
......@@ -39,7 +39,7 @@ function resid(period)
z = zeros(n,period);
fh = str2func([M_.fname '_dynamic']);
for it_=M_.maximum_lag+1:period+M_.maximum_lag
z(:,it_-M_.maximum_lag) = feval(fh,y(iyr0),oo_.exo_simul);
z(:,it_-M_.maximum_lag) = feval(fh,y(iyr0),oo_.exo_simul, M_.params, it_);
iyr0 = iyr0 + n;
end
......
......@@ -60,17 +60,17 @@ if options_.steadystate_flag
[oo_.exo_steady_state; oo_.exo_det_steady_state]);
else
% testing if ys isn't a steady state or if we aren't computing Ramsey policy
if max(abs(feval(fh,dr.ys,[oo_.exo_steady_state; oo_.exo_det_steady_state]))) ...
if max(abs(feval(fh,dr.ys,[oo_.exo_steady_state; oo_.exo_det_steady_state], M_.params))) ...
> options_.dynatol & options_.ramsey_policy == 0
if options_.linear == 0
% nonlinear models
[dr.ys,check1] = dynare_solve(fh,dr.ys,options_.jacobian_flag,...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]);
oo_.exo_det_steady_state], M_.params);
else
% linear models
[fvec,jacob] = feval(fh,dr.ys,[oo_.exo_steady_state;...
oo_.exo_det_steady_state]);
oo_.exo_det_steady_state], M_.params);
dr.ys = dr.ys-jacob\fvec;
end
end
......@@ -78,7 +78,7 @@ end
% testing for problem
if check1
info(1) = 20;
resid = feval(fh,ys,oo_.exo_steady_state);
resid = feval(fh,ys,oo_.exo_steady_state, M_.params);
info(2) = resid'*resid; % penalty...
return
end
......@@ -92,7 +92,7 @@ end
if options_.dr_algo == 1 & options_.order > 1
dr.ys = dynare_solve('dr2',ys,0,dr);
dr.fbias = 2*feval([M_.fname '_static'],dr.ys,oo_.exo_steady_state);
dr.fbias = 2*feval([M_.fname '_static'],dr.ys,oo_.exo_steady_state, M_.params);
[dr, info1] = dr1(dr,check_flag);
if info1(1)
info(1) = info(1)+10;
......
......@@ -8,7 +8,7 @@ function [sR,sG] = restricted_steadystate(y,x,indx)
ss(indx) = y;
eval(['[R,G] = ' M_.fname '_static(ss,x);']);
eval(['[R,G] = ' M_.fname '_static(ss, x, M_.params);']);
sR = R(inde);
sG = G(inde,indx);
\ No newline at end of file
......@@ -57,14 +57,14 @@ for iter = 1:options_.maxit
it_ = it_init ;
z = [oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)] ;
[d1,M_.jacobia] = feval([M_.fname '_dynamic'],z,oo_.exo_simul);
[d1,M_.jacobia] = feval([M_.fname '_dynamic'],z,oo_.exo_simul, M_.params, it_);
M_.jacobia = [M_.jacobia(:,iz) -d1] ;
ic = [1:ny] ;
icp = iyp ;
c (ic,:) = M_.jacobia(:,is)\M_.jacobia(:,isf1) ;
for it_ = it_init+(1:options_.periods-1)
z = [oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)] ;
[d1,M_.jacobia] = feval([M_.fname '_dynamic'],z,oo_.exo_simul);
[d1,M_.jacobia] = feval([M_.fname '_dynamic'],z,oo_.exo_simul, M_.params, it_);
M_.jacobia = [M_.jacobia(:,iz) -d1] ;
M_.jacobia(:,[isf nrs]) = M_.jacobia(:,[isf nrs])-M_.jacobia(:,isp)*c(icp,:) ;
ic = ic + ny ;
......
......@@ -22,7 +22,7 @@ function simk
global M_ options_ oo_
global it_ iyr0 ct_ broyden_
func_name = [M_.fname '_static'];
%func_name = [M_.fname '_static'];
nk = M_.maximum_endo_lag + M_.maximum_endo_lead + 1 ;
ny = size(M_.lead_lag_incidence,2) ;
icc1 = M_.lead_lag_incidence(nk,:) > 0;
......@@ -126,10 +126,10 @@ for iter = 1:options_.maxit
h3 = clock ;
if broyden_ & iter > 1
%d1_ = -feval(fh,oo_.endo_simul(iyr));
d1 = -feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),z,oo_.exo_simul);
d1 = -feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),z,oo_.exo_simul, M_.params, it_);
else
%jacob(func_name,oo_.endo_simul(iyr)) ;
[d1,M_.jacobia] = feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),oo_.exo_simul);
[d1,M_.jacobia] = feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),oo_.exo_simul, M_.params, it_);
d1 = -d1 ;
end
err_f = max(err_f,max(abs(d1)));
......@@ -219,10 +219,10 @@ for iter = 1:options_.maxit
while it_ <= options_.periods+M_.maximum_lag
if broyden_
%d1_ = -feval(fh,oo_.endo_simul(iyr));
d1 = -feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),z,oo_.exo_simul);
d1 = -feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),z,oo_.exo_simul, M_.params, it_);
else
%jacob(func_name,oo_.endo_simul(iyr)) ;
[d1,M_.jacobia] = feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),oo_.exo_simul);
[d1,M_.jacobia] = feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),oo_.exo_simul, M_.params, it_);
d1 = -d1 ;
end
err_f = max(err_f,max(abs(d1)));
......
......@@ -32,7 +32,7 @@ function steady_()
check1 = max(abs(feval([M_.fname '_static'],...
oo_.steady_state,...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]))) > options_.dynatol ;
oo_.exo_det_steady_state], M_.params))) > options_.dynatol ;
if check1
error(['The seadystate values returned by ' M_.fname ...
'_steadystate.m don''t solve the static model!' ])
......@@ -56,7 +56,7 @@ function steady_()
oo_.steady_state,...
options_.jacobian_flag, ...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]);
oo_.exo_det_steady_state], M_.params);
end
if check ~= 0
......
......@@ -1530,7 +1530,7 @@ ModelTree::writeStaticMFile(const string &static_basename) const
exit(-1);
}
// Writing comments and function definition command
mStaticModelFile << "function [residual, g1, g2] = " << static_basename << "( y, x )\n";
mStaticModelFile << "function [residual, g1, g2] = " << static_basename << "(y, x, params)" << endl;
mStaticModelFile << interfaces::comment()+"\n"+interfaces::comment();
mStaticModelFile << "Status : Computes static model for Dynare\n" << interfaces::comment() << "\n";
mStaticModelFile << interfaces::comment();
......@@ -1557,7 +1557,7 @@ ModelTree::writeDynamicMFile(const string &dynamic_basename) const
cerr << "Error: Can't open file " << filename << " for writing" << endl;
exit(-1);
}
mDynamicModelFile << "function [residual, g1, g2, g3] = " << dynamic_basename << "(y, x)\n";
mDynamicModelFile << "function [residual, g1, g2, g3] = " << dynamic_basename << "(y, x, params, it_)" << endl;
mDynamicModelFile << interfaces::comment()+"\n"+interfaces::comment();
mDynamicModelFile << "Status : Computes dynamic model for Dynare\n" << interfaces::comment() << "\n";
mDynamicModelFile << interfaces::comment();
......@@ -1590,9 +1590,7 @@ ModelTree::writeStaticCFile(const string &static_basename) const
<< endl
<< " */" << endl
<< "#include <math.h>" << endl
<< "#include \"mex.h\"" << endl
// A global variable for model parameters
<< "double *params;" << endl;
<< "#include \"mex.h\"" << endl;
// Writing the function Static
writeStaticModel(mStaticModelFile);
......@@ -1601,9 +1599,8 @@ ModelTree::writeStaticCFile(const string &static_basename) const
mStaticModelFile << "/* The gateway routine */" << endl
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
<< "{" << endl
<< " double *y, *x;" << endl
<< " double *y, *x, *params;" << endl
<< " double *residual, *g1;" << endl
<< " mxArray *M_;" << endl
<< endl
<< " /* Create a pointer to the input matrix y. */" << endl
<< " y = mxGetPr(prhs[0]);" << endl
......@@ -1611,6 +1608,9 @@ ModelTree::writeStaticCFile(const string &static_basename) const
<< " /* Create a pointer to the input matrix x. */" << endl
<< " x = mxGetPr(prhs[1]);" << endl
<< endl
<< " /* Create a pointer to the input matrix params. */" << endl
<< " params = mxGetPr(prhs[2]);" << endl
<< endl
<< " residual = NULL;" << endl
<< " if (nlhs >= 1)" << endl
<< " {" << endl
......@@ -1629,15 +1629,8 @@ ModelTree::writeStaticCFile(const string &static_basename) const
<< " g1 = mxGetPr(plhs[1]);" << endl
<< " }" << endl
<< endl
<< " /* Gets model parameters from global workspace of Matlab */" << endl
<< " M_ = mexGetVariable(\"global\",\"M_\");" << endl
<< " if (M_ == NULL ){" << endl
<< " mexPrintf(\"Global variable not found : \");" << endl
<< " mexErrMsgTxt(\"M_ \\n\");" << endl
<< " }" << endl
<< " params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"params\")));" << endl
<< " /* Call the C Static. */" << endl
<< " Static(y, x, residual, g1);" << endl
<< " Static(y, x, params, residual, g1);" << endl
<< "}" << endl;
mStaticModelFile.close();
......@@ -1663,12 +1656,7 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename) const
<< endl
<< " */" << endl
<< "#include <math.h>" << endl
<< "#include \"mex.h\"" << endl
// A global variable for model parameters
<< "double *params;" << endl
// A global variable for it_
<< "int it_;" << endl
<< "int nb_row_x;" << endl;
<< "#include \"mex.h\"" << endl;
// Writing the function body
writeDynamicModel(mDynamicModelFile);
......@@ -1677,15 +1665,22 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename) const
mDynamicModelFile << "/* The gateway routine */" << endl
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
<< "{" << endl
<< " double *y, *x;" << endl
<< " double *y, *x, *params;" << endl
<< " double *residual, *g1, *g2;" << endl
<< " mxArray *M_;" << endl
<< " int nb_row_x, it_;" << endl
<< endl
<< " /* Create a pointer to the input matrix y. */" << endl
<< " y = mxGetPr(prhs[0]);" << endl
<< endl
<< " /* Create a pointer to the input matrix x. */" << endl
<< " x = mxGetPr(prhs[1]);" << endl
<< endl
<< " /* Create a pointer to the input matrix params. */" << endl
<< " params = mxGetPr(prhs[2]);" << endl
<< endl
<< " /* Fetch time index */" << endl
<< " it_ = (int) mxGetScalar(prhs[3]) - 1;" << endl
<< endl
<< " /* Gets number of rows of matrix x. */" << endl
<< " nb_row_x = mxGetM(prhs[1]);" << endl
<< endl
......@@ -1721,18 +1716,8 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename) const
<< " g2 = mxGetPr(plhs[2]);" << endl
<< " }" << endl
<< endl
<< " /* Gets model parameters from global workspace of Matlab */" << endl
<< " M_ = mexGetVariable(\"global\",\"M_\");" << endl
<< " if (M_ == NULL)" << endl
<< " {" << endl
<< " mexPrintf(\"Global variable not found : \");" << endl
<< " mexErrMsgTxt(\"M_ \\n\");" << endl
<< " }" << endl
<< " params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"params\")));" << endl
<< " /* Gets it_ from global workspace of Matlab */" << endl
<< " it_ = (int) mxGetScalar(mexGetVariable(\"global\", \"it_\"))-1;" << endl
<< " /* Call the C subroutines. */" << endl
<< " Dynamic(y, x, residual, g1, g2);" << endl
<< " Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2);" << endl
<< "}" << endl;
mDynamicModelFile.close();
}
......@@ -1814,9 +1799,6 @@ ModelTree::writeStaticModel(ostream &StaticOutput) const
// Writing ouputs
if (mode != eDLLMode)
{
StaticOutput << "global M_ \n";
StaticOutput << "if M_.param_nbr > 0\n params = M_.params;\nend\n";
StaticOutput << " residual = zeros( " << equations.size() << ", 1);\n";
StaticOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
StaticOutput << "Model equations\n\t";
......@@ -1854,7 +1836,7 @@ ModelTree::writeStaticModel(ostream &StaticOutput) const
}
else
{
StaticOutput << "void Static(double *y, double *x, double *residual, double *g1)" << endl
StaticOutput << "void Static(double *y, double *x, double *params, double *residual, double *g1)" << endl
<< "{" << endl
<< " double lhs, rhs;" << endl
// Writing residual equations
......@@ -3419,10 +3401,8 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
if (mode == eStandardMode)
{
DynamicOutput << "global M_ it_\n";
DynamicOutput << "if M_.param_nbr > 0\n params = M_.params;\nend\n";
DynamicOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
DynamicOutput << "Model equations\n\t";
DynamicOutput << interfaces::comment() << endl << interfaces::comment();
DynamicOutput << "Model equations" << endl;
DynamicOutput << interfaces::comment() + "\n\n";
DynamicOutput << "residual = zeros(" << nrows << ", 1);\n";
......@@ -3466,7 +3446,7 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
}
else
{
DynamicOutput << "void Dynamic(double *y, double *x, double *residual, double *g1, double *g2)" << endl
DynamicOutput << "void Dynamic(double *y, double *x, int nb_row_x, double *params, int it_, double *residual, double *g1, double *g2)" << endl
<< "{" << endl
<< " double lhs, rhs;" << endl
<< endl
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment