1function myoutput = posterior_sampler_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
2% function myoutput = posterior_sampler_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
3% Contains the most computationally intensive portion of code in
4% posterior_sampler (the 'for xxx = fblck:nblck' loop). The branches in  that 'for'
5% cycle are completely independent to be suitable for parallel execution.
6%
7% INPUTS
8%   o myimput            [struc]     The mandatory variables for local/remote
9%                                    parallel computing obtained from posterior_sampler.m
10%                                    function.
11%   o fblck and nblck    [integer]   The Metropolis-Hastings chains.
12%   o whoiam             [integer]   In concurrent programming a modality to refer to the different threads running in parallel is needed.
13%                                    The integer whoaim is the integer that
14%                                    allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the
15%                                    cluster.
16%   o ThisMatlab         [integer]   Allows us to distinguish between the
17%                                    'main' Matlab, the slave Matlab worker, local Matlab, remote Matlab,
18%                                     ... Then it is the index number of this slave machine in the cluster.
19% OUTPUTS
20%   o myoutput  [struc]
21%               If executed without parallel, this is the original output of 'for b =
22%               fblck:nblck'. Otherwise, it's a portion of it computed on a specific core or
23%               remote machine. In this case:
24%                               record;
25%                               irun;
26%                               NewFile;
27%                               OutputFileName
28%
29% ALGORITHM
30%   Portion of Posterior Sampler.
31%
32% SPECIAL REQUIREMENTS.
33%   None.
34%
35% PARALLEL CONTEXT
36% See the comments in the posterior_sampler.m funtion.
37
38
39% Copyright (C) 2006-2017 Dynare Team
40%
41% This file is part of Dynare.
42%
43% Dynare is free software: you can redistribute it and/or modify
44% it under the terms of the GNU General Public License as published by
45% the Free Software Foundation, either version 3 of the License, or
46% (at your option) any later version.
47%
48% Dynare is distributed in the hope that it will be useful,
49% but WITHOUT ANY WARRANTY; without even the implied warranty of
50% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
51% GNU General Public License for more details.
52%
53% You should have received a copy of the GNU General Public License
54% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
55
56if nargin<4
57    whoiam=0;
58end
59
60% reshape 'myinputs' for local computation.
61% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
62
63TargetFun=myinputs.TargetFun;
64ProposalFun=myinputs.ProposalFun;
65xparam1=myinputs.xparam1;
66mh_bounds=myinputs.mh_bounds;
67last_draw=myinputs.ix2;
68last_posterior=myinputs.ilogpo2;
69fline=myinputs.fline;
70npar=myinputs.npar;
71nruns=myinputs.nruns;
72NewFile=myinputs.NewFile;
73MAX_nruns=myinputs.MAX_nruns;
74sampler_options=myinputs.sampler_options;
75d=myinputs.d;
76InitSizeArray=myinputs.InitSizeArray;
77record=myinputs.record;
78dataset_ = myinputs.dataset_;
79dataset_info = myinputs.dataset_info;
80bayestopt_ = myinputs.bayestopt_;
81estim_params_ = myinputs.estim_params_;
82options_ = myinputs.options_;
83M_ = myinputs.M_;
84oo_ = myinputs.oo_;
85% Necessary only for remote computing!
86if whoiam
87    % initialize persistent variables in priordens()
88    priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1);
89    % initialize persistent variables in prior_draw()
90    prior_draw(bayestopt_,options_.prior_trunc);
91end
92
93MetropolisFolder = CheckPath('metropolis',M_.dname);
94ModelName = M_.fname;
95BaseName = [MetropolisFolder filesep ModelName];
96save_tmp_file = sampler_options.save_tmp_file;
97
98options_.lik_algo = 1;
99OpenOldFile = ones(nblck,1);
100if strcmpi(ProposalFun,'rand_multivariate_normal')
101    sampler_options.n = npar;
102    sampler_options.ProposalDensity = 'multivariate_normal_pdf';
103elseif strcmpi(ProposalFun,'rand_multivariate_student')
104    sampler_options.n = sampler_options.student_degrees_of_freedom;
105    sampler_options.ProposalDensity = 'multivariate_student_pdf';
106end
107
108%
109% Now I run the (nblck-fblck+1) MCMC chains
110%
111
112sampler_options.xparam1 = xparam1;
113if ~isempty(d)
114    sampler_options.proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
115    %store information for load_mh_file
116    record.ProposalCovariance=d;
117    record.ProposalScaleVec=bayestopt_.jscale;
118end
119
120block_iter=0;
121
122for curr_block = fblck:nblck
123    LastSeeds=[];
124    block_iter=block_iter+1;
125    try
126        % This will not work if the master uses a random number generator not
127        % available in the slave (different Matlab version or
128        % Matlab/Octave cluster). Therefore the trap.
129        %
130        % Set the random number generator type (the seed is useless but needed by the function)
131        if ~isoctave
132            set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed);
133        else
134            set_dynare_seed(options_.DynareRandomStreams.seed+curr_block);
135        end
136        % Set the state of the RNG
137        set_dynare_random_generator_state(record.InitialSeeds(curr_block).Unifor, record.InitialSeeds(curr_block).Normal);
138    catch
139        % If the state set by master is incompatible with the slave, we only reseed
140        set_dynare_seed(options_.DynareRandomStreams.seed+curr_block);
141    end
142    mh_recover_flag=0;
143    if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
144        load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'])
145        x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)];
146        logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)];
147        OpenOldFile(curr_block) = 0;
148    else
149        if options_.mh_recover && exist([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'file')==2
150            load([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']);
151            draw_iter = size(neval_this_chain,2)+1;
152            draw_index_current_file = draw_iter;
153            feval_this_chain = sum(sum(neval_this_chain));
154            feval_this_file = sum(sum(neval_this_chain));
155            if feval_this_chain>draw_iter-1
156                % non Metropolis type of sampler
157                accepted_draws_this_chain = draw_iter-1;
158                accepted_draws_this_file = draw_iter-1;
159            else
160                accepted_draws_this_chain = 0;
161                accepted_draws_this_file = 0;
162            end
163            mh_recover_flag=1;
164            set_dynare_random_generator_state(LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal);
165            last_draw(curr_block,:)=x2(draw_iter-1,:);
166            last_posterior(curr_block)=logpo2(draw_iter-1);
167
168        else
169
170            x2 = zeros(InitSizeArray(curr_block),npar);
171            logpo2 = zeros(InitSizeArray(curr_block),1);
172        end
173    end
174    %Prepare waiting bars
175    if whoiam
176        refresh_rate = sampler_options.parallel_bar_refresh_rate;
177        bar_title = sampler_options.parallel_bar_title;
178        prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode);
179        hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
180    else
181        refresh_rate = sampler_options.serial_bar_refresh_rate;
182        bar_title = sampler_options.serial_bar_title;
183        hh = dyn_waitbar(0,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
184        set(hh,'Name',bar_title);
185    end
186    if mh_recover_flag==0
187        accepted_draws_this_chain = 0;
188        accepted_draws_this_file = 0;
189        feval_this_chain = 0;
190        feval_this_file = 0;
191        draw_iter = 1;
192        draw_index_current_file = fline(curr_block); %get location of first draw in current block
193    end
194    sampler_options.curr_block = curr_block;
195    while draw_iter <= nruns(curr_block)
196
197        [par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFun, last_draw(curr_block,:), last_posterior(curr_block), sampler_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
198
199        x2(draw_index_current_file,:) = par;
200        last_draw(curr_block,:) = par;
201        logpo2(draw_index_current_file) = logpost;
202        last_posterior(curr_block) = logpost;
203        neval_this_chain(:, draw_iter) = neval;
204        feval_this_chain = feval_this_chain + sum(neval);
205        feval_this_file = feval_this_file + sum(neval);
206        accepted_draws_this_chain = accepted_draws_this_chain + accepted;
207        accepted_draws_this_file = accepted_draws_this_file + accepted;
208
209        prtfrc = draw_iter/nruns(curr_block);
210        if mod(draw_iter, refresh_rate)==0
211            if accepted_draws_this_chain/draw_iter==1 && sum(neval)>1
212                dyn_waitbar(prtfrc,hh,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ') ' sprintf('Function eval per draw %4.3f', feval_this_chain/draw_iter)]);
213            else
214                dyn_waitbar(prtfrc,hh,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', accepted_draws_this_chain/draw_iter)]);
215            end
216            if save_tmp_file
217                [LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal] = get_dynare_random_generator_state();
218                save([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'x2','logpo2','LastSeeds','neval_this_chain','accepted_draws_this_chain','accepted_draws_this_file','feval_this_chain','feval_this_file');
219            end
220        end
221        if (draw_index_current_file == InitSizeArray(curr_block)) || (draw_iter == nruns(curr_block)) % Now I save the simulations, either because the current file is full or the chain is done
222            [LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal] = get_dynare_random_generator_state();
223            if save_tmp_file
224                delete([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']);
225            end
226            save([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'],'x2','logpo2','LastSeeds','accepted_draws_this_chain','accepted_draws_this_file','feval_this_chain','feval_this_file');
227            fidlog = fopen([MetropolisFolder '/metropolis.log'],'a');
228            fprintf(fidlog,['\n']);
229            fprintf(fidlog,['%% Mh' int2str(NewFile(curr_block)) 'Blck' int2str(curr_block) ' (' datestr(now,0) ')\n']);
230            fprintf(fidlog,' \n');
231            fprintf(fidlog,['  Number of simulations.: ' int2str(length(logpo2)) '\n']);
232            fprintf(fidlog,['  Acceptance ratio......: ' num2str(accepted_draws_this_file/length(logpo2)) '\n']);
233            fprintf(fidlog,['  Feval per iteration...: ' num2str(feval_this_file/length(logpo2)) '\n']);
234            fprintf(fidlog,['  Posterior mean........:\n']);
235            for i=1:length(x2(1,:))
236                fprintf(fidlog,['    params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']);
237            end
238            fprintf(fidlog,['    log2po:' num2str(mean(logpo2)) '\n']);
239            fprintf(fidlog,['  Minimum value.........:\n']);
240            for i=1:length(x2(1,:))
241                fprintf(fidlog,['    params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']);
242            end
243            fprintf(fidlog,['    log2po:' num2str(min(logpo2)) '\n']);
244            fprintf(fidlog,['  Maximum value.........:\n']);
245            for i=1:length(x2(1,:))
246                fprintf(fidlog,['    params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']);
247            end
248            fprintf(fidlog,['    log2po:' num2str(max(logpo2)) '\n']);
249            fprintf(fidlog,' \n');
250            fclose(fidlog);
251            accepted_draws_this_file = 0;
252            feval_this_file = 0;
253            if draw_iter == nruns(curr_block) % I record the last draw...
254                record.LastParameters(curr_block,:) = x2(end,:);
255                record.LastLogPost(curr_block) = logpo2(end);
256            end
257            % size of next file in chain curr_block
258            InitSizeArray(curr_block) = min(nruns(curr_block)-draw_iter,MAX_nruns);
259            % initialization of next file if necessary
260            if InitSizeArray(curr_block)
261                x2 = zeros(InitSizeArray(curr_block),npar);
262                logpo2 = zeros(InitSizeArray(curr_block),1);
263                NewFile(curr_block) = NewFile(curr_block) + 1;
264                draw_index_current_file = 0;
265            end
266        end
267        draw_iter=draw_iter+1;
268        draw_index_current_file = draw_index_current_file + 1;
269    end % End of the simulations for one mh-block.
270    dyn_waitbar_close(hh);
271    if nruns(curr_block)
272        record.AcceptanceRatio(curr_block) = accepted_draws_this_chain/(draw_iter-1);
273        record.FunctionEvalPerIteration(curr_block) = feval_this_chain/(draw_iter-1);
274        [record.LastSeeds(curr_block).Unifor, record.LastSeeds(curr_block).Normal] = get_dynare_random_generator_state();
275    end
276    OutputFileName(block_iter,:) = {[MetropolisFolder,filesep], [ModelName '_mh*_blck' int2str(curr_block) '.mat']};
277end % End of the loop over the mh-blocks.
278
279
280myoutput.record = record;
281myoutput.irun = draw_index_current_file;
282myoutput.NewFile = NewFile;
283myoutput.OutputFileName = OutputFileName;
284