Merge pull request #15 from tsipkens/developer

Developer

+invert/tikhonov_lpr.m

@@ -1,12 +1,12 @@
 
 % TIKHONOV_LPR Generates Tikhonov smoothing operators/matrix, L. 
 % Author:   Timothy Sipkens, 2020-02-05
-%-------------------------------------------------------------------------%
+% 
 % Inputs:
 %   order       Order of the Tikhonov operator
 %   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)
+%   x_length    Length of x vector
+%               (only used if a Grid is not specified for n_grid)
 %
 % Outputs:
 %   Lpr0        Tikhonov matrix
@@ -47,7 +47,7 @@
 
 %== GENL1 ================================================================%
 %   Generates Tikhonov matrix for 1st order Tikhonov regularization.
-%-------------------------------------------------------------------------%
+% 
 % Inputs:
 %   n           Length of first dimension of solution
 %   x_length    Length of x vector
@@ -84,7 +84,7 @@
 
 %== GENL2 ================================================================%
 %   Generates Tikhonov matrix for 2nd order Tikhonov regularization.
-%-------------------------------------------------------------------------%
+% 
 % Inputs:
 %   n           Length of first dimension of solution
 %   x_length    Length of x vector

+kernel/gen.m

@@ -107,7 +107,7 @@
 A = bsxfun(@times,K,dr_log'); % multiply kernel by element area
 A = sparse(A); % exploit sparse structure
 
-disp('Completed computing A matrix.');
+disp('Completed computing kernel matrix, <strong>A</strong>.');
 disp(' ');
 
 end

+kernel/gen_grid.m

@@ -81,19 +81,16 @@
 
 %== 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
+Lambda_mat = cell(1,n_z); % pre-allocate for speed, one cell entry per charge state
+sp = tfer_pma.get_setpoint(prop_pma,...
+    'm_star',grid_b.edges{1}.*1e-18,varargin{:}); % get PMA setpoints
+
 for kk=1:n_z % loop over the charge state
     Lambda_mat{kk} = sparse(n_b(1),N_i);% pre-allocate for speed
 
     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,pname,pval(ii));
         Lambda_mat{kk}(ii,:) = kernel.tfer_pma(...
             sp(ii),m.*1e-18,...
             d.*1e-9,z_vec(kk),prop_pma)';
@@ -124,7 +121,7 @@
 A = bsxfun(@times,K,dr_log'); % multiply kernel by element area
 A = sparse(A); % exploit sparse structure
 
-disp('Completed computing A matrix.');
+disp('Completed computing kernel matrix, <strong>A</strong>.');
 disp(' ');
 
 

+kernel/gen_grid_c2.m

@@ -85,19 +85,16 @@
 
 %== 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
+Lambda_mat = cell(1,n_z); % pre-allocate for speed, one cell entry per charge state
+sp = tfer_pma.get_setpoint(prop_pma,...
+    'm_star',grid_b.edges{2}.*1e-18,varargin{:}); % get PMA setpoints
+
 for kk=1:n_z % loop over the charge state
     Lambda_mat{kk} = sparse(n_b(1),N_i);% pre-allocate for speed
 
     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,pname,pval(ii));
         Lambda_mat{kk}(ii,:) = kernel.tfer_pma(...
             sp(ii),m.*1e-18,d.*1e-9,...
             z_vec(kk),prop_pma)';
@@ -129,7 +126,7 @@
 A = bsxfun(@times,K,dr_log'); % multiply kernel by element area
 A = sparse(A); % exploit sparse structure
 
-disp('Completed computing A matrix.');
+disp('Completed computing kernel matrix, <strong>A</strong>.');
 disp(' ');
 
 

+kernel/tfer_charge.m

@@ -3,7 +3,7 @@
 % Author: Timothy Sipkens, 2018-12-27
 % Based on code from Buckley et al., J. Aerosol Sci. (2008) and Olfert
 % laboratory at the University of Alberta.
-%-------------------------------------------------------------------------%
+% 
 % Inputs:
 %   d       Particle diameter [m]
 %   z       Integer particle charge state (optional, default = 0:6)

+optimize/bayesf.m

@@ -9,25 +9,18 @@
 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.*Lpr0*x)^2);
