Skip to content

Commit

Permalink
Merge pull request #10 from tsipkens/developer
Browse files Browse the repository at this point in the history 
Developer
  • Loading branch information
@tsipkens
tsipkens authored Oct 28, 2019
2 parents 862ff81 + b64f530 commit c2d4879
Show file tree
Hide file tree
Showing 64 changed files with 1,444 additions and 852 deletions.
45 changes: 25 additions & 20 deletions +invert/exp_dist.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,17 +3,17 @@
% Author: Timothy Sipkens, 2018-10-22
%=========================================================================%

function [x,D,Lx,Lpo] = exp_dist(A,b,d_vec,m_vec,lambda,Lex,x0,solver,sigma)
function [x,D,Lpr,Gpo_inv] = exp_dist(A,b,d_vec,m_vec,lambda,Lex,x0,solver,sigma)
%-------------------------------------------------------------------------%
% Inputs:
% A Model matrix
% b Data
% A Model matrix
% b Data
%
% Outputs:
% x Regularized estimate
% D Inverse operator (x = D*[b;0])
% Lx Cholesky factorization of prior covariance
% Lpo Cholesky factorization of posterior covariance
% x Regularized estimate
% D Inverse operator (x = D*[b;0])
% Lpr Cholesky factorization of prior covariance
% Gpo_inv Inverse of posterior covariance
%-------------------------------------------------------------------------%


Expand All @@ -22,7 +22,6 @@

%-- Parse inputs ---------------------------------------------%
if ~exist('solver','var'); solver = []; end
if isempty(solver); solver = 'interior-point'; end
% if computation method not specified

if ~exist('Lex','var'); Lex = []; end
Expand All @@ -39,34 +38,40 @@

d1 = log(vec_m1)-log(vec_m2);
d2 = log(vec_d1)-log(vec_d2);
d = sqrt((d1.*Lex(1,1)+d2.*Lex(1,2)).^2+(d1.*Lex(2,1)+d2.*Lex(2,2)).^2); % distance
d = sqrt((d1.*Lex(1,1)+d2.*Lex(1,2)).^2+...
(d1.*Lex(2,1)+d2.*Lex(2,2)).^2); % distance

%-- Generate prior covariance matrix --------------------------------------
Gx = exp(-d);
Gpr = exp(-d);
if exist('sigma','var') % incorporate structure into covariance, if specified
for ii=1:x_length
for jj=1:x_length
Gx(ii,jj) = Gx(ii,jj).*sigma(ii).*sigma(jj);
Gpr(ii,jj) = Gpr(ii,jj).*sigma(ii).*sigma(jj);
end
end
end
Gx(Gx<(0.05.*mean(mean(Gx)))) = 0; % remove any entries below thershold
Gx = Gx./max(max(Gx)); % normalize matrix structure
Gpr(Gpr<(0.05.*mean(mean(Gpr)))) = 0; % remove any entries below thershold
Gpr = Gpr./max(max(Gpr)); % normalize matrix structure

Gxi = inv(Gx);
Lx = chol(Gxi);
Lx = lambda.*Lx./max(max(Lx));
Lx(abs(Lx)<(0.01.*mean(mean(abs(Lx))))) = 0;
Lx = sparse(Lx);
Gpr_inv = inv(Gpr);
Lpr = chol(Gpr_inv);
clear Gpr_inv; % to save memory
Lpr = lambda.*Lpr./max(max(Lpr));
Lpr(abs(Lpr)<(0.005.*max(max(abs(Lpr))))) = 0;
[x,y] = meshgrid(1:size(Lpr,1),1:size(Lpr,2));
Lpr = Lpr.*(abs(x-y)<200); % crop distant pixels
Lpr = sparse(Lpr);


%-- Choose and execute solver --------------------------------------------%
[x,D] = invert.lsq(...
[A;Lx],[b;sparse(x_length,1)],solver,x0);
[A;Lpr],[b;sparse(x_length,1)],x0,solver);


%-- Uncertainty quantification -------------------------------------------%
Lpo = [];
if nargout>=4
Gpo_inv = A'*A+Lpr'*Lpr;
end

end

15 changes: 13 additions & 2 deletions +invert/lsq.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
% Author: Timothy Sipkens, 2019-07-17
%=========================================================================%

function [x,D] = lsq(A,b,solver,x0)
function [x,D] = lsq(A,b,x0,solver)
%-------------------------------------------------------------------------%
% Inputs:
% A Model matrix
Expand Down Expand Up @@ -35,6 +35,11 @@
end

switch solver
case 'non-neg' % constrained, iterative linear least squares
options = optimset('Display','off','TolX',eps*norm(A,1)*length(A));
x = lsqnonneg(A,b,options);
D = []; % not specified when using this method

