Merge pull request #27 from tsipkens/developer

This merges numerous updates to the README, adjustments removing some less-often used code, support for YML configuration files, and refactoring (e.g., updating function headers to be closer to the standard Matlab format).

+invert/tikhonov.m

@@ -1,29 +1,49 @@
 
-% TIKHONOV  Performs inversion using various order Tikhonov regularization in 2D.
-% Author:   Timothy Sipkens, 2018-11-21
-% 
-% Inputs:
-%   A           Model matrix
-%   b           Data
-%   lambda      Regularization parameter
-%   order_L     Order of regularization -OR- 
-%               pre-computed Tikhonov matrix structure
-%                   (OPTIONAL, default is set by tikhonov_lpr)
-%   n           The length of first dimension of solution -OR-
-%               a grid, so as to support partial grids
-%                   (ONLY REQUIRED when order is specified)
-%   xi          Initial guess for solver
-%                   (OPTIONAL, default is zeros)
-%   solver      Solver (OPTIONAL, default is interior-point)
-%
-% Outputs:
-%   x           Regularized estimate
-%   D           Inverse operator (x = D*[b;0])
-%   Lpr0        Tikhonov matrix structure
-%   Gpo_inv     Inverse of posterior covariance
-%=========================================================================%
+% TIKHONOV  Performs inversion using various order Tikhonov regularization. 
+%  Regularization takes place in 2D. The type of regularization or prior
+%  information added depends on the order of Tikhonov applied. 
+% 
+%  X = invert.tikhonov(A,B,LAMBDA) inverts the system A*X = B using 
+%  Tikhonov regularization, where A is a model matrix or kernel and B is a 
+%  data vector, and appling a regularization parameter LAMBDA. By default, 
+%  the order of Tikhonov and thus the constructed Tikhonov matrix are given 
+%  as the default of the invert.tikhonov_lpr(...) function. 
+% 
+%  X = invert.tikhonov(A,B,LAMBDA,L) uses the pre-computed Tikhonov matrix
+%  L during the inversion. This matrix can be computed using the
+%  invert.tikhonov_lpr(...) function. 
+% 
+%  X = invert.tikhonov(A,B,LAMBDA,ORDER,N) applies Tikhonov regularization
+%  of the order specified by ORDER, where ORDER is an integer in the range 
+%  [0,2]. This requires the quantity N, an integer representing the length 
+%  of the first dimension of the solution, to build the Tikhonov matrix. 
+%  This is relevant for full (not partial) grids. 
+% 
+%  X = invert.tikhonov(A,B,LAMBDA,ORDER,GRID) applies Tikhonov regularization
+%  of the order specified by ORDER, as above. In the case of a partial
+%  grid, an instance of the Grid class is required to build the Tikhonov
+%  matrix in the place of an integer dimension.
+% 
+%  X = invert.tikhonov(...,XI) applies Tikhonov regularization using an
+%  initial guess of XI. The default if not included is zeros. 
+% 
+%  X = invert.tikhonov(...,XI,SOLVER) applies Tikhonov regularization 
+%  using the solver specified by the string in SOLVER. 
+% 
+%  [X,D] = invert.tikhonov(...) applies Tikhonov regularization, outputting
+%  explicit inverse operator implied by this procedure, specifically such
+%  that X = D*([B;0]). 
+%  
+%  [X,D,LPR0] = invert.tikhonov(...) applies Tikhonov regularization, also
+%  outputting the unscaled Tikhonov matrix. 
+% 
+%  [X,D,LPR0,GPO_INV] = invert.tikhonov(...) applies Tikhonov
+%  regularization, also outputting the inverse posterior covariance. 
+% 
+%  AUTHOR: Timothy Sipkens, 2018-11-21
 
-function [x,D,Lpr0,Gpo_inv] = tikhonov(A,b,lambda,order_L,n_grid,xi,solver)
+function [x, D, Lpr0, Gpo_inv] = ... 
+    tikhonov(A, b, lambda, order_L, n_grid, xi, solver)
 
 x_length = size(A,2);
 

+invert/twomey.m

@@ -1,25 +1,27 @@
 
