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.
INPUTS
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
ALGORITHM
Portion of Metropolis-Hastings.
SPECIAL REQUIREMENTS.
None.0001 function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab) 0002 % PARALLEL CONTEXT 0003 % This function contain the most computationally intensive portion of code in 0004 % random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in 'for' 0005 % cycle and are completely independent than suitable to be executed in parallel way. 0006 % 0007 % INPUTS 0008 % o myimput [struc] The mandatory variables for local/remote 0009 % parallel computing obtained from random_walk_metropolis_hastings.m 0010 % function. 0011 % o fblck and nblck [integer] The Metropolis-Hastings chains. 0012 % o whoiam [integer] In concurrent programming a modality to refer to the differents thread running in parallel is needed. 0013 % The integer whoaim is the integer that 0014 % allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the 0015 % cluster. 0016 % o ThisMatlab [integer] Allows us to distinguish between the 0017 % 'main' matlab, the slave matlab worker, local matlab, remote matlab, 0018 % ... Then it is the index number of this slave machine in the cluster. 0019 % OUTPUTS 0020 % o myoutput [struc] 0021 % If executed without parallel is the original output of 'for b = 0022 % fblck:nblck' otherwise a portion of it computed on a specific core or 0023 % remote machine. In this case: 0024 % record; 0025 % irun; 0026 % NewFile; 0027 % OutputFileName 0028 % 0029 % ALGORITHM 0030 % Portion of Metropolis-Hastings. 0031 % 0032 % SPECIAL REQUIREMENTS. 0033 % None. 0034 0035 % PARALLEL CONTEXT 0036 % The most computationally intensive part of this function may be executed 0037 % in parallel. The code sutable to be executed in parallel on multi core or cluster machine, 0038 % is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion. 0039 % Then the DYNARE parallel package contain a set of pairs matlab functios that can be executed in 0040 % parallel and called name_function.m and name_function_core.m. 0041 % In addition in the parallel package we have second set of functions used 0042 % to manage the parallel computation. 0043 % 0044 % This function was the first function to be parallelized, later other 0045 % functions have been parallelized using the same methodology. 0046 % Then the comments write here can be used for all the other pairs of 0047 % parallel functions and also for management funtions. 0048 0049 0050 % Copyright (C) 2006-2011 Dynare Team 0051 % 0052 % This file is part of Dynare. 0053 % 0054 % Dynare is free software: you can redistribute it and/or modify 0055 % it under the terms of the GNU General Public License as published by 0056 % the Free Software Foundation, either version 3 of the License, or 0057 % (at your option) any later version. 0058 % 0059 % Dynare is distributed in the hope that it will be useful, 0060 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0061 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0062 % GNU General Public License for more details. 0063 % 0064 % You should have received a copy of the GNU General Public License 0065 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0066 0067 if nargin<4, 0068 whoiam=0; 0069 end 0070 0071 0072 global bayestopt_ estim_params_ options_ M_ oo_ 0073 0074 % reshape 'myinputs' for local computation. 0075 % In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by: 0076 0077 TargetFun=myinputs.TargetFun; 0078 ProposalFun=myinputs.ProposalFun; 0079 xparam1=myinputs.xparam1; 0080 vv=myinputs.vv; 0081 mh_bounds=myinputs.mh_bounds; 0082 ix2=myinputs.ix2; 0083 ilogpo2=myinputs.ilogpo2; 0084 ModelName=myinputs.ModelName; 0085 fline=myinputs.fline; 0086 npar=myinputs.npar; 0087 nruns=myinputs.nruns; 0088 NewFile=myinputs.NewFile; 0089 MAX_nruns=myinputs.MAX_nruns; 0090 d=myinputs.d; 0091 InitSizeArray=myinputs.InitSizeArray; 0092 record=myinputs.record; 0093 varargin=myinputs.varargin; 0094 0095 % Necessary only for remote computing! 0096 if whoiam 0097 Parallel=myinputs.Parallel; 0098 % initialize persistent variables in priordens() 0099 priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ... 0100 bayestopt_.p3,bayestopt_.p4,1); 0101 end 0102 0103 % (re)Set the penalty 0104 bayestopt_.penalty = Inf; 0105 0106 MhDirectoryName = CheckPath('metropolis',M_.dname); 0107 0108 options_.lik_algo = 1; 0109 OpenOldFile = ones(nblck,1); 0110 if strcmpi(ProposalFun,'rand_multivariate_normal') 0111 n = npar; 0112 elseif strcmpi(ProposalFun,'rand_multivariate_student') 0113 n = options_.student_degrees_of_freedom; 0114 end 0115 % load([MhDirectoryName '/' ModelName '_mh_history.mat'],'record'); 0116 %%%% 0117 %%%% NOW i run the (nblck-fblck+1) metropolis-hastings chains 0118 %%%% 0119 0120 0121 if any(isnan(bayestopt_.jscale)) 0122 if exist([ModelName '_optimal_mh_scale_parameter.mat'])% This file is created by mode_compute=6. 0123 load([ModelName '_optimal_mh_scale_parameter']) 0124 proposal_covariance_Cholesky_decomposition = d*Scale; 0125 else 0126 error('mh:: Something is wrong. I can''t figure out the value of the scale parameter.') 0127 end 0128 else 0129 proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale); 0130 end 0131 0132 0133 jloop=0; 0134 0135 JSUM = 0; 0136 for b = fblck:nblck, 0137 jloop=jloop+1; 0138 try % Trap in the case matlab slave is called from an octave master. 0139 randn('state',record.Seeds(b).Normal); 0140 rand('state',record.Seeds(b).Unifor); 0141 catch 0142 JSUM = JSUM + sum(100*clock); 0143 randn('state',JSUM); 0144 rand('state',JSUM); 0145 end 0146 if (options_.load_mh_file~=0) && (fline(b)>1) && OpenOldFile(b) 0147 load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ... 0148 '_blck' int2str(b) '.mat']) 0149 x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)]; 0150 logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)]; 0151 OpenOldFile(b) = 0; 0152 else 0153 x2 = zeros(InitSizeArray(b),npar); 0154 logpo2 = zeros(InitSizeArray(b),1); 0155 end 0156 if whoiam 0157 prc0=(b-fblck)/(nblck-fblck+1)*(exist('OCTAVE_VERSION') || options_.console_mode); 0158 hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},['MH (' int2str(b) '/' int2str(options_.mh_nblck) ')...']); 0159 else 0160 hh = dyn_waitbar(0,['Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']); 0161 set(hh,'Name','Metropolis-Hastings'); 0162 end 0163 isux = 0; 0164 jsux = 0; 0165 irun = fline(b); 0166 j = 1; 0167 while j <= nruns(b) 0168 par = feval(ProposalFun, ix2(b,:), proposal_covariance_Cholesky_decomposition, n); 0169 if all( par(:) > mh_bounds(:,1) ) && all( par(:) < mh_bounds(:,2) ) 0170 try 0171 logpost = - feval(TargetFun, par(:),varargin{:}); 0172 catch 0173 logpost = -inf; 0174 end 0175 else 0176 logpost = -inf; 0177 end 0178 if (logpost > -inf) && (log(rand) < logpost-ilogpo2(b)) 0179 x2(irun,:) = par; 0180 ix2(b,:) = par; 0181 logpo2(irun) = logpost; 0182 ilogpo2(b) = logpost; 0183 isux = isux + 1; 0184 jsux = jsux + 1; 0185 else 0186 x2(irun,:) = ix2(b,:); 0187 logpo2(irun) = ilogpo2(b); 0188 end 0189 prtfrc = j/nruns(b); 0190 % if exist('OCTAVE_VERSION') || options_.console_mode 0191 % if mod(j, 10) == 0 0192 % if exist('OCTAVE_VERSION') 0193 % if (whoiam==0) 0194 % printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j); 0195 % end 0196 % else 0197 % s0=repmat('\b',1,length(newString)); 0198 % newString=sprintf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acceptance rate: %3.f%%', b, nblck, 100 * prtfrc, 100 * isux / j); 0199 % fprintf([s0,'%s'],newString); 0200 % end 0201 % end 0202 % if mod(j,50)==0 && whoiam 0203 % % keyboard; 0204 % if (strcmp([options_.parallel(ThisMatlab).MatlabOctavePath], 'octave')) 0205 % waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%',100 * isux / j)]; 0206 % fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab)); 0207 % else 0208 % waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%', 100 * isux/j)]; 0209 % fMessageStatus((b-fblck)/(nblck-fblck+1)+prtfrc/(nblck-fblck+1),whoiam,waitbarString, '', options_.parallel(ThisMatlab)); 0210 % end 0211 % end 0212 % else 0213 % if mod(j, 3)==0 && ~whoiam 0214 % waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]); 0215 % elseif mod(j,50)==0 && whoiam, 0216 % % keyboard; 0217 % waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]; 0218 % fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab)); 0219 % end 0220 % end 0221 if (mod(j, 3)==0 && ~whoiam) || (mod(j,50)==0 && whoiam) 0222 dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('acceptation rate %4.3f', isux/j)]); 0223 end 0224 if (irun == InitSizeArray(b)) || (j == nruns(b)) % Now I save the simulations 0225 save([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2'); 0226 fidlog = fopen([MhDirectoryName '/metropolis.log'],'a'); 0227 fprintf(fidlog,['\n']); 0228 fprintf(fidlog,['%% Mh' int2str(NewFile(b)) 'Blck' int2str(b) ' (' datestr(now,0) ')\n']); 0229 fprintf(fidlog,' \n'); 0230 fprintf(fidlog,[' Number of simulations.: ' int2str(length(logpo2)) '\n']); 0231 fprintf(fidlog,[' Acceptation rate......: ' num2str(jsux/length(logpo2)) '\n']); 0232 fprintf(fidlog,[' Posterior mean........:\n']); 0233 for i=1:length(x2(1,:)) 0234 fprintf(fidlog,[' params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']); 0235 end 0236 fprintf(fidlog,[' log2po:' num2str(mean(logpo2)) '\n']); 0237 fprintf(fidlog,[' Minimum value.........:\n']); 0238 for i=1:length(x2(1,:)) 0239 fprintf(fidlog,[' params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']); 0240 end 0241 fprintf(fidlog,[' log2po:' num2str(min(logpo2)) '\n']); 0242 fprintf(fidlog,[' Maximum value.........:\n']); 0243 for i=1:length(x2(1,:)) 0244 fprintf(fidlog,[' params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']); 0245 end 0246 fprintf(fidlog,[' log2po:' num2str(max(logpo2)) '\n']); 0247 fprintf(fidlog,' \n'); 0248 fclose(fidlog); 0249 jsux = 0; 0250 if j == nruns(b) % I record the last draw... 0251 record.LastParameters(b,:) = x2(end,:); 0252 record.LastLogLiK(b) = logpo2(end); 0253 end 0254 % size of next file in chain b 0255 InitSizeArray(b) = min(nruns(b)-j,MAX_nruns); 0256 % initialization of next file if necessary 0257 if InitSizeArray(b) 0258 x2 = zeros(InitSizeArray(b),npar); 0259 logpo2 = zeros(InitSizeArray(b),1); 0260 NewFile(b) = NewFile(b) + 1; 0261 irun = 0; 0262 end 0263 end 0264 j=j+1; 0265 irun = irun + 1; 0266 end% End of the simulations for one mh-block. 0267 record.AcceptationRates(b) = isux/j; 0268 % if exist('OCTAVE_VERSION') || options_.console_mode || whoiam 0269 % if exist('OCTAVE_VERSION') 0270 % printf('\n'); 0271 % else 0272 % fprintf('\n'); 0273 % end 0274 % diary on; 0275 % else %if ~whoiam 0276 % close(hh); 0277 % end 0278 dyn_waitbar_close(hh); 0279 record.Seeds(b).Normal = randn('state'); 0280 record.Seeds(b).Unifor = rand('state'); 0281 OutputFileName(jloop,:) = {[MhDirectoryName,filesep], [ModelName '_mh*_blck' int2str(b) '.mat']}; 0282 end% End of the loop over the mh-blocks. 0283 0284 0285 myoutput.record = record; 0286 myoutput.irun = irun; 0287 myoutput.NewFile = NewFile; 0288 myoutput.OutputFileName = OutputFileName;