Merge pull request #14 from tsipkens/developer

Merge ahead of v2.0 release.

+invert/tikhonov_lpr.m

@@ -4,35 +4,35 @@
 %-------------------------------------------------------------------------%
 % Inputs:
 %   order       Order of the Tikhonov operator
-%   n           Length of first dimension of solution or Grid for x
+%   n_grid      Length of first dimension of solution or Grid for x
 %   x_length    Length of x vector (only used if a Grid is not 
 %               specified for n_grid)
 %
 % Outputs:
 %   Lpr0        Tikhonov matrix
 %=========================================================================%
 
-function Lpr0 = tikhonov_lpr(order,n,x_length)
+function Lpr0 = tikhonov_lpr(order,n_grid,x_length)
 
 %-- Generate Tikhonov smoothing matrix -----------------------------------%
 switch order
     case 0 % 0th order Tikhonov
-        if 	isa(n,'Grid') % use Grid method (for partial grid support)
-            Lpr0 = -speye(n.Ne);
+        if 	isa(n_grid,'Grid') % use Grid method (for partial grid support)
+            Lpr0 = -speye(n_grid.Ne);
         else
             Lpr0 = -speye(x_length);
         end
     case 1 % 1st order Tikhonov
-        if 	isa(n,'Grid') % use Grid method (for partial grid support)
-            Lpr0 = n.l1;
+        if 	isa(n_grid,'Grid') % use Grid method (for partial grid support)
+            Lpr0 = n_grid.l1;
         else
-            Lpr0 = genL1(n,x_length);
+            Lpr0 = genL1(n_grid,x_length);
         end
     case 2 % 2nd order Tikhonov
-        if 	isa(n,'Grid') % use Grid method (for partial grid support)
-            Lpr0 = n.l2;
+        if 	isa(n_grid,'Grid') % use Grid method (for partial grid support)
+            Lpr0 = n_grid.l2;
         else
-            Lpr0 = genL2(n,x_length);
+            Lpr0 = genL2(n_grid,x_length);
         end
     otherwise
         disp('The specified order of Tikhonov is not available.');

+kernel/gen.m

@@ -1,22 +1,21 @@
 
 % GEN Generate A matrix describing kernel/transfer functions for DMA-PMA
-% Author:  Timothy Sipkens, 2018-11-27
+% Author:  Timothy Sipkens, 2020-02-04
 % 
 % Notes:
 %   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).
 % 
-%-------------------------------------------------------------------------%
 % Inputs:
-%   grid_b      Grid on which the data exists
+%   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
-%   varargin    Name-value pairs used in evaluating the PMA tfer. fun.
 %=========================================================================%
 
-function [A,sp] = gen(sp,d_star,grid_i,prop_pma)
+function A = gen(sp,d_star,grid_i,prop_pma)
 
 if ~exist('prop_pma','var'); prop_pma = []; end
 if isempty(prop_pma); prop_pma = kernel.prop_pma; end
@@ -31,7 +30,7 @@
 
 %-- Generate grid for intergration ---------------------------------------%
 n_i = grid_i.ne;
-N_i = prod(n_i); % length of integration vector
+N_i = grid_i.Ne; % length of integration vector
 
 r = grid_i.elements;
 m = r(:,1);

+kernel/gen_grid.m

@@ -81,6 +81,10 @@
 
 %== STEP 2: Evaluate PMA transfer function ===============================%
 disp('Computing PMA contribution:');
+pname = varargin{1}; % name of field for additional PMA field
+pval = varargin{2}; % value of field for current setpoint
+if length(pval)==1; pval = pval.*ones(n_b(2),1); end % stretch if scalar
+
 tools.textbar(0); % initiate textbar
 Lambda_mat = cell(1,n_z); % pre-allocate for speed
     % one cell entry per charge state
@@ -89,7 +93,7 @@
 
     for ii=1:n_b(1) % loop over m_star
         sp(ii) = tfer_pma.get_setpoint(...
-            prop_pma,'m_star',grid_b.edges{1}(ii).*1e-18,varargin{:});
+            prop_pma,'m_star',grid_b.edges{1}(ii).*1e-18,pname,pval(ii));
         Lambda_mat{kk}(ii,:) = kernel.tfer_pma(...
             sp(ii),m.*1e-18,...
             d.*1e-9,z_vec(kk),prop_pma)';

+kernel/gen_grid_c2.m

@@ -10,7 +10,6 @@
 %   store information on higher resolutions grids
 %   (such as those used for phantoms).
 % 
-%-------------------------------------------------------------------------%
 % Inputs:
 %   grid_b      Grid on which the data exists
 %   grid_i      Grid on which to perform integration