-% TWOMEY   Performs inversion using the iterative Twomey approach.
-% Author:  Timothy Sipkens, 2018-11-21
-%-------------------------------------------------------------------------%
-% Inputs:
-%   A             Model matrix
-%   b             Data
-%   Lb            Cholesky factorization of inverse covariance matrix
-%   xi            Initial guess
-%   iter          Max. number of iterations
-%   sigma_fun     Function of x to be evaluated to determine convergence
-%                   (Optional, default: ignore)
-%   sigma         Value to which sigma_fun is compared to determine convergence
-%                   (Optional, default: ignore)
-%   f_bar         Boolean to determine whether or not to show textbar
-%                   (Optional, default: 0)
-%
-% Outputs:
-%   x             Twomey estimate
-%=========================================================================%
+% TWOMEY  Performs inversion using the iterative Twomey approach.
+% 
+%  X = invert.twomey(A,B) inverts the system A*X = B using Twomey 
+%  inversion, where A is a model matrix or kernel and B is a data vector.
+% 
+%  X = invert.twomey(A,B,XI) applies Twomey inversion using an initial 
+%  guess of XI. If excluded, a vector of ones is used. 
+% 
+%  X = invert.twomey(A,B,XI,ITER) applies Twomey inversion using ITER 
+%  number of iterations. If excluded, 100 iterations is used. 
+% 
+%  X = invert.twomey(A,B,XI,ITER,SIGMA_FUN,SIGMA) applies Twomey inversion 
+%  using a function of X, SIGMA_FUN, and threshold, SIGMA, to determine 
+%  convergence. If excluded, convergence is not checked and the full number
+%  of iterations is used. 
+% 
+%  X = invert.twomey(A,B,XI,ITER,SIGMA_FUN,SIGMA,F_BAR) applies Twomey 
+%  inversion, as above, with F_BAR flagging whether a progress bar is 
+%  output to the console. 
+% 
+%  AUTHOR:  Timothy Sipkens, 2018-11-21
 
-function [x] = twomey(A,b,xi,iter,sigma_fun,sigma,f_bar)
+function [x] = twomey(A, b, xi, iter, sigma_fun, sigma, f_bar)
 
 %-- Parse inputs ---------------------------------------------------------%
 if ~exist('iter','var'); iter = []; end

+io/load_config.m

@@ -0,0 +1,25 @@
+
+% LOAD_CONFIG  Loads settings from configuration file (YML). 
+%  Files are loaded in order supplied, overwriting properties where
+%  relevant. 
+%  
+%  AUTHOR: Timothy Sipkens, 2021-03-25
+
+function config = load_config(fnames)
+
+if ~iscell(fnames); fnames = {fnames}; end
+
+config = struct();
+for ii=1:length(fnames)
+
+    config0 = io.read_yml(fnames{ii});  % read new settings
+
+    % Copy (or overwrite) existing settings.
+    f = fieldnames(config0);
+    for jj = 1:length(f)
+        config.(f{jj}) = config0.(f{jj});
+    end
+
+end
+
+end

+io/read_yml.m

@@ -0,0 +1,53 @@
+
+% READ_YML  Simple little function to read simple little YAML format parameter files.
+%  
+%  AUTHOR: Lloyd Russell, 2017
+%  REPO: https://github.com/llerussell/ReadYAML
+
+function results = read_yml(filePath)
+
+% read file line by line
+fid = fopen(filePath, 'r');
+data = textscan(fid, '%s', 'delimiter', '\n', 'whitespace', '');
+fclose(fid);
+
+% remove empty lines
+data = deblank(data{1});
+data(cellfun('isempty', data)) = [];
+
+% prepare final results structure
+results = [];
+
+% parse the contents (line by line)
+for i = 1:numel(data)
+
+    % extract this line
+    thisLine = data{i};
+
+    % ignore if this line is a comment
+    if strcmpi(thisLine(1), '#')
+        continue
+    end
+
+    % find the seperator between key and value
+    sepIndex = find(thisLine==':', 1, 'first');
+
+    % get the key name (remove whitespace)
+    key = strtrim(thisLine(1:sepIndex-1));  
+
+    % get the value, ignoring any comments (remove whitespace)
+    value = strsplit(thisLine(sepIndex+1:end), '#');
+    value = strtrim(value{1}); 
+
+    % attempt to convert value to numeric type
+    [convertedValue, success] = str2num(value);
+    if success
+        value = convertedValue;
+    end
+
+    % store the key and value in the results
+    results.(key) = value;
+end
+
+end
+

+kernel/gen.m → +kernel/gen_pma_dma.m

@@ -1,31 +1,43 @@
 
-% GEN  Evaluate kernel/transfer functions for DMA-PMA to find A.
-% Author:  Timothy Sipkens, 2020-02-04
+% GEN_PMA_DMA  Evaluate kernel/transfer functions for PMA-DMA.
+%  This function applies to miscellaneous PMA-DMA setpoints. 
+%  The kernel.gen_pma_dma_grid function is preferred for computional
+%  reasions. As such, this function should only be used when the data does
+%  not adhere to a reasonable grid. 
 % 
