tauto indent - 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
---
(DIR) commit 2c2e64edc6a4ff20b7ced2fd89aac4228e276817
(DIR) parent acc981ec7c797ca4cf883e70a2dad4f9cecd1dcf
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 21 Oct 2015 11:01:49 +0200
auto indent
Diffstat:
M matlab/mcmc_inversion.m | 354 ++++++++++++++++----------------
1 file changed, 177 insertions(+), 177 deletions(-)
---
(DIR) diff --git a/matlab/mcmc_inversion.m b/matlab/mcmc_inversion.m
t@@ -16,7 +16,7 @@ function [Ss, save_file] = mcmc_inversion(matlab_scripts_folder, debug, ...
beeps = false;
-%clear; close all;
+%clear; close all;
format compact;
%Set path so that we can find other required m-files
t@@ -24,177 +24,177 @@ addpath(matlab_scripts_folder)
fs.g_case = 'CosmoLongsteps'; %must match a case in function gz = linspace(0,10,100);
-switch fs.g_case
- case 'CosmoLongsteps'
- %>......... The observations and observation errors:
- %fs.Nucleides = {'10Be','26Al','14C','21Ne'}; %We may switch nucleides on and off
- %fs.Nucleides = {'10Be','26Al'}; %We may switch nucleides on and off
- fs.Nucleides = {};
- fs.RelErrorObs = [];
- concentrations = [];
- if be_conc > 0.
- fs.Nucleides = [fs.Nucleides '10Be'];
- fs.RelErrorObs = [fs.RelErrorObs be_uncer];
- concentrations = [concentrations be_conc];
- end
- if al_conc > 0.
- fs.Nucleides = [fs.Nucleides '26Al'];
- fs.RelErrorObs = [fs.RelErrorObs al_uncer];
- concentrations = [concentrations al_conc];
- end
- if c_conc > 0.
- fs.Nucleides = [fs.Nucleides '14C'];
- fs.RelErrorObs = [fs.RelErrorObs c_uncer];
- concentrations = [concentrations c_conc];
- end
- if ne_conc > 0.
- fs.Nucleides = [fs.Nucleides '21Ne'];
- fs.RelErrorObs = [fs.RelErrorObs ne_uncer];
- concentrations = [concentrations ne_conc];
- end
-
-
-% fs.RelErrorObs = 0.02;%0.02; %0.02 means 2% observational error
-% fs.RelErrorObs = 0.01*[2.0;2.04]';%0.02; %0.02 means 2% observational error
-
- % one depth per nuclide or not?
- fs.zobs = [0]; %Depths where nucleides are observed
-% fs.zobs = [0,0.3,1,3,10]; %Depths where nucleides are observed
-
- if debug
- disp(fs.Nucleides)
- end
-
- %fs.dobsMode = 'SyntheticNoNoise'; %'SyntheticNoNoise','SyntheticAndNoise','ObservedData'
- fs.dobsMode = 'ObservedData'; %'SyntheticNoNoise','SyntheticAndNoise','ObservedData'
- if strcmp(fs.dobsMode,'ObservedData')
- fs.d_obs = [];
- %fs.d_obs = repmat([5.3e8;3.7e9]',length(fs.zobs),1); %<<<<<< put in values as best you can
- fs.d_obs = repmat(concentrations',length(fs.zobs),1);
- fs.d_obs = fs.d_obs(:);
- end
- if length(fs.RelErrorObs) == 1
- fs.RelErrorObs = fs.RelErrorObs*ones(size(fs.Nucleides)); %extend to vector
- elseif length(fs.RelErrorObs) == length(fs.Nucleides) %ok
- else error('fs.RelErrorObs must have length 1 or length of fs.Nucleides')
- end
-
- %>........ For advec-solution
- fs.testmode = 'fast'; %Means "Skip advective model". may change to 'linearized' below
- % fs.testmode = 'run_advec'; %Means "Run advective model for comparisson"
- D =100;
- z = linspace(0,10,100); %Note! modified below
- fs.D = D;
- fs.z = D*z.^3/1000;
- fs.dt = 100;
-
- %>........ Structure of glacial cycles
- fs.CycleMode = 'd18OTimes'; %'FixedC','FixedQuaternary','FixedTimes', 'd18OTimes'
- switch fs.CycleMode
- case 'FixedC'
- fs.C = 20; %If isempty, C=round(2.6e6/(round(dtGla*(1+dtIdtG))));
- case 'FixedQuaternary'
- fs.tQuaternary = 2.6e6; %time of first glaciation, adjust as desired
- %First glaciation and first interglacial adjusted accordingly
- case 'FixedTimes'
- fs.tStarts = NaN; %load or compute fixed times of more or less glaciated periods
- fs.relExpos = NaN; %load or compute degree of exposure in periods
- case 'd18OTimes'
- %fs.d18Ofn = 'lisiecki_triinterp_2p6Ma_5ky.mat';
- %fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_30ky.mat'; % zachos_triinterp_2p6Ma
- if strcmp(record, 'rec_5kyr')
- fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_5ky.mat'; % zachos_triinterp_2p6Ma
- elseif strcmp(record, 'rec_20kyr')
- fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_20ky.mat'; % zachos_triinterp_2p6Ma
- elseif strcmp(record, 'rec_30kyr')
- fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_30ky.mat'; % zachos_triinterp_2p6Ma
- else
- error(['record ' record ' not understood']);
- end
- fs.tStarts = NaN; %load or compute fixed times of more or less glaciated periods
- fs.relExpos = NaN; %load or compute degree of exposure in periods
- end
-
- %>........ Starting condition
- fs.Cstart = 'extend interglacial';
- % fs.Cstart = 'zeros';
-
- %>........ Model parameters.
- fs.mname{1} = 'ErateInt';
- fs.mname{2} = 'ErateGla';
- fs.mname{3} = 'tDegla';
-% fs.mname{4} = 'dtGla';
-% fs.mname{5} = 'dtIdtG';
- fs.mname{4} = 'd18Oth';
- %>........ Prior information
- % m = [ErateInt,ErateGla,tDegla,dtGla,dtIdtG];
- %fs.ErateIntminmax = [1e-7,1e-3]; %0.26m to 2600 m pr. Quaternary
- %fs.ErateGlaminmax = [1e-7,1e-3];
- fs.ErateIntminmax = [epsilon_int_min, epsilon_int_max]; %0.26m to 2600 m pr. Quaternary
- fs.ErateGlaminmax = [epsilon_gla_min, epsilon_gla_max];
- %fs.tDeglaminmax = [10e3,12e3]; %8000 to 10000 yr Holocene
- fs.tDeglaminmax = [t_degla - t_degla_uncer, t_degla + t_degla_uncer];
- % fs.dtGlaminmax = [40e3,200e3];
- % fs.dtIdtGminmax = [0,0.5];
- %fs.d18Othminmax = [3.6,4.4];
- fs.d18Othminmax = [record_threshold_min, record_threshold_max];
-
- fs.ErateIntDistr = 'logunif';
- fs.ErateGlaDistr = 'logunif';
- fs.tDeglaDistr = 'uniform';
- % fs.dtGlaDistr = 'uniform';
- % fs.dtIdtGDistr = 'uniform';
- fs.d18OthDistr = 'uniform';
-
- for im=1:length(fs.mname)
- fs.mminmax(im,:) = eval(['fs.',fs.mname{im},'minmax']);
- fs.mDistr(im,:) = eval(['fs.',fs.mname{im},'Distr']);
- end
- switch fs.dobsMode %'SyntheticNoNoise','SyntheticAndNoise','ObservedData'
- case 'ObservedData'
- %print out for checking:
- disp('>>>>>> Check that values match nucleides and depths:')
- id = 0;
- for iNucl=1:length(fs.Nucleides)
- disp(fs.Nucleides{iNucl})
- for iz=1:length(fs.zobs)
- id = id+1;
- disp(['>>> z=',num2str(fs.zobs(iz)),' m:',sprintf('%10g',fs.d_obs(id))])
- end
+switch fs.g_case
+ case 'CosmoLongsteps'
+ %>......... The observations and observation errors:
+ %fs.Nucleides = {'10Be','26Al','14C','21Ne'}; %We may switch nucleides on and off
+ %fs.Nucleides = {'10Be','26Al'}; %We may switch nucleides on and off
+ fs.Nucleides = {};
+ fs.RelErrorObs = [];
+ concentrations = [];
+ if be_conc > 0.
+ fs.Nucleides = [fs.Nucleides '10Be'];
+ fs.RelErrorObs = [fs.RelErrorObs be_uncer];
+ concentrations = [concentrations be_conc];
+ end
+ if al_conc > 0.
+ fs.Nucleides = [fs.Nucleides '26Al'];
+ fs.RelErrorObs = [fs.RelErrorObs al_uncer];
+ concentrations = [concentrations al_conc];
+ end
+ if c_conc > 0.
+ fs.Nucleides = [fs.Nucleides '14C'];
+ fs.RelErrorObs = [fs.RelErrorObs c_uncer];
+ concentrations = [concentrations c_conc];
+ end
+ if ne_conc > 0.
+ fs.Nucleides = [fs.Nucleides '21Ne'];
+ fs.RelErrorObs = [fs.RelErrorObs ne_uncer];
+ concentrations = [concentrations ne_conc];
+ end
+
+
+ % fs.RelErrorObs = 0.02;%0.02; %0.02 means 2% observational error
+ % fs.RelErrorObs = 0.01*[2.0;2.04]';%0.02; %0.02 means 2% observational error
+
+ % one depth per nuclide or not?
+ fs.zobs = [0]; %Depths where nucleides are observed
+ % fs.zobs = [0,0.3,1,3,10]; %Depths where nucleides are observed
+
+ if debug
+ disp(fs.Nucleides)
+ end
+
+ %fs.dobsMode = 'SyntheticNoNoise'; %'SyntheticNoNoise','SyntheticAndNoise','ObservedData'
+ fs.dobsMode = 'ObservedData'; %'SyntheticNoNoise','SyntheticAndNoise','ObservedData'
+ if strcmp(fs.dobsMode,'ObservedData')
+ fs.d_obs = [];
+ %fs.d_obs = repmat([5.3e8;3.7e9]',length(fs.zobs),1); %<<<<<< put in values as best you can
+ fs.d_obs = repmat(concentrations',length(fs.zobs),1);
+ fs.d_obs = fs.d_obs(:);
+ end
+ if length(fs.RelErrorObs) == 1
+ fs.RelErrorObs = fs.RelErrorObs*ones(size(fs.Nucleides)); %extend to vector
+ elseif length(fs.RelErrorObs) == length(fs.Nucleides) %ok
+ else error('fs.RelErrorObs must have length 1 or length of fs.Nucleides')
+ end
+
+ %>........ For advec-solution
+ fs.testmode = 'fast'; %Means "Skip advective model". may change to 'linearized' below
+ % fs.testmode = 'run_advec'; %Means "Run advective model for comparisson"
+ D =100;
+ z = linspace(0,10,100); %Note! modified below
+ fs.D = D;
+ fs.z = D*z.^3/1000;
+ fs.dt = 100;
+
+ %>........ Structure of glacial cycles
+ fs.CycleMode = 'd18OTimes'; %'FixedC','FixedQuaternary','FixedTimes', 'd18OTimes'
+ switch fs.CycleMode
+ case 'FixedC'
+ fs.C = 20; %If isempty, C=round(2.6e6/(round(dtGla*(1+dtIdtG))));
+ case 'FixedQuaternary'
+ fs.tQuaternary = 2.6e6; %time of first glaciation, adjust as desired
+ %First glaciation and first interglacial adjusted accordingly
+ case 'FixedTimes'
+ fs.tStarts = NaN; %load or compute fixed times of more or less glaciated periods
+ fs.relExpos = NaN; %load or compute degree of exposure in periods
+ case 'd18OTimes'
+ %fs.d18Ofn = 'lisiecki_triinterp_2p6Ma_5ky.mat';
+ %fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_30ky.mat'; % zachos_triinterp_2p6Ma
+ if strcmp(record, 'rec_5kyr')
+ fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_5ky.mat'; % zachos_triinterp_2p6Ma
+ elseif strcmp(record, 'rec_20kyr')
+ fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_20ky.mat'; % zachos_triinterp_2p6Ma
+ elseif strcmp(record, 'rec_30kyr')
+ fs.d18O_filename = 'lisiecki_triinterp_2p6Ma_30ky.mat'; % zachos_triinterp_2p6Ma
+ else
+ error(['record ' record ' not understood']);
+ end
+ fs.tStarts = NaN; %load or compute fixed times of more or less glaciated periods
+ fs.relExpos = NaN; %load or compute degree of exposure in periods
+ end
+
+ %>........ Starting condition
+ fs.Cstart = 'extend interglacial';
+ % fs.Cstart = 'zeros';
+
+ %>........ Model parameters.
+ fs.mname{1} = 'ErateInt';
+ fs.mname{2} = 'ErateGla';
+ fs.mname{3} = 'tDegla';
+ % fs.mname{4} = 'dtGla';
+ % fs.mname{5} = 'dtIdtG';
+ fs.mname{4} = 'd18Oth';
+ %>........ Prior information
+ % m = [ErateInt,ErateGla,tDegla,dtGla,dtIdtG];
+ %fs.ErateIntminmax = [1e-7,1e-3]; %0.26m to 2600 m pr. Quaternary
+ %fs.ErateGlaminmax = [1e-7,1e-3];
+ fs.ErateIntminmax = [epsilon_int_min, epsilon_int_max]; %0.26m to 2600 m pr. Quaternary
+ fs.ErateGlaminmax = [epsilon_gla_min, epsilon_gla_max];
+ %fs.tDeglaminmax = [10e3,12e3]; %8000 to 10000 yr Holocene
+ fs.tDeglaminmax = [t_degla - t_degla_uncer, t_degla + t_degla_uncer];
+ % fs.dtGlaminmax = [40e3,200e3];
+ % fs.dtIdtGminmax = [0,0.5];
+ %fs.d18Othminmax = [3.6,4.4];
+ fs.d18Othminmax = [record_threshold_min, record_threshold_max];
+
+ fs.ErateIntDistr = 'logunif';
+ fs.ErateGlaDistr = 'logunif';
+ fs.tDeglaDistr = 'uniform';
+ % fs.dtGlaDistr = 'uniform';
+ % fs.dtIdtGDistr = 'uniform';
+ fs.d18OthDistr = 'uniform';
+
+ for im=1:length(fs.mname)
+ fs.mminmax(im,:) = eval(['fs.',fs.mname{im},'minmax']);
+ fs.mDistr(im,:) = eval(['fs.',fs.mname{im},'Distr']);
+ end
+ switch fs.dobsMode %'SyntheticNoNoise','SyntheticAndNoise','ObservedData'
+ case 'ObservedData'
+ %print out for checking:
+ disp('>>>>>> Check that values match nucleides and depths:')
+ id = 0;
+ for iNucl=1:length(fs.Nucleides)
+ disp(fs.Nucleides{iNucl})
+ for iz=1:length(fs.zobs)
+ id = id+1;
+ disp(['>>> z=',num2str(fs.zobs(iz)),' m:',sprintf('%10g',fs.d_obs(id))])
+ end
+ end
+ case {'SyntheticNoNoise','SyntheticAndNoise'}
+ % fs.m_true = [...
+ % 1e-5;...
+ % 1e-6;...
+ % 10e3;...
+ % 100e3;...
+ % 10/100];
+ %fs.m_true = [...
+ %1e-6;...
+ %1e-4;...
+ %12e3;...
+ %4.0];
+ fs.m_true = [... % ANDERS: This is eps_int, eps_gla, t_degla, d180O_threshold
+ 5e-5;...
+ 1e-6;...
+ 11e3;...
+ 3.8];
+ fs.d_true= g(fs.m_true,fs);
+ fs.d_obs = fs.d_true + 0; %no noise, updated if dobsMode=='SyntheticAndNoise'
+ end
+ % >>>> finalizeing fixed_stuff with the observational error stds
+ % We compute errors relative to the larger surface value
+ % Otherwise the small concentrations at depth may be deemed too accurate.
+ Nz = length(fs.zobs);
+ Nnucl = length(fs.Nucleides);
+ for iNucl = 1:Nnucl
+ dtop = fs.d_obs(1+(iNucl-1)*Nz),
+ fs.ErrorStdObs((1:Nz) + (iNucl-1)*Nz,1) = ...
+ ones(Nz,1)*dtop*fs.RelErrorObs(iNucl);
+ end
+ if strcmp(fs.dobsMode,'SyntheticAndNoise')
+ fs.d_obs = fs.d_true + fs.ErrorStdObs.*randn(size(fs.d_true));
end
- case {'SyntheticNoNoise','SyntheticAndNoise'}
-% fs.m_true = [...
-% 1e-5;...
-% 1e-6;...
-% 10e3;...
-% 100e3;...
-% 10/100];
- %fs.m_true = [...
- %1e-6;...
- %1e-4;...
- %12e3;...
- %4.0];
- fs.m_true = [... % ANDERS: This is eps_int, eps_gla, t_degla, d180O_threshold
- 5e-5;...
- 1e-6;...
- 11e3;...
- 3.8];
- fs.d_true= g(fs.m_true,fs);
- fs.d_obs = fs.d_true + 0; %no noise, updated if dobsMode=='SyntheticAndNoise'
- end
- % >>>> finalizeing fixed_stuff with the observational error stds
- % We compute errors relative to the larger surface value
- % Otherwise the small concentrations at depth may be deemed too accurate.
- Nz = length(fs.zobs);
- Nnucl = length(fs.Nucleides);
- for iNucl = 1:Nnucl
- dtop = fs.d_obs(1+(iNucl-1)*Nz),
- fs.ErrorStdObs((1:Nz) + (iNucl-1)*Nz,1) = ...
- ones(Nz,1)*dtop*fs.RelErrorObs(iNucl);
- end
- if strcmp(fs.dobsMode,'SyntheticAndNoise')
- fs.d_obs = fs.d_true + fs.ErrorStdObs.*randn(size(fs.d_true));
- end
end %switch fs.g_case
% keyboard
%>........ For the MetHas algorithm
t@@ -229,8 +229,8 @@ for iwalk=1:fixed_stuff.Nwalkers
d_starts(:,iwalk) = g(m_starts(:,iwalk),fixed_stuff);
%>>>>> ......Burn in:
- seed = fixed_stuff.WalkerSeeds(iwalk);
- isBurnIn=1;
+ seed = fixed_stuff.WalkerSeeds(iwalk);
+ isBurnIn=1;
[S.msBurnIn,S.acceptsBurnIn,S.QsBurnIn,S.QdsBurnIn,S.lump_MetHas_BurnIn]=MetHasLongstep4(...
m_starts(:,iwalk),seed,isBurnIn,fixed_stuff, ...
statusfile); %%%% ADDED BY ANDERS
t@@ -238,7 +238,7 @@ for iwalk=1:fixed_stuff.Nwalkers
%<<<<< ... End Burn in
%>>>>> ......Sampling the posterior:
- seed = fixed_stuff.WalkerSeeds(iwalk);
+ seed = fixed_stuff.WalkerSeeds(iwalk);
isBurnIn=0; %<<<<<<<<<<<<<<<<<< ------------------- must be changed when Sampling
[S.ms,S.accepts,S.Qs,S.Qds,S.lump_MetHas]=MetHasLongstep4(...
mStartSampling,seed,isBurnIn,fixed_stuff, ...
t@@ -248,9 +248,9 @@ for iwalk=1:fixed_stuff.Nwalkers
S.m_start = m_starts(:,iwalk);
S.d_start = d_starts(:,iwalk);
Ss{iwalk}=S;
-% S.ms{iwalk}=ms
-% S.lump_MetHass{iwalk}=lump_MetHas;
- if beeps
+ % S.ms{iwalk}=ms
+ % S.lump_MetHass{iwalk}=lump_MetHas;
+ if beeps
sound(sin(1:0.5:500))
end
end