GopherProxy

	tmcmc_inversion.m - cosmo - front and backend for Markov-Chain Monte Carlo inve…
	git clone git://src.adamsgaard.dk/cosmo
	Log
	Files
	Refs
	README
	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 = linspa…
	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 nuc…
	47 %fs.Nucleides = {'10Be','26Al'}; %We may switch nucleides on and…
	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 …
	74 % fs.RelErrorObs = 0.01*[2.0;2.04]';%0.02; %0.02 means 2% obs…
	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 obser…
	80
	81 if debug
	82 disp(fs.Nucleides)
	83 end
	84
	85 %fs.dobsMode = 'SyntheticNoNoise'; %'SyntheticNoNoise','Syntheti…
	86 fs.dobsMode = 'ObservedData'; %'SyntheticNoNoise','SyntheticAndN…
	87 if strcmp(fs.dobsMode,'ObservedData')
	88 fs.d_obs = [];
	89 %fs.d_obs = repmat([5.3e8;3.7e9]',length(fs.zobs),1); %<<<<<…
	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)); %e…
	95 elseif length(fs.RelErrorObs) == length(fs.Nucleides) %ok
	96 else error('fs.RelErrorObs must have length 1 or length of fs.Nu…
	97 end
	98
	99 %>........ For advec-solution
	100 fs.testmode = 'fast'; %Means "Skip advective model". may change …
	101 % fs.testmode = 'run_advec'; %Means "Run advective model for com…
	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','FixedTi…
	110 switch fs.CycleMode
	111 case 'FixedC'
	112 fs.C = 20; %If isempty, C=round(2.6e6/(round(dtGla*(1+dt…
	113 case 'FixedQuaternary'
	114 fs.tQuaternary = 2.6e6; %time of first glaciation, adjus…
	115 %First glaciation and first interglacial adjusted accord…
	116 case 'FixedTimes'
	117 fs.tStarts = NaN; %load or compute fixed times of more o…
	118 fs.relExpos = NaN; %load or compute degree of exposure i…
	119 case 'd18OTimes'
	120 %fs.d18Ofn = 'lisiecki_triinterp_2p6Ma_5ky.mat';
	121 %fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_30ky.mat';…
	122 if strcmp(record, 'rec_5kyr')
	123 fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_5ky.mat…
	124 elseif strcmp(record, 'rec_20kyr')
	125 fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_20ky.ma…
	126 elseif strcmp(record, 'rec_30kyr')
	127 fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_30ky.ma…
	128 else
	129 error(['record ' record ' not understood']);
	130 end
	131 fs.tStarts = NaN; %load or compute fixed times of more o…
	132 fs.relExpos = NaN; %load or compute degree of exposure i…
	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 t…
	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','Obse…
	171 case 'ObservedData'
	172 %print out for checking:
	173 disp('>>>>>> Check that values match nucleides and depth…
	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:',sprint…
	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_…
	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…
	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 …
	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)dtopfs.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';'PriorCorne…
	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 t…
	241 % ANDERS: consider parfor for parallel computing. However, fixed_stuff a…
	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 progr…
	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_Bur…
	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 …
	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.5sin(1:0.5:500));pause(0.3);sound(0.5sin(1:0.75:750))
	282 pause(0.6)
	283 sound(0.5sin(1:0.5:500));pause(0.3);sound(0.5sin(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')