diff --git a/matlab/McMCDiagnostics.m b/matlab/McMCDiagnostics.m
index 98c0923cb39676450a5b332b219ab83e3c4c6f65..08301172d681acf7b94f7764fbd13dd721eed0d0 100644
--- a/matlab/McMCDiagnostics.m
+++ b/matlab/McMCDiagnostics.m
@@ -83,52 +83,78 @@ end
disp('MCMC Diagnostics: Univariate convergence diagnostic, Brooks and Gelman (1998):')
-for j=1:npar
- fprintf(' Parameter %d... ',j);
- for b = 1:nblck
- startline = 0;
- for n = 1:NumberOfMcFilesPerBlock
- %load([MhDirectoryName '/' mcfiles(n,1,b).name],'x2');
- load([MhDirectoryName '/' M_.fname '_mh',int2str(n),'_blck' int2str(b) ...
- '.mat'],'x2');
- nx2 = size(x2,1);
- tmp((b-1)*NumberOfDraws+startline+(1:nx2),1) = x2(:,j);
-% clear x2;
- startline = startline + nx2;
- end
-% $$$ %load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'x2');
-% $$$ load([MhDirectoryName '/' M_.fname '_mh',int2str(NumberOfMcFilesPerBlock),'_blck' int2str(b) '.mat'],'x2');
-% $$$ tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = x2(:,j);
-% $$$ clear x2;
-% $$$ startline = startline + LastLineNumber;
- end
- tmp(:,2) = kron(transpose(1:nblck),ones(NumberOfDraws,1));
- tmp(:,3) = kron(ones(nblck,1),time');
- tmp = sortrows(tmp,1);
- ligne = 0;
- for iter = Origin:StepSize:NumberOfDraws
- ligne = ligne+1;
- linea = ceil(0.5*iter);
- n = iter-linea+1;
- cinf = round(n*ALPHA/2);
- csup = round(n*(1-ALPHA/2));
- CINF = round(nblck*n*ALPHA/2);
- CSUP = round(nblck*n*(1-ALPHA/2));
- temp = tmp(find((tmp(:,3)>=linea) & (tmp(:,3)<=iter)),1:2);
- UDIAG(ligne,1,j) = temp(CSUP,1)-temp(CINF,1);
- moyenne = mean(temp(:,1));%% Pooled mean.
- UDIAG(ligne,3,j) = sum((temp(:,1)-moyenne).^2)/(nblck*n-1);
- UDIAG(ligne,5,j) = sum(abs(temp(:,1)-moyenne).^3)/(nblck*n-1);
- for i=1:nblck
- pmet = temp(find(temp(:,2)==i));
- UDIAG(ligne,2,j) = UDIAG(ligne,2,j) + pmet(csup,1)-pmet(cinf,1);
- moyenne = mean(pmet,1); %% Within mean.
- UDIAG(ligne,4,j) = UDIAG(ligne,4,j) + sum((pmet(:,1)-moyenne).^2)/(n-1);
- UDIAG(ligne,6,j) = UDIAG(ligne,6,j) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1);
+
+localVars.MhDirectoryName = MhDirectoryName;
+localVars.nblck = nblck;
+localVars.NumberOfMcFilesPerBlock = NumberOfMcFilesPerBlock;
+localVars.Origin = Origin;
+localVars.StepSize = StepSize;
+localVars.NumberOfDraws = NumberOfDraws;
+localVars.NumberOfLines = NumberOfLines;
+localVars.time = time;
+localVars.M_ = M_;
+
+if isnumeric(options_.parallel),% | isunix, % for the moment exclude unix platform from parallel implementation
+ fout = McMCDiagnostics_core(localVars,1,npar,0);
+ UDIAG = fout.UDIAG;
+ clear fout
+
+else
+
+ [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, 1, npar,{},'McMCDiagnostics_core', localVars);
+ UDIAG = [];
+ for j=1:totCPU,
+ UDIAG = cat(3,UDIAG ,fout(j).UDIAG);
end
- end
- fprintf('Done! \n');
+
end
+
+% for j=1:npar
+% fprintf(' Parameter %d... ',j);
+% for b = 1:nblck
+% startline = 0;
+% for n = 1:NumberOfMcFilesPerBlock
+% %load([MhDirectoryName '/' mcfiles(n,1,b).name],'x2');
+% load([MhDirectoryName '/' M_.fname '_mh',int2str(n),'_blck' int2str(b) ...
+% '.mat'],'x2');
+% nx2 = size(x2,1);
+% tmp((b-1)*NumberOfDraws+startline+(1:nx2),1) = x2(:,j);
+% % clear x2;
+% startline = startline + nx2;
+% end
+% % $$$ %load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'x2');
+% % $$$ load([MhDirectoryName '/' M_.fname '_mh',int2str(NumberOfMcFilesPerBlock),'_blck' int2str(b) '.mat'],'x2');
+% % $$$ tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = x2(:,j);
+% % $$$ clear x2;
+% % $$$ startline = startline + LastLineNumber;
+% end
+% tmp(:,2) = kron(transpose(1:nblck),ones(NumberOfDraws,1));
+% tmp(:,3) = kron(ones(nblck,1),time');
+% tmp = sortrows(tmp,1);
+% ligne = 0;
+% for iter = Origin:StepSize:NumberOfDraws
+% ligne = ligne+1;
+% linea = ceil(0.5*iter);
+% n = iter-linea+1;
+% cinf = round(n*ALPHA/2);
+% csup = round(n*(1-ALPHA/2));
+% CINF = round(nblck*n*ALPHA/2);
+% CSUP = round(nblck*n*(1-ALPHA/2));
+% temp = tmp(find((tmp(:,3)>=linea) & (tmp(:,3)<=iter)),1:2);
+% UDIAG(ligne,1,j) = temp(CSUP,1)-temp(CINF,1);
+% moyenne = mean(temp(:,1));%% Pooled mean.
+% UDIAG(ligne,3,j) = sum((temp(:,1)-moyenne).^2)/(nblck*n-1);
+% UDIAG(ligne,5,j) = sum(abs(temp(:,1)-moyenne).^3)/(nblck*n-1);
+% for i=1:nblck
+% pmet = temp(find(temp(:,2)==i));
+% UDIAG(ligne,2,j) = UDIAG(ligne,2,j) + pmet(csup,1)-pmet(cinf,1);
+% moyenne = mean(pmet,1); %% Within mean.
+% UDIAG(ligne,4,j) = UDIAG(ligne,4,j) + sum((pmet(:,1)-moyenne).^2)/(n-1);
+% UDIAG(ligne,6,j) = UDIAG(ligne,6,j) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1);
+% end
+% end
+% fprintf('Done! \n');
+% end
UDIAG(:,[2 4 6],:) = UDIAG(:,[2 4 6],:)/nblck;
disp(' ')
clear pmet temp moyenne CSUP CINF csup cinf n linea iter tmp;
diff --git a/matlab/McMCDiagnostics_core.m b/matlab/McMCDiagnostics_core.m
new file mode 100644
index 0000000000000000000000000000000000000000..9827d2056693f726755a4b161802d13eeddcc186
--- /dev/null
+++ b/matlab/McMCDiagnostics_core.m
@@ -0,0 +1,76 @@
+function myoutput = McMCDiagnostics_core(myinputs,fpar,npar,whoiam, ThisMatlab)
+
+if nargin<4,
+ whoiam=0;
+end
+struct2local(myinputs);
+
+if ~exist('MhDirectoryName'),
+ MhDirectoryName = CheckPath('metropolis');
+end
+
+ALPHA = 0.2; % increase too much with the number of simulations.
+tmp = zeros(NumberOfDraws*nblck,3);
+UDIAG = zeros(NumberOfLines,6,npar-fpar+1);
+% keyboard;
+
+ if whoiam
+% keyboard;
+ waitbarString = ['Please wait... McMCDiagnostics (' int2str(fpar) 'of' int2str(npar) ')...'];
+ waitbarTitle=['McMCDiagnostics ',Parallel(ThisMatlab).PcName];
+ fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo);
+
+ end
+for j=fpar:npar,
+ fprintf(' Parameter %d... ',j);
+ for b = 1:nblck
+ startline = 0;
+ for n = 1:NumberOfMcFilesPerBlock
+ %load([MhDirectoryName '/' mcfiles(n,1,b).name],'x2');
+ load([MhDirectoryName '/' M_.fname '_mh',int2str(n),'_blck' int2str(b) ...
+ '.mat'],'x2');
+ nx2 = size(x2,1);
+ tmp((b-1)*NumberOfDraws+startline+(1:nx2),1) = x2(:,j);
+% clear x2;
+ startline = startline + nx2;
+ end
+% $$$ %load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'x2');
+% $$$ load([MhDirectoryName '/' M_.fname '_mh',int2str(NumberOfMcFilesPerBlock),'_blck' int2str(b) '.mat'],'x2');
+% $$$ tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = x2(:,j);
+% $$$ clear x2;
+% $$$ startline = startline + LastLineNumber;
+ end
+ tmp(:,2) = kron(transpose(1:nblck),ones(NumberOfDraws,1));
+ tmp(:,3) = kron(ones(nblck,1),time');
+ tmp = sortrows(tmp,1);
+ ligne = 0;
+ for iter = Origin:StepSize:NumberOfDraws
+ ligne = ligne+1;
+ linea = ceil(0.5*iter);
+ n = iter-linea+1;
+ cinf = round(n*ALPHA/2);
+ csup = round(n*(1-ALPHA/2));
+ CINF = round(nblck*n*ALPHA/2);
+ CSUP = round(nblck*n*(1-ALPHA/2));
+ temp = tmp(find((tmp(:,3)>=linea) & (tmp(:,3)<=iter)),1:2);
+ UDIAG(ligne,1,j-fpar+1) = temp(CSUP,1)-temp(CINF,1);
+ moyenne = mean(temp(:,1));%% Pooled mean.
+ UDIAG(ligne,3,j-fpar+1) = sum((temp(:,1)-moyenne).^2)/(nblck*n-1);
+ UDIAG(ligne,5,j-fpar+1) = sum(abs(temp(:,1)-moyenne).^3)/(nblck*n-1);
+ for i=1:nblck
+ pmet = temp(find(temp(:,2)==i));
+ UDIAG(ligne,2,j-fpar+1) = UDIAG(ligne,2,j-fpar+1) + pmet(csup,1)-pmet(cinf,1);
+ moyenne = mean(pmet,1); %% Within mean.
+ UDIAG(ligne,4,j-fpar+1) = UDIAG(ligne,4,j-fpar+1) + sum((pmet(:,1)-moyenne).^2)/(n-1);
+ UDIAG(ligne,6,j-fpar+1) = UDIAG(ligne,6,j-fpar+1) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1);
+ end
+ end
+ fprintf('Done! \n');
+ if whoiam,
+% keyboard;
+ waitbarString = [ 'Parameter ' int2str(j) '/' int2str(npar) ' done.'];
+ fMessageStatus((j-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo)
+ end
+end
+
+myoutput.UDIAG = UDIAG;
\ No newline at end of file