2009-12-16 18:17:34 +01:00
function myoutput = random_walk_metropolis_hastings_core ( myinputs,fblck,nblck,whoiam, ThisMatlab)
2015-04-24 18:12:46 +02:00
% function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
% Contains the most computationally intensive portion of code in
% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in that 'for'
% cycle are completely independent to be suitable for parallel execution.
2011-02-02 14:13:11 +01:00
%
% INPUTS
2010-05-31 11:55:25 +02:00
% 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.
2015-04-24 18:12:46 +02:00
% o whoiam [integer] In concurrent programming a modality to refer to the different threads running in parallel is needed.
2010-05-31 11:55:25 +02:00
% 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
2015-04-24 18:12:46 +02:00
% 'main' Matlab, the slave Matlab worker, local Matlab, remote Matlab,
2010-05-31 11:55:25 +02:00
% ... Then it is the index number of this slave machine in the cluster.
% OUTPUTS
% o myoutput [struc]
2015-04-24 18:12:46 +02:00
% 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
2010-05-31 11:55:25 +02:00
% remote machine. In this case:
% record;
% irun;
% NewFile;
% OutputFileName
%
2011-02-02 14:13:11 +01:00
% ALGORITHM
% Portion of Metropolis-Hastings.
2010-05-31 11:55:25 +02:00
%
% SPECIAL REQUIREMENTS.
% None.
2015-04-24 18:12:46 +02:00
%
2010-05-31 11:55:25 +02:00
% PARALLEL CONTEXT
2015-04-24 18:12:46 +02:00
% See the comments in the random_walk_metropolis_hastings.m funtion.
2009-12-16 18:17:34 +01:00
2010-05-31 11:55:25 +02:00
2015-04-24 18:12:46 +02:00
% Copyright (C) 2006-2015 Dynare Team
2009-12-16 18:17:34 +01: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
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin < 4 ,
whoiam = 0 ;
end
2010-05-31 11:55:25 +02:00
% 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 ;
2015-04-24 18:12:46 +02:00
last_draw = myinputs . ix2 ;
last_posterior = myinputs . ilogpo2 ;
2010-05-31 11:55:25 +02:00
fline = myinputs . fline ;
npar = myinputs . npar ;
nruns = myinputs . nruns ;
NewFile = myinputs . NewFile ;
MAX_nruns = myinputs . MAX_nruns ;
d = myinputs . d ;
2011-02-02 14:13:11 +01:00
InitSizeArray = myinputs . InitSizeArray ;
2010-05-31 11:55:25 +02:00
record = myinputs . record ;
2012-08-05 15:10:21 +02:00
dataset_ = myinputs . dataset_ ;
2014-06-16 17:41:59 +02:00
dataset_info = myinputs . dataset_info ;
2012-08-05 15:10:21 +02:00
bayestopt_ = myinputs . bayestopt_ ;
estim_params_ = myinputs . estim_params_ ;
options_ = myinputs . options_ ;
M_ = myinputs . M_ ;
oo_ = myinputs . oo_ ;
2010-05-31 11:55:25 +02:00
% Necessary only for remote computing!
if whoiam
2011-02-02 14:13:11 +01:00
% initialize persistent variables in priordens()
2013-11-20 18:03:12 +01:00
priordens ( xparam1 , bayestopt_ . pshape , bayestopt_ . p6 , bayestopt_ . p7 , bayestopt_ . p3 , bayestopt_ . p4 , 1 ) ;
2010-05-31 11:55:25 +02:00
end
2009-12-16 18:17:34 +01:00
2013-11-20 18:03:12 +01:00
MetropolisFolder = CheckPath ( ' metropolis' , M_ . dname ) ;
ModelName = M_ . fname ;
BaseName = [ MetropolisFolder filesep ModelName ] ;
2009-12-16 18:17:34 +01:00
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
2012-10-06 16:45:19 +02:00
2013-11-20 18:03:12 +01:00
%
2015-04-24 18:12:46 +02:00
% Now I run the (nblck-fblck+1) Metropolis-Hastings chains
2013-11-20 18:03:12 +01:00
%
2010-09-15 08:59:51 +02:00
2012-10-06 16:45:19 +02:00
proposal_covariance_Cholesky_decomposition = d * diag ( bayestopt_ . jscale ) ;
2009-12-16 18:17:34 +01:00
2015-04-24 18:12:46 +02:00
block_iter = 0 ;
2009-12-16 18:17:34 +01:00
2015-04-24 18:12:46 +02:00
for curr_block = fblck : nblck ,
block_iter = block_iter + 1 ;
2012-08-30 12:24:05 +02:00
try
2015-04-24 18:12:46 +02:00
% This will not work if the master uses a random number generator not
2012-08-30 12:24:05 +02:00
% available in the slave (different Matlab version or
2015-04-24 18:12:46 +02:00
% Matlab/Octave cluster). Therefore the trap.
2013-11-20 18:03:12 +01:00
%
2015-04-24 18:12:46 +02:00
% Set the random number generator type (the seed is useless but needed by the function)
2013-11-20 18:03:12 +01:00
set_dynare_seed ( options_ . DynareRandomStreams . algo , options_ . DynareRandomStreams . seed ) ;
2015-04-24 18:12:46 +02:00
% Set the state of the RNG
set_dynare_random_generator_state ( record . InitialSeeds ( curr_block ) . Unifor , record . InitialSeeds ( curr_block ) . Normal ) ;
2012-08-30 12:24:05 +02:00
catch
2015-04-24 18:12:46 +02:00
% If the state set by master is incompatible with the slave, we only reseed
set_dynare_seed ( options_ . DynareRandomStreams . seed + curr_block ) ;
2012-08-30 12:24:05 +02:00
end
2015-04-24 18:12:46 +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 ;
2009-12-16 18:17:34 +01:00
else
2015-04-24 18:12:46 +02:00
x2 = zeros ( InitSizeArray ( curr_block ) , npar ) ;
logpo2 = zeros ( InitSizeArray ( curr_block ) , 1 ) ;
2009-12-16 18:17:34 +01:00
end
2015-04-24 18:12:46 +02:00
%Prepare waiting bars
2010-12-17 09:21:30 +01:00
if whoiam
2015-04-24 18:12:46 +02:00
prc0 = ( curr_block - fblck ) / ( nblck - fblck + 1 ) * ( isoctave || options_ . console_mode ) ;
hh = dyn_waitbar ( { prc0 , whoiam , options_ . parallel ( ThisMatlab ) } , [ ' MH (' int2str ( curr_block ) ' /' int2str ( options_ . mh_nblck ) ' )...' ] ) ;
2011-12-13 18:32:57 +01:00
else
2015-04-24 18:12:46 +02:00
hh = dyn_waitbar ( 0 , [ ' Metropolis-Hastings (' int2str ( curr_block ) ' /' int2str ( options_ . mh_nblck ) ' )...' ] ) ;
2011-12-13 18:32:57 +01:00
set ( hh , ' Name' , ' Metropolis-Hastings' ) ;
2009-12-16 18:17:34 +01:00
end
2015-04-24 18:12:46 +02:00
accepted_draws_this_chain = 0 ;
accepted_draws_this_file = 0 ;
draw_index_current_file = fline ( curr_block ) ; %get location of first draw in current block
draw_iter = 1 ;
while draw_iter < = nruns ( curr_block )
par = feval ( ProposalFun , last_draw ( curr_block , : ) , proposal_covariance_Cholesky_decomposition , n ) ;
2014-10-20 16:18:54 +02:00
if all ( par ( : ) > mh_bounds . lb ) && all ( par ( : ) < mh_bounds . ub )
2009-12-16 18:17:34 +01:00
try
2014-10-20 16:18:54 +02:00
logpost = - feval ( TargetFun , par ( : ) , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , mh_bounds , oo_ ) ;
2009-12-16 18:17:34 +01:00
catch
logpost = - inf ;
end
else
logpost = - inf ;
end
2015-04-24 18:12:46 +02:00
if ( logpost > - inf ) && ( log ( rand ) < logpost - last_posterior ( curr_block ) )
x2 ( draw_index_current_file , : ) = par ;
last_draw ( curr_block , : ) = par ;
logpo2 ( draw_index_current_file ) = logpost ;
last_posterior ( curr_block ) = logpost ;
accepted_draws_this_chain = accepted_draws_this_chain + 1 ;
accepted_draws_this_file = accepted_draws_this_file + 1 ;
2011-02-02 14:13:11 +01:00
else
2015-04-24 18:12:46 +02:00
x2 ( draw_index_current_file , : ) = last_draw ( curr_block , : ) ;
logpo2 ( draw_index_current_file ) = last_posterior ( curr_block ) ;
2009-12-16 18:17:34 +01:00
end
2015-04-24 18:12:46 +02:00
prtfrc = draw_iter / nruns ( curr_block ) ;
if ( mod ( draw_iter , 3 ) == 0 && ~ whoiam ) || ( mod ( draw_iter , 50 ) == 0 && whoiam )
dyn_waitbar ( prtfrc , hh , [ ' MH (' int2str ( curr_block ) ' /' int2str ( options_ . mh_nblck ) ' ) ' sprintf ( ' Current acceptance ratio %4.3f' , accepted_draws_this_chain / draw_iter ) ] ) ;
2009-12-16 18:17:34 +01:00
end
2015-04-24 18:12:46 +02:00
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
save ( [ BaseName ' _mh' int2str ( NewFile ( curr_block ) ) ' _blck' int2str ( curr_block ) ' .mat' ] , ' x2' , ' logpo2' ) ;
2013-11-20 18:03:12 +01:00
fidlog = fopen ( [ MetropolisFolder ' /metropolis.log' ] , ' a' ) ;
2009-12-16 18:17:34 +01:00
fprintf ( fidlog , [ ' \n' ] ) ;
2015-04-24 18:12:46 +02:00
fprintf ( fidlog , [ ' %% Mh' int2str ( NewFile ( curr_block ) ) ' Blck' int2str ( curr_block ) ' (' datestr ( now , 0 ) ' )\n' ] ) ;
2009-12-16 18:17:34 +01:00
fprintf ( fidlog , ' \n' ) ;
fprintf ( fidlog , [ ' Number of simulations.: ' int2str ( length ( logpo2 ) ) ' \n' ] ) ;
2015-04-24 18:12:46 +02:00
fprintf ( fidlog , [ ' Acceptance ratio......: ' num2str ( accepted_draws_this_file / length ( logpo2 ) ) ' \n' ] ) ;
2009-12-16 18:17:34 +01:00
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' ] ) ;
2011-02-10 15:54:23 +01:00
fprintf ( fidlog , [ ' Minimum value.........:\n' ] ) ;
2009-12-16 18:17:34 +01:00
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 ) ;
2015-04-24 18:12:46 +02:00
accepted_draws_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 ) ;
2009-12-16 18:17:34 +01:00
end
2015-04-24 18:12:46 +02:00
% size of next file in chain curr_block
InitSizeArray ( curr_block ) = min ( nruns ( curr_block ) - draw_iter , MAX_nruns ) ;
2009-12-16 18:17:34 +01:00
% initialization of next file if necessary
2015-04-24 18:12:46 +02:00
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 ;
2009-12-16 18:17:34 +01:00
end
end
2015-04-24 18:12:46 +02:00
draw_iter = draw_iter + 1 ;
draw_index_current_file = draw_index_current_file + 1 ;
2009-12-16 18:17:34 +01:00
end % End of the simulations for one mh-block.
2015-04-24 18:12:46 +02:00
record . AcceptanceRatio ( curr_block ) = accepted_draws_this_chain / draw_iter ;
2011-12-13 18:32:57 +01:00
dyn_waitbar_close ( hh ) ;
2015-04-24 18:12:46 +02:00
[ record . LastSeeds ( curr_block ) . Unifor , record . LastSeeds ( curr_block ) . Normal ] = get_dynare_random_generator_state ( ) ;
OutputFileName ( block_iter , : ) = { [ MetropolisFolder , filesep ] , [ ModelName ' _mh*_blck' int2str ( curr_block ) ' .mat' ] } ;
2009-12-16 18:17:34 +01:00
end % End of the loop over the mh-blocks.
myoutput . record = record ;
2015-04-24 18:12:46 +02:00
myoutput . irun = draw_index_current_file ;
2009-12-16 18:17:34 +01:00
myoutput . NewFile = NewFile ;
2010-06-25 14:58:53 +02:00
myoutput . OutputFileName = OutputFileName ;