tuse new variables, convert to matlab units after reading - 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 a93f27d9ebd13f7d8566ba0d712dfb453feda10d
(DIR) parent 2b9186117adfb7aabad5c46eeae74f0cf628f493
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 18 Nov 2015 13:35:59 +0100
use new variables, convert to matlab units after reading
Diffstat:
M matlab/import_php_file.m | 21 +++++++++++++++++++++
D matlab/import_php_file.m~ | 163 -------------------------------
M matlab/mcmc_inversion.m | 7 ++++---
3 files changed, 25 insertions(+), 166 deletions(-)
---
(DIR) diff --git a/matlab/import_php_file.m b/matlab/import_php_file.m
t@@ -162,3 +162,23 @@ record_threshold_max = cell2mat(rawNumericColumns(:, 29)); % 35
nwalkers = cell2mat(rawNumericColumns(:, 30)); % 36
+%% Change units
+be_conc = be_conc*1000.; % atoms/g to atoms/kg
+al_conc = al_conc*1000.; % atoms/g to atoms/kg
+c_conc = c_conc*1000.; % atoms/g to atoms/kg
+ne_conc = ne_conc*1000.; % atoms/g to atoms/kg
+
+be_prod_muon = be_prod_muon*1000.; % atoms/g/yr to atoms/kg/yr
+al_prod_muon = al_prod_muon*1000.; % atoms/g/yr to atoms/kg/yr
+c_prod_muon = c_prod_muon*1000.; % atoms/g/yr to atoms/kg/yr
+ne_prod_muon = ne_prod_muon*1000.; % atoms/g/yr to atoms/kg/yr
+
+be_prod_spall = be_prod_spall*1000.; % atoms/g/yr to atoms/kg/yr
+al_prod_spall = al_prod_spall*1000.; % atoms/g/yr to atoms/kg/yr
+c_prod_spall = c_prod_spall*1000.; % atoms/g/yr to atoms/kg/yr
+ne_prod_spall = ne_prod_spall*1000.; % atoms/g/yr to atoms/kg/yr
+
+epsilon_gla_min = epsilon_gla_min/1000.; % m/Myr to mm/yr
+epsilon_gla_max = epsilon_gla_max/1000.; % m/Myr to mm/yr
+epsilon_int_min = epsilon_int_min/1000.; % m/Myr to mm/yr
+epsilon_int_max = epsilon_int_max/1000.; % m/Myr to mm/yr
+\ No newline at end of file
(DIR) diff --git a/matlab/import_php_file.m~ b/matlab/import_php_file.m~
t@@ -1,163 +0,0 @@
-function [sample_id, name, email, ...
- lat, long, ...
- be_conc, al_conc, c_conc, ne_conc, ...
- be_uncer, al_uncer, c_uncer, ne_uncer, ...
- be_zobs, al_zobs, c_zobs, ne_zobs, ...
- be_prod_muon, al_prod_muon, c_prod_muon, ne_prod_muon, ...
- be_prod_spall, al_prod_spall, c_prod_spall, ne_prod_spall, ...
- rock_density, ...
- epsilon_gla_min, epsilon_gla_max, ...
- epsilon_int_min, epsilon_int_max, ...
- t_degla_min, t_degla_max, ...
- record, record_threshold_min, record_threshold_max, ...
- nwalkers] = ...
- import_php_file(filename, startRow, endRow)
-
-%% import_php_file.m
-% Automatically generated using the `uiimport` tool in Matlab.
-% If the output format in "uploadhistory.php" is changed, update this file
-% accordingly.
-% All columns are initially read as strings. Some of the columns (specified
-% by the col vector) are converted to numbers.
-
-%IMPORTFILE Import numeric data from a text file as column vectors.
-% [SAMPLEID,NAME,EMAIL,LAT,LONG,BE_CONC,AL_CONC,C_CONC,NE_CONC,BE_UNCER,AL_UNCER,C_UNCER,NE_UNCER,BE_ZOBS,AL_ZOBS,C_ZOBS,NE_ZOBS,BE_PROD,AL_PROD,C_PROD,NE_PROD,ROCK_DENSITY,EPSILON_GLA_MIN,EPSILON_GLA_MAX,EPSILON_INT_MIN,EPSILON_INT_MAX,T_DEGLA,T_DEGLA_UNCER,RECORD,RECORD_THRESHOLD_MIN,RECORD_THRESHOLD_MAX]
-% = IMPORTFILE(FILENAME) Reads data from text file FILENAME for the
-% default selection.
-%
-% [SAMPLEID,NAME,EMAIL,LAT,LONG,BE_CONC,AL_CONC,C_CONC,NE_CONC,BE_UNCER,AL_UNCER,C_UNCER,NE_UNCER,BE_ZOBS,AL_ZOBS,C_ZOBS,NE_ZOBS,BE_PROD,AL_PROD,C_PROD,NE_PROD,ROCK_DENSITY,EPSILON_GLA_MIN,EPSILON_GLA_MAX,EPSILON_INT_MIN,EPSILON_INT_MAX,T_DEGLA,T_DEGLA_UNCER,RECORD,RECORD_THRESHOLD_MIN,RECORD_THRESHOLD_MAX]
-% = IMPORTFILE(FILENAME, STARTROW, ENDROW) Reads data from rows STARTROW
-% through ENDROW of text file FILENAME.
-%
-% Example:
-% [sampleid,name,email,lat,long,be_conc,al_conc,c_conc,ne_conc,be_uncer,al_uncer,c_uncer,ne_uncer,be_zobs,al_zobs,c_zobs,ne_zobs,be_prod,al_prod,c_prod,ne_prod,rock_density,epsilon_gla_min,epsilon_gla_max,epsilon_int_min,epsilon_int_max,t_degla,t_degla_uncer,record,record_threshold_min,record_threshold_max]
-% = importfile('cosmo_pgpzvt',1, 1);
-%
-% See also TEXTSCAN.
-
-% Auto-generated by MATLAB on 2015/08/24 12:47:00
-
-%% Initialize variables.
-delimiter = '\t';
-if nargin<=2
- startRow = 1;
- endRow = inf;
-end
-
-%% Read columns of data as strings:
-% For more information, see the TEXTSCAN documentation.
-%formatSpec = '%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%[^\n\r]';
-formatSpec = '%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%[^\n\r]';
-
-%% Open the text file.
-fileID = fopen(filename,'r');
-
-%% Read columns of data according to format string.
-% This call is based on the structure of the file used to generate this
-% code. If an error occurs for a different file, try regenerating the code
-% from the Import Tool.
-dataArray = textscan(fileID, formatSpec, endRow(1)-startRow(1)+1, 'Delimiter', delimiter, 'HeaderLines', startRow(1)-1, 'ReturnOnError', false);
-for block=2:length(startRow)
- frewind(fileID);
- dataArrayBlock = textscan(fileID, formatSpec, endRow(block)-startRow(block)+1, 'Delimiter', delimiter, 'HeaderLines', startRow(block)-1, 'ReturnOnError', false);
- for col=1:length(dataArray)
- dataArray{col} = [dataArray{col};dataArrayBlock{col}];
- end
-end
-
-%% Close the text file.
-fclose(fileID);
-
-%% Convert the contents of columns containing numeric strings to numbers.
-% Replace non-numeric strings with NaN.
-raw = repmat({''},length(dataArray{1}),length(dataArray)-1);
-for col=1:length(dataArray)-1
- raw(1:length(dataArray{col}),col) = dataArray{col};
-end
-numericData = NaN(size(dataArray{1},1),size(dataArray,2));
-
-% the columns in col are numeric
-%for col=[6,7,8,9,10,11,12,13,18,19,20,21,22,23,24,26,27]
-for col=[6:32, 34:36]
- % Converts strings in the input cell array to numbers. Replaced non-numeric
- % strings with NaN.
- rawData = dataArray{col};
- for row=1:size(rawData, 1);
- % Create a regular expression to detect and remove non-numeric prefixes and
- % suffixes.
- regexstr = '(?<prefix>.*?)(?<numbers>([-]*(\d+[\,]*)+[\.]{0,1}\d*[eEdD]{0,1}[-+]*\d*[i]{0,1})|([-]*(\d+[\,]*)*[\.]{1,1}\d+[eEdD]{0,1}[-+]*\d*[i]{0,1}))(?<suffix>.*)';
- try
- result = regexp(rawData{row}, regexstr, 'names');
- numbers = result.numbers;
-
- % Detected commas in non-thousand locations.
- invalidThousandsSeparator = false;
- if any(numbers==',');
- thousandsRegExp = '^\d+?(\,\d{3})*\.{0,1}\d*$';
- if isempty(regexp(thousandsRegExp, ',', 'once'));
- numbers = NaN;
- invalidThousandsSeparator = true;
- end
- end
- % Convert numeric strings to numbers.
- if ~invalidThousandsSeparator;
- numbers = textscan(strrep(numbers, ',', ''), '%f');
- numericData(row, col) = numbers{1};
- raw{row, col} = numbers{1};
- end
- catch me
- end
- end
-end
-
-
-%% Split data into numeric and cell columns.
-
-% rows with numbers, check that range matches values in for loop l. 68 and
-% the list below
-rawNumericColumns = raw(:, [6:28, 30:31]);
-
-% rows with strings
-rawCellColumns = raw(:, [1:5, 29]);
-
-
-%% Allocate imported array to column variable names
-% use rawCellColumns(:, i) for text fields and
-% cell2mat(rawNumericColumns(:, i)) for numeric fields
-sample_id = rawCellColumns(:, 1); % 1
-name = rawCellColumns(:, 2); % 2
-email = rawCellColumns(:, 3); % 3
-lat = rawCellColumns(:, 4); % 4
-long = rawCellColumns(:, 5); % 5
-be_conc = cell2mat(rawNumericColumns(:, 1)); % 6
-al_conc = cell2mat(rawNumericColumns(:, 2)); % 7
-c_conc = cell2mat(rawNumericColumns(:, 3)); % 8
-ne_conc = cell2mat(rawNumericColumns(:, 4)); % 9
-be_uncer = cell2mat(rawNumericColumns(:, 5)); % 10
-al_uncer = cell2mat(rawNumericColumns(:, 6)); % 11
-c_uncer = cell2mat(rawNumericColumns(:, 7)); % 12
-ne_uncer = cell2mat(rawNumericColumns(:, 8)); % 13
-be_zobs = cell2mat(rawNumericColumns(:, 9)); % 14
-al_zobs = cell2mat(rawNumericColumns(:, 10)); % 15
-c_zobs = cell2mat(rawNumericColumns(:, 11)); % 16
-ne_zobs = cell2mat(rawNumericColumns(:, 12)); % 17
-be_prod_muon = cell2mat(rawNumericColumns(:, 13)); % 18
-al_prod_muon = cell2mat(rawNumericColumns(:, 14)); % 19
-c_prod_muon = cell2mat(rawNumericColumns(:, 15)); % 20
-ne_prod_muon = cell2mat(rawNumericColumns(:, 16)); % 21
-be_prod_spall = cell2mat(rawNumericColumns(:, 13)); % 18
-al_prod_spall = cell2mat(rawNumericColumns(:, 14)); % 19
-c_prod_spall = cell2mat(rawNumericColumns(:, 15)); % 20
-ne_prod_spall = cell2mat(rawNumericColumns(:, 16)); % 21
-rock_density = cell2mat(rawNumericColumns(:, 17)); % 22
-epsilon_gla_min = cell2mat(rawNumericColumns(:, 18)); % 23
-epsilon_gla_max = cell2mat(rawNumericColumns(:, 19)); % 24
-epsilon_int_min = cell2mat(rawNumericColumns(:, 20)); % 25
-epsilon_int_max = cell2mat(rawNumericColumns(:, 21)); % 26
-t_degla_min = cell2mat(rawNumericColumns(:, 22)); % 27
-t_degla_max = cell2mat(rawNumericColumns(:, 23)); % 28
-record = rawCellColumns(:, 6); % 29
-record_threshold_min = cell2mat(rawNumericColumns(:, 24)); % 30
-record_threshold_max = cell2mat(rawNumericColumns(:, 25)); % 31
-
-
(DIR) diff --git a/matlab/mcmc_inversion.m b/matlab/mcmc_inversion.m
t@@ -3,11 +3,12 @@ function [Ss, save_file] = mcmc_inversion(matlab_scripts_folder, debug, ...
be_conc, al_conc, c_conc, ne_conc, ...
be_uncer, al_uncer, c_uncer, ne_uncer, ...
be_zobs, al_zobs, c_zobs, ne_zobs, ...
- be_prod, al_prod, c_prod, ne_prod, ...
+ be_prod_muons, al_prod_muons, c_prod_muons, ne_prod_muons, ...
+ be_prod_spall, al_prod_spall, c_prod_spall, ne_prod_spall, ...
rock_density, ...
epsilon_gla_min, epsilon_gla_max, ...
epsilon_int_min, epsilon_int_max, ...
- t_degla, t_degla_uncer, ...
+ t_degla_min, t_degla_max, ...
record, ...
record_threshold_min, record_threshold_max, ...
statusfile, sim_id)
t@@ -137,7 +138,7 @@ switch fs.g_case
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.tDeglaminmax = [t_degla_min, t_degla_max];
% fs.dtGlaminmax = [40e3,200e3];
% fs.dtIdtGminmax = [0,0.5];
%fs.d18Othminmax = [3.6,4.4];