@@ -20,10 +19,17 @@
 
 function [A,sp] = gen_grid_c2(grid_b,grid_i,prop_pma,varargin)
 
+
+%-- Parse inputs ---------------------------------------------------------%
 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 or(isempty(varargin),length(varargin)~=2) % parse extra information for PMA
+    error('Invalid additional information for PMA setpoint.');
+end
+%-------------------------------------------------------------------------%
 
 
 %-- Parse measurement set points (b) -------------------------------------%
@@ -55,16 +61,17 @@
 n_z = length(z_vec);
 
 
-%== STEP 1: Evaluate DMA transfer function ===============================%
+%== STEP 1: Evaluate SP2 transfer function ===============================%
 %   Note: The SP2 contribution is 1D and does not depend on the charge
 %   state. It is boxcar function that takes into account discretization
 %   only. 
 disp('Computing SP2 contribution...');
 Omega_mat = sparse(n_b(1),n_i(1));% pre-allocate for speed
 for ii=1:n_b(1)
-    Omega_mat(ii,:) = ...
-        and(grid_b.nodes{1}(ii)<grid_i.edges{1},...
-        grid_b.nodes{1}(ii+1)>grid_i.edges{1});
+    Omega_mat(ii,:) = max(...
+        min(grid_i.nodes{1}(2:end),grid_b.nodes{1}(ii+1))-... % lower bound
+        max(grid_i.nodes{1}(1:(end-1)),grid_b.nodes{1}(ii))... % upper bound
+        ,0)./(grid_b.nodes{1}(ii+1)-grid_b.nodes{1}(ii)); % normalize by SP2 bin size
 end
 
 [~,jj] = max(mrbc==grid_i.edges{1},[],2);
@@ -78,6 +85,10 @@
 
 %== STEP 2: Evaluate PMA transfer function ===============================%
 disp('Computing PMA contribution:');
+pname = varargin{1}; % name of field for additional PMA field
+pval = varargin{2}; % value of field for current setpoint
+if length(pval)==1; pval = pval.*ones(n_b(2),1); end % stretch if scalar
+
 tools.textbar(0); % initiate textbar
 Lambda_mat = cell(1,n_z); % pre-allocate for speed
     % one cell entry per charge state
@@ -86,7 +97,7 @@
 
     for ii=1:n_b(2) % loop over m_star
         sp(ii) = tfer_pma.get_setpoint(...
-            prop_pma,'m_star',grid_b.edges{2}(ii).*1e-18,varargin{:});
+            prop_pma,'m_star',grid_b.edges{2}(ii).*1e-18,pname,pval(ii));
         Lambda_mat{kk}(ii,:) = kernel.tfer_pma(...
             sp(ii),m.*1e-18,d.*1e-9,...
             z_vec(kk),prop_pma)';
@@ -113,6 +124,7 @@
 disp('Completed kernel.');
 %=========================================================================%
 
+
 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
 A = sparse(A); % exploit sparse structure

+kernel/prop_pma.m

@@ -56,6 +56,8 @@
         prop.T = 293; % system temperature [K]
         prop.Q = 0.3/1000/60; % volume flow rate (m^3/s) (prev: ~1 lpm)
         prop.omega_hat = 32/33; % ratio of angular speed
+        prop.mass_mob_pref = 524;
+        prop.mass_mob_exp = 3;
 
     %-- APM parameters from Ehara et al. -------------%
     case 'Ehara'

+optimize/exp_dist_bayesf.m → +optimize/bayesf.m

@@ -1,20 +1,20 @@
 
-% EXP_DIST_BAYESF Evaluate the Bayes factor for the exponential distance prior.
+% BAYESF Evaluate the Bayes factor given Lpr0 and lambda.
 % Author: Timothy Sipkens, 2019-12-21
 %=========================================================================%
 
-function [B,F,C] = exp_dist_bayesf(A,b,lambda,Lpr,x)
+function [B,F,C] = bayesf(A,b,x,Lpr0,lambda)
 
 r = length(x);
 lx = size(A,2); % n in Thompson and Kay
 lb = size(A,1); % s in Thompson and Kay
 