case 'interior-point' % constrained, iterative linear least squares
options = optimoptions('lsqlin','Algorithm','interior-point','Display','none');
x = lsqlin(A,b,...
Expand All @@ -48,7 +53,13 @@
options = optimoptions('lsqlin','Algorithm','trust-region-reflective');
x = lsqlin(A,b,...
[],[],[],[],x_lb,[],x0,options);


case 'interior-point-neg'
options = optimoptions('lsqlin','Algorithm','interior-point','Display','none');
x = lsqlin(A,b,...
[],[],[],[],[],[],x0,options);
D = []; % not specified when using this method

case 'algebraic' % matrix multiplication least squares (not non-negative constrained)
D = (A'*A)\A'; % invert combined matrices
x = D*b;
Expand Down
85 changes: 42 additions & 43 deletions +invert/tikhonov.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,45 +3,44 @@
% Author: Timothy Sipkens, 2018-11-21
%=========================================================================%

function [x,D,Lx,Gpo] = tikhonov(A,b,n,lambda,order,x0,solver)
function [x,D,Lpr,Gpo_inv] = tikhonov(A,b,n,lambda,order,x0,solver)
%-------------------------------------------------------------------------%
% Inputs:
% A Model matrix
% b Data
% n Length of first dimension of solution
% lambda Regularization parameter
% order Order of regularization (Optional, default is 1)
% x0 Initial guess for solver (Optional, default is zeros)
% solver Solver (Optional, default is interior-point)
% A Model matrix
% b Data
% n Length of first dimension of solution
% lambda Regularization parameter
% order Order of regularization (Optional, default is 1)
% x0 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])
% Lx Tikhonov matrix
% Lpo Cholesky factorization of posterior covariance
% x Regularized estimate
% D Inverse operator (x = D*[b;0])
% Lpr Tikhonov matrix
% Gpo_inv Inverse of posterior covariance
%-------------------------------------------------------------------------%

x_length = length(A(1,:));

%-- Parse inputs ---------------------------------------------------------%
if ~exist('order','var'); order = []; end
if isempty(order); order = 1; end
% if order not specified

if ~exist('x0','var'); x0 = []; end % if initial guess is not specified
if ~exist('solver','var'); solver = []; end

if isempty(order); order = 1; end % if order not specified
% if isempty(x0); x0 = sparse(x_length,1); end % if initial guess is not specified
if isempty(solver); solver = 'interior-point'; end % if computation method not specified
%-------------------------------------------------------------------------%


%-- Generate Tikhonov smoothing matrix -----------------------------------%
switch order
case 0 % 0th order Tikhonov
Lx = -lambda.*speye(x_length);
Lpr = -lambda.*speye(x_length);
case 1 % 1st order Tikhonov
Lx = lambda.*genLx1(n,x_length);
Lpr = lambda.*genL1(n,x_length);
case 2 % 2nd order Tikhonov
Lx = lambda.*genLx2(n,x_length);
Lpr = lambda.*genL2(n,x_length);
otherwise
disp('The specified order of Tikhonov is not available.');
disp(' ');
Expand All @@ -51,92 +50,92 @@

%-- Choose and execute solver --------------------------------------------%
[x,D] = invert.lsq(...
[A;Lx],[b;sparse(x_length,1)],solver,x0);
[A;Lpr],[b;sparse(x_length,1)],x0,solver);


%-- Uncertainty quantification -------------------------------------------%
if nargout>=4
Gpo = A'*A+Lx'*Lx;
Gpo_inv = A'*A+Lpr'*Lpr;
end

end
%=========================================================================%


%== GENLX1 ===============================================================%
%== GENL1 ================================================================%
% Generates Tikhonov matrix for 1st order Tikhonov regularization.
function Lx = genLx1(n,x_length)
function L = genL1(n,x_length)
%-------------------------------------------------------------------------%
% Inputs:
% n Length of first dimension of solution
% x_length Length of x vector
%
% Outputs:
% Lx Tikhonov matrix
% L Tikhonov matrix
%-------------------------------------------------------------------------%

% Dx = speye(n);
% Dx = spdiag(-ones(n,1),1,Dx);
% Dx = kron(Dx,speye(x_length/n));

Lx = -eye(x_length);
L = -eye(x_length);
for jj=1:x_length
if ~(mod(jj,n)==0)
Lx(jj,jj+1) = 0.5;
L(jj,jj+1) = 0.5;
else % if on right edge
Lx(jj,jj) = Lx(jj,jj)+0.5;
L(jj,jj) = L(jj,jj)+0.5;
end

if jj<=(x_length-n)
Lx(jj,jj+n) = 0.5;
L(jj,jj+n) = 0.5;
else % if on bottom
Lx(jj,jj) = Lx(jj,jj)+0.5;
L(jj,jj) = L(jj,jj)+0.5;
end
end
Lx = sparse(Lx);
L = sparse(L);

end
%=========================================================================%


%== GENLX2 ===============================================================%
%== GENL2 ================================================================%
% Generates Tikhonov matrix for 2nd order Tikhonov regularization.
function Lx = genLx2(n,x_length)
function L = genL2(n,x_length)
% Inputs:
% n Length of first dimension of solution
% x_length Length of x vector
%
% Outputs:
% Lx Tikhonov matrix
% L Tikhonov matrix
%-------------------------------------------------------------------------%

