random_walk_metropolis_hastings_core.m 12.2 KB
Newer Older
1
function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
2
3
4
5
% PARALLEL CONTEXT
% This function contain the most computationally intensive portion of code in
% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in 'for'
% cycle and are completely independent than suitable to be executed in parallel way.
6
7
%
% INPUTS
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
%   o myimput            [struc]     The mandatory variables for local/remote
%                                    parallel computing obtained from random_walk_metropolis_hastings.m
%                                    function.
%   o fblck and nblck    [integer]   The Metropolis-Hastings chains.
%   o whoiam             [integer]   In concurrent programming a modality to refer to the differents thread running in parallel is needed.
%                                    The integer whoaim is the integer that
%                                    allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the
%                                    cluster.
%   o ThisMatlab         [integer]   Allows us to distinguish between the
%                                    'main' matlab, the slave matlab worker, local matlab, remote matlab,
%                                     ... Then it is the index number of this slave machine in the cluster.
% OUTPUTS
%   o myoutput  [struc]
%               If executed without parallel is the original output of 'for b =
%               fblck:nblck' otherwise a portion of it computed on a specific core or
%               remote machine. In this case:
%                               record;
%                               irun;
%                               NewFile;
%                               OutputFileName
%
29
30
% ALGORITHM
%   Portion of Metropolis-Hastings.
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
%
% SPECIAL REQUIREMENTS.
%   None.

% PARALLEL CONTEXT
% The most computationally intensive part of this function may be executed
% in parallel. The code sutable to be executed in parallel on multi core or cluster machine,
% is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion.
% Then the DYNARE parallel package contain a set of pairs matlab functios that can be executed in
% parallel and called name_function.m and name_function_core.m.
% In addition in the parallel package we have second set of functions used
% to manage the parallel computation.
%
% This function was the first function to be parallelized, later other
% functions have been parallelized using the same methodology.
% Then the comments write here can be used for all the other pairs of
% parallel functions and also for management funtions.
48

49

50
% Copyright (C) 2006-2011 Dynare Team
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
%
% 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/>.

if nargin<4,
    whoiam=0;
end

71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
% reshape 'myinputs' for local computation.
% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:

TargetFun=myinputs.TargetFun;
ProposalFun=myinputs.ProposalFun;
xparam1=myinputs.xparam1;
vv=myinputs.vv;
mh_bounds=myinputs.mh_bounds;
ix2=myinputs.ix2;
ilogpo2=myinputs.ilogpo2;
ModelName=myinputs.ModelName;
fline=myinputs.fline;
npar=myinputs.npar;
nruns=myinputs.nruns;
NewFile=myinputs.NewFile;
MAX_nruns=myinputs.MAX_nruns;
d=myinputs.d;
88
InitSizeArray=myinputs.InitSizeArray;
89
record=myinputs.record;
90
91
92
93
94
95
dataset_ = myinputs.dataset_;
bayestopt_ = myinputs.bayestopt_;
estim_params_ = myinputs.estim_params_;
options_ = myinputs.options_;
M_ = myinputs.M_;
oo_ = myinputs.oo_;
96
97
98
99
varargin=myinputs.varargin;

% Necessary only for remote computing!
if whoiam
100
101
102
    Parallel=myinputs.Parallel;
    % initialize persistent variables in priordens()
    priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ...
103
        bayestopt_.p3,bayestopt_.p4,1);
104
end
105

106
MhDirectoryName = CheckPath('metropolis',M_.dname);
107
108
109
110
111
112
113
114

options_.lik_algo = 1;
OpenOldFile = ones(nblck,1);
if strcmpi(ProposalFun,'rand_multivariate_normal')
    n = npar;
elseif strcmpi(ProposalFun,'rand_multivariate_student')
    n = options_.student_degrees_of_freedom;
end
115

116
117
118
%%%%
%%%% NOW i run the (nblck-fblck+1) metropolis-hastings chains
%%%%
119

120
proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
121
122
123

jloop=0;

124
JSUM = 0;
125
126
for b = fblck:nblck,
    jloop=jloop+1;
127
128
129
130
    try
        % this will not work if the master uses a random generator not
        % available in the slave (different Matlab version or
        % Matlab/Octave cluster). Therefor the trap.
131
132
133
134
135
136
        
        % this set the random generator type (the seed is useless but
        % needed by the function)
        set_dynare_seed(options_.DynareRandomStreams.algo,...
                        options_.DynareRandomStreams.seed);
        % this set the state 
137
138
139
        set_dynare_random_generator_state(record.Seeds(b).Unifor, ...
                                          record.Seeds(b).Normal);
    catch
140
141
        % if the state set by master is incompatible with the slave, we
        % only reseed 
142
143
        set_dynare_seed(options_.DynareRandomStreams.seed+b);
    end
144
    if (options_.load_mh_file~=0) && (fline(b)>1) && OpenOldFile(b)
145
        load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ...
146
            '_blck' int2str(b) '.mat'])
147
148
149
150
151
152
153
        x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)];
        logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)];
        OpenOldFile(b) = 0;
    else
        x2 = zeros(InitSizeArray(b),npar);
        logpo2 = zeros(InitSizeArray(b),1);
    end
154
155
    if whoiam
        prc0=(b-fblck)/(nblck-fblck+1)*(exist('OCTAVE_VERSION') || options_.console_mode);
156
157
158
159
        hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},['MH (' int2str(b) '/' int2str(options_.mh_nblck) ')...']);
    else
        hh = dyn_waitbar(0,['Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']);
        set(hh,'Name','Metropolis-Hastings');
