1function [x,info,fvec,fjac] = dynare_solve(func,x,options,varargin)
2% function [x,info,fvec,fjac] = dynare_solve(func,x,options,varargin)
3% proposes different solvers
4%
5% INPUTS
6%    func:             name of the function to be solved
7%    x:                guess values
8%    options:          struct of Dynare options
9%    varargin:         list of arguments following jacobian_flag
10%
11% OUTPUTS
12%    x:                solution
13%    info=1:           the model can not be solved
14%    fvec:             Function value (used for debugging when check=1)
15%    fjac:             Jacobian (used for debugging when check=1)
16%
17% SPECIAL REQUIREMENTS
18%    none
19
20% Copyright (C) 2001-2017 Dynare Team
21%
22% This file is part of Dynare.
23%
24% Dynare is free software: you can redistribute it and/or modify
25% it under the terms of the GNU General Public License as published by
26% the Free Software Foundation, either version 3 of the License, or
27% (at your option) any later version.
28%
29% Dynare is distributed in the hope that it will be useful,
30% but WITHOUT ANY WARRANTY; without even the implied warranty of
31% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
32% GNU General Public License for more details.
33%
34% You should have received a copy of the GNU General Public License
35% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
36
37% jacobian_flag=true:  jacobian given by the 'func' function
38% jacobian_flag=false:  jacobian obtained numerically
39jacobian_flag = options.jacobian_flag;
40
41% Set tolerance parameter depending the the caller function.
42stack = dbstack;
43if isoctave
44    [pathstr,name,ext]=fileparts(stack(2).file);
45    caller_file_name=[name,ext];
46else
47    caller_file_name=stack(2).file;
48end
49if strcmp(caller_file_name,'solve_stacked_problem.m')
50    tolf = options.dynatol.f;
51    tolx = options.dynatol.x;
52else
53    tolf = options.solve_tolf;
54    tolx = options.solve_tolx;
55end
56
57if strcmp(caller_file_name,'dyn_ramsey_static.m')
58    maxit = options.ramsey.maxit;
59else
60    maxit = options.steady.maxit;
61end
62
63info = 0;
64nn = size(x,1);
65
66% Get status of the initial guess (default values?)
67if any(x)
68    % The current initial guess is not the default for all the variables.
69    idx = find(x);      % Indices of the variables with default initial guess values.
70    in0 = length(idx);
71else
72    % The current initial guess is the default for all the variables.
73    idx = transpose(1:nn);
74    in0 = nn;
75end
76
77% checking initial values
78if jacobian_flag
79    [fvec, fjac] = feval(func, x, varargin{:});
80    wrong_initial_guess_flag = false;
81    if ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) ...
82            || any(~isreal(fvec)) || any(~isreal(fjac(:)))
83        if max(abs(fvec)) < tolf %return if initial value solves problem
84            info = 0;
85            return;
86        end
87        disp('Randomize initial guess...')
88        % Let's try random numbers for the variables initialized with the default value.
89        wrong_initial_guess_flag = true;
90        % First try with positive numbers.
91        tentative_number = 0;
92        while wrong_initial_guess_flag && tentative_number<=in0*10
93            tentative_number = tentative_number+1;
94            x(idx) = rand(in0, 1)*10;
95            [fvec, fjac] = feval(func, x, varargin{:});
96            wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
97        end
98        % If all previous attempts failed, try with real numbers.
99        tentative_number = 0;
100        while wrong_initial_guess_flag && tentative_number<=in0*10
101            tentative_number = tentative_number+1;
102            x(idx) = randn(in0, 1)*10;
103            [fvec, fjac] = feval(func, x, varargin{:});
104            wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
105        end
106        % Last tentative, ff all previous attempts failed, try with negative numbers.
107        tentative_number = 0;
108        while wrong_initial_guess_flag && tentative_number<=in0*10
109            tentative_number = tentative_number+1;
110            x(idx) = -rand(in0, 1)*10;
111            [fvec, fjac] = feval(func, x, varargin{:});
112            wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
113        end
114    end
115else
116    fvec = feval(func,x,varargin{:});
117    fjac = zeros(nn,nn);
118    wrong_initial_guess_flag = false;
119    if ~all(isfinite(fvec))
120        % Let's try random numbers for the variables initialized with the default value.
121        wrong_initial_guess_flag = true;
122        % First try with positive numbers.
123        tentative_number = 0;
124        while wrong_initial_guess_flag && tentative_number<=in0*10
125            tentative_number = tentative_number+1;
126            x(idx) = rand(in0, 1)*10;
127            fvec = feval(func, x, varargin{:});
128            wrong_initial_guess_flag = ~all(isfinite(fvec));
129        end
130        % If all previous attempts failed, try with real numbers.
131        tentative_number = 0;
132        while wrong_initial_guess_flag && tentative_number<=in0*10
133            tentative_number = tentative_number+1;
134            x(idx) = randn(in0, 1)*10;
135            fvec = feval(func, x, varargin{:});
136            wrong_initial_guess_flag = ~all(isfinite(fvec));
137        end
138        % Last tentative, ff all previous attempts failed, try with negative numbers.
139        tentative_number = 0;
140        while wrong_initial_guess_flag && tentative_number<=in0*10
141            tentative_number = tentative_number+1;
142            x(idx) = -rand(in0, 1)*10;
143            fvec = feval(func, x, varargin{:});
144            wrong_initial_guess_flag = ~all(isfinite(fvec));
145        end
146    end
147end
148
149% Exit with error if no initial guess has been found.
150if wrong_initial_guess_flag
151    info=1;
152    x = NaN(size(fvec));
153    return
154end
155
156% this test doesn't check complementarity conditions and is not used for
157% mixed complementarity problems
158if (~ismember(options.solve_algo,[10,11])) && (max(abs(fvec)) < tolf)
159    return ;
160end
161
162if options.solve_algo == 0
163    if ~isoctave
164        if ~user_has_matlab_license('optimization_toolbox')
165            error('You can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
166        end
167    end
168    options4fsolve=optimset('fsolve');
169    options4fsolve.MaxFunEvals = 50000;
170    options4fsolve.MaxIter = maxit;
171    options4fsolve.TolFun = tolf;
172    if options.debug==1
173        options4fsolve.Display = 'final';
174    else
175        options4fsolve.Display = 'off';
176    end
177    if jacobian_flag
178        options4fsolve.Jacobian = 'on';
179    else
180        options4fsolve.Jacobian = 'off';
181    end
182    if ~isoctave
183        [x,fval,exitval,output] = fsolve(func,x,options4fsolve,varargin{:});
184    else
185        % Under Octave, use a wrapper, since fsolve() does not have a 4th arg
186        if ischar(func)
187            func2 = str2func(func);
188        else
189            func2 = func;
190        end
191        func = @(x) func2(x, varargin{:});
192        % The Octave version of fsolve does not converge when it starts from the solution
193        fvec = feval(func,x);
194        if max(abs(fvec)) >= tolf
195            [x,fval,exitval,output] = fsolve(func,x,options4fsolve);
196        else
197            exitval = 3;
198        end
199    end
200
201    if exitval == 1
202        info = 0;
203    elseif exitval > 1
204        if ischar(func)
205            func2 = str2func(func);
206        else
207            func2 = func;
208        end
209        func = @(x) func2(x, varargin{:});
210        fvec = feval(func,x);
211        if max(abs(fvec)) >= tolf
212            info = 1;
213        else
214            info = 0;
215        end
216    else
217        info = 1;
218    end
219elseif options.solve_algo == 1
220    [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,options.gstep, ...
221                    tolf,tolx, ...
222                    maxit,options.debug,varargin{:});
223elseif options.solve_algo == 9
224    [x,info]=trust_region(func,x,1:nn,1:nn,jacobian_flag,options.gstep, ...
225                          tolf,tolx, ...
226                          maxit,options.debug,varargin{:});
227elseif options.solve_algo == 2 || options.solve_algo == 4
228
229    if options.solve_algo == 2
230        solver = @solve1;
231    else
232        solver = @trust_region;
233    end
234
235    if ~jacobian_flag
236        fjac = zeros(nn,nn) ;
237        dh = max(abs(x),options.gstep(1)*ones(nn,1))*eps^(1/3);
238        for j = 1:nn
239            xdh = x ;
240            xdh(j) = xdh(j)+dh(j) ;
241            fjac(:,j) = (feval(func,xdh,varargin{:}) - fvec)./dh(j) ;
242        end
243    end
244
245    [j1,j2,r,s] = dmperm(fjac);
246
247    if options.debug
248        disp(['DYNARE_SOLVE (solve_algo=2|4): number of blocks = ' num2str(length(r))]);
249    end
250
251    for i=length(r)-1:-1:1
252        if options.debug
253            disp(['DYNARE_SOLVE (solve_algo=2|4): solving block ' num2str(i) ', of size ' num2str(r(i+1)-r(i)) ]);
254        end
255        [x,info]=solver(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag, ...
256                        options.gstep, ...
257                        tolf,tolx, ...
258                        maxit,options.debug,varargin{:});
259        if info
260            return
261        end
262    end
263    fvec = feval(func,x,varargin{:});
264    if max(abs(fvec)) > tolf
265        [x,info]=solver(func,x,1:nn,1:nn,jacobian_flag, ...
266                        options.gstep, tolf,tolx, ...
267                        maxit,options.debug,varargin{:});
268    end
269elseif options.solve_algo == 3
270    if jacobian_flag
271        [x,info] = csolve(func,x,func,1e-6,500,varargin{:});
272    else
273        [x,info] = csolve(func,x,[],1e-6,500,varargin{:});
274    end
275    [fvec, fjac] = feval(func, x, varargin{:});
276elseif options.solve_algo == 10
277    % LMMCP
278    olmmcp = options.lmmcp;
279
280    [x,fval,exitflag] = lmmcp(func,x,olmmcp.lb,olmmcp.ub,olmmcp,varargin{:});
281    if exitflag == 1
282        info = 0;
283    else
284        info = 1;
285    end
286elseif options.solve_algo == 11
287    % PATH mixed complementary problem
288    % PATH linear mixed complementary problem
289    if ~exist('mcppath')
290        error(['PATH can''t be provided with Dynare. You need to install it ' ...
291               'yourself and add its location to Matlab/Octave path before ' ...
292               'running Dynare'])
293    end
294    omcppath = options.mcppath;
295    global mcp_data
296    mcp_data.func = func;
297    mcp_data.args = varargin;
298    info=0;
299    try
300        [x,fval,jac,mu] = pathmcp(x,omcppath.lb,omcppath.ub,'mcp_func',omcppath.A,omcppath.b,omcppath.t,omcppath.mu0);
301    catch
302        info = 1;
303    end
304else
305    error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11]')
306end
307