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