-% Notes:
-%   1. Cell arrays are used for Omega_mat and Lambda_mat in order to 
-%      allow for the use of sparse matrices, which is necessary to 
-%      store information on higher resolutions grids
-%      (such as those used for phantoms).
-%   2. By default, this function uses the analytical PMA transfer function 
-%      corresponding to Case 1C from Sipkens et al. (Aerosol Sci. Technol. 2020b).
-%   3. Unlike gen_grid, requires a pre-computed setpoint structure (and,
-%      as such, a varargin argument is not included).
+%  A = kernel.gen_pma_dma(SP,D_STAR,GRID_I) evaluates the transfer function
+%  for the PMA setpoints specified by SP and the DMA setpoints specified by 
+%  D_STAR. The kernel is evaluated by integrating the transfer function
+%  over the elements in GRID_I. This require the same numder of entries in
+%  SP and D_STAR. 
+%  Please refer to the get_setpoint(...) function in tfer_pma folder for
+%  more details on generating the SP struture.
 % 
-% Inputs:
-%   sp          PMA setpoint structure
-%   d_star      DMA setpoints
-%   grid_i      Grid on which to perform integration
-%   prop_pma    Structure defining the properties of the PMA
-%=========================================================================%
+%  A = kernel.gen_pma_dma(SP,D_STAR,GRID_I,PROP_PMA) specifies a
+%  pre-computed PMA property data structure. If not given, the function
+%  uses the defaults of kernel.prop_pma(...).
+% 
+%  A = kernel.gen_pma_dma(SP,D_STAR,GRID_I,PROP_PMA,PROP_DMA) specifies a
+%  pre-computed DMA property data structure. 
+%  
+%  ------------------------------------------------------------------------
+% 
+%  NOTE: Cell arrays are used for Omega_mat and Lambda_mat in order to 
+%  allow for the use of sparse matrices, which is necessary to 
+%  store information on higher resolutions grids (such as those 
+%  used for phantoms).
+% 
+%  AUTHOR: Timothy Sipkens, 2020-02-04
 
-function A = gen(sp,d_star,grid_i,prop_dma)
+function A = gen_pma_dma(sp, d_star, grid_i, prop_pma, prop_dma)
 
+% If not given, import default properties of PMA, 
+% as selected by prop_pma function.
 if ~exist('prop_pma','var'); prop_pma = []; end
 if isempty(prop_pma); prop_pma = kernel.prop_pma; end
-    % import properties of PMA
-    % use default properties selected by prop_pma function
+
 if length(sp)~=length(d_star); error('Setpoint / d_star mismatch.'); end
+
 if ~exist('prop_dma','var'); prop_dma = []; end
 
 
@@ -55,7 +67,7 @@
 %   Note: The DMA transfer function is 1D (only a function of mobility),
 %   which is exploited to speed evaluation. The results is 1 by 3 cell, 
 %   with one entry per charge state.
-disp('Computing DMA contribution...');
+disp(' Computing DMA contribution...');
 Omega_mat = cell(1,n_z); % pre-allocate for speed, one cell entry per charge state
 for kk=1:n_z
     Omega_mat{kk} = sparse(N_b,n_i(2)); % pre-allocate for speed
@@ -75,12 +87,12 @@
     Omega_mat{kk} = Omega_mat{kk}(:,jj);
         % repeat transfer function for repeated mass setpoint
 end
-disp('Completed DMA contribution.');
+disp(' Completed DMA contribution.');
 disp(' ');
 
 
 %== STEP 2: Evaluate PMA transfer function ===============================%
-disp('Computing PMA contribution:');
+disp(' Computing PMA contribution:');
 tools.textbar(0); % initiate textbar
 Lambda_mat = cell(1,n_z); % pre-allocate for speed
     % one cell entry per charge state
@@ -96,18 +108,19 @@
         tools.textbar((N_b*(kk-1)+ii)/(n_z*N_b));
     end
 end
+disp(' Complete.');
 disp(' ');
 
 
 %== SETP 3: Combine to compile kernel ====================================%
-disp('Compiling kernel...');
+disp(' Compiling kernel...');
 K = sparse(N_b,N_i);
 for kk=1:n_z
     K = K+f_z(z_vec(kk),:).*... % charging contribution
         Lambda_mat{kk}(:,:).*... % PMA contribution
         Omega_mat{kk}(:,:); % DMA contribution
 end
-disp('Kernel compiled.');
+disp(' Kernel compiled.');
 
 dr_log = grid_i.dr; % area of integral elements in [logm,logd]T space
 A = bsxfun(@times,K,dr_log'); % multiply kernel by element area