Lx = -eye(x_length);
L = -eye(x_length);
for jj=1:x_length
if ~(mod(jj,n)==0)
Lx(jj,jj+1) = 0.25;
L(jj,jj+1) = 0.25;
else
Lx(jj,jj) = Lx(jj,jj)+0.25;
L(jj,jj) = L(jj,jj)+0.25;
end

if ~(mod(jj-1,n)==0)
Lx(jj,jj-1) = 0.25;
L(jj,jj-1) = 0.25;
else
Lx(jj,jj) = Lx(jj,jj)+0.25;
L(jj,jj) = L(jj,jj)+0.25;
end

if jj>n
Lx(jj,jj-n) = 0.25;
L(jj,jj-n) = 0.25;
else
Lx(jj,jj) = Lx(jj,jj)+0.25;
L(jj,jj) = L(jj,jj)+0.25;
end

if jj<=(x_length-n)
Lx(jj,jj+n) = 0.25;
L(jj,jj+n) = 0.25;
else
Lx(jj,jj) = Lx(jj,jj)+0.25;
L(jj,jj) = L(jj,jj)+0.25;
end
end
Lx = sparse(Lx);
L = sparse(L);

end
%=========================================================================%
Expand Down
28 changes: 15 additions & 13 deletions +kernel/gen_A.m
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@

% GEN_A Generate A matrix describing kernel/transfer functions
% GEN_A Generate A matrix describing kernel/transfer functions for DMA-PMA
% Author: Timothy Sipkens, 2018-11-27
%-------------------------------------------------------------------------%
% Notes:
Expand All @@ -9,14 +9,17 @@
% (such as those used for phantoms).
%=========================================================================%

function A = gen_A(grid_b,grid_i,varargin)
function A = gen_A(grid_b,grid_i,prop_PMA,varargin)
%-------------------------------------------------------------------------%
% Inputs:
% grid_b Grid on which the data exists
% grid_i Grid on which to perform integration
% varargin Name-value pairs used in evaluating the CPMA tfer. fun.
% varargin Name-value pairs used in evaluating the PMA tfer. fun.
%-------------------------------------------------------------------------%

if ~exist('prop_PMA','var'); prop_PMA = []; end
if isempty(prop_PMA); prop_PMA = kernel.prop_PMA('Olfert'); end

%-- Parse measurement set points (b) -------------------------------------%
r_star = grid_b.elements;
m_star = r_star(:,1);
Expand All @@ -37,13 +40,13 @@
%-- Start evaluate kernel ------------------------------------------------%
disp('Evaluating kernel...');

%-- Evaluate particle charging fractions --%
%== Evaluate particle charging fractions =================================%
z_vec = (1:3)';
f_z = sparse(kernel.tfer_charge(d.*1e-9,z_vec)); % get fraction charged for d
n_z = length(z_vec);


%-- Evaluate DMA transfer function ---------------------------------------%
%== Evaluate DMA transfer function =======================================%
%-- Note: The DMA transfer function is 1D, speeding evaluation.
Omega_mat = cell(1,n_z); % pre-allocate for speed
for kk=1:n_z
Expand All @@ -63,18 +66,17 @@
end


%-- Evaluate CPMA transfer function --------------------------------------%
prop_CPMA = kernel.prop_CPMA('Olfert');
disp('Evaluating CPMA contribution:');
%-- Evaluate PMA transfer function ---------------------------------------%
disp('Evaluating PMA contribution:');
tools.textbar(0); % initiate textbar
Lambda_mat = cell(1,n_z); % pre-allocate for speed
for kk=1:n_z
Lambda_mat{kk} = sparse(n_b(1),N_i);% pre-allocate for speed
for ii=1:n_b(1)
Lambda_mat{kk}(ii,:) = kernel.tfer_CPMA(...
Lambda_mat{kk}(ii,:) = kernel.tfer_PMA(...
grid_b.edges{1}(ii).*1e-18,m.*1e-18,...
d.*1e-9,z_vec(kk),prop_CPMA,[],varargin{:})';
% CPMA transfer function
d.*1e-9,z_vec(kk),prop_PMA,[],varargin{:})';
% PMA transfer function

if or(max(Lambda_mat{kk}(ii,:))>(1+1e-9),any(sum(Lambda_mat{kk}(ii,:))<0))
disp(' ');
Expand All @@ -90,11 +92,11 @@
disp('Compiling kernel...');
K = sparse(N_b,N_i);
for kk=1:n_z
[~,i1] = max(m_star==grid_b.edges{1},[],2); % index corresponding to CPMA setpoint
[~,i1] = max(m_star==grid_b.edges{1},[],2); % index corresponding to PMA setpoint
[~,i2] = max(d_star==grid_b.edges{2},[],2); % index correspondng to DMA setpoint

K = K+f_z(z_vec(kk),:).*... % charging contribution
Lambda_mat{kk}(i1,:).*... % CPMA contribution
Lambda_mat{kk}(i1,:).*... % PMA contribution
Omega_mat{kk}(i2,:); % DMA contribution
end
disp('Completed kernel evaluation.');
Expand Down
Loading

0 comments on commit c2d4879

Please sign in to comment.