Skip to content
Snippets Groups Projects
Commit ea91ef25 authored by mingf2's avatar mingf2
Browse files

remove dependency on SignalProc. & Optim. Toolbox

parent f4b44ca1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
LMFsolve/LMFsolve.m 0 → 100644
+ 293
− 0
View file @ ea91ef25
function [xf, S, cnt] = LMFsolve(varargin)
% LMFSOLVE Solve a Set of Nonlinear Equations in Least-Squares Sense.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A solution is obtained by a shortened Fletcher version of the
% Levenberg-Maquardt algoritm for minimization of a sum of squares
% of equation residuals.
%
% [Xf, Ssq, CNT] = LMFsolve(FUN,Xo,Options)
% FUN is a function handle or a function M-file name that evaluates
% m-vector of equation residuals,
% Xo is n-vector of initial guesses of solution,
% Options is an optional set of Name/Value pairs of control parameters
% of the algorithm. It may be also preset by calling:
% Options = LMFsolve('default'), or by a set of Name/Value pairs:
% Options = LMFsolve('Name',Value, ... ), or updating the Options
% set by calling
% Options = LMFsolve(Options,'Name',Value, ...).
%
% Name Values {default} Description
% 'Display' integer Display iteration information
% {0} no display
% k display initial and every k-th iteration;
% 'FunTol' {1e-7} norm(FUN(x),1) stopping tolerance;
% 'XTol' {1e-7} norm(x-xold,1) stopping tolerance;
% 'MaxIter' {100} Maximum number of iterations;
% 'ScaleD' Scale control:
% value D = eye(m)*value;
% vector D = diag(vector);
% {[]} D(k,k) = JJ(k,k) for JJ(k,k)>0, or
% = 1 otherwise,
% where JJ = J.'*J
% Not defined fields of the Options structure are filled by default values.
%
% Output Arguments:
% Xf final solution approximation
% Ssq sum of squares of residuals
% Cnt >0 count of iterations
% -MaxIter, did not converge in MaxIter iterations
% Example: Rosenbrock valey inside circle with unit diameter
% R = @(x) sqrt(x'*x)-.5; % A distance from the radius r=0.5
% ros= @(x) [ 10*(x(2)-x(1)^2); 1-x(1); (R(x)>0)*R(x)*1000];
% [x,ssq,cnt]=LMFsolve(ros,[-1.2,1],'Display',1,'MaxIter',50)
% returns x = [0.4556; 0.2059], ssq = 0.2966, cnt = 18.
%
% Note: Users with old MATLAB versions (<7), which have no anonymous
% functions implemented, should call LMFsolve with named function for
% residuals. For above example it is
% [x,ssq,cnt]=LMFsolve('rosen',[-1.2,1]);
% where the function rosen.m is of the form
% function r = rosen(x)
%% Rosenbrock valey with a constraint
% R = sqrt(x(1)^2+x(2)^2)-.5;
%% Residuals:
% r = [ 10*(x(2)-x(1)^2) % first part
% 1-x(1) % second part
% (R>0)*R*1000. % penalty
% ];
% Reference:
% Fletcher, R., (1971): A Modified Marquardt Subroutine for Nonlinear Least
% Squares. Rpt. AERE-R 6799, Harwell
% Miroslav Balda,
% balda AT cdm DOT cas DOT cz
% 2007-07-02 v 1.0
% 2008-12-22 v 1.1 * Changed name of the function in LMFsolv
% * Removed part with wrong code for use of analytical
% form for assembling of Jacobian matrix
% 2009-01-08 v 1.2 * Changed subfunction printit.m for better one, and
% modified its calling from inside LMFsolve.
% * Repaired a bug, which caused an inclination to
% istability, in charge of slower convergence.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OPTIONS
%%%%%%%
% Default Options
if nargin==1 && strcmpi('default',varargin(1))
xf.Display = 0; % no print of iterations
xf.MaxIter = 100; % maximum number of iterations allowed
xf.ScaleD = []; % automatic scaling by D = diag(diag(J'*J))
xf.FunTol = 1e-7; % tolerace for final function value
xf.XTol = 1e-4; % tolerance on difference of x-solutions
return
% Updating Options
elseif isstruct(varargin{1}) % Options=LMFsolve(Options,'Name','Value',...)
if ~isfield(varargin{1},'Display')
error('Options Structure not correct for LMFsolve.')
end
xf=varargin{1}; % Options
for i=2:2:nargin-1
name=varargin{i}; % Option to be updated
if ~ischar(name)
error('Parameter Names Must be Strings.')
end
name=lower(name(isletter(name)));
value=varargin{i+1}; % value of the option
if strncmp(name,'d',1), xf.Display = value;
elseif strncmp(name,'f',1), xf.FunTol = value(1);
elseif strncmp(name,'x',1), xf.XTol = value(1);
elseif strncmp(name,'m',1), xf.MaxIter = value(1);
elseif strncmp(name,'s',1), xf.ScaleD = value;
else disp(['Unknown Parameter Name --> ' name])
end
end
return
% Pairs of Options
elseif ischar(varargin{1}) % check for Options=LMFSOLVE('Name',Value,...)
Pnames=char('display','funtol','xtol','maxiter','scaled');
if strncmpi(varargin{1},Pnames,length(varargin{1}))
xf=LMFsolve('default'); % get default values
xf=LMFsolve(xf,varargin{:});
return
end
end
% LMFSOLVE(FUN,Xo,Options)
%%%%%%%%%%%%%%%%%%%%%%%%
FUN=varargin{1}; % function handle
if ~(isvarname(FUN) || isa(FUN,'function_handle'))
error('FUN Must be a Function Handle or M-file Name.')
end
xc=varargin{2}; % Xo
if nargin>2 % OPTIONS
if isstruct(varargin{3})
options=varargin{3};
else
if ~exist('options','var')
options = LMFsolve('default');
end
for i=3:2:size(varargin,2)-1
options=LMFsolve(options, varargin{i},varargin{i+1});
end
end
else
if ~exist('options','var')
options = LMFsolve('default');
end
end
x = xc(:);
lx = length(x);
r = feval(FUN,x); % Residuals at starting point
%~~~~~~~~~~~~~~~~~
S = r'*r;
epsx = options.XTol(:);
epsf = options.FunTol(:);
if length(epsx)<lx, epsx=epsx*ones(lx,1); end
J = finjac(FUN,r,x,epsx);
%~~~~~~~~~~~~~~~~~~~~~~~
nfJ = 2;
A = J.'*J; % System matrix
v = J.'*r;
D = options.ScaleD;
if isempty(D)
D = diag(diag(A)); % automatic scaling
for i = 1:lx
if D(i,i)==0, D(i,i)=1; end
end
else
if numel(D)>1
D = diag(sqrt(abs(D(1:lx)))); % vector of individual scaling
else
D = sqrt(abs(D))*eye(lx); % scalar of unique scaling
end
end
Rlo = 0.25;
Rhi = 0.75;
l=1; lc=.75; is=0;
cnt = 0;
ipr = options.Display;
printit(ipr,-1); % Table header
d = options.XTol; % vector for the first cycle
maxit = options.MaxIter; % maximum permitted number of iterations
while cnt<maxit && ... % MAIN ITERATION CYCLE
any(abs(d) >= epsx) && ... %%%%%%%%%%%%%%%%%%%%
any(abs(r) >= epsf)
d = (A+l*D)\v; % negative solution increment
xd = x-d;
rd = feval(FUN,xd);
% ~~~~~~~~~~~~~~~~~~~
nfJ = nfJ+1;
Sd = rd.'*rd;
dS = d.'*(2*v-A*d); % predicted reduction
R = (S-Sd)/dS;
if R>Rhi % halve lambda if R too high
l = l/2;
if l<lc, l=0; end
elseif R<Rlo % find new nu if R too low
nu = (Sd-S)/(d.'*v)+2;
if nu<2
nu = 2;
elseif nu>10
nu = 10;
end
if l==0
lc = 1/max(abs(diag(inv(A))));
l = lc;
nu = nu/2;
end
l = nu*l;
end
cnt = cnt+1;
if ipr~=0 && (rem(cnt,ipr)==0 || cnt==1) % print iteration?
printit(ipr,cnt,nfJ,S,x,d,l,lc)
end
if Sd<S
S = Sd;
x = xd;
r = rd;
J = finjac(FUN,r,x,epsx);
% ~~~~~~~~~~~~~~~~~~~~~~~~~
nfJ = nfJ+1;
A = J'*J;
v = J'*r;
end
end % while
xf = x; % final solution
if cnt==maxit
cnt = -cnt;
end % maxit reached
rd = feval(FUN,xf);
nfJ = nfJ+1;
Sd = rd.'*rd;
if ipr, disp(' '), end
printit(ipr,cnt,nfJ,Sd,xf,d,l,lc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FINJAC numerical approximation to Jacobi matrix
% %%%%%%
function J = finjac(FUN,r,x,epsx)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
lx=length(x);
J=zeros(length(r),lx);
for k=1:lx
dx=.25*epsx(k);
xd=x;
xd(k)=xd(k)+dx;
rd=feval(FUN,xd);
% ~~~~~~~~~~~~~~~~
J(:,k)=((rd-r)/dx);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function printit(ipr,cnt,res,SS,x,dx,l,lc)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Printing of intermediate results
% ipr < 0 do not print lambda columns
% = 0 do not print at all
% > 0 print every (ipr)th iteration
% cnt = -1 print out the header
% 0 print out second row of results
% >0 print out first row of results
if ipr~=0
if cnt<0 % table header
disp('')
disp(char('*'*ones(1,75)))
fprintf(' itr nfJ SUM(r^2) x dx');
if ipr>0
fprintf(' l lc');
end
fprintf('\n');
disp(char('*'*ones(1,75)))
disp('')
else % iteration output
if rem(cnt,ipr)==0
f='%12.4e ';
if ipr>0
fprintf(['%4.0f %4.0f ' f f f f f '\n'],...
cnt,res,SS, x(1),dx(1),l,lc);
else
fprintf(['%4.0f %4.0f ' f f f '\n'],...
cnt,res,SS, x(1),dx(1));
end
for k=2:length(x)
fprintf([blanks(23) f f '\n'],x(k),dx(k));
end
end
end
end
LMFsolve/LMFsolveOLD.m 0 → 100644
+ 301
− 0
View file @ ea91ef25
function [xf, S, cnt] = LMFsolveOLD(varargin)
% Solve a Set of Overdetermined Nonlinear Equations in Least-Squares Sense.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A solution is obtained by a Fletcher version of the Levenberg-Maquardt
% algoritm for minimization of a sum of squares of equation residuals.
%
% [Xf, Ssq, CNT] = LMFsolveOLD(FUN,Xo,Options)
% FUN is a function handle or a function M-file name that evaluates
% m-vector of equation residuals,
% Xo is n-vector of initial guesses of solution,
% Options is an optional set of Name/Value pairs of control parameters
% of the algorithm. It may be also preset by calling:
% Options = LMFsolveOLD('default'), or by a set of Name/Value pairs:
% Options = LMFsolveOLD('Name',Value, ... ), or updating the Options
% set by calling
% Options = LMFsolveOLD(Options,'Name',Value, ...).
%
% Name Values {default} Description
% 'Display' integer Display iteration information
% {0} no display
% k display initial and every k-th iteration;
% 'Jacobian' handle Jacobian matrix function handle; {@finjac}
% 'FunTol' {1e-7} norm(FUN(x),1) stopping tolerance;
% 'XTol' {1e-7} norm(x-xold,1) stopping tolerance;
% 'MaxIter' {100} Maximum number of iterations;
% 'ScaleD' Scale control:
% value D = eye(m)*value;
% vector D = diag(vector);
% {[]} D(k,k) = JJ(k,k) for JJ(k,k)>0, or
% = 1 otherwise,
% where JJ = J.'*J
% Not defined fields of the Options structure are filled by default values.
%
% Output Arguments:
% Xf final solution approximation
% Ssq sum of squares of residuals
% Cnt >0 count of iterations
% -MaxIter, did not converge in MaxIter iterations
% Example:
% The general Rosenbrock's function has the form
% f(x) = 100(x(1)-x(2)^2)^2 + (1-x(1))^2
% Optimum solution gives f(x)=0 for x(1)=x(2)=1. Function f(x) can be
% expressed in the form
% f(x) = f1(x)^2 =f2(x)^2,
% where f1(x) = 10(x(1)-x(2)^2) and f2(x) = 1-x(1).
% Values of the functions f1(x) and f2(x) can be used as residuals.
% LMFsolveOLD finds the solution of this problem in 5 iterations. The more
% complicated problem sounds:
% Find the least squares solution of the Rosenbrock valey inside a circle
% of the unit diameter centered at the origin. It is necessary to build
% third function, which is zero inside the circle and increasing outside it.
% This property has, say, the next function:
% f3(x) = sqrt(x(1)^2 + x(2)^2) - r, where r is a radius of the circle.
% Its implementation using anonymous functions has the form
% R = @(x) sqrt(x'*x)-.5; % A distance from the radius r=0.5
% ros= @(x) [10*(x(2)-x(1)^2); 1-x(1); (R(x)>0)*R(x)*1000];
% [x,ssq,cnt]=LMFsolveOLD(ros,[-1.2,1],'Display',1,'MaxIter',50)
% Solution: x = [0.4556; 0.2059], |x| = 0.5000
% sum of squares: ssq = 0.2966,
% number of iterations: cnt = 18.
%
% Note:
% Users with older MATLAB versions, which have no anonymous functions
% implemented, have to call LMFsolveOLD with named function for residuals.
% For above example it is
%
% [x,ssq,cnt]=LMFsolveOLD('rosen',[-1.2,1]);
%
% where the function rosen.m is of the form
%
% function r = rosen(x)
%% Rosenbrock valey with a constraint
% R = sqrt(x(1)^2+x(2)^2)-.5;
%% Residuals:
% r = [ 10*(x(2)-x(1)^2) % first part
% 1-x(1) % second part
% (R>0)*R*1000. % penalty
% ];
% Reference:
% Fletcher, R., (1971): A Modified Marquardt Subroutine for Nonlinear Least
% Squares. Rpt. AERE-R 6799, Harwell
% M. Balda,
% Institute of Thermomechanics,
% Academy of Sciences of The Czech Republic,
% balda AT cdm DOT cas DOT cz
% 2007-07-02
% 2007-10-08 formal changes, improved description
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OPTIONS
%%%%%%%
% Default Options
if nargin==1 && strcmpi('default',varargin(1))
xf.Display = 0; % no print of iterations
xf.Jacobian = @finjac; % finite difference Jacobian approximation
xf.MaxIter = 100; % maximum number of iterations allowed
xf.ScaleD = []; % automatic scaling by D = diag(diag(J'*J))
xf.FunTol = 1e-7; % tolerace for final function value
xf.XTol = 1e-4; % tolerance on difference of x-solutions
return
% Updating Options
elseif isstruct(varargin{1}) % Options=LMFsolveOLD(Options,'Name','Value',...)
if ~isfield(varargin{1},'Jacobian')
error('Options Structure not Correct for LMFsolveOLD.')
end
xf=varargin{1}; % Options
for i=2:2:nargin-1
name=varargin{i}; % option to be updated
if ~ischar(name)
error('Parameter Names Must be Strings.')
end
name=lower(name(isletter(name)));
value=varargin{i+1}; % value of the option
if strncmp(name,'d',1), xf.Display = value;
elseif strncmp(name,'f',1), xf.FunTol = value(1);
elseif strncmp(name,'x',1), xf.XTol = value(1);
elseif strncmp(name,'j',1), xf.Jacobian = value;
elseif strncmp(name,'m',1), xf.MaxIter = value(1);
elseif strncmp(name,'s',1), xf.ScaleD = value;
else disp(['Unknown Parameter Name --> ' name])
end
end
return
% Pairs of Options
elseif ischar(varargin{1}) % check for Options=LMFsolveOLD('Name',Value,...)
Pnames=char('display','funtol','xtol','jacobian','maxiter','scaled');
if strncmpi(varargin{1},Pnames,length(varargin{1}))
xf=LMFsolveOLD('default'); % get default values
xf=LMFsolveOLD(xf,varargin{:});
return
end
end
% LMFSOLVEOLD(FUN,Xo,Options)
%%%%%%%%%%%%%%%%%%%%%%%%%%%
FUN=varargin{1}; % function handle
if ~(isvarname(FUN) || isa(FUN,'function_handle'))
error('FUN Must be a Function Handle or M-file Name.')
end
xc=varargin{2}; % Xo
if nargin>2 % OPTIONS
if isstruct(varargin{3})
options=varargin{3};
else
if ~exist('options','var')
options = LMFsolveOLD('default');
end
for i=3:2:size(varargin,2)-1
options=LMFsolveOLD(options, varargin{i},varargin{i+1});
end
end
else
if ~exist('options','var')
options = LMFsolveOLD('default');
end
end
x=xc(:);
lx=length(x);
r=feval(FUN,x); % Residuals at starting point
%~~~~~~~~~~~~~~
S=r'*r;
epsx=options.XTol(:);
epsf=options.FunTol(:);
if length(epsx)<lx, epsx=epsx*ones(lx,1); end
J=options.Jacobian(FUN,r,x,epsx);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
A=J.'*J; % System matrix
v=J.'*r;
D = options.ScaleD;
if isempty(D)
D=diag(diag(A)); % automatic scaling
for i=1:lx
if D(i,i)==0, D(i,i)=1; end
end
else
if numel(D)>1
D=diag(sqrt(abs(D(1:lx)))); % vector of individual scaling
else
D=sqrt(abs(D))*eye(lx); % scalar of unique scaling
end
end
Rlo=0.25; Rhi=0.75;
l=1; lc=.75; is=0;
cnt=0;
ipr=options.Display;
printit(-1); % Table header
d=options.XTol; % vector for the first cycle
maxit = options.MaxIter; % maximum permitted number of iterations
while cnt<maxit && ... % MAIN ITERATION CYCLE
any(abs(d)>=epsx) && ... %%%%%%%%%%%%%%%%%%%%
any(abs(r)>=epsf)
d=(A+l*D)\v; % negative solution increment
xd=x-d;
rd=feval(FUN,xd);
% ~~~~~~~~~~~~~~~~~
Sd=rd.'*rd;
dS=d.'*(2*v-A*d); % predicted reduction
R=(S-Sd)/dS;
if R>Rhi
l=l/2;
if l<lc, l=0; end
elseif R<Rlo
nu=(Sd-S)/(d.'*v)+2;
if nu<2
nu=2;
elseif nu>10
nu=10;
end
if l==0
lc=1/max(abs(diag(inv(A))));
l=lc;
nu=nu/2;
end
l=nu*l;
end
cnt=cnt+1;
if ipr>0 && (rem(cnt,ipr)==0 || cnt==1)
printit(cnt,[S,l,R,x(:).'])
printit(0,[lc,d(:).'])
end
if Sd>S && is==0
is=1;
St=S; xt=x; rt=r; Jt=J; At=A; vt=v;
end
if Sd<S || is>0
S=Sd; x=xd; r=rd;
J=options.Jacobian(FUN,r,x,epsx);
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
A=J.'*J; v=J.'*r;
else
S=St; x=xt; r=rt; J=Jt; A=At; v=vt;
end
if Sd<S, is=0; end
end
xf = x; % finat solution
if cnt==maxit, cnt=-cnt; end % maxit reached
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PRINTIT Printing of intermediate results
% %%%%%%%
function printit(cnt,mx)
% ~~~~~~~~~~~~~~~
% cnt = -1 print out the header
% 0 print out second row of results
% >0 print out first row of results
if ipr>0
if cnt<0 % table header
disp('')
disp(char('*'*ones(1,100)))
disp([' cnt SUM(r^2) l R' blanks(21) 'x(i) ...'])
disp([blanks(24) 'lc' blanks(32) 'dx(i) ...'])
disp(char('*'*ones(1,100)))
disp('')
else % iteration output
if cnt>0 || rem(cnt,ipr)==0
f='%12.4e ';
form=[f f f f '\n' blanks(49)];
if cnt>0
fprintf(['%4.0f ' f f f ' x = '],[cnt,mx(1:3)])
fprintf(form,mx(4:length(mx)))
else
fprintf([blanks(18) f blanks(13) 'dx = '], mx(1))
fprintf(form,mx(2:length(mx)))
end
fprintf('\n')
end
end
end
end
% FINJAC numerical approximation to Jacobi matrix
% %%%%%%
function J = finjac(FUN,r,x,epsx)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
lx=length(x);
J=zeros(length(r),lx);
for k=1:lx
dx=.25*epsx(k);
xd=x;
xd(k)=xd(k)+dx;
rd=feval(FUN,xd);
% ~~~~~~~~~~~~~~~~
J(:,k)=((rd-r)/dx);
end
end
end
\ No newline at end of file
LMFsolve/LMFsolvetest.m 0 → 100644
+ 40
− 0
View file @ ea91ef25
% LMFsolvetest Test of the function LMFtest
% ~~~~~~~~~~~~ { default = Rosenbrock } 2009-01-08
% x = inp('Starting point [xo]',[-1.2; 1]);
x = input('Starting point [xo]') ;
if isempty(x)
x = [-1.2; 1];
end
x = x(:);
lx = length(x);
% fun = inp('Function name ','rosen');
fun = input('Function name ') ;
if isempty(fun)
fun = 'rosen';
end
% D = eval(inp('vektor vah [D]','1'));
D = input('vektor vah [D]') ;
if isempty(D)
D = 1;
end
% ipr = inp('Print control ipr',1);
ipr = input('Print control ipr') ;
if isempty(ipr)
ipr = 1;
end
resid = eval(['@' fun]); % function handle
options = LMFsolve('default');
options = LMFsolve...
(options,...
'XTol',1e-6,...
'FTol',1e-12,...
'ScaleD',D,...
'Display',ipr...
);
[x,S,cnt]=LMFsolve(resid,x,options);
% ~~~~~~~~~~~~~~~~~~~~~~~~~
fprintf('\n radius = %15.8e\n', sqrt(x'*x));
LMFsolve/license.txt 0 → 100644
+ 24
− 0
View file @ ea91ef25
Copyright (c) 2007, Miroslav Balda
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
LMFsolve/rosen.m 0 → 100644
+ 12
− 0
View file @ ea91ef25
function r = rosen(x)
%%%%%%%%%%%%%%%%%%%%% Rosenbrock valey with a constraint
R = sqrt(x(1)^2+x(2)^2)-.5;
%R = sqrt(x(1)^2+x(2)^2)-1;
%R = sqrt(x(1)^2+x(2)^2)-2;
% Residuals:
r = [ 10*(x(2)-x(1)^2) % first part
1-x(1) % second part
(R>0)*R*500. % penalty
];
extractTimeStamp.m
+ 5
− 1
View file @ ea91ef25
...@@ -45,7 +45,11 @@ function pulseHeights = get_pulse_heights(voltagePulses, method, window_len, win ...@@ -45,7 +45,11 @@ function pulseHeights = get_pulse_heights(voltagePulses, method, window_len, win
if strcmp(window, 'flat') if strcmp(window, 'flat')
w = ones(window_len, 1); w = ones(window_len, 1);
else else
w = hann(window_len); % w = hann(window_len);
w = zeros(window_len, 1);
for n = 1:window_len
w(n) = 0.5 * (1 - cos(2*pi*(n-1)/(window_len-1)));
end
end end
pulseHeights = get_filtered_pulse_heights(voltagePulses, window_len, w / sum(w)); pulseHeights = get_filtered_pulse_heights(voltagePulses, window_len, w / sum(w));
end end
......
getRossiAlpha.m
+ 62
− 4
View file @ ea91ef25
...@@ -13,7 +13,7 @@ TME = mergedTimestamps_us(end) / 1e6; % seconds ...@@ -13,7 +13,7 @@ TME = mergedTimestamps_us(end) / 1e6; % seconds
[counts, bins] = getTimeDist(mergedTimestamps_us); [counts, bins] = getTimeDist(mergedTimestamps_us);
centers = (bins(2:end) + bins(1:end-1)) / 2; centers = (bins(2:end) + bins(1:end-1)) / 2;
fit_Rossi_alpha(centers(:), counts(:)./TME); fit_Rossi_alpha(@model, centers(:), counts(:)./TME);
% Define functions % Define functions
function [counts, bins] = getTimeDist(lst) function [counts, bins] = getTimeDist(lst)
...@@ -45,9 +45,11 @@ function y = model(p, xdata) ...@@ -45,9 +45,11 @@ function y = model(p, xdata)
y = p(1) * exp(-xdata ./ p(2)) + p(3); y = p(1) * exp(-xdata ./ p(2)) + p(3);
end end
function fit_Rossi_alpha(xdata, ydata) function fit_Rossi_alpha(model, xdata, ydata)
% Curve fitting % Curve fitting requring Optimization Toolbox
popt = lsqcurvefit(@model, [1, 20, 0.1], xdata, ydata./max(ydata)); % popt = lsqcurvefit(@model, [1, 20, 0.1], xdata, ydata./max(ydata));
popt_init = my_init_fit(xdata, ydata./max(ydata));
popt = LMF_fit(@model, xdata, ydata./max(ydata), popt_init);
disp(['best fit parameters = ', num2str(popt)]); disp(['best fit parameters = ', num2str(popt)]);
fig = figure; fig = figure;
...@@ -62,3 +64,59 @@ function fit_Rossi_alpha(xdata, ydata) ...@@ -62,3 +64,59 @@ function fit_Rossi_alpha(xdata, ydata)
legend(ax, 'Simulated data', sprintf('Fit: die-away time=%5.2f µs', popt(2))); legend(ax, 'Simulated data', sprintf('Fit: die-away time=%5.2f µs', popt(2)));
xlim(ax, [0, 300]); xlim(ax, [0, 300]);
end end
function popt = my_init_fit(xdata, ydata)
% Perform regression of functions of the form y = a + b exp(c x).
%
% Syntax: popt = my_exponential_fit(xdata, ydata)
%
% Ref: https://scikit-guess.readthedocs.io/en/latest/appendices/reei/translation.html#reei2
n = length(xdata);
S = zeros(n, 1);
for k = 2:n
S(k) = S(k-1) + 0.5 * (ydata(k)+ydata(k-1)) * (xdata(k) - xdata(k-1));
end
M11 = dot(xdata-xdata(1), xdata-xdata(1));
M12 = dot(xdata-xdata(1), S);
M21 = M12;
M22 = dot(S, S);
M_inv = 1/(M11*M22-M21*M12) .* [M22 -M12; -M21 M11];
V1 = dot(ydata-ydata(1), xdata-xdata(1));
V2 = dot(ydata-ydata(1), S);
V = [V1; V2];
U = M_inv * V;
c2 = U(2);
theta = exp(c2 * xdata);
N11 = n;
N12 = sum(theta);
N21 = N12;
N22 = dot(theta, theta);
N_inv = 1/(N11*N22-N21*N12) .* [N22 -N12; -N21 N11];
K = [sum(ydata); dot(ydata, theta)];
P = N_inv * K;
popt = [P(2) -1/c2 P(1)];
end
function popt = LMF_fit(model, xdata, ydata, init_guess)
% Perform non-linear least-square fit using the Levenberg-Marquardt-Fletcher algorithm.
%
% Syntax: popt = LMF_fit(model, xdata, ydata, init_guess)
%
addpath('./LMFsolve');
D = 1;
ipr = 0; % print control of interation results
% resid = eval(['@' fun]); % function handle
resid = @(p) model(p, xdata) - ydata;
options = LMFsolve('default');
options = LMFsolve...
(options,...
'XTol',1e-6,...
'FTol',1e-12,...
'ScaleD',D,...
'Display',ipr...
);
[x,S,cnt]=LMFsolve(resid, init_guess, options);
disp(['sum of residual squares = ', num2str(S)]);
% ~~~~~~~~~~~~~~~~~~~~~~~~~
popt = x';
end
python_code/get_rossi_alpha.py
+ 24
− 1
View file @ ea91ef25
...@@ -3,7 +3,7 @@ Description: Plot Ross-alpha based on list of time stamps ...@@ -3,7 +3,7 @@ Description: Plot Ross-alpha based on list of time stamps
Author: Ming Fang Author: Ming Fang
Date: 2021-09-11 23:53:54 Date: 2021-09-11 23:53:54
LastEditors: Ming Fang LastEditors: Ming Fang
LastEditTime: 2023-06-07 19:47:07 LastEditTime: 2023-06-08 10:52:53
''' '''
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
...@@ -52,14 +52,37 @@ def model(x, a, b, c): ...@@ -52,14 +52,37 @@ def model(x, a, b, c):
return a * np.exp(-x / b) + c return a * np.exp(-x / b) + c
def my_exponential_fit(xdata, ydata):
n = len(xdata)
S = np.zeros(n)
for k in range(1, n):
S[k] = S[k-1] + 0.5 * (ydata[k]+ydata[k-1]) * (xdata[k]-xdata[k-1])
# solve the regression system: V = M U
K = np.array([xdata-xdata[0], S])
M = np.matmul(K, K.T)
V = np.matmul(K, ydata - ydata[0])
U = np.linalg.solve(M, V)
c2 = U[1]
c2 = -1/22.32
theta = np.exp(c2 * xdata)
K = np.array([np.ones(n), theta])
M = np.matmul(K, K.T)
V = np.matmul(K, ydata)
P = np.linalg.solve(M, V)
return [P[1], -1/c2, P[0]]
def fit_Rossi_alpha(xdata, ydata): def fit_Rossi_alpha(xdata, ydata):
popt, pcov = curve_fit(model, xdata, ydata/np.max(ydata)) popt, pcov = curve_fit(model, xdata, ydata/np.max(ydata))
# popt = my_exponential_fit(xdata, ydata/np.max(ydata))
print("best fit parameters = ", popt) print("best fit parameters = ", popt)
fig, ax = plt.subplots(1, 1, figsize=(6, 6)) fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.scatter(xdata, ydata, s=4, label='Simulated data') ax.scatter(xdata, ydata, s=4, label='Simulated data')
ax.plot(xdata, np.max(ydata)*model(xdata, *popt), 'r', ax.plot(xdata, np.max(ydata)*model(xdata, *popt), 'r',
label=r'Fit: die-away time=%5.2f$\pm$%5.2f$\mu$s' % tuple([popt[1], np.sqrt(np.diag(pcov))[1]])) label=r'Fit: die-away time=%5.2f$\pm$%5.2f$\mu$s' % tuple([popt[1], np.sqrt(np.diag(pcov))[1]]))
# ax.plot(xdata, np.max(ydata)*model(xdata, *popt), 'r',
# label=r'Fit: die-away time=%5.2f $\mu$s' % tuple([popt[1]]))
ax.set_xlim([0, 300]) ax.set_xlim([0, 300])
# ax.set_ylim([0.5, 5.8]) # ax.set_ylim([0.5, 5.8])
......
readme.md
+ 6
− 1
View file @ ea91ef25
...@@ -11,14 +11,19 @@ ...@@ -11,14 +11,19 @@
│ │ │ └── DataR_CH*@DT5730_1770_ananlog_100MB.BIN: Data acquired from anaglog channel in a Cf-252 measurement. Truncated to reduce size to 100MB. │ │ │ └── DataR_CH*@DT5730_1770_ananlog_100MB.BIN: Data acquired from anaglog channel in a Cf-252 measurement. Truncated to reduce size to 100MB.
| | ├── ttl | | ├── ttl
│ │ │ └── DataR_CH*@DT5730_1770_ttl_100MB.BIN: Data acquired from ttl channel in a Cf-252 measurement. Truncated to reduce size to 100MB. │ │ │ └── DataR_CH*@DT5730_1770_ttl_100MB.BIN: Data acquired from ttl channel in a Cf-252 measurement. Truncated to reduce size to 100MB.
├── LMFsolve: external library for performing non-linear least-square fit.
└── unittest: scripts and data for unit testing. └── unittest: scripts and data for unit testing.
## Dependencies
No Matalab Toolbox is needed.
## How to run the test: ## How to run the test:
1. Run `extractTimeStamp.m` 1. Run `extractTimeStamp.m`
A mat file `mergedtimestamps.mat` is generated in folder `output_matlab/analog`. A mat file `mergedtimestamps.mat` is generated in folder `output_matlab/analog`.
2. Run `getRossiAlpha.m` 2. Run `getRossiAlpha.m`
The Rossi-alpha plot and the fit show up on your screen. The die-away time found thorugh the fit is 22.3 us. The Rossi-alpha plot and the fit show up on your screen. The die-away time found through data-fitting is 22.3 us.
3. Run `getSDT.m` 3. Run `getSDT.m`
The S,D,T count rates and their associated uncertainties are printed, as shown below: The S,D,T count rates and their associated uncertainties are printed, as shown below:
``` ```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment