2015-09-25 12:54:52 +02:00
function myoutput = posterior_sampler_core ( myinputs,fblck,nblck,whoiam, ThisMatlab)
% function myoutput = posterior_sampler_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
% Contains the most computationally intensive portion of code in
2016-05-19 17:13:39 +02:00
% posterior_sampler (the 'for xxx = fblck:nblck' loop). The branches in that 'for'
2015-09-25 12:54:52 +02:00
% cycle are completely independent to be suitable for parallel execution.
%
% INPUTS
% o myimput [struc] The mandatory variables for local/remote
2016-05-19 17:13:39 +02:00
% parallel computing obtained from posterior_sampler.m
2015-09-25 12:54:52 +02:00
% function.
% o fblck and nblck [integer] The Metropolis-Hastings chains.
% o whoiam [integer] In concurrent programming a modality to refer to the different threads 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, this is the original output of 'for b =
% fblck:nblck'. Otherwise, it's a portion of it computed on a specific core or
% remote machine. In this case:
% record;
% irun;
% NewFile;
% OutputFileName
%
% ALGORITHM
2016-05-19 17:13:39 +02:00
% Portion of Posterior Sampler.
2015-09-25 12:54:52 +02:00
%
% SPECIAL REQUIREMENTS.
% None.
2017-05-16 15:10:20 +02:00
%
2015-09-25 12:54:52 +02:00
% PARALLEL CONTEXT
2016-05-19 17:13:39 +02:00
% See the comments in the posterior_sampler.m funtion.
2015-09-25 12:54:52 +02:00
2023-04-26 10:34:25 +02:00
% Copyright © 2006-2023 Dynare Team
2015-09-25 12:54:52 +02:00
%
% 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
2021-06-09 17:33:48 +02:00
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
2015-09-25 12:54:52 +02:00
2017-05-16 12:42:01 +02:00
if nargin < 4
2015-09-25 12:54:52 +02:00
whoiam = 0 ;
end
% 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 ;
mh_bounds = myinputs . mh_bounds ;
last_draw = myinputs . ix2 ;
last_posterior = myinputs . ilogpo2 ;
fline = myinputs . fline ;
npar = myinputs . npar ;
nruns = myinputs . nruns ;
NewFile = myinputs . NewFile ;
MAX_nruns = myinputs . MAX_nruns ;
sampler_options = myinputs . sampler_options ;
d = myinputs . d ;
InitSizeArray = myinputs . InitSizeArray ;
record = myinputs . record ;
dataset_ = myinputs . dataset_ ;
dataset_info = myinputs . dataset_info ;
bayestopt_ = myinputs . bayestopt_ ;
estim_params_ = myinputs . estim_params_ ;
options_ = myinputs . options_ ;
M_ = myinputs . M_ ;
oo_ = myinputs . oo_ ;
% Necessary only for remote computing!
if whoiam
% initialize persistent variables in priordens()
priordens ( xparam1 , bayestopt_ . pshape , bayestopt_ . p6 , bayestopt_ . p7 , bayestopt_ . p3 , bayestopt_ . p4 , 1 ) ;
end
MetropolisFolder = CheckPath ( ' metropolis' , M_ . dname ) ;
ModelName = M_ . fname ;
BaseName = [ MetropolisFolder filesep ModelName ] ;
save_tmp_file = sampler_options . save_tmp_file ;
options_ . lik_algo = 1 ;
OpenOldFile = ones ( nblck , 1 ) ;
if strcmpi ( ProposalFun , ' rand_multivariate_normal' )
sampler_options . n = npar ;
sampler_options . ProposalDensity = ' multivariate_normal_pdf' ;
elseif strcmpi ( ProposalFun , ' rand_multivariate_student' )
2016-05-15 14:18:03 +02:00
sampler_options . n = sampler_options . student_degrees_of_freedom ;
2015-09-25 12:54:52 +02:00
sampler_options . ProposalDensity = ' multivariate_student_pdf' ;
end
%
% Now I run the (nblck-fblck+1) MCMC chains
%
sampler_options . xparam1 = xparam1 ;
2017-05-16 12:42:01 +02:00
if ~ isempty ( d )
2015-09-25 17:25:41 +02:00
sampler_options . proposal_covariance_Cholesky_decomposition = d * diag ( bayestopt_ . jscale ) ;
2016-06-14 15:09:05 +02:00
%store information for load_mh_file
record . ProposalCovariance = d ;
record . ProposalScaleVec = bayestopt_ . jscale ;
2015-09-25 17:25:41 +02:00
end
2015-09-25 12:54:52 +02:00
block_iter = 0 ;
2017-05-16 12:42:01 +02:00
for curr_block = fblck : nblck
2016-04-20 10:00:26 +02:00
LastSeeds = [ ] ;
2015-09-25 12:54:52 +02:00
block_iter = block_iter + 1 ;
try
% This will not work if the master uses a random number generator not
% available in the slave (different Matlab version or
% Matlab/Octave cluster). Therefore the trap.
%
% Set the random number generator type (the seed is useless but needed by the function)
2016-10-03 22:11:54 +02:00
if ~ isoctave
set_dynare_seed ( options_ . DynareRandomStreams . algo , options_ . DynareRandomStreams . seed ) ;
else
set_dynare_seed ( options_ . DynareRandomStreams . seed + curr_block ) ;
end
2015-09-25 12:54:52 +02:00
% Set the state of the RNG
2017-05-16 15:10:20 +02:00
set_dynare_random_generator_state ( record . InitialSeeds ( curr_block ) . Unifor , record . InitialSeeds ( curr_block ) . Normal ) ;
2015-09-25 12:54:52 +02:00
catch
2017-05-16 15:10:20 +02:00
% If the state set by master is incompatible with the slave, we only reseed
2015-09-25 12:54:52 +02:00
set_dynare_seed ( options_ . DynareRandomStreams . seed + curr_block ) ;
end
2019-07-24 16:16:03 +02:00
mh_recover_flag = 0 ;
2015-09-25 12:54:52 +02:00
if ( options_ . load_mh_file ~= 0 ) && ( fline ( curr_block ) > 1 ) && OpenOldFile ( curr_block ) %load previous draws and likelihood
load ( [ BaseName ' _mh' int2str ( NewFile ( curr_block ) ) ' _blck' int2str ( curr_block ) ' .mat' ] )
x2 = [ x2 ; zeros ( InitSizeArray ( curr_block ) - fline ( curr_block ) + 1 , npar ) ] ;
logpo2 = [ logpo2 ; zeros ( InitSizeArray ( curr_block ) - fline ( curr_block ) + 1 , 1 ) ] ;
OpenOldFile ( curr_block ) = 0 ;
else
2019-07-24 16:16:03 +02:00
if options_ . mh_recover && exist ( [ BaseName ' _mh_tmp_blck' int2str ( curr_block ) ' .mat' ] , ' file' ) == 2
load ( [ BaseName ' _mh_tmp_blck' int2str ( curr_block ) ' .mat' ] ) ;
draw_iter = size ( neval_this_chain , 2 ) + 1 ;
draw_index_current_file = draw_iter ;
feval_this_chain = sum ( sum ( neval_this_chain ) ) ;
feval_this_file = sum ( sum ( neval_this_chain ) ) ;
if feval_this_chain > draw_iter - 1
% non Metropolis type of sampler
accepted_draws_this_chain = draw_iter - 1 ;
accepted_draws_this_file = draw_iter - 1 ;
else
accepted_draws_this_chain = 0 ;
accepted_draws_this_file = 0 ;
end
mh_recover_flag = 1 ;
set_dynare_random_generator_state ( LastSeeds . ( [ ' file' int2str ( NewFile ( curr_block ) ) ] ) . Unifor , LastSeeds . ( [ ' file' int2str ( NewFile ( curr_block ) ) ] ) . Normal ) ;
last_draw ( curr_block , : ) = x2 ( draw_iter - 1 , : ) ;
last_posterior ( curr_block ) = logpo2 ( draw_iter - 1 ) ;
2019-12-20 16:28:06 +01:00
2019-07-24 16:16:03 +02:00
else
2019-12-20 16:28:06 +01:00
2019-07-24 16:16:03 +02:00
x2 = zeros ( InitSizeArray ( curr_block ) , npar ) ;
logpo2 = zeros ( InitSizeArray ( curr_block ) , 1 ) ;
end
2015-09-25 12:54:52 +02:00
end
%Prepare waiting bars
if whoiam
refresh_rate = sampler_options . parallel_bar_refresh_rate ;
bar_title = sampler_options . parallel_bar_title ;
prc0 = ( curr_block - fblck ) / ( nblck - fblck + 1 ) * ( isoctave || options_ . console_mode ) ;
hh = dyn_waitbar ( { prc0 , whoiam , options_ . parallel ( ThisMatlab ) } , [ bar_title ' (' int2str ( curr_block ) ' /' int2str ( options_ . mh_nblck ) ' )...' ] ) ;
else
refresh_rate = sampler_options . serial_bar_refresh_rate ;
bar_title = sampler_options . serial_bar_title ;
hh = dyn_waitbar ( 0 , [ bar_title ' (' int2str ( curr_block ) ' /' int2str ( options_ . mh_nblck ) ' )...' ] ) ;
set ( hh , ' Name' , bar_title ) ;
end
2019-07-24 16:16:03 +02:00
if mh_recover_flag == 0
accepted_draws_this_chain = 0 ;
accepted_draws_this_file = 0 ;
feval_this_chain = 0 ;
feval_this_file = 0 ;
draw_iter = 1 ;
draw_index_current_file = fline ( curr_block ) ; %get location of first draw in current block
end
sampler_options . curr_block = curr_block ;
2015-09-25 12:54:52 +02:00
while draw_iter < = nruns ( curr_block )
2017-05-16 15:10:20 +02:00
2015-09-25 12:54:52 +02:00
[ 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_ ) ;
x2 ( draw_index_current_file , : ) = par ;
last_draw ( curr_block , : ) = par ;
logpo2 ( draw_index_current_file ) = logpost ;
last_posterior ( curr_block ) = logpost ;
2019-12-20 16:28:06 +01:00
neval_this_chain ( : , draw_iter ) = neval ;
2015-09-25 12:54:52 +02:00
feval_this_chain = feval_this_chain + sum ( neval ) ;
feval_this_file = feval_this_file + sum ( neval ) ;
accepted_draws_this_chain = accepted_draws_this_chain + accepted ;
accepted_draws_this_file = accepted_draws_this_file + accepted ;
prtfrc = draw_iter / nruns ( curr_block ) ;
if mod ( draw_iter , refresh_rate ) == 0
if accepted_draws_this_chain / draw_iter == 1 && sum ( neval ) > 1
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 ) ] ) ;
else
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 ) ] ) ;
end
if save_tmp_file
2016-04-20 10:00:26 +02:00
[ LastSeeds . ( [ ' file' int2str ( NewFile ( curr_block ) ) ] ) . Unifor , LastSeeds . ( [ ' file' int2str ( NewFile ( curr_block ) ) ] ) . Normal ] = get_dynare_random_generator_state ( ) ;
2019-07-24 16:16:03 +02:00
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' ) ;
2015-09-25 12:54:52 +02:00
end
end
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
2016-04-20 10:00:26 +02:00
[ LastSeeds . ( [ ' file' int2str ( NewFile ( curr_block ) ) ] ) . Unifor , LastSeeds . ( [ ' file' int2str ( NewFile ( curr_block ) ) ] ) . Normal ] = get_dynare_random_generator_state ( ) ;
2017-05-16 12:42:01 +02:00
if save_tmp_file
2015-09-25 12:54:52 +02:00
delete ( [ BaseName ' _mh_tmp_blck' int2str ( curr_block ) ' .mat' ] ) ;
end
2019-07-24 16:16:03 +02:00
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' ) ;
2015-09-25 12:54:52 +02:00
fidlog = fopen ( [ MetropolisFolder ' /metropolis.log' ] , ' a' ) ;
fprintf ( fidlog , [ ' \n' ] ) ;
fprintf ( fidlog , [ ' %% Mh' int2str ( NewFile ( curr_block ) ) ' Blck' int2str ( curr_block ) ' (' datestr ( now , 0 ) ' )\n' ] ) ;
fprintf ( fidlog , ' \n' ) ;
fprintf ( fidlog , [ ' Number of simulations.: ' int2str ( length ( logpo2 ) ) ' \n' ] ) ;
fprintf ( fidlog , [ ' Acceptance ratio......: ' num2str ( accepted_draws_this_file / length ( logpo2 ) ) ' \n' ] ) ;
fprintf ( fidlog , [ ' Feval per iteration...: ' num2str ( feval_this_file / 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' ] ) ;
fprintf ( fidlog , [ ' Minimum value.........:\n' ] ) ;
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 ) ;
accepted_draws_this_file = 0 ;
feval_this_file = 0 ;
if draw_iter == nruns ( curr_block ) % I record the last draw...
record . LastParameters ( curr_block , : ) = x2 ( end , : ) ;
record . LastLogPost ( curr_block ) = logpo2 ( end ) ;
end
% size of next file in chain curr_block
InitSizeArray ( curr_block ) = min ( nruns ( curr_block ) - draw_iter , MAX_nruns ) ;
% initialization of next file if necessary
if InitSizeArray ( curr_block )
x2 = zeros ( InitSizeArray ( curr_block ) , npar ) ;
logpo2 = zeros ( InitSizeArray ( curr_block ) , 1 ) ;
NewFile ( curr_block ) = NewFile ( curr_block ) + 1 ;
draw_index_current_file = 0 ;
end
end
draw_iter = draw_iter + 1 ;
draw_index_current_file = draw_index_current_file + 1 ;
2019-12-20 16:28:06 +01:00
end % End of the simulations for one mh-block.
2015-09-25 12:54:52 +02:00
dyn_waitbar_close ( hh ) ;
2019-07-24 16:16:03 +02:00
if nruns ( curr_block )
record . AcceptanceRatio ( curr_block ) = accepted_draws_this_chain / ( draw_iter - 1 ) ;
record . FunctionEvalPerIteration ( curr_block ) = feval_this_chain / ( draw_iter - 1 ) ;
[ record . LastSeeds ( curr_block ) . Unifor , record . LastSeeds ( curr_block ) . Normal ] = get_dynare_random_generator_state ( ) ;
end
2015-09-25 12:54:52 +02:00
OutputFileName ( block_iter , : ) = { [ MetropolisFolder , filesep ] , [ ModelName ' _mh*_blck' int2str ( curr_block ) ' .mat' ] } ;
2019-12-20 16:28:06 +01:00
end % End of the loop over the mh-blocks.
2015-09-25 12:54:52 +02:00
myoutput . record = record ;
myoutput . irun = draw_index_current_file ;
myoutput . NewFile = NewFile ;
2019-07-24 16:16:03 +02:00
myoutput . OutputFileName = OutputFileName ;