160
161
162
163
164
165
    end
    isux = 0;
    jsux = 0;
    irun = fline(b);
    j = 1;
    while j <= nruns(b)
166
        par = feval(ProposalFun, ix2(b,:), proposal_covariance_Cholesky_decomposition, n);
167
        if all( par(:) > mh_bounds(:,1) ) && all( par(:) < mh_bounds(:,2) )
168
            try
169
                logpost = - feval(TargetFun, par(:),dataset_,options_,M_,estim_params_,bayestopt_,oo_);
170
171
172
173
174
175
            catch
                logpost = -inf;
            end
        else
            logpost = -inf;
        end
176
        if (logpost > -inf) && (log(rand) < logpost-ilogpo2(b))
177
178
            x2(irun,:) = par;
            ix2(b,:) = par;
179
            logpo2(irun) = logpost;
180
181
182
            ilogpo2(b) = logpost;
            isux = isux + 1;
            jsux = jsux + 1;
183
        else
184
185
186
187
            x2(irun,:) = ix2(b,:);
            logpo2(irun) = ilogpo2(b);
        end
        prtfrc = j/nruns(b);
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
%         if exist('OCTAVE_VERSION') || options_.console_mode
%             if mod(j, 10) == 0
%                 if exist('OCTAVE_VERSION')
%                     if (whoiam==0)
%                         printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
%                     end
%                 else
%                     s0=repmat('\b',1,length(newString));
%                     newString=sprintf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acceptance rate: %3.f%%', b, nblck, 100 * prtfrc, 100 * isux / j);
%                     fprintf([s0,'%s'],newString);
%                 end
%             end
%             if mod(j,50)==0 && whoiam
%                 %             keyboard;
%                 if (strcmp([options_.parallel(ThisMatlab).MatlabOctavePath], 'octave'))
%                     waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%',100 * isux / j)];
%                     fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab));
%                 else
%                     waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%', 100 * isux/j)];
%                     fMessageStatus((b-fblck)/(nblck-fblck+1)+prtfrc/(nblck-fblck+1),whoiam,waitbarString, '', options_.parallel(ThisMatlab));
%                 end
%             end
%         else
%             if mod(j, 3)==0 && ~whoiam
%                 waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]);
%             elseif mod(j,50)==0 && whoiam,
%                 %             keyboard;
%                 waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)];
%                 fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab));
%             end
%         end
        if (mod(j, 3)==0 && ~whoiam) || (mod(j,50)==0 && whoiam)
            dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('acceptation rate %4.3f', isux/j)]);
221
        end
222
        if (irun == InitSizeArray(b)) || (j == nruns(b)) % Now I save the simulations
223
224
225
226
227
228
229
230
231
232
233
234
            save([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2');
            fidlog = fopen([MhDirectoryName '/metropolis.log'],'a');
            fprintf(fidlog,['\n']);
            fprintf(fidlog,['%% Mh' int2str(NewFile(b)) 'Blck' int2str(b) ' (' datestr(now,0) ')\n']);
            fprintf(fidlog,' \n');
            fprintf(fidlog,['  Number of simulations.: ' int2str(length(logpo2)) '\n']);
            fprintf(fidlog,['  Acceptation rate......: ' num2str(jsux/length(logpo2)) '\n']);
            fprintf(fidlog,['  Posterior mean........:\n']);
            for i=1:length(x2(1,:))
                fprintf(fidlog,['    params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']);
            end
            fprintf(fidlog,['    log2po:' num2str(mean(logpo2)) '\n']);
235
            fprintf(fidlog,['  Minimum value.........:\n']);
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
            for i=1:length(x2(1,:))
                fprintf(fidlog,['    params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']);
            end
            fprintf(fidlog,['    log2po:' num2str(min(logpo2)) '\n']);
            fprintf(fidlog,['  Maximum value.........:\n']);
            for i=1:length(x2(1,:))
                fprintf(fidlog,['    params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']);
            end
            fprintf(fidlog,['    log2po:' num2str(max(logpo2)) '\n']);
            fprintf(fidlog,' \n');
            fclose(fidlog);
            jsux = 0;
            if j == nruns(b) % I record the last draw...
                record.LastParameters(b,:) = x2(end,:);
                record.LastLogLiK(b) = logpo2(end);
            end
            % size of next file in chain b
            InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
            % initialization of next file if necessary
            if InitSizeArray(b)
                x2 = zeros(InitSizeArray(b),npar);
                logpo2 = zeros(InitSizeArray(b),1);
                NewFile(b) = NewFile(b) + 1;
                irun = 0;
            end
        end
        j=j+1;
        irun = irun + 1;
    end% End of the simulations for one mh-block.
    record.AcceptationRates(b) = isux/j;
266
267
268
269
270
271
272
273
274
275
276
%     if exist('OCTAVE_VERSION') || options_.console_mode || whoiam
%         if exist('OCTAVE_VERSION')
%             printf('\n');
%         else
%             fprintf('\n');
%         end
%         diary on;
%     else %if ~whoiam
%         close(hh);
%     end
    dyn_waitbar_close(hh);
277
    [record.Seeds(b).Unifor, record.Seeds(b).Normal] = get_dynare_random_generator_state();
278
279
280
281
282
283
284
    OutputFileName(jloop,:) = {[MhDirectoryName,filesep], [ModelName '_mh*_blck' int2str(b) '.mat']};
end% End of the loop over the mh-blocks.


myoutput.record = record;
myoutput.irun = irun;
myoutput.NewFile = NewFile;
285
myoutput.OutputFileName = OutputFileName;