-F = -1/2*(norm(A*x-b)^2 + norm(lambda.*Lpr*x)^2);
+F = -1/2*(norm(A*x-b)^2 + norm(lambda.*Lpr0*x)^2);
 
 % Gpo_inv = (A')*A+(Lpr')*Lpr;
 % C = tools.logdet(Lpr'*Lpr)/2 - tools.logdet(Gpo_inv)/2;
 
-[~,~,~,S1,S2] = gsvd(full(A),full(Lpr));
+[~,~,~,S1,S2] = gsvd(full(A),full(Lpr0));
 h = max(S1,[],2); % picks off diagonal elements of each row
 the = diag(S2);
 

+optimize/tikhonov_bayesf.m → +optimize/bayesf_precomp.m

@@ -1,9 +1,15 @@
 
-% TIKHONOV_BAYESF Evaluate the Bayes factor for the Tikhonov prior.
+% BAYESF_PRECOMP Evaluate the Bayes factor with pre-comptued GSVD.
+% This is useful in the context that Lpr0 does not change, but lambda does.
 % Author: Timothy Sipkens, 2019-12-21
 %=========================================================================%
 
-function [B,F,C] = tikhonov_bayesf(A,b,x,Lpr,lambda,order,S1,S2)
+function [B,F,C] = bayesf_precomp(A,b,x,Lpr0,lambda,S1,S2,order)
+
+if ~exist('order','var'); order = []; end
+if isempty(order); order = 0; end
+    % if not Tikhonov and Lpr0 is full rank, use order = 0
+
 
 r = length(x)-order;
 lx = size(A,2); % n in Thompson and Kay
@@ -12,10 +18,10 @@
 %-- If the gsvd was not pre-computed --%
 if ~exist('S1','var'); S1 = []; S2 = []; end
 if isempty(S1)
-    [~,~,~,S1,S2] = gsvd(full(A),full(Lpr));
+    [~,~,~,S1,S2] = gsvd(full(A),full(Lpr0));
 end
 
-F = -1/2*(norm(A*x-b)^2 + norm(lambda.*Lpr*x)^2);
+F = -1/2*(norm(A*x-b)^2 + norm(lambda.*Lpr0*x)^2);
 
 % Gpo_inv = (A')*A+lambda^2.*(Lpr')*Lpr;
 % C = (length(x)-order)*log(lambda) - tools.logdet(Gpo_inv)/2;

+optimize/exp_dist_op.m

@@ -17,7 +17,7 @@
 %   x       Regularized estimate
 %=========================================================================%
 
-function [x,lambda,out] = exp_dist_op(A,b,span,Gd,d_vec,m_vec,x_ex,xi,solver)
+function [x,lambda,out] = exp_dist_op(A,b,span,Gd,d_vec,m_vec,x_ex,xi,solver,n)
 
 
 %-- Parse inputs ---------------------------------------------%
@@ -30,10 +30,16 @@
 
 if ~exist('xi','var'); xi = []; end % if no initial x is given
 if ~exist('x_ex','var'); x_ex = []; end
+
+if ~exist('n','var'); n = []; end
+if isempty(n); n = 30; end % default number of lambda entries to consider
 %--------------------------------------------------------------%
 
 
-lambda = logspace(log10(span(1)),log10(span(2)),30);
+lambda = logspace(log10(span(1)),log10(span(2)),n);
+
+Lpr = invert.exp_dist_lpr(Gd,d_vec,m_vec); % same structure throughout
+[~,~,~,S1,S2] = gsvd(full(A),full(Lpr)); % pre-compute GSVD
 
 disp('Optimizing exponential distance regularization:');
 tools.textbar(0);
@@ -45,7 +51,7 @@
     out(ii).R12 = Gd(1,2)/sqrt(Gd(1,1)*Gd(2,2));
 
     %-- Perform inversion --------------------------%
-    [out(ii).x,~,Lpr] = invert.exp_dist(...
+    out(ii).x = invert.exp_dist(...
         A,b,lambda(ii),Gd,d_vec,m_vec,xi,solver);
 
     %-- Store ||Ax-b|| and Euclidean error ---------%
@@ -54,15 +60,17 @@
 
     %-- Compute credence, fit, and Bayes factor ----%
     [out(ii).B,out(ii).F,out(ii).C] = ...
-        optimize.exp_dist_bayesf(A,b,lambda(ii),Lpr,out(ii).x);
+        optimize.bayesf_precomp(A,b,out(ii).x,Lpr,out(ii).lambda,S1,S2,0);
+            % bypass exp_dist_bayesf, because GSVD is pre-computed
+        % optimize.exp_dist_bayesf(A,b,lambda(ii),Lpr,out(ii).x);
 
     tools.textbar((length(lambda)-ii+1)/length(lambda));
 end
 
 if ~isempty(x_ex)
     [~,ind_min] = min([out.chi]);
 else
-    ind_min = [];
+    [~,ind_min] = max([out.B]);
 end
 lambda = out(ind_min).lambda;
 x = out(ind_min).x;

+optimize/exp_dist_op1d.m

@@ -2,7 +2,7 @@
 % EXP_DIST_OP1D  Single parameter sensitivity study for exponential distance regularization.
 %=========================================================================%
 
-function [out] = exp_dist_op1d(A,b,lambda,Gd,d_vec,m_vec,x_ex,xi,solver)
+function [out] = exp_dist_op1d(A,b,lambda,Gd,d_vec,m_vec,x_ex,xi,solver,type)
 
 %-- Parse inputs ---------------------------------------------%
 if ~exist('solver','var'); solver = []; end
@@ -14,27 +14,40 @@
 
 if ~exist('xi','var'); xi = []; end % if no initial x is given
 if ~exist('x_ex','var'); x_ex = []; end
+
+if ~exist('type','var'); type = []; end % if parameter to study is not given
+if isempty(type); type = 'lmld'; end
 %--------------------------------------------------------------%
 
 
-%-- For lm/ld scaling --%
-% param = logspace(log10(0.01),log10(100),42);
+%-- Set up parameter vector ----------%
+if strcmp(type,'lmld') % lm/ld scaling
+    beta_vec = logspace(log10(0.01),log10(100),42);
+
+elseif strcmp(type,'corr') % corr. scaling
+    beta_vec = 1-[logspace(log10(1e-3),log10(1.5),26),1.85,1.95,1.97,1.99,1.999];
+
+else % if invalid parameter is specified
+    error('Invalid parameter name.')
+end
+%-------------------------------------%
 
-%-- For corr. scaling --%
-param = 1-[logspace(log10(1e-3),log10(1.5),26),1.85,1.95,1.97,1.99,1.999];
 
 disp('Optimizing exponential distance regularization:');
 tools.textbar(0);
-for ii=length(param):-1:1
+for ii=length(beta_vec):-1:1
 
-    %-- For lm/ld scaling --%
-    % out(ii).alpha = param(ii);
-    % Gd_alt = Gd.*param(ii);
+    %-- Evaluate Gd/parameters for current vector entry ------%
+    if strcmp(type,'lmld') % lm/ld scaling
+        out(ii).alpha = beta_vec(ii);
+        Gd_alt = Gd.*beta_vec(ii);
 
-    %-- For corr. scaling --%
-    out(ii).R12 = 1-param(ii);
-    Gd_12_alt = sqrt(Gd(1,1)*Gd(2,2))*param(ii);
-    Gd_alt = diag(diag(Gd))+[0,Gd_12_alt;Gd_12_alt,0];
+    elseif strcmp(type,'corr') % corr. scaling
+        out(ii).R12 = 1-beta_vec(ii);
+        Gd_12_alt = sqrt(Gd(1,1)*Gd(2,2))*beta_vec(ii);
+        Gd_alt = diag(diag(Gd))+[0,Gd_12_alt;Gd_12_alt,0];
+    end
+    %---------------------------------------------------------%
 
     %-- Store other case parameters --%
     out(ii).lambda = lambda;
@@ -52,9 +65,9 @@
 
     %-- Compute credence, fit, and Bayes factor --%
     [out(ii).B,out(ii).F,out(ii).C] = ...
-        optimize.exp_dist_bayesf(A,b,lambda,Lpr,out(ii).x);
+        optimize.bayesf(A,b,out(ii).x,Lpr,lambda);
 
-    tools.textbar((length(param)-ii+1)/length(param));
+    tools.textbar((length(beta_vec)-ii+1)/length(beta_vec));
 end
 
 

+optimize/exp_dist_opbf.m

@@ -58,7 +58,7 @@
     out(ii).corr = vec_corr(ii);
 
     [out(ii).B,out(ii).F,out(ii).C] = ...
-        optimize.exp_dist_bayesf(A,b,out(ii).lambda,Lpr,out(ii).x);
+        optimize.bayesf(A,b,out(ii).x,Lpr,out(ii).lambda);
 
     tools.textbar(ii/length(vec_lambda));
 end

+optimize/tikhonov_op.m

@@ -54,8 +54,8 @@
 
     %-- Compute credence, fit, and Bayes factor --%
     [out(ii).B,out(ii).F,out(ii).C] = ...
-        optimize.tikhonov_bayesf(...
-        A,b,out(ii).x,Lpr0,lambda(ii),order,S1,S2);
+        optimize.bayesf_precomp(...
+        A,b,out(ii).x,Lpr0,lambda(ii),S1,S2,order);
 
     tools.textbar((length(lambda)-ii+1)/length(lambda));
 end