diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m index 4b4ba596457c8eed3be97f190fd047d6e4390063..fb4bf8f992a4212bbd772784770015b19f906faf 100644 --- a/matlab/posterior_sampler_core.m +++ b/matlab/posterior_sampler_core.m @@ -141,31 +141,32 @@ for curr_block = fblck:nblck options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+curr_block); end mh_recover_flag=0; - 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)]; + if options_.mh_recover && exist([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'file')==2 && OpenOldFile(curr_block) + % this should be done whatever value of load_mh_file + load([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']); + draw_iter = size(neval_this_chain,2)+1; + draw_index_current_file = draw_iter+fline(curr_block)-1; + feval_this_chain = sum(sum(neval_this_chain)); + feval_this_file = sum(sum(neval_this_chain)); + if feval_this_chain>draw_index_current_file-fline(curr_block) + % non Metropolis type of sampler + accepted_draws_this_chain = draw_index_current_file-fline(curr_block); + accepted_draws_this_file = draw_index_current_file-fline(curr_block); + 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_index_current_file-1,:); + last_posterior(curr_block)=logpo2(draw_index_current_file-1); OpenOldFile(curr_block) = 0; else - 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); - + 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 x2 = zeros(InitSizeArray(curr_block),npar); diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m index f18ee5558d6de29bb971d07cdb29d315cdd093f6..d80644b1b5dafedf4b6c5afbfe9f8da8a3578dd5 100644 --- a/matlab/posterior_sampler_initialization.m +++ b/matlab/posterior_sampler_initialization.m @@ -465,7 +465,7 @@ elseif options_.mh_recover AllMhFiles = dir([BaseName '_mh*_blck*.mat']); TotalNumberOfMhFiles = length(AllMhFiles)-length(dir([BaseName '_mh_tmp*_blck*.mat'])); % Quit if no crashed mcmc chain can be found as there are as many files as expected - if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles) + if (ExpectedNumberOfMhFilesPerBlock>LastFileNumberInThePreviousMh) && (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles) if isnumeric(options_.parallel) fprintf('%s: It appears that you don''t need to use the mh_recover option!\n',dispString); fprintf(' You have to edit the mod file and remove the mh_recover option\n'); @@ -476,15 +476,23 @@ elseif options_.mh_recover % 2. Something needs to be done; find out what % Count the number of saved mh files per block. NumberOfMhFilesPerBlock = zeros(NumberOfBlocks,1); + is_chain_complete = true(NumberOfBlocks,1); for b = 1:NumberOfBlocks BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']); NumberOfMhFilesPerBlock(b) = length(BlckMhFiles)-length(dir([BaseName '_mh_tmp*_blck' int2str(b) '.mat'])); + if (ExpectedNumberOfMhFilesPerBlock==LastFileNumberInThePreviousMh) + % here we need to check for the number of lines + tmpdata = load([BaseName '_mh' int2str(LastFileNumberInThePreviousMh) '_blck' int2str(b) '.mat'],'logpo2'); + if length(tmpdata.logpo2) == LastLineNumberInThePreviousMh + is_chain_complete(b) = false; + end + end end % Find FirstBlock (First block), an integer targeting the crashed mcmc chain. FirstBlock = 1; %initialize FBlock = zeros(NumberOfBlocks,1); while FirstBlock <= NumberOfBlocks - if NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock + if (NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock) || not(is_chain_complete(FirstBlock)) fprintf('%s: Chain %u is not complete!\n', dispString,FirstBlock); FBlock(FirstBlock)=1; else