tgenerate and display per-walker data - cosmo - front and backend for Markov-Chain Monte Carlo inversion of cosmogenic nuclide concentrations
git clone git://src.adamsgaard.dk/cosmo
Log
Files
Refs
README
LICENSE
---
commit 95859a1a81f7efd051b16dc279538eaec8c94c7f
parent 12eb539c6c9b5f7db4845203856aa03806071042
Author: Anders Damsgaard 
Date:   Mon,  7 Dec 2015 16:31:51 +0100

generate and display per-walker data

Diffstat:
  M index.php                           |       2 ++
  M matlab/generate_plots.m             |      53 ++++++++++++++++++++++++++++++

2 files changed, 55 insertions(+), 0 deletions(-)
---
diff --git a/index.php b/index.php
t@@ -253,6 +253,8 @@ if (isset($_GET['wait_id']) && !empty($_GET['wait_id'])) {
                   Link to FIG
+                  
+
diff --git a/matlab/generate_plots.m b/matlab/generate_plots.m
t@@ -227,6 +227,11 @@ if titles
 end
 
 %% Plot one histogram per model parameter per walker
+epsilon_int_data_w = cell(1, Nwalkers);
+epsilon_gla_data_w = cell(1, Nwalkers);
+t_degla_data_w = cell(1, Nwalkers);
+record_threshold_data_w = cell(1, Nwalkers);
+
 disp('One histogram per model parameter per walker')
 fh = [fh;figure('visible', show_figures)];
 for i1 = 1:M % for each model parameter
t@@ -257,6 +262,7 @@ for i1 = 1:M % for each model parameter
         epsilon_int_25(iwalk) = prctile(Ss{iwalk}.ms(i1,:)*1.0e6, 25);
         epsilon_int_50(iwalk) = prctile(Ss{iwalk}.ms(i1,:)*1.0e6, 50);
         epsilon_int_75(iwalk) = prctile(Ss{iwalk}.ms(i1,:)*1.0e6, 75);
+        epsilon_int_data_w{iwalk} = Ss{iwalk}.ms(i1,:)*1.0e6;
         
     elseif i1 == 2
         %xlabel('Glacial erosion rate [m/yr]')
t@@ -267,12 +273,14 @@ for i1 = 1:M % for each model parameter
         epsilon_gla_25(iwalk) = prctile(Ss{iwalk}.ms(i1,:)*1.0e6, 25);
         epsilon_gla_50(iwalk) = prctile(Ss{iwalk}.ms(i1,:)*1.0e6, 50);
         epsilon_gla_75(iwalk) = prctile(Ss{iwalk}.ms(i1,:)*1.0e6, 75);
+        epsilon_gla_data_w{iwalk} = Ss{iwalk}.ms(i1,:)*1.0e6;
         
     elseif i1 == 3
         %xlabel('Timing of last deglaciation [yr]')
         xlabel('t_{degla} [yr]')
         text(0.02,0.98,'c', 'Units', ...
             'Normalized', 'VerticalAlignment', 'Top')
+        t_degla_data_w{iwalk} = Ss{iwalk}.ms(i1,:);
         
     elseif i1 == 4
         %xlabel('$\delta^{18}$O$_\mathrm{threshold}$ [$^o/_{oo}$]', ...
t@@ -283,6 +291,7 @@ for i1 = 1:M % for each model parameter
         record_threshold_25(iwalk) = prctile(Ss{iwalk}.ms(i1,:), 25);
         record_threshold_50(iwalk) = prctile(Ss{iwalk}.ms(i1,:), 50);
         record_threshold_75(iwalk) = prctile(Ss{iwalk}.ms(i1,:), 75);
+        record_threshold_data_w{iwalk} = Ss{iwalk}.ms(i1,:);
         
     else
         disp(['Using mname for i1 = ' i1])
t@@ -1100,3 +1109,47 @@ dlmwrite([save_file, '-eps_int.txt'], epsilon_int_data);
 dlmwrite([save_file, '-eps_gla.txt'], epsilon_gla_data);
 dlmwrite([save_file, '-t_degla.txt'], t_degla_data);
 dlmwrite([save_file, '-d18O_threshold.txt'], record_threshold_data);
+
+%% Save per-walker data
+for i=1:Nwalkers
+    dlmwrite([save_file, '-eps_int-w' num2str(i) '.txt'],...
+        epsilon_int_data_w{i});
+    dlmwrite([save_file, '-eps_gla-w' num2str(i) '.txt'],...
+        epsilon_gla_data_w{i});
+    dlmwrite([save_file, '-t_degla-w' num2str(i) '.txt'],...
+        t_degla_data_w{i});
+    dlmwrite([save_file, '-d18O_threshold-w' num2str(i) '.txt'],...
+        record_threshold_data_w{i});
+end
+
+% create HTML snippet with results
+html = '\n';
+for i=1:Nwalkers
+    html = [html, ...
+        '                  Walker ', num2str(i), ...
+        ' εint ',...
+        'data\n'];
+    html = [html, ...
+        '                  Walker ', num2str(i), ...
+        ' εgla ',...
+        'data\n'];
+    html = [html, ...
+        '                  Walker ', num2str(i), ...
+        ' tdegla ',...
+        'data\n'];
+    html = [html, ...
+        '                  Walker ', num2str(i), ...
+        ' δ18O', ...
+        'threshold data\n'];
+end
+fileID = fopen([save_file, '-walker-data.html','w');
+fprintf(fileID, html);
+fclose(fileID);