Commit 12c4a52e authored by Houtan Bastani's avatar Houtan Bastani
Browse files

Add support for external functions

parent 38b8ccf5
function d=hess_element(func,element1,element2,args)
% function d=hess_element(func,element1,element2,args)
% returns an entry of the finite differences approximation to the hessian of func
%
% INPUTS
% func [function handle] associated with the function
% element1 [int] the indices showing the element within the hessian that should be returned
% element2 [int]
% args [cell array] arguments provided to func
%
% OUTPUTS
% d [double] the (element1,element2) entry of the hessian
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2010 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
h=10e-6;
p10 = args;
p01 = args;
m10 = args;
m01 = args;
p11 = args;
m11 = args;
for i=1:size(args,2)
if i==element1
p10{i} = p10{i} + h;
m10{i} = m10{i} - h;
p11{i} = p11{i} + h;
m11{i} = m11{i} - h;
end
if i==element2
p01{i} = p01{i} + h;
m01{i} = m01{i} - h;
p11{i} = p11{i} + h;
m11{i} = m11{i} - h;
end
end
% From Abramowitz and Stegun. Handbook of Mathematical Functions (1965)
% formulas 25.3.24 and 25.3.27 p. 884
if element1==element2
d = (16*func(p10{:})...
+16*func(m10{:})...
-30*func(args{:})...
-func(p11{:})...
-func(m11{:}))/(12*h^2);
else
d = (func(p10{:})...
+func(m10{:})...
+func(p01{:})...
+func(m01{:})...
-2*func(args{:})...
-func(p11{:})...
-func(m11{:}))/(-2*h^2);
end
function d=jacob_element(func,element,args)
% function d=jacob_element(func,element,args)
% returns an entry of the finite differences approximation to the jacobian of func
%
% INPUTS
% func [function handle] associated with the function
% element [int] the index showing the element within the jacobian that should be returned
% args [cell array] arguments provided to func
%
% OUTPUTS
% d [double] jacobian[element]
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2010 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
h=10e-6;
pargs=args;
margs=args;
for i=1:size(args,2)
if i==element
pargs{i} = pargs{i} + h;
margs{i} = margs{i} - h;
end
end
d=(func(pargs{:})...
-func(margs{:}))/(2*h);
function deriv = subscript_get(nargsout, func, args, varargin)
% function deriv = subscript_get(nargsout, func, args, varargin)
% returns the appropriate entry of the return argument from a user-defined function
% which returns either the jacobian or hessian (or both)
%
% INPUTS
% nargsout [int] integer indicating the number of the return argument containing the jac/hess
% func [function handle] associated with the function
% args [cell array] arguments provided to func
% varargin [cell array] arguments showing the index (or indices) of the element to be returned
%
% OUTPUTS
% deriv [double] the (element1,element2) entry of the hessian
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2010 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
switch size(varargin,2)
case 1 %first deriv
switch nargsout
case 1
[outargs{1}] = func(args{:});
deriv = outargs{1}(varargin{1});
case 2
[outargs{1} outargs{2}] = func(args{:});
deriv = outargs{2}(varargin{1});
otherwise
error('Wrong number of output arguments (%d) passed to subscript_get().',nargsout);
end
case 2 %second deriv
switch nargsout
case 1
[outargs{1}] = func(args{:});
deriv = outargs{1}(varargin{1},varargin{2});
case 3
[outargs{1} outargs{2} outargs{3}] = func(args{:});
deriv = outargs{3}(varargin{1},varargin{2});
otherwise
error('Wrong number of output arguments (%d) passed to subscript_get().',nargsout);
end
otherwise
error('Wrong number of indices (%d) was passed to subscript_get().',size(varargin,2));
end
end
/*
* Copyright (C) 2007-2009 Dynare Team
* Copyright (C) 2007-2010 Dynare Team
*
* This file is part of Dynare.
*
......@@ -139,7 +139,7 @@ enum SymbolType
eParameter = 4, //!< Parameter
eModelLocalVariable = 10, //!< Local variable whose scope is model (pound expression)
eModFileLocalVariable = 11, //!< Local variable whose scope is mod file (model excluded)
eUnknownFunction = 12 //!< Function unknown to the preprocessor
eExternalFunction = 12 //!< External (user-defined) function
};
enum ExpressionType
......
/*
* Copyright (C) 2003-2009 Dynare Team
* Copyright (C) 2003-2010 Dynare Team
*
* This file is part of Dynare.
*
......@@ -23,9 +23,12 @@
#include "DataTree.hh"
DataTree::DataTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg) :
DataTree::DataTree(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg,
ExternalFunctionsTable &external_functions_table_arg) :
symbol_table(symbol_table_arg),
num_constants(num_constants_arg),
external_functions_table(external_functions_table_arg),
node_counter(0)
{
Zero = AddNumConstant("0");
......@@ -454,13 +457,30 @@ DataTree::AddLocalVariable(const string &name, NodeID value) throw (LocalVariabl
}
NodeID
DataTree::AddUnknownFunction(const string &function_name, const vector<NodeID> &arguments)
DataTree::AddExternalFunction(const string &function_name, const vector<NodeID> &arguments)
{
int id = symbol_table.getID(function_name);
return AddExternalFunction(symbol_table.getID(function_name), arguments);
}
NodeID
DataTree::AddExternalFunction(int symb_id, const vector<NodeID> &arguments)
{
assert(symbol_table.getType(symb_id) == eExternalFunction);
return new ExternalFunctionNode(*this, symb_id, arguments);
}
assert(symbol_table.getType(id) == eUnknownFunction);
NodeID
DataTree::AddFirstDerivExternalFunctionNode(int top_level_symb_id, const vector<NodeID> &arguments, int input_index)
{
assert(symbol_table.getType(top_level_symb_id) == eExternalFunction);
return new FirstDerivExternalFunctionNode(*this, top_level_symb_id, arguments, input_index);
}
return new UnknownFunctionNode(*this, id, arguments);
NodeID
DataTree::AddSecondDerivExternalFunctionNode(int top_level_symb_id, const vector<NodeID> &arguments, int input_index1, int input_index2)
{
assert(symbol_table.getType(top_level_symb_id) == eExternalFunction);
return new SecondDerivExternalFunctionNode(*this, top_level_symb_id, arguments, input_index1, input_index2);
}
bool
......
/*
* Copyright (C) 2003-2009 Dynare Team
* Copyright (C) 2003-2010 Dynare Team
*
* This file is part of Dynare.
*
......@@ -30,6 +30,7 @@ using namespace std;
#include "SymbolTable.hh"
#include "NumericalConstants.hh"
#include "ExternalFunctionsTable.hh"
#include "ExprNode.hh"
#define CONSTANTS_PRECISION 16
......@@ -42,12 +43,16 @@ class DataTree
friend class UnaryOpNode;
friend class BinaryOpNode;
friend class TrinaryOpNode;
friend class UnknownFunctionNode;
friend class ExternalFunctionNode;
friend class FirstDerivExternalFunctionNode;
friend class SecondDerivExternalFunctionNode;
protected:
//! A reference to the symbol table
SymbolTable &symbol_table;
//! Reference to numerical constants table
NumericalConstants &num_constants;
//! A reference to the external functions table
ExternalFunctionsTable &external_functions_table;
typedef map<int, NumConstNode *> num_const_node_map_type;
num_const_node_map_type num_const_node_map;
......@@ -80,7 +85,7 @@ private:
inline NodeID AddTrinaryOp(NodeID arg1, TrinaryOpcode op_code, NodeID arg2, NodeID arg3);
public:
DataTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg);
DataTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg, ExternalFunctionsTable &external_functions_table_arg);
virtual ~DataTree();
//! Some predefined constants
......@@ -175,9 +180,14 @@ public:
NodeID AddEqual(NodeID iArg1, NodeID iArg2);
//! Adds a model local variable with its value
void AddLocalVariable(const string &name, NodeID value) throw (LocalVariableException);
//! Adds an unknown function node
//! Adds an external function node
/*! \todo Use a map to share identical nodes */
NodeID AddUnknownFunction(const string &function_name, const vector<NodeID> &arguments);
NodeID AddExternalFunction(const string &function_name, const vector<NodeID> &arguments);
NodeID AddExternalFunction(int symb_id, const vector<NodeID> &arguments);
//! Adds an external function node for the first derivative of an external function
NodeID AddFirstDerivExternalFunctionNode(int top_level_symb_id, const vector<NodeID> &arguments, int input_index);
//! Adds an external function node for the second derivative of an external function
NodeID AddSecondDerivExternalFunctionNode(int top_level_symb_id, const vector<NodeID> &arguments, int input_index1, int input_index2);
//! Checks if a given symbol is used somewhere in the data tree
bool isSymbolUsed(int symb_id) const;
//! Checks if a given unary op is used somewhere in the data tree
......
/*
* Copyright (C) 2003-2009 Dynare Team
* Copyright (C) 2003-2010 Dynare Team
*
* This file is part of Dynare.
*
......@@ -36,8 +36,9 @@
#endif
DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg) :
ModelTree(symbol_table_arg, num_constants_arg),
NumericalConstants &num_constants_arg,
ExternalFunctionsTable &external_functions_table_arg) :
ModelTree(symbol_table_arg, num_constants_arg, external_functions_table_arg),
max_lag(0), max_lead(0),
max_endo_lag(0), max_endo_lead(0),
max_exo_lag(0), max_exo_lead(0),
......
/*
* Copyright (C) 2003-2009 Dynare Team
* Copyright (C) 2003-2010 Dynare Team
*
* This file is part of Dynare.
*
......@@ -199,7 +199,7 @@ private:
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;
public:
DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg, ExternalFunctionsTable &external_functions_table_argx);
//! Adds a variable node
/*! This implementation allows for non-zero lag */
virtual VariableNode *AddVariable(int symb_id, int lag = 0);
......
/*
* Copyright (C) 2003-2009 Dynare Team
* Copyright (C) 2003-2010 Dynare Team
*
* This file is part of Dynare.
*
......@@ -158,9 +158,10 @@ class ParsingDriver;
%token SVAR_IDENTIFICATION EQUATION EXCLUSION LAG UPPER_CHOLESKY LOWER_CHOLESKY
%token MARKOV_SWITCHING CHAIN STATE DURATION NUMBER_OF_STATES
%token SVAR COEFFICIENTS VARIANCES CONSTANTS EQUATIONS
%token EXTERNAL_FUNCTION EXT_FUNC_NAME EXT_FUNC_NARGS FIRST_DERIV_PROVIDED SECOND_DERIV_PROVIDED
%type <node_val> expression expression_or_empty
%type <node_val> equation hand_side model_var
%type <node_val> equation hand_side
%type <string_val> signed_float signed_integer prior
%type <string_val> filename symbol expectation_input
%type <string_val> value value1
......@@ -239,6 +240,7 @@ statement : parameters
| svar_identification
| markov_switching
| svar
| external_function
;
dsample : DSAMPLE INT_NUMBER ';'
......@@ -413,8 +415,8 @@ expression : '(' expression ')'
{ $$ = driver.add_max($3, $5); }
| MIN '(' expression COMMA expression ')'
{ $$ = driver.add_min($3, $5); }
| symbol '(' comma_expression ')'
{ $$ = driver.add_unknown_function($1); }
| symbol { driver.push_external_function_arg_vector_onto_stack(); } '(' comma_expression ')'
{ $$ = driver.add_model_var_or_external_function($1); }
| NORMCDF '(' expression COMMA expression COMMA expression ')'
{ $$ = driver.add_normcdf($3, $5, $7); }
| NORMCDF '(' expression ')'
......@@ -426,14 +428,14 @@ expression : '(' expression ')'
;
comma_expression : expression
{ driver.add_unknown_function_arg($1); }
| comma_expression COMMA expression
{ driver.add_unknown_function_arg($3); }
{ driver.add_external_function_arg($1); }
| comma_expression COMMA expression
{ driver.add_external_function_arg($3); }
;
expression_or_empty : {$$ = driver.add_nan_constant();}
| expression
;
| expression
;
initval : INITVAL ';' initval_list END
{ driver.end_initval(); }
......@@ -503,7 +505,8 @@ tag_pair : NAME EQUAL QUOTED_STRING
hand_side : '(' hand_side ')'
{ $$ = $2;}
| model_var
| symbol
{ $$ = driver.add_model_variable($1); }
| FLOAT_NUMBER
{ $$ = driver.add_constant($1); }
| INT_NUMBER
......@@ -562,6 +565,8 @@ hand_side : '(' hand_side ')'
{ $$ = driver.add_max($3, $5); }
| MIN '(' hand_side COMMA hand_side ')'
{ $$ = driver.add_min($3, $5); }
| symbol { driver.push_external_function_arg_vector_onto_stack(); } '(' comma_hand_side ')'
{ $$ = driver.add_model_var_or_external_function($1); }
| NORMCDF '(' hand_side COMMA hand_side COMMA hand_side ')'
{ $$ = driver.add_normcdf($3, $5, $7); }
| NORMCDF '(' hand_side ')'
......@@ -570,6 +575,11 @@ hand_side : '(' hand_side ')'
{ $$ = driver.add_steady_state($3); }
;
comma_hand_side : hand_side
{ driver.add_external_function_arg($1); }
| comma_hand_side COMMA hand_side
{ driver.add_external_function_arg($3); }
expectation_input : signed_integer
| VAROBS { string *varobs = new string("varobs"); $$ = varobs; }
| FULL { string *full = new string("full"); $$ = full; }
......@@ -578,12 +588,6 @@ expectation_input : signed_integer
pound_expression: '#' symbol EQUAL hand_side ';'
{ driver.declare_and_init_model_local_variable($2, $4); };
model_var : symbol
{ $$ = driver.add_model_variable($1); }
| symbol '(' signed_integer ')'
{ $$ = driver.add_model_variable($1, $3); }
;
shocks : SHOCKS ';' shock_list END { driver.end_shocks(); };
shock_list : shock_list shock_elem
......@@ -772,6 +776,20 @@ simul_options : o_periods
| o_minimal_solving_periods
;
external_function : EXTERNAL_FUNCTION '(' external_function_options_list ')' ';'
{ driver.external_function(); }
;
external_function_options_list : external_function_options_list COMMA external_function_options
| external_function_options
;
external_function_options : o_ext_func_name
| o_ext_func_nargs
| o_first_deriv_provided
| o_second_deriv_provided
;
stoch_simul : STOCH_SIMUL ';'
{ driver.stoch_simul(); }
| STOCH_SIMUL '(' stoch_simul_options_list ')' ';'
......@@ -1888,6 +1906,19 @@ o_equations : EQUATIONS EQUAL vec_int
o_instruments : INSTRUMENTS EQUAL '(' symbol_list ')' {driver.option_symbol_list("instruments"); };
o_ext_func_name : EXT_FUNC_NAME EQUAL filename { driver.external_function_option("name", $3); };
o_ext_func_nargs : EXT_FUNC_NARGS EQUAL INT_NUMBER { driver.external_function_option("nargs",$3); };
o_first_deriv_provided : FIRST_DERIV_PROVIDED EQUAL filename
{ driver.external_function_option("first_deriv_provided", $3); }
| FIRST_DERIV_PROVIDED
{ driver.external_function_option("first_deriv_provided", ""); }
;
o_second_deriv_provided : SECOND_DERIV_PROVIDED EQUAL filename
{ driver.external_function_option("second_deriv_provided", $3); }
| SECOND_DERIV_PROVIDED
{ driver.external_function_option("second_deriv_provided", ""); }
;
range : symbol ':' symbol
{
$1->append(":");
......
/*
* Copyright (C) 2003-2009 Dynare Team
* Copyright (C) 2003-2010 Dynare Team
*
* This file is part of Dynare.
*
......@@ -144,6 +144,7 @@ int sigma_e = 0;
<INITIAL>markov_switching {BEGIN DYNARE_STATEMENT; return token::MARKOV_SWITCHING;}
<INITIAL>svar {BEGIN DYNARE_STATEMENT; return token::SVAR;}
<INITIAL>external_function {BEGIN DYNARE_STATEMENT; return token::EXTERNAL_FUNCTION;}
/* End of a Dynare statement */
<DYNARE_STATEMENT>; {
......@@ -229,6 +230,10 @@ int sigma_e = 0;
<DYNARE_STATEMENT>aim_solver {return token::AIM_SOLVER;}
<DYNARE_STATEMENT>partial_information {return token::PARTIAL_INFORMATION;}
<DYNARE_STATEMENT>conditional_variance_decomposition {return token::CONDITIONAL_VARIANCE_DECOMPOSITION;}
<DYNARE_STATEMENT>name {return token::EXT_FUNC_NAME;}
<DYNARE_STATEMENT>nargs {return token::EXT_FUNC_NARGS;}
<DYNARE_STATEMENT>first_deriv_provided {return token::FIRST_DERIV_PROVIDED;}
<DYNARE_STATEMENT>second_deriv_provided {return token::SECOND_DERIV_PROVIDED;}
<DYNARE_STATEMENT>freq {return token::FREQ;}
<DYNARE_STATEMENT>initial_year {return token::INITIAL_YEAR;}
......@@ -531,16 +536,16 @@ int sigma_e = 0;
}
/* An instruction starting with a recognized symbol (which is not a modfile local
or an unknown function) is passed as NAME, otherwise it is a native statement
or an external function) is passed as NAME, otherwise it is a native statement
until the end of the line.
We exclude modfile local vars because the user may want to modify their value
using a Matlab assignment statement.
We also exclude unknown functions because the user may have used a Matlab matrix
element in initval (in which case Dynare recognizes the matrix name as an unknown
We also exclude external functions because the user may have used a Matlab matrix
element in initval (in which case Dynare recognizes the matrix name as an external
function symbol), and may want to modify the matrix later with Matlab statements.
*/
<INITIAL>[A-Za-z_][A-Za-z0-9_]* {
if (driver.symbol_exists_and_is_not_modfile_local_or_unknown_function(yytext))
if (driver.symbol_exists_and_is_not_modfile_local_or_external_function(yytext))
{
BEGIN DYNARE_STATEMENT;
yylval->string_val = new string(yytext);
......
/*
* Copyright (C) 2007-2009 Dynare Team
* Copyright (C) 2007-2010 Dynare Team
*
* This file is part of Dynare.
*
......@@ -359,7 +359,7 @@ VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg)
datatree.variable_node_map[make_pair(symb_id, lag)] = this;
// It makes sense to allow a lead/lag on parameters: during steady state calibration, endogenous and parameters can be swapped
assert(type != eUnknownFunction
assert(type != eExternalFunction
&& (lag == 0 || (type != eModelLocalVariable && type != eModFileLocalVariable)));
}
......@@ -389,7 +389,7 @@ VariableNode::prepareForDerivation()
case eModFileLocalVariable:
// Such a variable is never derived
break;
case eUnknownFunction:
case eExternalFunction:
cerr << "VariableNode::prepareForDerivation: impossible case" << endl;
exit(EXIT_FAILURE);
}
......@@ -413,7 +413,7 @@ VariableNode::computeDerivative(int deriv_id)
case eModFileLocalVariable:
cerr << "ModFileLocalVariable is not derivable" << endl;
exit(EXIT_FAILURE);
case eUnknownFunction:
case eExternalFunction:
cerr << "Impossible case!" << endl;
exit(EXIT_FAILURE);
}
......@@ -601,7 +601,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
}
break;
case eUnknownFunction:
case eExternalFunction:
cerr << "Impossible case" << endl;
exit(EXIT_FAILURE);
}
......@@ -767,7 +767,7 @@ VariableNode::getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recur
case eModFileLocalVariable:
cerr << "ModFileLocalVariable is not derivable" << endl;
exit(EXIT_FAILURE);
case eUnknownFunction:
case eExternalFunction:
cerr << "Impossible case!" << endl;
exit(EXIT_FAILURE);
}
......@@ -3174,7 +3174,7 @@ TrinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOp
return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
}
UnknownFunctionNode::UnknownFunctionNode(DataTree &datatree_arg,
ExternalFunctionNode::ExternalFunctionNode(DataTree &datatree_arg,
int symb_id_arg,
const vector<NodeID> &arguments_arg) :
ExprNode(datatree_arg),
......@@ -3184,37 +3184,65 @@ UnknownFunctionNode::UnknownFunctionNode(DataTree &datatree_arg,
}
void
UnknownFunctionNode::prepareForDerivation()
ExternalFunctionNode::prepareForDerivation()
{
cerr << "UnknownFunctionNode::prepareForDerivation: operation impossible!" << endl;
exit(EXIT_FAILURE);
if (preparedForDerivation)
return;
for (vector<NodeID>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
(*it)->prepareForDerivation();
non_null_derivatives = arguments.at(0)->non_null_derivatives;
for (int i = 1; i < (int)arguments.size(); i++)
set_union(non_null_derivatives.begin(),
non_null_derivatives.end(),
arguments.at(i)->non_null_derivatives.begin(),
arguments.at(i)->non_null_derivatives.end(),
inserter(non_null_derivatives, non_null_derivatives.begin()));
preparedForDerivation = true;
}
NodeID
UnknownFunctionNode::computeDerivative(int deriv_id)
ExternalFunctionNode::computeDerivative(int deriv_id)
{
cerr << "UnknownFunctionNode::computeDerivative: operation impossible!" << endl;
exit(EXIT_FAILURE);
vector<NodeID> dargs;
for (vector<NodeID>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
dargs.push_back((*it)->getDerivative(deriv_id));
return composeDerivatives(dargs);
}