1function SampleAddress = selec_posterior_draws(SampleSize,drsize)
2% Selects a sample of draws from the posterior distribution and if nargin>1
3% saves the draws in _pdraws mat files (metropolis folder). If drsize>0
4% the dr structure, associated to the parameters, is also saved in _pdraws.
5% This routine is more efficient than metropolis_draw.m because here an
6% _mh file cannot be opened twice.
7%
8% INPUTS
9%   o SampleSize     [integer]  Size of the sample to build.
10%   o drsize         [double]   structure dr is drsize megaoctets.
11%
12% OUTPUTS
13%   o SampleAddress  [integer]  A (SampleSize*4) array, each line specifies the
14%                               location of a posterior draw:
15%                                  Column 2 --> Chain number
16%                                  Column 3 --> (mh) File number
17%                                  Column 4 --> (mh) line number
18%
19% SPECIAL REQUIREMENTS
20%   None.
21%
22
23% Copyright (C) 2006-2017 Dynare Team
24%
25% This file is part of Dynare.
26%
27% Dynare is free software: you can redistribute it and/or modify
28% it under the terms of the GNU General Public License as published by
29% the Free Software Foundation, either version 3 of the License, or
30% (at your option) any later version.
31%
32% Dynare is distributed in the hope that it will be useful,
33% but WITHOUT ANY WARRANTY; without even the implied warranty of
34% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
35% GNU General Public License for more details.
36%
37% You should have received a copy of the GNU General Public License
38% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
39
40global M_ options_ estim_params_ oo_
41
42% Number of parameters:
43npar = estim_params_.nvx;
44npar = npar + estim_params_.nvn;
45npar = npar + estim_params_.ncx;
46npar = npar + estim_params_.ncn;
47npar = npar + estim_params_.np;
48
49% Select one task:
50switch nargin
51  case 1
52    info = 0;
53  case 2
54    MAX_mega_bytes = 10;% Should be an option...
55    if drsize>0
56        info=2;
57    else
58        info=1;
59    end
60    drawsize = drsize+npar*8/1048576;
61  otherwise
62    error(['selec_posterior_draws:: Unexpected number of input arguments!'])
63end
64
65MetropolisFolder = CheckPath('metropolis',M_.dname);
66ModelName = M_.fname;
67BaseName = [MetropolisFolder filesep ModelName];
68
69% Get informations about the mcmc:
70load_last_mh_history_file(MetropolisFolder, ModelName);
71FirstMhFile = record.KeepedDraws.FirstMhFile;
72FirstLine = record.KeepedDraws.FirstLine;
73TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
74LastMhFile = TotalNumberOfMhFiles;
75TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
76NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
77MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
78mh_nblck = options_.mh_nblck;
79
80% Randomly select draws in the posterior distribution:
81SampleAddress = zeros(SampleSize,4);
82for i = 1:SampleSize
83    ChainNumber = ceil(rand*mh_nblck);
84    DrawNumber  = ceil(rand*NumberOfDraws);
85    SampleAddress(i,1) = DrawNumber;
86    SampleAddress(i,2) = ChainNumber;
87    if DrawNumber <= MAX_nruns-FirstLine+1
88        MhFileNumber = FirstMhFile;
89        MhLineNumber = FirstLine+DrawNumber-1;
90    else
91        DrawNumber  = DrawNumber-(MAX_nruns-FirstLine+1);
92        MhFileNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
93        MhLineNumber = DrawNumber-(MhFileNumber-FirstMhFile-1)*MAX_nruns;
94    end
95    SampleAddress(i,3) = MhFileNumber;
96    SampleAddress(i,4) = MhLineNumber;
97end
98SampleAddress = sortrows(SampleAddress,[3 2]);
99
100% Selected draws in the posterior distribution, and if drsize>0
101% reduced form solutions, are saved on disk.
102if info
103    %delete old stale files before creating new ones
104    delete_stale_file([BaseName '_posterior_draws*.mat'])
105    if  SampleSize*drawsize <= MAX_mega_bytes% The posterior draws are saved in one file.
106        pdraws = cell(SampleSize,info);
107        old_mhfile = 0;
108        old_mhblck = 0;
109        for i = 1:SampleSize
110            mhfile = SampleAddress(i,3);
111            mhblck = SampleAddress(i,2);
112            if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
113                load([BaseName '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2')
114            end
115            pdraws(i,1) = {x2(SampleAddress(i,4),:)};
116            if info-1
117                set_parameters(pdraws{i,1});
118                [dr,info,M_,options_,oo_] =compute_decision_rules(M_,options_,oo_);
119                pdraws(i,2) = { dr };
120            end
121            old_mhfile = mhfile;
122            old_mhblck = mhblck;
123        end
124        clear('x2')
125        save([BaseName '_posterior_draws1.mat'],'pdraws','estim_params_')
126    else% The posterior draws are saved in xx files.
127        NumberOfDrawsPerFile = fix(MAX_mega_bytes/drawsize);
128        NumberOfFiles = ceil(SampleSize*drawsize/MAX_mega_bytes);
129        NumberOfLines = SampleSize - (NumberOfFiles-1)*NumberOfDrawsPerFile;
130        linee = 0;
131        fnum  = 1;
132        pdraws = cell(NumberOfDrawsPerFile,info);
133        old_mhfile = 0;
134        old_mhblck = 0;
135        for i=1:SampleSize
136            linee = linee+1;
137            mhfile = SampleAddress(i,3);
138            mhblck = SampleAddress(i,2);
139            if (mhfile ~= old_mhfile) || (mhblck ~= old_mhblck)
140                load([BaseName '_mh' num2str(mhfile) '_blck' num2str(mhblck) '.mat'],'x2')
141            end
142            pdraws(linee,1) = {x2(SampleAddress(i,4),:)};
143            if info-1
144                set_parameters(pdraws{linee,1});
145                [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
146                pdraws(linee,2) = { dr };
147            end
148            old_mhfile = mhfile;
149            old_mhblck = mhblck;
150            if fnum < NumberOfFiles && linee == NumberOfDrawsPerFile
151                linee = 0;
152                save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_')
153                fnum = fnum+1;
154                if fnum < NumberOfFiles
155                    pdraws = cell(NumberOfDrawsPerFile,info);
156                else
157                    pdraws = cell(NumberOfLines,info);
158                end
159            end
160        end
161        save([BaseName '_posterior_draws' num2str(fnum) '.mat'],'pdraws','estim_params_')
162    end
163end