-
-% Gpo_inv = (A')*A+(Lpr')*Lpr;
-% C = tools.logdet(Lpr'*Lpr)/2 - tools.logdet(Gpo_inv)/2;
+F = -1/2*(norm(A*x-b)^2 + norm(lambda.*Lpr0*x)^2); % fit
 
 [~,~,~,S1,S2] = gsvd(full(A),full(Lpr0));
 h = max(S1,[],2); % picks off diagonal elements of each row
 the = diag(S2);
 
-% ct = 2*sum(log(the((lx-r+1):lb)));
-% % cp = 2*tools.logdet(P); % term cancels out between |Lpr| and |Gpo|
-% det_po = sum(log(h((lx-r+1):lb).^2 + the((lx-r+1):lb).^2));
-% C = det_pr/2 - det_po/2;
-
 ct = 2*sum(log(the((lx-r+1):lb))); % contributes to the constant term
-% cp = 2*tools.logdet(P); % term cancels out between |Lpr| and |Gpo|
 det_po = sum(log(h((lx-r+1):lb).^2 + lambda^2.*the((lx-r+1):lb).^2));
 
 C = (r+lb-lx)*log(lambda) - det_po/2 + ct/2;
+    % measurement credence
+    % ignores det_b contributions
 
 B = F+C;
 

+optimize/bayesf_b.m

@@ -0,0 +1,42 @@
+
+% BAYESF_B Evaluate the Bayes factor given Lpr0, lambda, and evolving A matrix.
+% Author: Timothy Sipkens, 2020-02-27
+%=========================================================================%
+
+function [B,F,C] = bayesf_b(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.*Lpr0*x)^2); % fit
+
+%-- Depreciated code ----------------------%
+% Gpo_inv = (A')*A+(Lpr')*Lpr;
+% C = tools.logdet(Lpr'*Lpr)/2 - tools.logdet(Gpo_inv)/2;
+%------------------------------------------%
+
+[~,~,P,S1,S2] = gsvd(full(A),full(Lpr0));
+h = max(S1,[],2); % picks off diagonal elements of each row
+the = diag(S2);
+
+%-- Depreciated code ----------------------%
+% ct = 2*sum(log(the((lx-r+1):lb)));
+% % cp = 2*tools.logdet(P); % term cancels out between |Gpr| and |Gpo|
+% det_po = sum(log(h((lx-r+1):lb).^2 + the((lx-r+1):lb).^2));
+% C = det_pr/2 - det_po/2;
+%------------------------------------------%
+
+ct = 2*sum(log(the((lx-r+1):lb))); % contributes to the constant term
+det_po = sum(log(h((lx-r+1):lb).^2 + lambda^2.*the((lx-r+1):lb).^2));
+
+cp = 2*tools.logdet(P); % contributes due to presence of |Gb|
+det_b = cp + 2.*sum(log(h(1:lb)));
+
+C = (r+lb-lx)*log(lambda) - det_po/2 + ct/2 + det_b/2;
+    % measurement credence, including |Gb| contribution
+
+B = F+C;
+
+end
+

+optimize/bayesf_precomp.m

@@ -21,19 +21,17 @@
     [~,~,~,S1,S2] = gsvd(full(A),full(Lpr0));
 end
 
-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;
+F = -1/2*(norm(A*x-b)^2 + norm(lambda.*Lpr0*x)^2); % fit
 
 h = max(S1,[],2); % picks off diagonal elements of each row
 the = diag(S2);
 
 ct = 2*sum(log(the((lx-r+1):lb))); % contributes to the constant term
-% cp = 2*tools.logdet(P); % term cancels out between |Lpr| and |Gpo|
 det_po = sum(log(h((lx-r+1):lb).^2 + lambda^2.*the((lx-r+1):lb).^2));
 
 C = (r+lb-lx)*log(lambda) - det_po/2 + ct/2;
+    % measurement credence
+    % ignores det_b contributions
 
 B = F+C;
 

+optimize/exp_dist_op.m

@@ -15,9 +15,12 @@
 %
 % Outputs:
 %   x       Regularized estimate
+%   lambda  Semi-optimal regularization parameter
+%           (against exact solution if x_ex is specified or using Bayes factor)
+%   output  Output structure with information for a range of the regularization parameter
 %=========================================================================%
 
-function [x,lambda,out] = exp_dist_op(A,b,span,Gd,d_vec,m_vec,x_ex,xi,solver,n)
+function [x,lambda,output] = exp_dist_op(A,b,span,Gd,d_vec,m_vec,x_ex,xi,solver,n)
 
 
 %-- Parse inputs ---------------------------------------------%
@@ -45,35 +48,35 @@
 tools.textbar(0);
 for ii=length(lambda):-1:1
     %-- Store case parameters ----------------------%
-    out(ii).lambda = lambda(ii);
-    out(ii).lm = sqrt(Gd(1,1));
-    out(ii).ld = sqrt(Gd(2,2));
-    out(ii).R12 = Gd(1,2)/sqrt(Gd(1,1)*Gd(2,2));
+    output(ii).lambda = lambda(ii);
+    output(ii).lm = sqrt(Gd(1,1));
+    output(ii).ld = sqrt(Gd(2,2));
+    output(ii).R12 = Gd(1,2)/sqrt(Gd(1,1)*Gd(2,2));
 
     %-- Perform inversion --------------------------%
-    out(ii).x = invert.exp_dist(...
+    output(ii).x = invert.exp_dist(...
         A,b,lambda(ii),Gd,d_vec,m_vec,xi,solver);
 
     %-- Store ||Ax-b|| and Euclidean error ---------%
-    if ~isempty(x_ex); out(ii).chi = norm(out(ii).x-x_ex); end
-    out(ii).Axb = norm(A*out(ii).x-b);
+    if ~isempty(x_ex); output(ii).chi = norm(output(ii).x-x_ex); end
+    output(ii).Axb = norm(A*output(ii).x-b);
 
     %-- Compute credence, fit, and Bayes factor ----%
-    [out(ii).B,out(ii).F,out(ii).C] = ...
-        optimize.bayesf_precomp(A,b,out(ii).x,Lpr,out(ii).lambda,S1,S2,0);
+    [output(ii).B,output(ii).F,output(ii).C] = ...
+        optimize.bayesf_precomp(A,b,output(ii).x,Lpr,output(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]);
+    [~,ind_min] = min([output.chi]);
 else
-    [~,ind_min] = max([out.B]);
+    [~,ind_min] = max([output.B]);
 end
-lambda = out(ind_min).lambda;
-x = out(ind_min).x;
+lambda = output(ind_min).lambda;
+x = output(ind_min).x;
 
 end
 

+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,type)
+function [output] = exp_dist_op1d(A,b,lambda,Gd,d_vec,m_vec,x_ex,xi,solver,type)
 
 %-- Parse inputs ---------------------------------------------%
 if ~exist('solver','var'); solver = []; end
@@ -39,33 +39,33 @@
 
     %-- Evaluate Gd/parameters for current vector entry ------%
     if strcmp(type,'lmld') % lm/ld scaling
-        out(ii).alpha = beta_vec(ii);
+        output(ii).alpha = beta_vec(ii);
         Gd_alt = Gd.*beta_vec(ii);
 
     elseif strcmp(type,'corr') % corr. scaling
-        out(ii).R12 = 1-beta_vec(ii);
+        output(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;
-    out(ii).lm = sqrt(Gd_alt(1,1));
-    out(ii).ld = sqrt(Gd_alt(2,2));
-    out(ii).R12 = Gd_alt(1,2)/sqrt(Gd_alt(1,1)*Gd_alt(2,2));
+    output(ii).lambda = lambda;
+    output(ii).lm = sqrt(Gd_alt(1,1));
+    output(ii).ld = sqrt(Gd_alt(2,2));
+    output(ii).R12 = Gd_alt(1,2)/sqrt(Gd_alt(1,1)*Gd_alt(2,2));
 
     %-- Perform inversion --%
-    [out(ii).x,~,Lpr] = invert.exp_dist(...
+    [output(ii).x,~,Lpr] = invert.exp_dist(...
         A,b,lambda,Gd_alt,d_vec,m_vec,xi,solver);
 
     %-- Store ||Ax-b|| and Euclidean error --%
-    if ~isempty(x_ex); out(ii).chi = norm(out(ii).x-x_ex); end
-    out(ii).Axb = norm(A*out(ii).x-b);
+    if ~isempty(x_ex); output(ii).chi = norm(output(ii).x-x_ex); end
+    output(ii).Axb = norm(A*output(ii).x-b);
 
     %-- Compute credence, fit, and Bayes factor --%
-    [out(ii).B,out(ii).F,out(ii).C] = ...
-        optimize.bayesf(A,b,out(ii).x,Lpr,lambda);
+    [output(ii).B,output(ii).F,output(ii).C] = ...
+        optimize.bayesf(A,b,output(ii).x,Lpr,lambda);
 
     tools.textbar((length(beta_vec)-ii+1)/length(beta_vec));
 end

+optimize/exp_dist_opbf.m

@@ -1,7 +1,7 @@
 
-% EXP_DIST_OPBF  Approximates optimal lambda for exponential distance solver by brute force method.
+% EXP_DIST_OPBF  Approximates optimal prior parameter set using the brute force method.
 % Author: Timothy Sipkens, 2019-12
-%-------------------------------------------------------------------------%
+% 
 % Inputs:
 %   A       Model matrix
 %   b       Data
@@ -15,9 +15,12 @@
 %
 % Outputs:
 %   x       Regularized estimate
+%   lambda  Semi-optimal regularization parameter
+%           (against exact solution if x_ex is specified or using Bayes factor)
+%   output  Output structure with information for a range of the prior parameters
 %=========================================================================%
 
-function [x,lambda,out] = exp_dist_opbf(A,b,d_vec,m_vec,x_ex,xi,solver)
+function [x,lambda,output] = exp_dist_opbf(A,b,d_vec,m_vec,x_ex,xi,solver)
 
 
 %-- Parse inputs ---------------------------------------------%
@@ -44,30 +47,30 @@
 vec_corr = vec_corr(:);
 
 tools.textbar(0);
-out(length(vec_lambda)).chi = [];
+output(length(vec_lambda)).chi = [];
 for ii=1:length(vec_lambda)
     y = [vec_lambda(ii),vec_ratio(ii),vec_ld(ii),vec_corr(ii)];
 
-    out(ii).x = invert.exp_dist(...
+    output(ii).x = invert.exp_dist(...
         A,b,y(1),Gd_fun(y),d_vec,m_vec,xi,solver);
-    out(ii).chi = norm(out(ii).x-x_ex);
+    output(ii).chi = norm(output(ii).x-x_ex);
 
-    out(ii).lambda = vec_lambda(ii);
-    out(ii).ratio = vec_ratio(ii);
-    out(ii).ld = vec_ld(ii);
-    out(ii).corr = vec_corr(ii);
+    output(ii).lambda = vec_lambda(ii);
+    output(ii).ratio = vec_ratio(ii);
+    output(ii).ld = vec_ld(ii);
+    output(ii).corr = vec_corr(ii);
 
-    [out(ii).B,out(ii).F,out(ii).C] = ...
-        optimize.bayesf(A,b,out(ii).x,Lpr,out(ii).lambda);
+    [output(ii).B,output(ii).F,output(ii).C] = ...
+        optimize.bayesf(A,b,output(ii).x,Lpr,output(ii).lambda);
 
     tools.textbar(ii/length(vec_lambda));
 end
 
 tools.textbar(1);
 
-[~,ind_min] = min([out(ii).chi]);
-x = [];%out(ind_min).x;
-lambda = [];%out(ind_min).lambda;
+[~,ind_min] = min([output(ii).chi]);
+x = out(ind_min).x;
+lambda = out(ind_min).lambda;
 
 end