From c8329915b2114c20c6958bb6874854999cd7fdb7 Mon Sep 17 00:00:00 2001
From: ratto <ratto@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Fri, 15 May 2009 16:33:20 +0000
Subject: [PATCH] Seeds are always initialized and stored for each chain
 independently

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2675 ac1d8469-bf42-47a9-8791-bf33cf982152
---
 matlab/independent_metropolis_hastings.m    |  6 +++--
 matlab/metropolis_hastings_initialization.m | 29 ++++++++++++++-------
 2 files changed, 23 insertions(+), 12 deletions(-)

diff --git a/matlab/independent_metropolis_hastings.m b/matlab/independent_metropolis_hastings.m
index 7f74a813a..a39b2dc5e 100644
--- a/matlab/independent_metropolis_hastings.m
+++ b/matlab/independent_metropolis_hastings.m
@@ -57,6 +57,8 @@ load([MhDirectoryName '/' ModelName '_mh_history'],'record');
 InitSizeArray = min([MAX_nruns*ones(nblck) nruns],[],2);
 jscale = diag(bayestopt_.jscale);
 for b = fblck:nblck
+   randn('state',record.Seeds(b).Normal);
+   rand('state',record.Seeds(b).Unifor);
     if (options_.load_mh_file~=0)  & (fline(b)>1) & OpenOldFile(b)
         load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ...
               '_blck' int2str(b) '.mat'])
@@ -141,9 +143,9 @@ for b = fblck:nblck
     end% End of the simulations for one mh-block.
     record.AcceptationRates(b) = isux/j;
     close(hh);
+    record.Seeds(b).Normal = randn('state');
+    record.Seeds(b).Unifor = rand('state');
 end% End of the loop over the mh-blocks.
-record.Seeds.Normal = randn('state');
-record.Seeds.Unifor = rand('state');
 save([MhDirectoryName '/' ModelName '_mh_history'],'record');
 disp(['MH: Number of mh files			: ' int2str(NewFile(1)) ' per block.'])
 disp(['MH: Total number of generated files	: ' int2str(NewFile(1)*nblck) '.'])
diff --git a/matlab/metropolis_hastings_initialization.m b/matlab/metropolis_hastings_initialization.m
index d8ad0e4de..9d65fbb6b 100644
--- a/matlab/metropolis_hastings_initialization.m
+++ b/matlab/metropolis_hastings_initialization.m
@@ -160,8 +160,15 @@ if ~options_.load_mh_file & ~options_.mh_recover
     record.MhDraws(1,2) = AnticipatedNumberOfFiles;
     record.MhDraws(1,3) = AnticipatedNumberOfLinesInTheLastFile;
     record.AcceptationRates = zeros(1,nblck);
-    record.Seeds.Normal = randn('state');
-    record.Seeds.Unifor = rand('state');
+    % separate initializaton for each chain
+    JSUM = 0;
+    for j=1:nblck,
+      JSUM  = JSUM + sum(100*clock);
+      randn('state',JSUM);
+      rand('state',JSUM);
+      record.Seeds(j).Normal = randn('state');
+      record.Seeds(j).Unifor = rand('state');
+    end
     record.InitialParameters = ix2;
     record.InitialLogLiK = ilogpo2;
     record.LastParameters = zeros(nblck,npar);
@@ -174,14 +181,16 @@ if ~options_.load_mh_file & ~options_.mh_recover
     fprintf(fidlog,['    Expected number of lines in the last file: ' ...
                     int2str(AnticipatedNumberOfLinesInTheLastFile) '.\n']);
     fprintf(fidlog,['\n']);
-    fprintf(fidlog,['    Initial seed (randn):\n']);
-    for i=1:length(record.Seeds.Normal)
-        fprintf(fidlog,['      ' num2str(record.Seeds.Normal(i)') '\n']);
+    for j = 1:nblck,
+    fprintf(fidlog,['    Initial seed (randn) for chain number ',int2str(j),':\n']);
+    for i=1:length(record.Seeds(j).Normal)
+        fprintf(fidlog,['      ' num2str(record.Seeds(j).Normal(i)') '\n']);
     end
-    fprintf(fidlog,['    Initial seed (rand).:\n']);
-    for i=1:length(record.Seeds.Unifor)
-        fprintf(fidlog,['      ' num2str(record.Seeds.Unifor(i)') '\n']);
+    fprintf(fidlog,['    Initial seed (rand) for chain number ',int2str(j),':\n']);
+    for i=1:length(record.Seeds(j).Unifor)
+        fprintf(fidlog,['      ' num2str(record.Seeds(j).Unifor(i)') '\n']);
     end
+    end,
     fprintf(fidlog,' \n');
     fclose(fidlog);
 elseif options_.load_mh_file & ~options_.mh_recover
@@ -239,8 +248,8 @@ elseif options_.load_mh_file & ~options_.mh_recover
     record.MhDraws(end,1) = nruns(1);
     record.MhDraws(end,2) = AnticipatedNumberOfFiles;
     record.MhDraws(end,3) = AnticipatedNumberOfLinesInTheLastFile;
-    randn('state',record.Seeds.Normal);
-    rand('state',record.Seeds.Unifor);
+%     randn('state',record.Seeds.Normal);
+%     rand('state',record.Seeds.Unifor);
     save([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
     disp(['MH: ... It''s done. I''ve loaded ' int2str(NumberOfPreviousSimulations) ' simulations.'])
     disp(' ')