tadd preliminary plot generator - 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 975cbd537b90d0449110c9a21a0a7effc32341ba
(DIR) parent 18f347b2531e003200eb970faa5c465cb5d0cfa8
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 26 Aug 2015 16:19:20 +0200
add preliminary plot generator
Diffstat:
M matlab/file_scanner_mcmc_starter.m | 6 +++++-
A matlab/generate_plots.m | 206 +++++++++++++++++++++++++++++++
M matlab/mcmc_inversion.m | 12 ++++++++----
3 files changed, 219 insertions(+), 5 deletions(-)
---
(DIR) diff --git a/matlab/file_scanner_mcmc_starter.m b/matlab/file_scanner_mcmc_starter.m
t@@ -24,6 +24,9 @@ matlab_scripts_folder = 'm_pakke2014maj11/';
%% general settings
debug = true; % show debugging output to matlab console
+% output graphics formats
+graphics_formats = ['fig', 'png', 'pdf'];
+
%% initialization
addpath(matlab_scripts_folder);
t@@ -73,7 +76,8 @@ while 1
record_threshold_min, record_threshold_max);
% generate plots
- CompareWalks2(Ss, save_file)
+ %CompareWalks2(Ss, save_file)
+ generate_plots(Ss, save_file, graphics_formats)
% delete or archive the file so it is not processed again
%delete(infile)
(DIR) diff --git a/matlab/generate_plots.m b/matlab/generate_plots.m
t@@ -0,0 +1,206 @@
+function generate_plots(Ss, save_file, formats)
+
+%% Copied from m_pakke2014maj11/CompareWalks2.m
+% Generates and saves all relevant figures
+
+S = Ss{1};
+fixed_stuff = S.fs;
+Nwalkers = fixed_stuff.Nwalkers;
+M = size(fixed_stuff.mminmax,1);
+
+%Compare BurnIn
+fh(1)=figure;
+for iwalk=1:min(4,Nwalkers) %only up to the first four walks
+ QsBurnIns(:,iwalk)=Ss{iwalk}.QsBurnIn;
+ Qss(:,iwalk) =Ss{iwalk}.Qs;
+ acceptss(:,iwalk)=Ss{iwalk}.accepts;
+ normdmpropspr(:,iwalk) = sqrt(sum(Ss{iwalk}.lump_MetHas.dmpropspr.^2)/M);
+ normdmcurspr(:,iwalk) = sqrt(sum(Ss{iwalk}.lump_MetHas.dmcurspr.^2)/M);
+end
+subplot(2,4,1)
+semilogy((1:size(QsBurnIns,1)),QsBurnIns,'-');
+ylim([0.1,1e4])
+ylabel('St. sum of sq.')
+title('Burn in')
+grid on
+
+subplot(2,4,2:4)
+semilogy(size(QsBurnIns,1)+(1:size(Qss,1)),Qss,'.')
+legend('1','2','3','4','5','6','7','8','9','10','location','NorthEast')
+ylim([0.1,1e4])
+% ylabel('St. sum of sq.')
+title(['Sampling. Result file =',save_file],'interp','none')
+grid on
+
+subplot(2,4,6)
+hist(normdmcurspr)
+xlabel('||dm_{cur}||')
+ylimit = ylim;
+title('Cur2cur step lengths (Stdz)')
+
+subplot(2,4,5)
+hist(normdmpropspr)
+xlabel('||dm_{prop}||')
+ylim(ylimit)
+title('Proposed step lengths (Stdz)')
+
+subplot(2,4,7)
+textfontsize=8;
+fs = Ss{1}.fs;
+title('Data and Model')
+%text out nucleides, RelError, depths
+axis([-0.03,1,-1,9])
+str = ['Nucleides'];
+for i=1:length(fs.Nucleides)
+ str = [str,' - ',fs.Nucleides{i}];
+end
+str = {str,['>> Rel.Error=',num2str(fs.RelErrorObs)],...
+ ['>> zobs = ',sprintf('%g, ',fs.zobs)]};
+text(0,1,str,'fontsize',textfontsize,'interp','none')
+%text out prior intervals
+for im=1:size(fs.mminmax,1)
+ %str=[fs.mname{im},'_minmax=[',sprintf('%g,%g]. True=%g',fs.mminmax(im,:),fs.m_true(im))];
+ text(0,im+2,str,'fontsize',textfontsize,'interp','none')
+end
+
+axis ij
+box on
+
+subplot(2,4,8)
+title('MCMC-analysis')
+str = ['Nwalkers:',num2str(fs.Nwalkers),'. ',fs.WalkerStartMode];
+text(0,1,str,'fontsize',textfontsize,'interp','none')
+
+fsamp = fs.BurnIn;
+str = {['BurnIn: ',num2str(fsamp.Nsamp),' skip ',num2str(fsamp.Nskip),' of ',fsamp.ProposerType],...
+ ['>> StepFactor ',fsamp.StepFactorMode]},
+switch fsamp.StepFactorMode
+ case 'Fixed', str{2} = [str{2},'=',num2str(fsamp.StepFactor)];
+ case 'Adaptive', str{2} = [str{2},', AccRat=',num2str(fsamp.TargetAcceptanceRatio)];
+end
+text(0,2,str,'fontsize',textfontsize,'interp','none')
+
+fsamp = fs.Sampling;
+str = {['Sampling: ',num2str(fsamp.Nsamp),' skip ',num2str(fsamp.Nskip),' of ',fsamp.ProposerType],...
+ ['>> StepFactor ',fsamp.StepFactorMode]},
+switch fsamp.StepFactorMode
+ case 'Fixed', str{2} = [str{2},'=',num2str(fsamp.StepFactor)];
+ case 'Adaptive', str{2} = [str{2},', AccRat=',num2str(fsamp.TargetAcceptanceRatio)];
+end
+text(0,3.5,str,'fontsize',textfontsize,'interp','none')
+
+axis([-0.03,1,0,10])
+axis ij
+box on
+
+%Compare sampling cross-plots
+fh = [fh;figure];
+Nbin = 50;
+
+mminmax = fixed_stuff.mminmax; %bounds of uniform prior intervals
+mDistr = fixed_stuff.mDistr; %'uniform' or 'logunif' for each coordinate
+% M = size(ms,1);
+
+isub1 = 0;
+for i1=1:M
+ mminmax_i1 = mminmax(i1,:);
+ % mDistr_i1 = mDistr(i1,:);
+ % xi = ms(i1,:);
+ switch mDistr(i1,:)
+ case 'uniform',
+ dxbin = diff(mminmax_i1)/Nbin; %we want outside bins to check that we do not sample outside
+ xbins{i1} = linspace(mminmax_i1(1)-dxbin,mminmax_i1(2)+2*dxbin,Nbin+4); %one extra to make pcolor work
+ case 'logunif'
+ dfacxbin = (mminmax_i1(2)/mminmax_i1(1))^(1/(Nbin)); %we want outside bins to check that we do not sample outside
+ xbins{i1} = logspace(log10(mminmax_i1(1)/dfacxbin),log10(mminmax_i1(2)*dfacxbin^2),Nbin+4);
+ end
+
+ % for i2 = 1:M
+ for i2 = (i1+1):M
+ % yi = ms(i2,:);
+
+ mminmax_i2 = mminmax(i2,:);
+ switch mDistr(i2,:)
+ case 'uniform',
+ dybin = diff(mminmax_i2)/Nbin; %we want outside bins to check that we do not sample outside
+ ybins = linspace(mminmax_i2(1)-dybin,mminmax_i2(2)+2*dybin,Nbin+4);
+ case 'logunif'
+ dfacybin = (mminmax_i2(2)/mminmax_i2(1))^(1/Nbin); %we want outside bins to check that we do not sample outside
+ ybins = logspace(log10(mminmax_i2(1)/dfacybin),log10(mminmax_i2(2)*dfacybin^2),Nbin+4);
+ end
+ W = ones(2); W = conv2(W,W); W = W/sum(W(:));
+ %loop over walks
+ isub1 = isub1+1;
+ for iwalk = 1:min(4,Nwalkers)
+ xi = Ss{iwalk}.ms(i1,:);yi = Ss{iwalk}.ms(i2,:);
+ [smoothgrid,histgrid] = smoothdens(xi,yi,xbins{i1},ybins,W);
+ isub2 = iwalk; isub = isub2 + (isub1-1)*4;
+ if isub>4*5
+ isub1 = 1; isub=1;
+ fh = [fh;figure];
+ end
+ subplot(5,4,isub)
+ pcolor(xbins{i1},ybins,smoothgrid);
+ xlabel(fixed_stuff.mname{i1})
+ ylabel(fixed_stuff.mname{i2})
+ grid on; shading flat; axis tight; set(gca,'fontsize',8);
+ switch mDistr(i1,:)
+ case 'uniform', set(gca,'xscale','lin')
+ case 'logunif', set(gca,'xscale','log')
+ end
+ switch mDistr(i2,:)
+ case 'uniform', set(gca,'yscale','lin')
+ case 'logunif', set(gca,'yscale','log')
+ end
+ axis([mminmax_i1,mminmax_i2])
+ end
+ end
+end
+
+fh = [fh;figure];
+for i1 = 1:M
+ for iwalk=1:min(4,Nwalkers)
+ isub = (i1-1)*4 + iwalk;
+ subplot(5,4,isub)
+ Nhistc=histc(Ss{iwalk}.ms(i1,:),xbins{i1});
+ bar(xbins{i1},Nhistc,'histc')
+ xlabel(fixed_stuff.mname{i1})
+ %set (gca,'xtick',[1e-7:1e-3]);
+
+ switch mDistr(i1,:)
+ case 'uniform', set(gca,'xscale','lin','xlim',mminmax(i1,:))
+ case 'logunif', set(gca,'xscale','log','xlim',mminmax(i1,:))
+ end
+ end
+end
+
+%Putting in titles over figure 2:4
+figure(fh(2))
+subplot(5,4,2)
+title(['Density cross-plots A. Result file =',save_file],'interp','none')
+figure(fh(3))
+subplot(5,4,2)
+title(['Density cross-plots B. Result file =',save_file],'interp','none')
+figure(fh(4))
+subplot(5,4,2)
+title(['Histograms. Result file =',save_file],'interp','none')
+
+figpos1 = [6 474 1910 504];
+figpos2 =[ 12 94 645 887];
+figpos3 =[ 610 94 645 887];
+figpos4 =[ 1207 94 740 887];
+set(fh(2),'pos',figpos2)
+set(fh(3),'pos',figpos3)
+set(fh(4),'pos',figpos4)
+set(fh(1),'pos',figpos1)
+figure(fh(1))
+
+% for i1 = 1:M
+% hold off
+% subplot(M,M,(M-1)*M+i1)
+% xlabel(fixed_stuff.mname{i1})
+% subplot(M,M,(i1-1)*M+1)
+% ylabel(fixed_stuff.mname{i1})
+% end
+% i=1;
+
(DIR) diff --git a/matlab/mcmc_inversion.m b/matlab/mcmc_inversion.m
t@@ -126,12 +126,16 @@ switch fs.g_case
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.tDeglaminmax = [10e3,12e3]; %8000 to 10000 yr Holocene
+ %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 = [3.6,4.4];
+ fs.d18Othminmax = [record_threshold_min, record_threshold_max];
fs.ErateIntDistr = 'logunif';
fs.ErateGlaDistr = 'logunif';