diff --git a/matlab/AnalyseComputationalEnviroment.m b/matlab/AnalyseComputationalEnviroment.m index 0028bcbba..3453e1cc0 100644 --- a/matlab/AnalyseComputationalEnviroment.m +++ b/matlab/AnalyseComputationalEnviroment.m @@ -1,5 +1,12 @@ function [ErrorCode] = AnalyseComputationalEnviroment(DataInput) -% Input/Output description: + +% DESCRIPTION + +% This function is used to check the user computational request. +% If no error happen the function return 0. + + +% INPUT/OUTPUT description: % % DataInput is the strcture option_.parallel, with the follow fields: % @@ -37,10 +44,17 @@ function [ErrorCode] = AnalyseComputationalEnviroment(DataInput) % % Value 3: The remote computer is unreachable!!! % -% Value 4: The user name and/or password is/are incorrect on the -% remote computer! +% Value 4: The fields user name and/or password are/is empty! +% +% Value 5: Remote Drive and/or Remote Folder not exist! +% +% Value 6: It is impossible write/read file on remote computer. +% +% Value 7: The valu user and/or passwd are incorret or the user have +% no permission to execute a Matlab section. +% +% % -% Value 5: It is impossible write/read file on remote computer. % % Then at the point call of this function it is possible react in a best way, in accord % with the ErrorCode. @@ -95,13 +109,13 @@ if (DataInput.Local == 1) % We look for the information on local computer hardware. - si=[]; - de=[]; + si0=[]; + de0=[]; - [si de]=system(['psinfo \\']); + [si0 de0]=system(['psinfo \\']); RealNumCPU=-1; - RealNumCPU=GiveCPUnumber(de); + RealNumCPU=GiveCPUnumber(de0); % Trasforming the input data provided in a form [n1:n2] in a single numerical % value. @@ -125,15 +139,15 @@ end % In this case we need more sophisticated check. - + if (DataInput.Local == 0) - si=[]; - de=[]; + si1=[]; + de1=[]; - [si de]=system(['ping ', DataInput.PcName]); + [si1 de1]=system(['ping ', DataInput.PcName]); - if si==1 + if si1==1 % It is impossiblie to be connected to the % remote computer. @@ -141,56 +155,93 @@ if (DataInput.Local == 0) return; end - - % -> IL CODICE SEGUENTE E' DA CONTROLLARE E VERIFICARE! - % The Local Machine can be connetted with Remote Computer. + % Now we verify if user name and password are correct and if remote % drive and remote folder exist on the remote computer and it is % possible to exchange data with them. + + + if (isempty(DataInput.user)) || (isempty(DataInput.passwd)) + + % The fields user name and/or password are/is empty! + + ErrorCode=4; + return - si=[]; - de=[]; - - [si de]=system(['psinfo \\', DataInput.PcName, ' -u ',DataInput.user, ' -p ',DataInput.passwd ]); - - if si<0 - % It is possible to be connected to the remote computer but it is not usable because the user - % name and/or password is/are incorrect. - - ErrorCodeComputer=4; - return; - else - % Username and Password are correct! end - % Now we verify if it possible to exchange data with the remote - % computer: + + % Now we very if RemoteDrive and/or RemoteFolder exist on remote + % computer + StartPwd=pwd; - fid = fopen('Tracing.txt','w+'); - fclose (fid); + try + cd(['\\',DataInput.PcName,'\',DataInput.RemoteDrive,'$\',DataInput.RemoteFolder]); + catch + + cd ([StartPwd]); + disp 'Remote Drive and/or Remote Folder not exist!'; + + ErrorCode=5; + return + + end - % ATTENZIONE: verificare perche sembra funzionare anche se il RemoteFolder non - % esiste. + cd ([StartPwd]); - Status=movefile('Tracing.txt', ['\\',DataInput.PcName,'\',DataInput.RemoteDrive,'$\',DataInput.RemoteFolder]); + + % Now we verify if it possible to exchange data with the remote + % computer: + + + Status=copyfile('Tracing.m', ['\\',DataInput.PcName,'\',DataInput.RemoteDrive,'$\',DataInput.RemoteFolder]); if Status==1 % Remote Drive/Folder exist on Remote computer and % it is possible to exchange data with him. else - % Move file error! - ErrorCodeComputer=5; + ErrorCode=6; return; end + + % Now we verify if it is possible execute a matlab section on remote + % machine when the user is .user with password .passwd + + si2=[]; + de2=[]; + + [si2 de2]=system(['start /B /WAIT psexec \\',DataInput.PcName,' -e -u ',DataInput.user,' -p ',DataInput.passwd,' -W ',DataInput.RemoteDrive,':\',DataInput.RemoteFolder, ' -low matlab -nosplash -nodesktop -minimize -r Tracing']); + + NoError='error code 0'; + + StrError= findstr(NoError,de2); + + + if isempty (StrError) + + % Bad user and/or passwd! + ErrorCode=7; + return; + + else + + % No error it is possible execute a matlab section on remote + % machine when the user is .user with password .passwd + end + + % At this point we can to analyze the remote computer hardware. - RealNumCPU=-1; - RealNumCPU=GiveCPUnumber(de); + + + [si0 de0]=system(['psinfo \\']); + RealNumCPU=GiveCPUnumber(de0); + % Trasforming the input data provided in a form [n1:n2] in a single numerical % value. diff --git a/matlab/DiffuseKalmanSmoother1.m b/matlab/DiffuseKalmanSmoother1.m index ef64185df..def1b080c 100644 --- a/matlab/DiffuseKalmanSmoother1.m +++ b/matlab/DiffuseKalmanSmoother1.m @@ -131,7 +131,7 @@ while notsteady & t. + + +nCPU=-1; + +OffSet=27; + +SringPosition=strfind(ComputerInformations, 'Processors:'); +nCPU=ComputerInformations(SringPosition+OffSet); + +% We check if there are Processors/Cores more than 9. + + + t0=ComputerInformations(SringPosition+OffSet+1); + t1=str2num(t0); + t1=isempty(t1); + + % if t1 is 0 the machine have more than 9 CPU. + + if t1==0 + nCPU=strcat(nCPU,t0); + end + + nCPU=str2num(nCPU); + + return + \ No newline at end of file diff --git a/matlab/McMCDiagnostics.m b/matlab/McMCDiagnostics.m index c449d9158..c90b67464 100644 --- a/matlab/McMCDiagnostics.m +++ b/matlab/McMCDiagnostics.m @@ -101,7 +101,7 @@ if isnumeric(options_.parallel), UDIAG = fout.UDIAG; clear fout else - [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, 1, npar,{},'McMCDiagnostics_core', localVars); + [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, 1, npar,{},'McMCDiagnostics_core', localVars, [], options_.parallel_info); UDIAG = fout(1).UDIAG; for j=2:totCPU, UDIAG = cat(3,UDIAG ,fout(j).UDIAG); diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m index c9e75fe13..47548110d 100644 --- a/matlab/PosteriorIRF.m +++ b/matlab/PosteriorIRF.m @@ -28,6 +28,7 @@ function PosteriorIRF(type) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . + global options_ estim_params_ oo_ M_ bayestopt_ % Set the number of periods if isempty(options_.irf) | ~options_.irf @@ -71,32 +72,32 @@ else end DirectoryName = CheckPath('Output'); if strcmpi(type,'posterior') - MhDirectoryName = CheckPath('metropolis'); + MhDirectoryName = CheckPath('metropolis'); elseif strcmpi(type,'gsa') - MhDirectoryName = CheckPath('GSA'); + MhDirectoryName = CheckPath('GSA'); else - MhDirectoryName = CheckPath('prior'); + MhDirectoryName = CheckPath('prior'); end if strcmpi(type,'posterior') - load([ MhDirectoryName '/' M_.fname '_mh_history.mat']) - TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); - NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); + load([ MhDirectoryName filesep M_.fname '_mh_history.mat']) + TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); + NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); elseif strcmpi(type,'gsa') - load([ MhDirectoryName '/' M_.fname '_prior.mat'],'lpmat0','lpmat','istable') - x=[lpmat0(istable,:) lpmat(istable,:)]; - clear lpmat istable - NumberOfDraws=size(x,1); - B=NumberOfDraws; options_.B = B; + load([ MhDirectoryName filesep M_.fname '_prior.mat'],'lpmat0','lpmat','istable') + x=[lpmat0(istable,:) lpmat(istable,:)]; + clear lpmat istable + NumberOfDraws=size(x,1); + B=NumberOfDraws; options_.B = B; else% type = 'prior' - NumberOfDraws = 500; + NumberOfDraws = 500; end if ~strcmpi(type,'gsa') - B = min([round(.5*NumberOfDraws),500]); options_.B = B; + B = min([round(.5*NumberOfDraws),500]); options_.B = B; end -try delete([MhDirectoryName '/' M_.fname '_irf_dsge*.mat']) +try delete([MhDirectoryName filesep M_.fname '_irf_dsge*.mat']) catch disp('No _IRFs (dsge) files to be deleted!') end -try delete([MhDirectoryName '/' M_.fname '_irf_bvardsge*.mat']) +try delete([MhDirectoryName filesep M_.fname '_irf_bvardsge*.mat']) catch disp('No _IRFs (bvar-dsge) files to be deleted!') end irun = 0; @@ -105,23 +106,16 @@ irun2 = 0; NumberOfIRFfiles_dsge = 1; NumberOfIRFfiles_dsgevar = 1; ifil2 = 1; -if strcmpi(type,'posterior') - h = waitbar(0,'Bayesian (posterior) IRFs...'); -elseif strcmpi(type,'gsa') - h = waitbar(0,'GSA (prior) IRFs...'); -else - h = waitbar(0,'Bayesian (prior) IRFs...'); -end % Create arrays if B <= MAX_nruns - stock_param = zeros(B, npar); + stock_param = zeros(B, npar); else - stock_param = zeros(MAX_nruns, npar); + stock_param = zeros(MAX_nruns, npar); end if B >= MAX_nirfs_dsge - stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge); + stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge); else - stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,B); + stock_irf_dsge = zeros(options_.irf,nvar,M_.exo_nbr,B); end if MAX_nirfs_dsgevar if B >= MAX_nirfs_dsgevar @@ -142,149 +136,99 @@ if MAX_nirfs_dsgevar end Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs); end + +%%%%%%%%% START the FIRST BLOCK of CODE EXECUTED in PARALLEL! %%%%%%%%% +% +% This portion of code is execute in parallel by PosteriorIRF_core1.m +% function. + b = 0; nosaddle = 0; -while b<=B + +localVars=[]; + +% Save the local variables. + +localVars.IRUN = IRUN; +localVars.irun = irun; +localVars.irun2=irun2; +localVars.nosaddle=nosaddle; + +% It is necessary to rename 'type' to avoid conflict with +% a native matlab funtion. + +localVars.typee=type; + +if strcmpi(type,'posterior'), +while b<=B b = b + 1; - irun = irun+1; - irun2 = irun2+1; - if ~strcmpi(type,'gsa') - deep = GetOneDraw(type); - else - deep = x(b,:); - end - stock_param(irun2,:) = deep; - set_parameters(deep); - [dr,info] = resol(oo_.steady_state,0); - if info(1) - nosaddle = nosaddle + 1; - b = b - 1; - irun = irun-1; - irun2 = irun2-1; - if info(1) == 1 - errordef = 'Static variables are not uniquely defined'; - elseif info(1) == 2 - errordef = 'Dll problem'; - elseif info(1) == 3 - errordef = 'No stable trajectory'; - elseif info(1) == 4 - errordef = 'Indeterminacy'; - elseif info(1) == 5 - errordef = 'Rank condition is not satisfied'; - end - disp(['PosteriorIRF :: Dynare is unable to solve the model (' errordef ')']) - continue - end - SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr); - SS = transpose(chol(SS)); - for i = 1:M_.exo_nbr - if SS(i,i) > 1e-13 - y=irf(dr,SS(M_.exo_names_orig_ord,i), options_.irf, options_.drop,options_.replic,options_.order); - if options_.relative_irf - y = 100*y/cs(i,i); - end - for j = 1:nvar - if max(y(IndxVariables(j),:)) - min(y(IndxVariables(j),:)) > 1e-12 - stock_irf_dsge(:,j,i,irun) = transpose(y(IndxVariables(j),:)); - end - end - end - end - if MAX_nirfs_dsgevar - IRUN = IRUN+1; - %tmp_dsgevar = zeros(options_.irf,nvobs*M_.exo_nbr); - [fval,cost_flag,info,PHI,SIGMAu,iXX] = DsgeVarLikelihood(deep',gend); - dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names)); - DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight)); - SIGMA_inv_upper_chol = chol(inv(SIGMAu*gend*(dsge_prior_weight+1))); - explosive_var = 1; - while explosive_var - % draw from the marginal posterior of SIGMA - SIGMAu_draw = rand_inverse_wishart(nvobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ... - SIGMA_inv_upper_chol); - % draw from the conditional posterior of PHI - PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,nvobs, PHI, ... - chol(SIGMAu_draw)', chol(iXX)'); - Companion_matrix(1:nvobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:)); - % Check for stationarity - explosive_var = any(abs(eig(Companion_matrix))>1.000000001); - end - % Get the mean -% $$$ if options_.noconstant - mu = zeros(1,nvobs); -% $$$ else -% $$$ AA = eye(nvobs); -% $$$ for lag=1:NumberOfLags -% $$$ AA = AA-PHI_draw((lag-1)*nvobs+1:lag*nvobs,:); -% $$$ end -% $$$ mu = transpose(AA\transpose(PHI_draw(end,:))); -% $$$ end -% Get rotation - if dsge_prior_weight > 0 - Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e); - A0 = Atheta(bayestopt_.mfys,:); - [OMEGAstar,SIGMAtr] = qr2(A0'); - end - SIGMAu_chol = chol(SIGMAu_draw)'; - SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar'; - PHIpower = eye(NumberOfLagsTimesNvobs); - irfs = zeros (options_.irf,nvobs*M_.exo_nbr); - tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA; - irfs(1,:) = tmp3(:)'; - for t = 2:options_.irf - PHIpower = Companion_matrix*PHIpower; - tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA; - irfs(t,:) = tmp3(:)'+kron(ones(1,M_.exo_nbr),mu); - end - tmp_dsgevar = kron(ones(options_.irf,1),mu); - for j = 1:(nvobs*M_.exo_nbr) - if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10 - tmp_dsgevar(:,j) = (irfs(:,j)); - end - end - if IRUN < MAX_nirfs_dsgevar - stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr); - else - stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr); - instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ... - int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];, - eval(['save ' instr]); - NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1; - IRUN =0; - stock_irf_dsgevar = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar); - end - end - if irun == MAX_nirfs_dsge | irun == B | b == B - if b == B - stock_irf_dsge = stock_irf_dsge(:,:,:,1:irun); - if MAX_nirfs_dsgevar & (b == B | IRUN == B) - stock_irf_bvardsge = stock_irf_bvardsge(:,:,:,1:IRUN); - instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ... - int2str(NumberOfIRFfiles_dsgevar) ' stock_irf_bvardsge;'];, - eval(['save ' instr]); - NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1; - irun = 0; - end - end - save([MhDirectoryName '/' M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat'],'stock_irf_dsge'); - NumberOfIRFfiles_dsge = NumberOfIRFfiles_dsge+1; - irun = 0; - end - if irun2 == MAX_nruns | b == B - if b == B - stock_param = stock_param(1:irun2,:); - end - stock = stock_param; - save([MhDirectoryName '/' M_.fname '_param_irf' int2str(ifil2) '.mat'],'stock'); - ifil2 = ifil2 + 1; - irun2 = 0; - end - waitbar(b/B,h); + x(b,:) = GetOneDraw(type); end +end + +if ~strcmpi(type,'prior'), + localVars.x=x; +end + +b=0; + +localVars.nvar=nvar; +localVars.IndxVariables=IndxVariables; +localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar; +localVars.MAX_nirfs_dsge=MAX_nirfs_dsge; +localVars.MAX_nruns=MAX_nruns; + +localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge; +localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar; +localVars.ifil2=ifil2; + + + +if isnumeric(options_.parallel),% | isunix, % for the moment exclude unix platform from parallel implementation + [fout] = PosteriorIRF_core1(localVars,1,B,0); +else + + [nCPU, totCPU, nBlockPerCPU] = distributeJobs(options_.parallel, 1, B); + for j=1:totCPU-1, + nfiles = ceil(nBlockPerCPU(j)/MAX_nirfs_dsge); + NumberOfIRFfiles_dsge(j+1) =NumberOfIRFfiles_dsge(j)+nfiles; + nfiles = ceil(nBlockPerCPU(j)/MAX_nirfs_dsgevar); + NumberOfIRFfiles_dsgevar(j+1) =NumberOfIRFfiles_dsgevar(j)+nfiles; + nfiles = ceil(nBlockPerCPU(j)/MAX_nruns); + ifil2(j+1) =ifil2(j)+nfiles; + end + localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge; + localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar; + localVars.ifil2=ifil2; + + globalVars = struct('M_',M_, ... + 'options_', options_, ... + 'bayestopt_', bayestopt_, ... + 'estim_params_', estim_params_, ... + 'oo_', oo_); + + % which files have to be copied to run remotely + NamFileInput(1,:) = {'',[M_.fname '_static.m']}; + NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']}; + if options_.steadystate_flag, + NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']}; + end + [fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'PosteriorIRF_core1', localVars, globalVars, options_.parallel_info); + +end + + +% END parallel code! +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + if nosaddle - disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))]) -end -close(h); + disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))]) +end + +% if isnumeric(options_.parallel) +% close(h); +% end ReshapeMatFiles('irf_dsge') if MAX_nirfs_dsgevar @@ -292,13 +236,13 @@ if MAX_nirfs_dsgevar end if strcmpi(type,'gsa') - return + return end -IRF_DSGEs = dir([MhDirectoryName '/' M_.fname '_IRF_DSGEs*.mat']); +IRF_DSGEs = dir([MhDirectoryName filesep M_.fname '_IRF_DSGEs*.mat']); NumberOfIRFfiles_dsge = length(IRF_DSGEs); -IRF_BVARDSGEs = dir([MhDirectoryName '/' M_.fname '_IRF_BVARDSGEs*.mat']); +IRF_BVARDSGEs = dir([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs*.mat']); NumberOfIRFfiles_dsgevar = length(IRF_BVARDSGEs); @@ -310,10 +254,10 @@ DistribIRF = zeros(options_.irf,9,nvar,M_.exo_nbr); HPDIRF = zeros(options_.irf,2,nvar,M_.exo_nbr); if options_.TeX - varlist_TeX = []; - for i=1:nvar - varlist_TeX = strvcat(varlist_TeX,M_.endo_names_tex(IndxVariables(i),:)); - end + varlist_TeX = []; + for i=1:nvar + varlist_TeX = strvcat(varlist_TeX,M_.endo_names_tex(IndxVariables(i),:)); + end end fprintf('MH: Posterior (dsge) IRFs...\n'); @@ -321,31 +265,32 @@ tit(M_.exo_names_orig_ord,:) = M_.exo_names; kdx = 0; for file = 1:NumberOfIRFfiles_dsge - load([MhDirectoryName '/' M_.fname '_IRF_DSGEs' int2str(file) '.mat']); - for i = 1:M_.exo_nbr - for j = 1:nvar - for k = 1:size(STOCK_IRF_DSGE,1) - kk = k+kdx; - [MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),... - DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(STOCK_IRF_DSGE(k,j,i,:)),0,options_.mh_conf_sig); - end + load([MhDirectoryName filesep M_.fname '_IRF_DSGEs' int2str(file) '.mat']); + for i = 1:M_.exo_nbr + for j = 1:nvar + for k = 1:size(STOCK_IRF_DSGE,1) + kk = k+kdx; + [MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),... + DistribIRF(kk,:,j,i)] = posterior_moments(squeeze(STOCK_IRF_DSGE(k,j,i,:)),0,options_.mh_conf_sig); end end - kdx = kdx + size(STOCK_IRF_DSGE,1); + end + kdx = kdx + size(STOCK_IRF_DSGE,1); + end clear STOCK_IRF_DSGE; for i = 1:M_.exo_nbr - for j = 1:nvar - name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))]; - eval(['oo_.PosteriorIRF.dsge.Mean.' name ' = MeanIRF(:,j,i);']); - eval(['oo_.PosteriorIRF.dsge.Median.' name ' = MedianIRF(:,j,i);']); - eval(['oo_.PosteriorIRF.dsge.Var.' name ' = VarIRF(:,j,i);']); - eval(['oo_.PosteriorIRF.dsge.Distribution.' name ' = DistribIRF(:,:,j,i);']); - eval(['oo_.PosteriorIRF.dsge.HPDinf.' name ' = HPDIRF(:,1,j,i);']); - eval(['oo_.PosteriorIRF.dsge.HPDsup.' name ' = HPDIRF(:,2,j,i);']); - end + for j = 1:nvar + name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))]; + eval(['oo_.PosteriorIRF.dsge.Mean.' name ' = MeanIRF(:,j,i);']); + eval(['oo_.PosteriorIRF.dsge.Median.' name ' = MedianIRF(:,j,i);']); + eval(['oo_.PosteriorIRF.dsge.Var.' name ' = VarIRF(:,j,i);']); + eval(['oo_.PosteriorIRF.dsge.Distribution.' name ' = DistribIRF(:,:,j,i);']); + eval(['oo_.PosteriorIRF.dsge.HPDinf.' name ' = HPDIRF(:,1,j,i);']); + eval(['oo_.PosteriorIRF.dsge.HPDsup.' name ' = HPDIRF(:,2,j,i);']); + end end @@ -359,7 +304,7 @@ if MAX_nirfs_dsgevar tit(M_.exo_names_orig_ord,:) = M_.exo_names; kdx = 0; for file = 1:NumberOfIRFfiles_dsgevar - load([MhDirectoryName '/' M_.fname '_IRF_BVARDSGEs' int2str(file) '.mat']); + load([MhDirectoryName filesep M_.fname '_IRF_BVARDSGEs' int2str(file) '.mat']); for i = 1:M_.exo_nbr for j = 1:nvar for k = 1:size(STOCK_IRF_BVARDSGE,1) @@ -386,114 +331,101 @@ if MAX_nirfs_dsgevar end end %% -%% Finally I build the plots. +%% Finally I build the plots. %% + + + +%%%%%%%%% START the SECOND BLOCK of CODE EXECUTED in PARALLEL! %%%%%%%%% +% +% This portion of code is execute in parallel by PosteriorIRF_core2.m +% function. + +% Save the local variables. +localVars=[]; + + Check=options_.TeX + if (Check) + localVars.varlist_TeX=varlist_TeX; + end + + +localVars.nvar=nvar; +localVars.MeanIRF=MeanIRF; +localVars.tit=tit; +localVars.nn=nn; +localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar; +localVars.HPDIRF=HPDIRF; +localVars.varlist=varlist; +localVars.MaxNumberOfPlotPerFigure=MaxNumberOfPlotPerFigure; + + +%%% The files .TeX are genereted in sequential way! + if options_.TeX - fidTeX = fopen([DirectoryName '/' M_.fname '_BayesianIRF.TeX'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by PosteriorIRF.m (Dynare).\n'); - fprintf(fidTeX,['%% ' datestr(now,0) '\n']); - fprintf(fidTeX,' \n'); - titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex; -end -%% -subplotnum = 0; -for i=1:M_.exo_nbr - NAMES = []; - if options_.TeX + fidTeX = fopen([DirectoryName filesep M_.fname '_BayesianIRF.TeX'],'w'); + fprintf(fidTeX,'%% TeX eps-loader file generated by PosteriorIRF.m (Dynare).\n'); + fprintf(fidTeX,['%% ' datestr(now,0) '\n']); + fprintf(fidTeX,' \n'); + titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex; + + for i=1:M_.exo_nbr + NAMES = []; TEXNAMES = []; - end - figunumber = 0; - for j=1:nvar - if max(abs(MeanIRF(:,j,i))) > 10^(-6) - subplotnum = subplotnum+1; - if options_.nograph - if subplotnum == 1 & options_.relative_irf - hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)],'Visible','off'); - elseif subplotnum == 1 & ~options_.relative_irf - hh = figure('Name',['Orthogonalized shock to ' tit(i,:)],'Visible','off'); - end - else - if subplotnum == 1 & options_.relative_irf - hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)]); - elseif subplotnum == 1 & ~options_.relative_irf - hh = figure('Name',['Orthogonalized shock to ' tit(i,:)]); - end - end - set(0,'CurrentFigure',hh) - subplot(nn,nn,subplotnum); - if ~MAX_nirfs_dsgevar - h1 = area(1:options_.irf,HPDIRF(:,2,j,i)); - set(h1,'FaceColor',[.9 .9 .9]); - set(h1,'BaseValue',min(HPDIRF(:,1,j,i))); - hold on - h2 = area(1:options_.irf,HPDIRF(:,1,j,i),'FaceColor',[1 1 1],'BaseValue',min(HPDIRF(:,1,j,i))); - set(h2,'FaceColor',[1 1 1]); - set(h2,'BaseValue',min(HPDIRF(:,1,j,i))); - plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3) - % plot([1 options_.irf],[0 0],'-r','linewidth',0.5); - box on - axis tight - xlim([1 options_.irf]); - hold off - else - h1 = area(1:options_.irf,HPDIRF(:,2,j,i)); - set(h1,'FaceColor',[.9 .9 .9]); - set(h1,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))])); - hold on; - h2 = area(1:options_.irf,HPDIRF(:,1,j,i)); - set(h2,'FaceColor',[1 1 1]); - set(h2,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))])); - plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3) - % plot([1 options_.irf],[0 0],'-r','linewidth',0.5); - plot(1:options_.irf,MeanIRFdsgevar(:,j,i),'--k','linewidth',2) - plot(1:options_.irf,HPDIRFdsgevar(:,1,j,i),'--k','linewidth',1) - plot(1:options_.irf,HPDIRFdsgevar(:,2,j,i),'--k','linewidth',1) - box on - axis tight - xlim([1 options_.irf]); - hold off - end + + for j=1:nvar + if max(abs(MeanIRF(:,j,i))) > 10^(-6) + name = deblank(varlist(j,:)); NAMES = strvcat(NAMES,name); - if options_.TeX - texname = deblank(varlist_TeX(j,:)); - TEXNAMES = strvcat(TEXNAMES,['$' texname '$']); + + texname = deblank(varlist_TeX(j,:)); + TEXNAMES = strvcat(TEXNAMES,['$' texname '$']); + end + + end + fprintf(fidTeX,'\\begin{figure}[H]\n'); + for jj = 1:size(TEXNAMES,1) + fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:))); end - title(name,'Interpreter','none') - end - if subplotnum == MaxNumberOfPlotPerFigure | (j == nvar & subplotnum> 0) - figunumber = figunumber+1; - set(hh,'visible','on') - eval(['print -depsc2 ' DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.eps']); - if ~exist('OCTAVE_VERSION') - eval(['print -dpdf ' DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)]); - saveas(hh,[DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.fig']); + fprintf(fidTeX,'\\centering \n'); + fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s}\n',M_.fname,deblank(tit(i,:))); + if options_.relative_irf + fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']); + else + fprintf(fidTeX,'\\caption{Bayesian IRF.}'); end - set(hh,'visible','off') - if options_.nograph, close(hh), end - if options_.TeX - fprintf(fidTeX,'\\begin{figure}[H]\n'); - for jj = 1:size(TEXNAMES,1) - fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:))); - end - fprintf(fidTeX,'\\centering \n'); - fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s}\n',M_.fname,deblank(tit(i,:))); - if options_.relative_irf - fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']); - else - fprintf(fidTeX,'\\caption{Bayesian IRF.}'); - end - fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s}\n',deblank(tit(i,:))); - fprintf(fidTeX,'\\end{figure}\n'); - fprintf(fidTeX,' \n'); - end - subplotnum = 0; - end - end% loop over selected endo_var -end% loop over exo_var -%% -if options_.TeX - fprintf(fidTeX,'%% End of TeX file.\n'); - fclose(fidTeX); + fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s}\n',deblank(tit(i,:))); + fprintf(fidTeX,'\\end{figure}\n'); + fprintf(fidTeX,' \n'); + end + + fprintf(fidTeX,'%% End of TeX file.\n'); + fclose(fidTeX); + end -fprintf('MH: Posterior IRFs, done!\n'); \ No newline at end of file + +% The others format in parallel by PosteriorIRF_core2! + + +if isnumeric(options_.parallel) || (M_.exo_nbr*ceil(size(varlist,1)/MaxNumberOfPlotPerFigure))<8,% | isunix, % for the moment exclude unix platform from parallel implementation + [fout] = PosteriorIRF_core2(localVars,1,M_.exo_nbr,0); +else + globalVars = struct('M_',M_, ... + 'options_', options_); + + [fout] = masterParallel(options_.parallel, 1, M_.exo_nbr,NamFileInput,'PosteriorIRF_core2', localVars, globalVars, options_.parallel_info); + +end + +% END parallel code! +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +fprintf('MH: Posterior IRFs, done!\n'); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + + diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m new file mode 100644 index 000000000..27637e636 --- /dev/null +++ b/matlab/PosteriorIRF_core1.m @@ -0,0 +1,278 @@ +function myoutput=PosteriorIRF_core1(myinputs,fpar,npar,whoiam, ThisMatlab) + +% Perfome in parallel a portion of PosteriorIRF +% +% INPUTS +% +% ... +% +% OUTPUTS +% ... +% +% SPECIAL REQUIREMENTS +% None + +% Copyright (C) 2006-2008 Dynare Team +% +% 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 . + + +global options_ estim_params_ oo_ M_ bayestopt_ + + + +if nargin<4, + whoiam=0; +end +struct2local(myinputs); + +MhDirectoryName = CheckPath('metropolis'); + +RemoteFlag = 0; + +if whoiam + waitbarString = ['Please wait... Bayesian (posterior) IRFs computing. (' int2str(fpar) 'of' int2str(npar) ')...']; + if Parallel(ThisMatlab).Local, + waitbarTitle=['Local ']; + else + waitbarTitle=[Parallel(ThisMatlab).PcName]; + RemoteFlag =1; + end + fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo); +else + if exist('OCTAVE_VERSION') + diary off; + printf('\n') + else + if strcmpi(typee,'posterior') + h = waitbar(0,'Bayesian (posterior) IRFs...'); + set(h,'Name','Bayesian (posterior) IRFs.'); + elseif strcmpi(typee,'gsa') + h = waitbar(0,'GSA (prior) IRFs...'); + else + h = waitbar(0,'Bayesian (prior) IRFs...'); + end + end + +end + +OutputFileName_bvardsge = {}; +OutputFileName_dsge = {}; +OutputFileName_param = {}; + + +fpar0=fpar; +fpar = fpar-1; + +if whoiam + ifil2=ifil2(whoiam); + NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge(whoiam); + NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar(whoiam); +end + while fpar<=npar + + fpar = fpar + 1; + irun = irun+1; + irun2 = irun2+1; + if strcmpi(typee,'prior') + + deep = GetOneDraw(typee); + + else + deep = x(fpar,:); + end + stock_param(irun2,:) = deep; + set_parameters(deep); + [dr,info] = resol(oo_.steady_state,0); + + + if info(1) + nosaddle = nosaddle + 1; + fpar = fpar - 1; + irun = irun-1; + irun2 = irun2-1; + if info(1) == 1 + errordef = 'Static variables are not uniquely defined'; + elseif info(1) == 2 + errordef = 'Dll problem'; + elseif info(1) == 3 + errordef = 'No stable trajectory'; + elseif info(1) == 4 + errordef = 'Indeterminacy'; + elseif info(1) == 5 + errordef = 'Rank condition is not satisfied'; + end + disp(['PosteriorIRF :: Dynare is unable to solve the model (' errordef ')']) + continue + end + SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr); + SS = transpose(chol(SS)); + for i = 1:M_.exo_nbr + if SS(i,i) > 1e-13 + y=irf(dr,SS(M_.exo_names_orig_ord,i), options_.irf, options_.drop,options_.replic,options_.order); + if options_.relative_irf + y = 100*y/cs(i,i); + end + for j = 1:nvar + if max(y(IndxVariables(j),:)) - min(y(IndxVariables(j),:)) > 1e-12 + stock_irf_dsge(:,j,i,irun) = transpose(y(IndxVariables(j),:)); + end + end + end + end + if MAX_nirfs_dsgevar + IRUN = IRUN+1; + %tmp_dsgevar = zeros(options_.irf,nvobs*M_.exo_nbr); + [fval,cost_flag,info,PHI,SIGMAu,iXX] = DsgeVarLikelihood(deep(fpar,:)',gend); + dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names)); + DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight)); + SIGMA_inv_upper_chol = chol(inv(SIGMAu*gend*(dsge_prior_weight+1))); + explosive_var = 1; + while explosive_var + % draw from the marginal posterior of SIGMA + SIGMAu_draw = rand_inverse_wishart(nvobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ... + SIGMA_inv_upper_chol); + % draw from the conditional posterior of PHI + PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,nvobs, PHI, ... + chol(SIGMAu_draw)', chol(iXX)'); + Companion_matrix(1:nvobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:)); + % Check for stationarity + explosive_var = any(abs(eig(Companion_matrix))>1.000000001); + end + % Get the mean +% $$$ if options_.noconstant + mu = zeros(1,nvobs); +% $$$ else +% $$$ AA = eye(nvobs); +% $$$ for lag=1:NumberOfLags +% $$$ AA = AA-PHI_draw((lag-1)*nvobs+1:lag*nvobs,:); +% $$$ end +% $$$ mu = transpose(AA\transpose(PHI_draw(end,:))); +% $$$ end + % Get rotation + if dsge_prior_weight > 0 + Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e); + A0 = Atheta(bayestopt_.mfys,:); + [OMEGAstar,SIGMAtr] = qr2(A0'); + end + SIGMAu_chol = chol(SIGMAu_draw)'; + SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar'; + PHIpower = eye(NumberOfLagsTimesNvobs); + irfs = zeros (options_.irf,nvobs*M_.exo_nbr); + tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA; + irfs(1,:) = tmp3(:)'; + for t = 2:options_.irf + PHIpower = Companion_matrix*PHIpower; + tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA; + irfs(t,:) = tmp3(:)'+kron(ones(1,M_.exo_nbr),mu); + end + tmp_dsgevar = kron(ones(options_.irf,1),mu); + for j = 1:(nvobs*M_.exo_nbr) + if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10 + tmp_dsgevar(:,j) = (irfs(:,j)); + end + end + if IRUN < MAX_nirfs_dsgevar + stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr); + else + stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr); + instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ... + int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];, + eval(['save ' instr]); + if RemoteFlag==1, + OutputFileName_bvardsge = [OutputFileName_bvardsge; {[MhDirectoryName filesep], [M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat']}]; + end + NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1; + IRUN =0; + stock_irf_dsgevar = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar); + end + end + if irun == MAX_nirfs_dsge | irun == npar | fpar == npar + if fpar == npar + stock_irf_dsge = stock_irf_dsge(:,:,:,1:irun); + if MAX_nirfs_dsgevar & (fpar == npar | IRUN == npar) + stock_irf_bvardsge = stock_irf_bvardsge(:,:,:,1:IRUN); + instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ... + int2str(NumberOfIRFfiles_dsgevar) ' stock_irf_bvardsge;'];, + eval(['save ' instr]); + NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1; + if RemoteFlag==1, + OutputFileName_bvardsge = [OutputFileName_bvardsge; {[MhDirectoryName filesep], [M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) '.mat']}]; + end + irun = 0; + end + end + save([MhDirectoryName '/' M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat'],'stock_irf_dsge'); + if RemoteFlag==1, + OutputFileName_dsge = [OutputFileName_dsge; {[MhDirectoryName filesep], [M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge) '.mat']}]; + end + NumberOfIRFfiles_dsge = NumberOfIRFfiles_dsge+1; + irun = 0; + end + if irun2 == MAX_nruns | fpar == npar + if fpar == npar + stock_param = stock_param(1:irun2,:); + end + stock = stock_param; + save([MhDirectoryName '/' M_.fname '_param_irf' int2str(ifil2) '.mat'],'stock'); + if RemoteFlag==1, + OutputFileName_param = [OutputFileName_param; {[MhDirectoryName filesep], [M_.fname '_param_irf' int2str(ifil2) '.mat']}]; + end + + ifil2 = ifil2 + 1; + irun2 = 0; + end + if exist('OCTAVE_VERSION') + printf(['Posterior IRF %3.f%% done\r'],(fpar/npar*100)); + elseif ~whoiam + waitbar(fpar/npar,h); + end + if mod(fpar,10)==0 & whoiam, + fprintf('Done! \n'); + waitbarString = [ 'Subdraw ' int2str(fpar) '/' int2str(npar) ' done.']; + fMessageStatus((fpar-fpar0)/(npar-fpar0),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo) + end + + end + +if whoiam==0 + if nosaddle + disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(npar+nosaddle))]) + end + if exist('h') + close(h); + end + if exist('OCTAVE_VERSION') + printf('\n'); + diary on; + end +end + + + + + +% Copy the rusults of computation on the call machine (specifically in the +% dyrectory on call machine that contain the model). + +myoutput.OutputFileName = [OutputFileName_dsge; + OutputFileName_param; + OutputFileName_bvardsge]; + +% MhDirectoryName=TempPath; + + + diff --git a/matlab/PosteriorIRF_core2.m b/matlab/PosteriorIRF_core2.m new file mode 100644 index 000000000..2a8846b63 --- /dev/null +++ b/matlab/PosteriorIRF_core2.m @@ -0,0 +1,121 @@ +function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab) + +global options_ M_ + +if nargin<4, + whoiam=0; +end +struct2local(myinputs); + +% To save the figures where the function is computed! + +DirectoryName = CheckPath('Output'); + +RemoteFlag = 0; +if whoiam, + waitbarString = ['Please wait... PosteriorIRF Plots (exog. shocks ' int2str(fpar) 'of' int2str(npar) ')...']; + if Parallel(ThisMatlab).Local, + waitbarTitle=['Local ']; + else + waitbarTitle=[Parallel(ThisMatlab).PcName]; + RemoteFlag = 1; + end + fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo); +end + +OutputFileName={}; + +subplotnum = 0; +for i=fpar:npar, + NAMES = []; + figunumber = 0; + + for j=1:nvar + if max(abs(MeanIRF(:,j,i))) > 10^(-6) + subplotnum = subplotnum+1; + if options_.nograph + if subplotnum == 1 & options_.relative_irf + hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)],'Visible','off'); + elseif subplotnum == 1 & ~options_.relative_irf + hh = figure('Name',['Orthogonalized shock to ' tit(i,:)],'Visible','off'); + end + else + if subplotnum == 1 & options_.relative_irf + hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)]); + elseif subplotnum == 1 & ~options_.relative_irf + hh = figure('Name',['Orthogonalized shock to ' tit(i,:)]); + end + end + set(0,'CurrentFigure',hh) + subplot(nn,nn,subplotnum); + if ~MAX_nirfs_dsgevar + h1 = area(1:options_.irf,HPDIRF(:,2,j,i)); + set(h1,'FaceColor',[.9 .9 .9]); + set(h1,'BaseValue',min(HPDIRF(:,1,j,i))); + hold on + h2 = area(1:options_.irf,HPDIRF(:,1,j,i),'FaceColor',[1 1 1],'BaseValue',min(HPDIRF(:,1,j,i))); + set(h2,'FaceColor',[1 1 1]); + set(h2,'BaseValue',min(HPDIRF(:,1,j,i))); + plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3) + % plot([1 options_.irf],[0 0],'-r','linewidth',0.5); + box on + axis tight + xlim([1 options_.irf]); + hold off + else + h1 = area(1:options_.irf,HPDIRF(:,2,j,i)); + set(h1,'FaceColor',[.9 .9 .9]); + set(h1,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))])); + hold on; + h2 = area(1:options_.irf,HPDIRF(:,1,j,i)); + set(h2,'FaceColor',[1 1 1]); + set(h2,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))])); + plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3) + % plot([1 options_.irf],[0 0],'-r','linewidth',0.5); + plot(1:options_.irf,MeanIRFdsgevar(:,j,i),'--k','linewidth',2) + plot(1:options_.irf,HPDIRFdsgevar(:,1,j,i),'--k','linewidth',1) + plot(1:options_.irf,HPDIRFdsgevar(:,2,j,i),'--k','linewidth',1) + box on + axis tight + xlim([1 options_.irf]); + hold off + end + name = deblank(varlist(j,:)); + NAMES = strvcat(NAMES,name); + title(name,'Interpreter','none') + end + + +% TempPath=DirectoryName; +% DirectoryNamePar='C:\dynare_calcs\ModelTest\ls2003\metropolis'; +% DirectoryName=DirectoryNamePar; + + + if subplotnum == MaxNumberOfPlotPerFigure | (j == nvar & subplotnum> 0) + figunumber = figunumber+1; + set(hh,'visible','on') + eval(['print -depsc2 ' DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.eps']); + if ~exist('OCTAVE_VERSION') + eval(['print -dpdf ' DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)]); + saveas(hh,[DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.fig']); + end + if RemoteFlag==1, + OutputFileName = [OutputFileName; {[DirectoryName,filesep], [M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber) '.*']}]; + end + set(hh,'visible','off') + if options_.nograph, close(hh), end + subplotnum = 0; + end + end% loop over selected endo_var + if whoiam, + fprintf('Done! \n'); + waitbarString = [ 'Exog. shocks ' int2str(i) '/' int2str(npar) ' done.']; + fMessageStatus((i-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo) + end +end% loop over exo_var + + +% DirectoryName=TempPath; + +myoutput.OutputFileName = OutputFileName; + diff --git a/matlab/ReshapeMatFiles.m b/matlab/ReshapeMatFiles.m index 4ecbbaeb6..950958e02 100644 --- a/matlab/ReshapeMatFiles.m +++ b/matlab/ReshapeMatFiles.m @@ -45,20 +45,20 @@ function ReshapeMatFiles(type, type2) global M_ options_ if nargin==1, - MhDirectoryName = [ CheckPath('metropolis') '/' ]; + MhDirectoryName = [ CheckPath('metropolis') filesep ]; else if strcmpi(type2,'posterior') - MhDirectoryName = [CheckPath('metropolis') '/' ]; + MhDirectoryName = [CheckPath('metropolis') filesep ]; elseif strcmpi(type2,'gsa') if options_.opt_gsa.morris==1, - MhDirectoryName = [CheckPath('GSA\SCREEN') '/' ]; + MhDirectoryName = [CheckPath('GSA\SCREEN') filesep ]; elseif options_.opt_gsa.morris==2, - MhDirectoryName = [CheckPath('GSA\IDENTIF') '/' ]; + MhDirectoryName = [CheckPath('GSA\IDENTIF') filesep ]; else - MhDirectoryName = [CheckPath('GSA') '/' ]; + MhDirectoryName = [CheckPath('GSA') filesep ]; end else - MhDirectoryName = [CheckPath('prior') '/' ]; + MhDirectoryName = [CheckPath('prior') filesep ]; end end switch type diff --git a/matlab/Tracing.m b/matlab/Tracing.m new file mode 100644 index 000000000..e3814552f --- /dev/null +++ b/matlab/Tracing.m @@ -0,0 +1,40 @@ +function [] = Tracing() + +% DESCRIPTION +% This function is used to test the correct execution of a matlab section +% on remote machine. +% +% If no error happen the function simply create a file. + +% INPUTS +% ... +% +% OUTPUTS +% ... +% +% SPECIAL REQUIREMENTS +% none + +% Copyright (C) 2005-2009 Dynare Team +% +% 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 . + + + fid = fopen('Tracing.txt','w+'); + fclose (fid); + +exit + diff --git a/matlab/closeSlave.m b/matlab/closeSlave.m new file mode 100644 index 000000000..8615606b6 --- /dev/null +++ b/matlab/closeSlave.m @@ -0,0 +1,14 @@ +function closeSlave(Parallel), +% In a parallelc context, this utility closes all remote matlab instances +% called by masteParallelMan (which leaves open remote matlab instances) + +delete( 'slaveParallel_input*.mat'); +for indPC=1:length(Parallel), + if Parallel(indPC).Local==0, + if isunix || (~matlab_ver_less_than('7.4') && ismac), + system(['ssh ',Parallel(indPC).user,'@',Parallel(indPC).PcName,' rm -fr ',Parallel(indPC).RemoteFolder,'/slaveParallel_input*.mat']); + else + mydelete('slaveParallel_input*.mat',['\\',Parallel(indPC).PcName,'\',Parallel(indPC).RemoteDrive,'$\',Parallel(indPC).RemoteFolder,'\']); + end + end +end diff --git a/matlab/distributeJobs.m b/matlab/distributeJobs.m new file mode 100644 index 000000000..492f47cd8 --- /dev/null +++ b/matlab/distributeJobs.m @@ -0,0 +1,37 @@ +function [nCPU, totCPU, nBlockPerCPU] = distributeJobs(Parallel, fBlock, nBlock) + +% Determine the total number of available CPUs, and the number of threads to run on each CPU + +% Copyright (C) 2009 Dynare Team +% +% 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 . + +totCPU=0; +for j=1:length(Parallel), + nCPU(j)=length(Parallel(j).NumCPU); + totCPU=totCPU+nCPU(j); +end + +nCPU=cumsum(nCPU); +offset0 = fBlock-1; +if (nBlock-offset0)>totCPU, + diff = mod((nBlock-offset0),totCPU); + nBlockPerCPU(1:diff) = ceil((nBlock-offset0)/totCPU); + nBlockPerCPU(diff+1:totCPU) = floor((nBlock-offset0)/totCPU); +else + nBlockPerCPU(1:nBlock-offset0)=1; + totCPU = nBlock-offset0; +end diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 14fb4c590..bff8b85c3 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -191,6 +191,7 @@ options_.mh_autocorrelation_function_size = 30; options_.plot_priors = 1; options_.cova_compute = 1; options_.parallel = 0; +options_.parallel_info.leaveSlaveOpen = 0; options_.number_of_grid_points_for_kde = 2^9; quarter = 1; years = [1 2 3 4 5 10 20 30 40 50]; diff --git a/matlab/independent_metropolis_hastings.m b/matlab/independent_metropolis_hastings.m index 05fa09579..4aa6916ab 100644 --- a/matlab/independent_metropolis_hastings.m +++ b/matlab/independent_metropolis_hastings.m @@ -95,7 +95,7 @@ else % from where to get back results % NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'}; - [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'independent_metropolis_hastings_core', localVars, globalVars); + [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'independent_metropolis_hastings_core', localVars, globalVars, options_.parallel_info); for j=1:totCPU, offset = sum(nBlockPerCPU(1:j-1))+fblck-1; record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j))); diff --git a/matlab/masterParallel.m b/matlab/masterParallel.m index 0709e98fb..9b4e698ef 100644 --- a/matlab/masterParallel.m +++ b/matlab/masterParallel.m @@ -1,4 +1,4 @@ -function [fOutVar,nBlockPerCPU, totCPU] = masterParallel(Parallel,fBlock,nBlock,NamFileInput,fname,fInputVar,fGlobalVar) +function [fOutVar,nBlockPerCPU, totCPU] = masterParallel(Parallel,fBlock,nBlock,NamFileInput,fname,fInputVar,fGlobalVar,Parallel_info) % Top-level function called on the master computer when parallelizing a task. % % The number of parallelized threads will be equal to (nBlock-fBlock+1). @@ -46,9 +46,13 @@ function [fOutVar,nBlockPerCPU, totCPU] = masterParallel(Parallel,fBlock,nBlock, % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -totCPU=0; - % Determine my hostname and my working directory +if ~isempty(Parallel_info) + if Parallel_info.leaveSlaveOpen, + [fOutVar,nBlockPerCPU, totCPU] = masterParallelMan(Parallel,fBlock,nBlock,NamFileInput,fname,fInputVar,fGlobalVar); + return + end +end DyMo=pwd; fInputVar.DyMo=DyMo; if isunix || (~matlab_ver_less_than('7.4') && ismac) , @@ -69,21 +73,24 @@ end save([fname,'_input.mat'],'Parallel','-append') % Determine the total number of available CPUs, and the number of threads to run on each CPU -for j=1:length(Parallel), - nCPU(j)=length(Parallel(j).NumCPU); - totCPU=totCPU+nCPU(j); -end - -nCPU=cumsum(nCPU); +[nCPU, totCPU, nBlockPerCPU] = distributeJobs(Parallel, fBlock, nBlock); offset0 = fBlock-1; -if (nBlock-offset0)>totCPU, - diff = mod((nBlock-offset0),totCPU); - nBlockPerCPU(1:diff) = ceil((nBlock-offset0)/totCPU); - nBlockPerCPU(diff+1:totCPU) = floor((nBlock-offset0)/totCPU); -else - nBlockPerCPU(1:nBlock-offset0)=1; - totCPU = nBlock-offset0; -end + +% for j=1:length(Parallel), +% nCPU(j)=length(Parallel(j).NumCPU); +% totCPU=totCPU+nCPU(j); +% end +% +% nCPU=cumsum(nCPU); +% offset0 = fBlock-1; +% if (nBlock-offset0)>totCPU, +% diff = mod((nBlock-offset0),totCPU); +% nBlockPerCPU(1:diff) = ceil((nBlock-offset0)/totCPU); +% nBlockPerCPU(diff+1:totCPU) = floor((nBlock-offset0)/totCPU); +% else +% nBlockPerCPU(1:nBlock-offset0)=1; +% totCPU = nBlock-offset0; +% end % Clean up remnants of previous runs mydelete(['comp_status_',fname,'*.mat']) @@ -214,6 +221,7 @@ t00=cputime; hh=NaN(1,nBlock); if exist('OCTAVE_VERSION'), diary off; + printf('\n'); else hfigstatus = figure('name',['Parallel ',fname],... 'MenuBar', 'none', ... @@ -241,7 +249,7 @@ while (1) load(stax(j).name) pcerdone(j) = prtfrc; if exist('OCTAVE_VERSION'), - statusString = [statusString, waitbarString, ', %3.f%% done! ']; + statusString = [statusString, int2str(j), ' %3.f%% done! ']; else status_String{j} = waitbarString; status_Title{j} = waitbarTitle; @@ -257,6 +265,7 @@ while (1) else figure(hfigstatus), for j=1:length(stax), + try axes(hstatus(idCPU(j))), hpat = findobj(hstatus(idCPU(j)),'Type','patch'); if ~isempty(hpat), @@ -265,6 +274,9 @@ while (1) patch([0 0 pcerdone(j) pcerdone(j)],[0 1 1 0],'r','EdgeColor','r') end title([status_Title{j},' - ',status_String{j}]); + catch + + end end end if isempty(dir(['P_',fname,'_*End.txt'])) diff --git a/matlab/masterParallelMan.m b/matlab/masterParallelMan.m index 855c9c47e..d4bd67b23 100644 --- a/matlab/masterParallelMan.m +++ b/matlab/masterParallelMan.m @@ -52,6 +52,7 @@ persistent initialize if isempty(initialize), mydelete(['P_slave_*End.txt']); mydelete(['slaveParallel_input*.mat']); + initialize = 0; pause(1), end @@ -101,6 +102,7 @@ mydelete(['P_',fname,'*End.txt']); % Create a shell script containing the commands to launch the required tasks on the slaves fid = fopen('ConcurrentCommand1.bat','w+'); for j=1:totCPU, + command1 = ' '; indPC=min(find(nCPU>=j)); @@ -156,7 +158,7 @@ for j=1:totCPU, command1=['start /B psexec -W ',DyMo, ' -a ',int2str(Parallel(indPC).NumCPU(j-nCPU0)),' -low matlab -nosplash -nodesktop -minimize -r slaveParallel(',int2str(j),',',int2str(indPC),')']; end end - else + elseif Parallel(indPC).Local==0, if isunix || (~matlab_ver_less_than('7.4') && ismac), % [tempo, RemoteName]=system(['ssh ',Parallel(indPC).user,'@',Parallel(indPC).PcName,' "ifconfig | grep \''inet addr:\''| grep -v \''127.0.0.1\'' | cut -d: -f2 | awk \''{ print $1}\''"']); [tempo, RemoteName]=system(['ssh ',Parallel(indPC).user,'@',Parallel(indPC).PcName,' "hostname --fqdn"']); @@ -263,6 +265,7 @@ t00=cputime; hh=NaN(1,nBlock); if exist('OCTAVE_VERSION'), diary off; + printf('\n'); else hfigstatus = figure('name',['Parallel ',fname],... 'MenuBar', 'none', ... @@ -290,7 +293,7 @@ while (1) load(stax(j).name) pcerdone(j) = prtfrc; if exist('OCTAVE_VERSION'), - statusString = [statusString, waitbarString, ', %3.f%% done! ']; + statusString = [statusString, int2str(j), ' %3.f%% done! ']; else status_String{j} = waitbarString; status_Title{j} = waitbarTitle; diff --git a/matlab/pm3.m b/matlab/pm3.m index d843c2587..7c682e4db 100644 --- a/matlab/pm3.m +++ b/matlab/pm3.m @@ -77,77 +77,80 @@ for i = 1:nvar eval(['oo_.' name3 '.HPDsup.' name ' = HPD(2,:,i);']); end %% -%% Finally I build the plots. +%% Finally I build the plots. %% + + +% %%%%%%%%% PARALLEL BLOCK % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% %%% The file .TeX! are not saved in parallel. + + +subplotnum = 0; + if options_.TeX fidTeX = fopen([M_.dname '/Output/' M_.fname '_' name3 '.TeX'],'w'); fprintf(fidTeX,'%% TeX eps-loader file generated by Dynare.\n'); fprintf(fidTeX,['%% ' datestr(now,0) '\n']); fprintf(fidTeX,' \n'); -end -%% -figunumber = 0; -subplotnum = 0; -hh = figure('Name',[tit1 ' ' int2str(figunumber+1)]); -for i=1:nvar - NAMES = []; - if options_.TeX - TEXNAMES = []; - end - if max(abs(Mean(:,i))) > 10^(-6) - subplotnum = subplotnum+1; - set(0,'CurrentFigure',hh) - subplot(nn,nn,subplotnum); - plot([1 n2],[0 0],'-r','linewidth',0.5); - hold on - for k = 1:9 - plot(1:n2,squeeze(Distrib(k,:,i)),'-g','linewidth',0.5) - end - plot(1:n2,Mean(:,i),'-k','linewidth',1) - xlim([1 n2]); - hold off - name = deblank(varlist(i,:)); - NAMES = strvcat(NAMES,name); - if options_.TeX + + for i=1:nvar + NAMES = []; + TEXNAMES = []; + if max(abs(Mean(:,i))) > 10^(-6) + subplotnum = subplotnum+1; + name = deblank(varlist(i,:)); + NAMES = strvcat(NAMES,name); texname = deblank(varlist_TeX(i,:)); TEXNAMES = strvcat(TEXNAMES,['$' texname '$']); end - title(name,'Interpreter','none') - end - if subplotnum == MaxNumberOfPlotsPerFigure | i == nvar - eval(['print -depsc2 ' M_.dname '/Output/' M_.fname '_' name3 '_' deblank(tit3(i,:)) '.eps' ]); - if ~exist('OCTAVE_VERSION') - eval(['print -dpdf ' M_.dname '/Output/' M_.fname '_' name3 '_' deblank(tit3(i,:))]); - saveas(hh,[M_.dname '/Output/' M_.fname '_' name3 '_' deblank(tit3(i,:)) '.fig']); - end - if options_.nograph, close(hh), end - if options_.TeX + if subplotnum == MaxNumberOfPlotsPerFigure | i == nvar fprintf(fidTeX,'\\begin{figure}[H]\n'); for jj = 1:size(TEXNAMES,1) fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:))); - end + end fprintf(fidTeX,'\\centering \n'); fprintf(fidTeX,['\\includegraphics[scale=0.5]{%s_' name3 '_%s}\n'],M_.fname,deblank(tit3(i,:))); - if options_.relative_irf - fprintf(fidTeX,['\\caption{' caption '.}']); - else - fprintf(fidTeX,['\\caption{' caption '.}']); - end fprintf(fidTeX,'\\label{Fig:%s:%s}\n',name3,deblank(tit3(i,:))); fprintf(fidTeX,'\\end{figure}\n'); fprintf(fidTeX,' \n'); - end - subplotnum = 0; - figunumber = figunumber+1; - if (i ~= nvar) - hh = figure('Name',[name3 ' ' int2str(figunumber+1)]); + subplotnum = 0; end end - -end -%% -if options_.TeX fprintf(fidTeX,'%% End of TeX file.\n'); fclose(fidTeX); end + +% Store the variable mandatory for local/remote parallel computing. + +localVars=[]; + +localVars.tit1=tit1; +localVars.nn=nn; +localVars.n2=n2; +localVars.Distrib=Distrib; +localVars.varlist=varlist; +localVars.MaxNumberOfPlotsPerFigure=MaxNumberOfPlotsPerFigure; +localVars.name3=name3; +localVars.tit3=tit3; +localVars.Mean=Mean; + + if isnumeric(options_.parallel) || ceil(size(varlist,1)/MaxNumberOfPlotsPerFigure)<4, + fout = pm3_core(localVars,1,nvar,0); + + else + globalVars = struct('M_',M_, ... + 'options_', options_, ... + 'oo_', oo_); + [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, 1, nvar, [],'pm3_core', localVars,globalVars, options_.parallel_info); + end + +%%%%%%%%% END PARALLEL BLOCK %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + fprintf(['MH: ' tit1 ', done!\n']); + + + + + + diff --git a/matlab/pm3_core.m b/matlab/pm3_core.m new file mode 100644 index 000000000..4ab6b40d0 --- /dev/null +++ b/matlab/pm3_core.m @@ -0,0 +1,81 @@ +function myoutput=pm3_core(myinputs,fpar,nvar,whoiam, ThisMatlab) + + +if nargin<4, + whoiam=0; +end +struct2local(myinputs); + + +global options_ M_ oo_ + + +if whoiam + waitbarString = ['Parallel plots pm3 ...']; + if Parallel(ThisMatlab).Local, + waitbarTitle=['Local ']; + else + waitbarTitle=[Parallel(ThisMatlab).PcName]; + end + fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo); + end + + + +figunumber = 0; +subplotnum = 0; +hh = figure('Name',[tit1 ' ' int2str(figunumber+1)]); +RemoteFlag = 0; +if whoiam, + if Parallel(ThisMatlab).Local ==0 + RemoteFlag=1; + end +end + +OutputFileName = {}; + +for i=fpar:nvar + NAMES = []; + if max(abs(Mean(:,i))) > 10^(-6) + subplotnum = subplotnum+1; + set(0,'CurrentFigure',hh) + subplot(nn,nn,subplotnum); + plot([1 n2],[0 0],'-r','linewidth',0.5); + hold on + for k = 1:9 + plot(1:n2,squeeze(Distrib(k,:,i)),'-g','linewidth',0.5) + end + plot(1:n2,Mean(:,i),'-k','linewidth',1) + xlim([1 n2]); + hold off + name = deblank(varlist(i,:)); + NAMES = strvcat(NAMES,name); + title(name,'Interpreter','none') + end + if subplotnum == MaxNumberOfPlotsPerFigure | i == nvar + eval(['print -depsc2 ' M_.dname '/Output/' M_.fname '_' name3 '_' deblank(tit3(i,:)) '.eps' ]); + if ~exist('OCTAVE_VERSION') + eval(['print -dpdf ' M_.dname '/Output/' M_.fname '_' name3 '_' deblank(tit3(i,:))]); + saveas(hh,[M_.dname '/Output/' M_.fname '_' name3 '_' deblank(tit3(i,:)) '.fig']); + end + if RemoteFlag==1, + OutputFileName = [OutputFileName; {[M_.dname, filesep, 'Output',filesep], [M_.fname '_' name3 '_' deblank(tit3(i,:)) '.*']}]; + end + if options_.nograph, close(hh), end + subplotnum = 0; + figunumber = figunumber+1; + if (i ~= nvar) + hh = figure('Name',[name3 ' ' int2str(figunumber+1)]); + end + end + + if whoiam, + waitbarString = [ 'Variable ' int2str(i) '/' int2str(nvar) ' done.']; + fMessageStatus((i-fpar+1)/(nvar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo) + end + + +end + + +myoutput.OutputFileName=OutputFileName; diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m index affbe3de6..7c60f7d17 100644 --- a/matlab/prior_posterior_statistics.m +++ b/matlab/prior_posterior_statistics.m @@ -1,4 +1,5 @@ function prior_posterior_statistics(type,Y,gend,data_index,missing_value) + % function PosteriorFilterSmootherAndForecast(Y,gend, type) % Computes posterior filter smoother and forecasts % @@ -7,13 +8,13 @@ function prior_posterior_statistics(type,Y,gend,data_index,missing_value) % prior % gsa % Y: data -% gend: number of observations +% gend: number of observations % data_index [cell] 1*smpl cell of column vectors of indices. % missing_value 1 if missing values, 0 otherwise -% +% % OUTPUTS % none -% +% % SPECIAL REQUIREMENTS % none @@ -36,6 +37,8 @@ function prior_posterior_statistics(type,Y,gend,data_index,missing_value) global options_ estim_params_ oo_ M_ bayestopt_ +localVars=[]; + nvx = estim_params_.nvx; nvn = estim_params_.nvn; ncx = estim_params_.ncx; @@ -61,8 +64,8 @@ maxlag = M_.maximum_endo_lag; DirectoryName = CheckPath('metropolis'); load([ DirectoryName '/' M_.fname '_mh_history']) FirstMhFile = record.KeepedDraws.FirstMhFile; -FirstLine = record.KeepedDraws.FirstLine; -TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles; +FirstLine = record.KeepedDraws.FirstLine; +TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles; TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); clear record; @@ -82,14 +85,14 @@ MAX_ninno = min(B,ceil(MaxNumberOfBytes/(exo_nbr*gend)/8)); MAX_nerro = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)*gend)/8)); if naK MAX_naK = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)* ... - length(options_.filter_step_ahead)*gend)/8)); + length(options_.filter_step_ahead)*gend)/8)); end if horizon MAX_nforc1 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8)); MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/ ... - 8)); + 8)); IdObs = bayestopt_.mfys; - + end MAX_momentsno = min(B,ceil(MaxNumberOfBytes/(get_moments_size(options_)*8))); %% @@ -108,11 +111,6 @@ end irun = ones(7,1); ifil = zeros(7,1); -if exist('OCTAVE_VERSION') - diary off; -else - h = waitbar(0,'Taking subdraws...'); -end stock_param = zeros(MAX_nruns, npar); stock_logpo = zeros(MAX_nruns,1); @@ -128,7 +126,7 @@ end if options_.filter_step_ahead stock_filter_step_ahead = zeros(naK,endo_nbr,gend+ ... - options_.filter_step_ahead(end),MAX_naK); + options_.filter_step_ahead(end),MAX_naK); run_smoother = 1; end if options_.forecast @@ -139,153 +137,128 @@ end %if moments_varendo % stock_moments = cell(MAX_momentsno,1); %end -for b=1:B - [deep, logpo] = GetOneDraw(type); - set_all_parameters(deep); - dr = resol(oo_.steady_state,0); - %if moments_varendo - % stock_moments{irun(8)} = compute_model_moments(dr,M_,options_); - %end - if run_smoother - [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK] = ... - DsgeSmoother(deep,gend,Y,data_index,missing_value); - if options_.loglinear - stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ... - repmat(log(dr.ys(dr.order_var)),1,gend); - stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ... - repmat(log(dr.ys(dr.order_var)),1,gend); - else - stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ... - repmat(dr.ys(dr.order_var),1,gend); - stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ... - repmat(dr.ys(dr.order_var),1,gend); - end - stock_innov(:,:,irun(2)) = etahat; - if nvn - stock_error(:,:,irun(3)) = epsilonhat; - end - if naK - stock_filter_step_ahead(:,dr.order_var,:,irun(4)) = aK(options_.filter_step_ahead,1:endo_nbr,:); - end - if horizon - yyyy = alphahat(iendo,i_last_obs); - yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr)); - if options_.prefilter == 1 - yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ... - horizon+maxlag,1); - end - yf(:,IdObs) = yf(:,IdObs)+(gend+[1-maxlag:horizon]')*trend_coeff'; - if options_.loglinear == 1 - yf = yf+repmat(log(SteadyState'),horizon+maxlag,1); - else - yf = yf+repmat(SteadyState',horizon+maxlag,1); - end - yf1 = forcst2(yyyy,horizon,dr,1); - if options_.prefilter == 1 - yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ... - repmat(bayestopt_.mean_varobs',[horizon+maxlag,1,1]); - end - yf1(:,IdObs,:) = yf1(:,IdObs,:)+repmat((gend+[1-maxlag:horizon]')* ... - trend_coeff',[1,1,1]); - if options_.loglinear == 1 - yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]); - else - yf1 = yf1 + repmat(SteadyState',[horizon+maxlag,1,1]); - end +% Store the variable mandatory for local/remote parallel computing. - stock_forcst_mean(:,:,irun(6)) = yf'; - stock_forcst_point(:,:,irun(7)) = yf1'; - end - - end - stock_param(irun(5),:) = deep; - stock_logpo(irun(5),1) = logpo; - stock_ys(irun(5),:) = SteadyState'; +localVars.typee=type; +localVars.run_smoother=run_smoother; +localVars.gend=gend; +localVars.Y=Y; +localVars.data_index=data_index; +localVars.missing_value=missing_value; +localVars.options_.varobs=options_.varobs; +localVars.irun=irun; +localVars.endo_nbr=endo_nbr; +localVars.nvn=nvn; +localVars.naK=naK; +localVars.horizon=horizon; +localVars.iendo=iendo; +if horizon + localVars.i_last_obs=i_last_obs; + localVars.IdObs=IdObs; +end +localVars.exo_nbr=exo_nbr; +localVars.maxlag=maxlag; +localVars.MAX_nsmoo=MAX_nsmoo; +localVars.MAX_ninno=MAX_ninno; +localVars.MAX_nerro = MAX_nerro; +if naK + localVars.MAX_naK=MAX_naK; +end +localVars.MAX_nruns=MAX_nruns; +if horizon + localVars.MAX_nforc1=MAX_nforc1; + localVars.MAX_nforc2=MAX_nforc2; +end +localVars.MAX_momentsno = MAX_momentsno; +localVars.ifil=ifil; - irun = irun + ones(7,1); - if irun(1) > MAX_nsmoo || b == B - stock = stock_smooth(:,:,1:irun(1)-1); - ifil(1) = ifil(1) + 1; - save([DirectoryName '/' M_.fname '_smooth' int2str(ifil(1)) '.mat'],'stock'); - stock = stock_update(:,:,1:irun(1)-1); - save([DirectoryName '/' M_.fname '_update' int2str(ifil(1)) '.mat'],'stock'); - irun(1) = 1; - end - - if irun(2) > MAX_ninno || b == B - stock = stock_innov(:,:,1:irun(2)-1); - ifil(2) = ifil(2) + 1; - save([DirectoryName '/' M_.fname '_inno' int2str(ifil(2)) '.mat'],'stock'); - irun(2) = 1; - end - - if nvn && (irun(3) > MAX_nerro || b == B) - stock = stock_error(:,:,1:irun(3)-1); - ifil(3) = ifil(3) + 1; - save([DirectoryName '/' M_.fname '_error' int2str(ifil(3)) '.mat'],'stock'); - irun(3) = 1; - end - - if naK && (irun(4) > MAX_naK || b == B) - stock = stock_filter_step_ahead(:,:,:,1:irun(4)-1); - ifil(4) = ifil(4) + 1; - save([DirectoryName '/' M_.fname '_filter_step_ahead' int2str(ifil(4)) '.mat'],'stock'); - irun(4) = 1; - end - - if irun(5) > MAX_nruns || b == B - stock = stock_param(1:irun(5)-1,:); - ifil(5) = ifil(5) + 1; - save([DirectoryName '/' M_.fname '_param' int2str(ifil(5)) '.mat'],'stock','stock_logpo','stock_ys'); - irun(5) = 1; - end - - if horizon && (irun(6) > MAX_nforc1 || b == B) - stock = stock_forcst_mean(:,:,1:irun(6)-1); - ifil(6) = ifil(6) + 1; - save([DirectoryName '/' M_.fname '_forc_mean' int2str(ifil(6)) '.mat'],'stock'); - irun(6) = 1; - end - - if horizon && (irun(7) > MAX_nforc2 || b == B) - stock = stock_forcst_point(:,:,1:irun(7)-1); - ifil(7) = ifil(7) + 1; - save([DirectoryName '/' M_.fname '_forc_point' int2str(ifil(7)) '.mat'],'stock'); - irun(7) = 1; - end - - % if moments_varendo && (irun(8) > MAX_momentsno || b == B) - % stock = stock_moments(1:irun(8)-1); - % ifil(8) = ifil(8) + 1; - % save([DirectoryName '/' M_.fname '_moments' int2str(ifil(8)) '.mat'],'stock'); - % irun(8) = 1; - % end - - if exist('OCTAVE_VERSION') - printf('Taking subdraws: %3.f%% done\r', b/B); - else - waitbar(b/B,h); +if strcmpi(type,'posterior'), + b=0; + while b<=B + b = b + 1; + [x(b,:), logpost(b,1)] = GetOneDraw(type); end end +if ~strcmpi(type,'prior'), + localVars.x=x; + localVars.logpost=logpost; +end + +b=0; + + + +if isnumeric(options_.parallel),% | isunix, % for the moment exclude unix platform from parallel implementation + [fout] = prior_posterior_statistics_core(localVars,1,B,0); +else + [nCPU, totCPU, nBlockPerCPU] = distributeJobs(options_.parallel, 1, B); + for j=1:totCPU-1, + nfiles = ceil(nBlockPerCPU(j)/MAX_nsmoo); + ifil(1,j+1) =ifil(1,j)+nfiles; + nfiles = ceil(nBlockPerCPU(j)/MAX_ninno); + ifil(2,j+1) =ifil(2,j)+nfiles; + nfiles = ceil(nBlockPerCPU(j)/MAX_nerro); + ifil(3,j+1) =ifil(3,j)+nfiles; + nfiles = ceil(nBlockPerCPU(j)/MAX_naK); + ifil(4,j+1) =ifil(4,j)+nfiles; + nfiles = ceil(nBlockPerCPU(j)/MAX_nruns); + ifil(5,j+1) =ifil(5,j)+nfiles; + nfiles = ceil(nBlockPerCPU(j)/MAX_nforc1); + ifil(6,j+1) =ifil(6,j)+nfiles; + nfiles = ceil(nBlockPerCPU(j)/MAX_nforc2); + ifil(7,j+1) =ifil(7,j)+nfiles; +% nfiles = ceil(nBlockPerCPU(j)/MAX_momentsno); +% ifil(8,j+1) =ifil(8,j)+nfiles; + end + localVars.ifil = ifil; + globalVars = struct('M_',M_, ... + 'options_', options_, ... + 'bayestopt_', bayestopt_, ... + 'estim_params_', estim_params_, ... + 'oo_', oo_); + + % which files have to be copied to run remotely + NamFileInput(1,:) = {'',[M_.fname '_static.m']}; + NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']}; + if options_.steadystate_flag, + NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']}; + end + [fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'prior_posterior_statistics_core', localVars,globalVars, options_.parallel_info); + +end +ifil = fout(end).ifil; + + + if exist('OCTAVE_VERSION') printf('\n'); diary on; else - close(h) + if exist('h') + close(h) + end + end stock_gend=gend; stock_data=Y; save([DirectoryName '/' M_.fname '_data.mat'],'stock_gend','stock_data'); +if ~isnumeric(options_.parallel), + leaveSlaveOpen = options_.parallel_info.leaveSlaveOpen; + if options_.parallel_info.leaveSlaveOpen == 0, + options_.parallel_info.leaveSlaveOpen = 1; % force locally to leave open remote matlab sessions (repeated pm3 calls) + end +end + if options_.smoother pm3(endo_nbr,gend,ifil(1),B,'Smoothed variables',... - '',varlist,'tit_tex',M_.endo_names,... + '',M_.endo_names(1:M_.orig_endo_nbr, :),'tit_tex',M_.endo_names,... varlist,'SmoothedVariables',[M_.dname '/metropolis'],'_smooth'); pm3(exo_nbr,gend,ifil(2),B,'Smoothed shocks',... '',M_.exo_names,'tit_tex',M_.exo_names,... @@ -315,4 +288,12 @@ if options_.forecast pm3(endo_nbr,horizon+maxlag,ifil(6),B,'Forecasted variables (point)',... '',varlist,'tit_tex',M_.endo_names,... varlist,'PointForecast',[M_.dname '/metropolis'],'_forc_point'); -end \ No newline at end of file +end + + +if ~isnumeric(options_.parallel), + options_.parallel_info.leaveSlaveOpen = leaveSlaveOpen; + if leaveSlaveOpen == 0, + closeSlave(options_.parallel), + end +end diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m new file mode 100644 index 000000000..8be66c9df --- /dev/null +++ b/matlab/prior_posterior_statistics_core.m @@ -0,0 +1,238 @@ +function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMatlab) + + +if nargin<4, + whoiam=0; +end +struct2local(myinputs); + + +global options_ oo_ M_ bayestopt_ estim_params_ + +DirectoryName = CheckPath('metropolis'); + +RemoteFlag = 0; +if whoiam + ifil=ifil(:,whoiam); + waitbarString = ['Please wait... Bayesian (posterior) subdraws (' int2str(fpar) 'of' int2str(B) ')...']; + if Parallel(ThisMatlab).Local, + waitbarTitle=['Local ']; + else + waitbarTitle=[Parallel(ThisMatlab).PcName]; + RemoteFlag = 1; + end + fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo); +else + if exist('OCTAVE_VERSION') + diary off; + else + h = waitbar(0,'Taking subdraws...'); + end + +end + +if RemoteFlag==1, +OutputFileName_smooth = {}; +OutputFileName_update = {}; +OutputFileName_inno = {}; +OutputFileName_error = {}; +OutputFileName_filter_step_ahead = {}; +OutputFileName_param = {}; +OutputFileName_forc_mean = {}; +OutputFileName_forc_point = {}; +% OutputFileName_moments = {}; +end + +for b=fpar:B + + % [deep, logpo] = GetOneDraw(typee); + % set_all_parameters(deep); + % dr = resol(oo_.steady_state,0); + if strcmpi(typee,'prior') + + [deep, logpo] = GetOneDraw(typee); + + else + deep = x(b,:); + logpo = logpost(b); + end + set_all_parameters(deep); + [dr,info] = resol(oo_.steady_state,0); + + if run_smoother + [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK] = ... + DsgeSmoother(deep,gend,Y,data_index,missing_value); + if options_.loglinear + stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ... + repmat(log(dr.ys(dr.order_var)),1,gend); + stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ... + repmat(log(dr.ys(dr.order_var)),1,gend); + else + stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ... + repmat(dr.ys(dr.order_var),1,gend); + stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ... + repmat(dr.ys(dr.order_var),1,gend); + end + stock_innov(:,:,irun(2)) = etahat; + if nvn + stock_error(:,:,irun(3)) = epsilonhat; + end + if naK + stock_filter_step_ahead(:,dr.order_var,:,irun(4)) = aK(options_.filter_step_ahead,1:endo_nbr,:); + end + + if horizon + yyyy = alphahat(iendo,i_last_obs); + yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr)); + if options_.prefilter == 1 + yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ... + horizon+maxlag,1); + end + yf(:,IdObs) = yf(:,IdObs)+(gend+[1-maxlag:horizon]')*trend_coeff'; + if options_.loglinear == 1 + yf = yf+repmat(log(SteadyState'),horizon+maxlag,1); + else + yf = yf+repmat(SteadyState',horizon+maxlag,1); + end + yf1 = forcst2(yyyy,horizon,dr,1); + if options_.prefilter == 1 + yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ... + repmat(bayestopt_.mean_varobs',[horizon+maxlag,1,1]); + end + yf1(:,IdObs,:) = yf1(:,IdObs,:)+repmat((gend+[1-maxlag:horizon]')* ... + trend_coeff',[1,1,1]); + if options_.loglinear == 1 + yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]); + else + yf1 = yf1 + repmat(SteadyState',[horizon+maxlag,1,1]); + end + + stock_forcst_mean(:,:,irun(6)) = yf'; + stock_forcst_point(:,:,irun(7)) = yf1'; + end + + end + stock_param(irun(5),:) = deep; + stock_logpo(irun(5),1) = logpo; + stock_ys(irun(5),:) = SteadyState'; + + irun = irun + ones(7,1); + +% +% TempPath=DirectoryName; +% DirectoryNamePar='C:\dynare_calcs\ModelTest\ls2003\metropolis' +% DirectoryName=DirectoryNamePar; + + + + if irun(1) > MAX_nsmoo || b == B + stock = stock_smooth(:,:,1:irun(1)-1); + ifil(1) = ifil(1) + 1; + save([DirectoryName '/' M_.fname '_smooth' int2str(ifil(1)) '.mat'],'stock'); + + stock = stock_update(:,:,1:irun(1)-1); + save([DirectoryName '/' M_.fname '_update' int2str(ifil(1)) '.mat'],'stock'); + if RemoteFlag==1, + OutputFileName_smooth = [OutputFileName_smooth; {[DirectoryName filesep], [M_.fname '_smooth' int2str(ifil(1)) '.mat']}]; + OutputFileName_update = [OutputFileName_update; {[DirectoryName filesep], [M_.fname '_update' int2str(ifil(1)) '.mat']}]; + end + irun(1) = 1; + end + + if irun(2) > MAX_ninno || b == B + stock = stock_innov(:,:,1:irun(2)-1); + ifil(2) = ifil(2) + 1; + save([DirectoryName '/' M_.fname '_inno' int2str(ifil(2)) '.mat'],'stock'); + if RemoteFlag==1, + OutputFileName_inno = [OutputFileName_inno; {[DirectoryName filesep], [M_.fname '_inno' int2str(ifil(2)) '.mat']}]; + end + irun(2) = 1; + end + + if nvn && (irun(3) > MAX_nerro || b == B) + stock = stock_error(:,:,1:irun(3)-1); + ifil(3) = ifil(3) + 1; + save([DirectoryName '/' M_.fname '_error' int2str(ifil(3)) '.mat'],'stock'); + if RemoteFlag==1, + OutputFileName_error = [OutputFileName_error; {[DirectoryName filesep], [M_.fname '_error' int2str(ifil(3)) '.mat']}]; + end + irun(3) = 1; + end + + if naK && (irun(4) > MAX_naK || b == B) + stock = stock_filter_step_ahead(:,:,:,1:irun(4)-1); + ifil(4) = ifil(4) + 1; + save([DirectoryName '/' M_.fname '_filter_step_ahead' int2str(ifil(4)) '.mat'],'stock'); + if RemoteFlag==1, + OutputFileName_filter_step_ahead = [OutputFileName_filter_step_ahead; {[DirectoryName filesep], [M_.fname '_filter_step_ahead' int2str(ifil(4)) '.mat']}]; + end + irun(4) = 1; + end + + if irun(5) > MAX_nruns || b == B + stock = stock_param(1:irun(5)-1,:); + ifil(5) = ifil(5) + 1; + save([DirectoryName '/' M_.fname '_param' int2str(ifil(5)) '.mat'],'stock','stock_logpo','stock_ys'); + if RemoteFlag==1, + OutputFileName_param = [OutputFileName_param; {[DirectoryName filesep], [M_.fname '_param' int2str(ifil(5)) '.mat']}]; + end + irun(5) = 1; + end + + if horizon && (irun(6) > MAX_nforc1 || b == B) + stock = stock_forcst_mean(:,:,1:irun(6)-1); + ifil(6) = ifil(6) + 1; + save([DirectoryName '/' M_.fname '_forc_mean' int2str(ifil(6)) '.mat'],'stock'); + if RemoteFlag==1, + OutputFileName_forc_mean = [OutputFileName_forc_mean; {[DirectoryName filesep], [M_.fname '_forc_mean' int2str(ifil(6)) '.mat']}]; + end + irun(6) = 1; + end + + if horizon && (irun(7) > MAX_nforc2 || b == B) + stock = stock_forcst_point(:,:,1:irun(7)-1); + ifil(7) = ifil(7) + 1; + save([DirectoryName '/' M_.fname '_forc_point' int2str(ifil(7)) '.mat'],'stock'); + if RemoteFlag==1, + OutputFileName_forc_point = [OutputFileName_forc_point; {[DirectoryName filesep], [M_.fname '_forc_point' int2str(ifil(7)) '.mat']}]; + end + irun(7) = 1; + end + + % if moments_varendo && (irun(8) > MAX_momentsno || b == B) + % stock = stock_moments(1:irun(8)-1); + % ifil(8) = ifil(8) + 1; + % save([DirectoryName '/' M_.fname '_moments' int2str(ifil(8)) '.mat'],'stock'); + % if RemoteFlag==1, + % OutputFileName_moments = [OutputFileName_moments; {[DirectoryName filesep], [M_.fname '_moments' int2str(ifil(8)) '.mat']}]; + % end + % irun(8) = 1; + % end + +% DirectoryName=TempPath; + + + if exist('OCTAVE_VERSION') + printf('Taking subdraws: %3.f%% done\r', b/B*100); + elseif ~whoiam, + waitbar(b/B,h); + end + if mod(b,10)==0 & whoiam, + fprintf('Done! \n'); + waitbarString = [ 'Subdraw ' int2str(b) '/' int2str(B) ' done.']; + fMessageStatus((b-fpar+1)/(B-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo) + end +end + +myoutput.ifil=ifil; +if RemoteFlag==1, +myoutput.OutputFileName = [OutputFileName_smooth; +OutputFileName_update; +OutputFileName_inno; +OutputFileName_error; +OutputFileName_filter_step_ahead; +OutputFileName_param; +OutputFileName_forc_mean; +OutputFileName_forc_point]; +% OutputFileName_moments]; +end diff --git a/matlab/random_walk_metropolis_hastings.m b/matlab/random_walk_metropolis_hastings.m index e566f9c10..825998d41 100644 --- a/matlab/random_walk_metropolis_hastings.m +++ b/matlab/random_walk_metropolis_hastings.m @@ -69,7 +69,7 @@ localVars.varargin=varargin; % tic, -if isnumeric(options_.parallel), +if isnumeric(options_.parallel) || (nblck-fblck)==0, fout = random_walk_metropolis_hastings_core(localVars, fblck, nblck, 0); record = fout.record; @@ -94,7 +94,7 @@ else % from where to get back results % NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'}; - [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'random_walk_metropolis_hastings_core', localVars, globalVars); + [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'random_walk_metropolis_hastings_core', localVars, globalVars, options_.parallel_info); for j=1:totCPU, offset = sum(nBlockPerCPU(1:j-1))+fblck-1; record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)));