tmcmc_inversion.m - cosmo - front and backend for Markov-Chain Monte Carlo inversion of cosmogenic nuclide concentrations
(HTM) git clone git://src.adamsgaard.dk/cosmo
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tmcmc_inversion.m (11936B)
---
1 function [Ss, save_file] = mcmc_inversion(matlab_scripts_folder, debug, ...
2 n_walkers, outfolder, ...
3 be_conc, al_conc, c_conc, ne_conc, ...
4 be_uncer, al_uncer, c_uncer, ne_uncer, ...
5 zobs, ...
6 be_prod_spall, al_prod_spall, c_prod_spall, ne_prod_spall, ...
7 be_prod_muons, al_prod_muons, c_prod_muons, ne_prod_muons, ...
8 rock_density, ...
9 epsilon_gla_min, epsilon_gla_max, ...
10 epsilon_int_min, epsilon_int_max, ...
11 t_degla_min, t_degla_max, ...
12 record, ...
13 record_threshold_min, record_threshold_max, ...
14 statusfile, sim_id)
15
16 %% mcmc_inversion.m
17 % function is called from `file_scanner_mcmc_starter.m`
18
19 beeps = false;
20
21 %clear; close all;
22 format compact;
23
24 %Set path so that we can find other required m-files
25 addpath(matlab_scripts_folder)
26
27 % save density for later use in subfunctions
28 fs.rho = rock_density;
29
30 % save production rates for later use in subfunctions
31 fs.be_prod_spall = be_prod_spall;
32 fs.al_prod_spall = al_prod_spall;
33 fs.c_prod_spall = c_prod_spall;
34 fs.ne_prod_spall = ne_prod_spall;
35
36 fs.be_prod_muons = be_prod_muons;
37 fs.al_prod_muons = al_prod_muons;
38 fs.c_prod_muons = c_prod_muons;
39 fs.ne_prod_muons = ne_prod_muons;
40
41 fs.g_case = 'CosmoLongsteps'; %must match a case in function gz = linspace(0,10,100);
42
43 switch fs.g_case
44 case 'CosmoLongsteps'
45 %>......... The observations and observation errors:
46 %fs.Nucleides = {'10Be','26Al','14C','21Ne'}; %We may switch nucleides on and off
47 %fs.Nucleides = {'10Be','26Al'}; %We may switch nucleides on and off
48 fs.Nucleides = {};
49 fs.RelErrorObs = [];
50 concentrations = [];
51 if be_conc > 0.
52 fs.Nucleides = [fs.Nucleides '10Be'];
53 fs.RelErrorObs = [fs.RelErrorObs be_uncer];
54 concentrations = [concentrations be_conc];
55 end
56 if al_conc > 0.
57 fs.Nucleides = [fs.Nucleides '26Al'];
58 fs.RelErrorObs = [fs.RelErrorObs al_uncer];
59 concentrations = [concentrations al_conc];
60 end
61 if c_conc > 0.
62 fs.Nucleides = [fs.Nucleides '14C'];
63 fs.RelErrorObs = [fs.RelErrorObs c_uncer];
64 concentrations = [concentrations c_conc];
65 end
66 if ne_conc > 0.
67 fs.Nucleides = [fs.Nucleides '21Ne'];
68 fs.RelErrorObs = [fs.RelErrorObs ne_uncer];
69 concentrations = [concentrations ne_conc];
70 end
71
72
73 % fs.RelErrorObs = 0.02;%0.02; %0.02 means 2% observational error
74 % fs.RelErrorObs = 0.01*[2.0;2.04]';%0.02; %0.02 means 2% observational error
75
76 % one depth per nuclide or not?
77 %fs.zobs = [0]; %Depths where nucleides are observed
78 fs.zobs = zobs;
79 % fs.zobs = [0,0.3,1,3,10]; %Depths where nucleides are observed
80
81 if debug
82 disp(fs.Nucleides)
83 end
84
85 %fs.dobsMode = 'SyntheticNoNoise'; %'SyntheticNoNoise','SyntheticAndNoise','ObservedData'
86 fs.dobsMode = 'ObservedData'; %'SyntheticNoNoise','SyntheticAndNoise','ObservedData'
87 if strcmp(fs.dobsMode,'ObservedData')
88 fs.d_obs = [];
89 %fs.d_obs = repmat([5.3e8;3.7e9]',length(fs.zobs),1); %<<<<<< put in values as best you can
90 fs.d_obs = repmat(concentrations',length(fs.zobs),1);
91 fs.d_obs = fs.d_obs(:);
92 end
93 if length(fs.RelErrorObs) == 1
94 fs.RelErrorObs = fs.RelErrorObs*ones(size(fs.Nucleides)); %extend to vector
95 elseif length(fs.RelErrorObs) == length(fs.Nucleides) %ok
96 else error('fs.RelErrorObs must have length 1 or length of fs.Nucleides')
97 end
98
99 %>........ For advec-solution
100 fs.testmode = 'fast'; %Means "Skip advective model". may change to 'linearized' below
101 % fs.testmode = 'run_advec'; %Means "Run advective model for comparisson"
102 D =100;
103 z = linspace(0,10,100); %Note! modified below
104 fs.D = D;
105 fs.z = D*z.^3/1000;
106 fs.dt = 100;
107
108 %>........ Structure of glacial cycles
109 fs.CycleMode = 'd18OTimes'; %'FixedC','FixedQuaternary','FixedTimes', 'd18OTimes'
110 switch fs.CycleMode
111 case 'FixedC'
112 fs.C = 20; %If isempty, C=round(2.6e6/(round(dtGla*(1+dtIdtG))));
113 case 'FixedQuaternary'
114 fs.tQuaternary = 2.6e6; %time of first glaciation, adjust as desired
115 %First glaciation and first interglacial adjusted accordingly
116 case 'FixedTimes'
117 fs.tStarts = NaN; %load or compute fixed times of more or less glaciated periods
118 fs.relExpos = NaN; %load or compute degree of exposure in periods
119 case 'd18OTimes'
120 %fs.d18Ofn = 'lisiecki_triinterp_2p6Ma_5ky.mat';
121 %fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_30ky.mat'; % zachos_triinterp_2p6Ma
122 if strcmp(record, 'rec_5kyr')
123 fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_5ky.mat'; % zachos_triinterp_2p6Ma
124 elseif strcmp(record, 'rec_20kyr')
125 fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_20ky.mat'; % zachos_triinterp_2p6Ma
126 elseif strcmp(record, 'rec_30kyr')
127 fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_30ky.mat'; % zachos_triinterp_2p6Ma
128 else
129 error(['record ' record ' not understood']);
130 end
131 fs.tStarts = NaN; %load or compute fixed times of more or less glaciated periods
132 fs.relExpos = NaN; %load or compute degree of exposure in periods
133 end
134
135 %>........ Starting condition
136 fs.Cstart = 'extend interglacial';
137 % fs.Cstart = 'zeros';
138
139 %>........ Model parameters.
140 fs.mname{1} = 'ErateInt';
141 fs.mname{2} = 'ErateGla';
142 fs.mname{3} = 'tDegla';
143 % fs.mname{4} = 'dtGla';
144 % fs.mname{5} = 'dtIdtG';
145 fs.mname{4} = 'd18Oth';
146 %>........ Prior information
147 % m = [ErateInt,ErateGla,tDegla,dtGla,dtIdtG];
148 %fs.ErateIntminmax = [1e-7,1e-3]; %0.26m to 2600 m pr. Quaternary
149 %fs.ErateGlaminmax = [1e-7,1e-3];
150 fs.ErateIntminmax = [epsilon_int_min, epsilon_int_max]; %0.26m to 2600 m pr. Quaternary
151 fs.ErateGlaminmax = [epsilon_gla_min, epsilon_gla_max];
152 %fs.tDeglaminmax = [10e3,12e3]; %8000 to 10000 yr Holocene
153 fs.tDeglaminmax = [t_degla_min, t_degla_max];
154 % fs.dtGlaminmax = [40e3,200e3];
155 % fs.dtIdtGminmax = [0,0.5];
156 %fs.d18Othminmax = [3.6,4.4];
157 fs.d18Othminmax = [record_threshold_min, record_threshold_max];
158
159 fs.ErateIntDistr = 'logunif';
160 fs.ErateGlaDistr = 'logunif';
161 fs.tDeglaDistr = 'uniform';
162 % fs.dtGlaDistr = 'uniform';
163 % fs.dtIdtGDistr = 'uniform';
164 fs.d18OthDistr = 'uniform';
165
166 for im=1:length(fs.mname)
167 fs.mminmax(im,:) = eval(['fs.',fs.mname{im},'minmax']);
168 fs.mDistr(im,:) = eval(['fs.',fs.mname{im},'Distr']);
169 end
170 switch fs.dobsMode %'SyntheticNoNoise','SyntheticAndNoise','ObservedData'
171 case 'ObservedData'
172 %print out for checking:
173 disp('>>>>>> Check that values match nucleides and depths:')
174 id = 0;
175 for iNucl=1:length(fs.Nucleides)
176 disp(fs.Nucleides{iNucl})
177 for iz=1:length(fs.zobs)
178 id = id+1;
179 disp(['>>> z=',num2str(fs.zobs(iz)),' m:',sprintf('%10g',fs.d_obs(id))])
180 end
181 end
182 case {'SyntheticNoNoise','SyntheticAndNoise'}
183 % fs.m_true = [...
184 % 1e-5;...
185 % 1e-6;...
186 % 10e3;...
187 % 100e3;...
188 % 10/100];
189 %fs.m_true = [...
190 %1e-6;...
191 %1e-4;...
192 %12e3;...
193 %4.0];
194 fs.m_true = [... % ANDERS: This is eps_int, eps_gla, t_degla, d180O_threshold
195 5e-5;...
196 1e-6;...
197 11e3;...
198 3.8];
199 fs.d_true= g(fs.m_true,fs);
200 fs.d_obs = fs.d_true + 0; %no noise, updated if dobsMode=='SyntheticAndNoise'
201 end
202 % >>>> finalizeing fixed_stuff with the observational error stds
203 % We compute errors relative to the larger surface value
204 % Otherwise the small concentrations at depth may be deemed too accurate.
205 Nz = length(fs.zobs);
206 Nnucl = length(fs.Nucleides);
207 for iNucl = 1:Nnucl
208 dtop = fs.d_obs(1+(iNucl-1)*Nz),
209 fs.ErrorStdObs((1:Nz) + (iNucl-1)*Nz,1) = ...
210 ones(Nz,1)*dtop*fs.RelErrorObs(iNucl);
211 end
212 if strcmp(fs.dobsMode,'SyntheticAndNoise')
213 fs.d_obs = fs.d_true + fs.ErrorStdObs.*randn(size(fs.d_true));
214 end
215 end %switch fs.g_case
216 % keyboard
217 %>........ For the MetHas algorithm
218 fs.Nwalkers = n_walkers; %Number of random walks
219 %fs.Nwalkers = 4; %Number of random walks
220 fs.WalkerStartMode = 'PriorEdge';%'PriorSample'; 'PriorMean';'PriorCorner';'PriorEdge'
221 fs.WalkerSeeds = 1:fs.Nwalkers; %must be at least fs.Nwalkers!
222
223 %%>... fs.BurnIn: Controlling the BurnIn phase:
224 fs.BurnIn.Nsamp = 1000; %number of samples in burn in
225 fs.BurnIn.Nskip = 1; %number of samples between samples kept
226 fs.BurnIn.ProposerType = 'Prior'; %'Native';'Prior';'ApproxPosterior'
227 fs.BurnIn.StepFactorMode = 'Fixed'; %'Fixed', 'Adaptive'
228 fs.BurnIn = CompleteFsSampling(fs.BurnIn);
229
230 %%>... fs.Sampling: Copy of fs.BurnIn but with different values set
231 fs.Sampling = fs.BurnIn;
232 fs.Sampling.Nsamp = 1e4;
233 fs.Sampling.Nskip = 1;
234 fs.Sampling.ProposerType = 'ApproxPosterior';
235 fs.Sampling.StepFactorMode = 'Adaptive';
236 fs.Sampling = CompleteFsSampling(fs.Sampling);
237
238 fixed_stuff = fs;
239
240 fixed_stuff.StartTime = now; %This should allow the program to predict time of finish
241 % ANDERS: consider parfor for parallel computing. However, fixed_stuff and
242 % S structures are incompatible with parfor.
243 for iwalk=1:fixed_stuff.Nwalkers
244
245 disp(' ')
246 disp(['#### iwalk = ' num2str(iwalk) '/' ...
247 num2str(fixed_stuff.Nwalkers) ' ####'])
248
249 fixed_stuff.iwalk = iwalk; %Helps program keep user updated on progress.
250 m_starts(:,iwalk) = WalkerStarter(iwalk,fixed_stuff);
251 d_starts(:,iwalk) = g(m_starts(:,iwalk),fixed_stuff);
252
253 %>>>>> ......Burn in:
254 seed = fixed_stuff.WalkerSeeds(iwalk);
255 isBurnIn=1;
256 [S.msBurnIn,S.acceptsBurnIn,S.QsBurnIn,S.QdsBurnIn,S.lump_MetHas_BurnIn]=MetHasLongstep4(...
257 m_starts(:,iwalk),seed,isBurnIn,fixed_stuff, ...
258 statusfile); %%%% ADDED BY ANDERS
259 mStartSampling = S.msBurnIn(:,end);
260 %<<<<< ... End Burn in
261
262 %>>>>> ......Sampling the posterior:
263 seed = fixed_stuff.WalkerSeeds(iwalk);
264 isBurnIn=0; %<<<<<<<<<<<<<<<<<< ------------------- must be changed when Sampling
265 [S.ms,S.accepts,S.Qs,S.Qds,S.lump_MetHas]=MetHasLongstep4(...
266 mStartSampling,seed,isBurnIn,fixed_stuff, ...
267 statusfile); %%%% ADDED BY ANDERS
268 %<<<<< ... End Sampling the posterior
269 S.fs = fixed_stuff;
270 S.m_start = m_starts(:,iwalk);
271 S.d_start = d_starts(:,iwalk);
272 Ss{iwalk}=S;
273 % S.ms{iwalk}=ms
274 % S.lump_MetHass{iwalk}=lump_MetHas;
275 if beeps
276 sound(sin(1:0.5:500))
277 end
278 end
279
280 if beeps
281 sound(0.5*sin(1:0.5:500));pause(0.3);sound(0.5*sin(1:0.75:750))
282 pause(0.6)
283 sound(0.5*sin(1:0.5:500));pause(0.3);sound(0.5*sin(1:0.75:750))
284 end
285
286 %save_file = [outfolder, '/', id, '_Walks_',datestr(now,'yyyymmdd_HHMMSS')];
287 save_file = strcat(outfolder, '/', char(sim_id), '_Walks');
288 save(save_file,'Ss','save_file')