?? lst_sq.m
字號:
%function [x,OPTIONS] = lst_sq(FUN,x,OPTIONS,GRADFUN,time_out,P1,P2,P3,P4,P5,P6,P7,P8,P9,
P10,P11,P12,P13,P14,P15,P16,P17,P18,P19,P20,P21)
% See leastsq.m for usage. Modified to time-out after "time_out" seconds.
%
% X=LST_SQ('FUN',X,OPTIONS,GRADFUN,time_out,P1,P2,..) allows
% coefficients, P1, P2, ... to be passed directly to FUN:
% F=FUN(X,P1,P2,...). Empty arguments ([]) are ignored.
% Author Date Predecessor Modification
% ====== ==== =========== ============
% B.McKay M.Willis 25/5/95 leastsq.m use etime instead of tic/toc
% J. Elsey 31/7/95 lst_sq.m increase the number of allowable inputs
% B.McKay 6/10/95 - returns latest x in case of time out (rather than initial)
%
% Function Calls:
% ===============
%
% Last Modification: 6/10/95
% =================
%
function [x,OPTIONS] = lst_sq(FUN,x,OPTIONS,GRADFUN,time_out,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,
P12,P13,P14,P15,P16,P17,P18,P19,P20,P21)
global optimea optimeb optimec;
% ------------Initialization----------------
t1=clock;
x_initial = x;
XOUT=x(:);
[nvars]=length(XOUT);
evalstr = [FUN];
if ~any(FUN<48)
evalstr=[evalstr, '(x'];
optimes=clock;
for i=1:nargin - 4
evalstr = [evalstr,',P',num2str(i)];
end
optimef=clock;
optimea=optimea+etime(optimef,optimes);
evalstr = [evalstr, ')'];
end
if nargin < 3, OPTIONS=[]; end
if nargin < 4, GRADFUN=[]; end
if length(GRADFUN)
evalstr2 = [GRADFUN];
if ~any(GRADFUN<48)
evalstr2 = [evalstr2, '(x'];
optimes=clock;
for i=1:nargin - 4
evalstr2 = [evalstr2,',P',num2str(i)];
end
optimef=clock;
optimeb=optimeb+etime(optimef,optimes);
evalstr2 = [evalstr2, ')'];
end
end
f = eval(evalstr);
f=f(:);
nfun=length(f);
GRAD=zeros(length(XOUT),nfun);
OLDX=XOUT;
MATX=zeros(3,1);
MATL=[f'*f;0;0];
OLDF=f'*f;
FIRSTF=f'*f;
[OLDX,OLDF,OPTIONS]=lsint(XOUT,f,OPTIONS);
EstSum=0.5;
GradFactor=0;
CHG = 1e-7*abs(XOUT)+1e-7*ones(nvars,1);
OPTIONS(10)=1;
status=-1;
while status~=1
if etime(clock,t1) > time_out; x(:)=OLDX; return; end
% Work Out Gradients
if ~length(GRADFUN) | OPTIONS(9)
OLDF=f;
CHG = sign(CHG+eps).*min(max(abs(CHG),OPTIONS(16)),OPTIONS(17));
optimes=clock;
for gcnt=1:nvars
if XOUT==[] | any(isnan(XOUT(:)))==1 | any(XOUT(:)==inf)==1 | any(...
XOUT(:)==-inf)==1,
x=x_initial;return, end;
temp = XOUT(gcnt);
XOUT(gcnt) = temp +CHG(gcnt);
x(:) = XOUT;
f(:) = eval(evalstr);
GRAD(gcnt,:)=(f-OLDF)'/(CHG(gcnt));
XOUT(gcnt) = temp;
end
optimef=clock;
optimec=optimec+etime(optimef,optimes);
f = OLDF;
OPTIONS(10)=OPTIONS(10)+nvars;
% Gradient check
if OPTIONS(9) == 1
GRADFD = GRAD;
x(:)=XOUT; GRAD = eval(evalstr2);
graderr(GRADFD, GRAD, evalstr2);
OPTIONS(9) = 0;
end
else
x(:) = XOUT;
OPTIONS(11)=OPTIONS(11)+1;
GRAD = eval(evalstr2);
end
% Try to set difference to 1e-8 for next iteration
if nfun==1
CHG = nfun*1e-8./GRAD;
else
CHG = nfun*1e-8./sum(abs(GRAD)')';
end
gradf = 2*GRAD*f;
fnew = f'*f;
%---------------Initialization of Search Direction------------------
if status==-1
if GRAD==[] | any(isnan(GRAD(:)))==1 | any(GRAD(:)==inf)==1 | any(...
GRAD(:)==-inf)==1,
x=x_initial;return, end;
if cond(GRAD)>1e8
SD=-(GRAD*GRAD'+(norm(GRAD)+1)*(eye(nvars)))\(GRAD*f);
if OPTIONS(5)==0, GradFactor=norm(GRAD)+1; end
how='COND';
else
% SD=GRAD'\(GRAD'*X-F)-X;
SD=-(GRAD*GRAD'+GradFactor*(eye(nvars)))\(GRAD*f);
end
FIRSTF=fnew;
OLDG=GRAD;
GDOLD=gradf'*SD;
OPTIONS(18)=1;
if OPTIONS(1)>0
disp([sprintf('%5.0f %12.3g %12.3g ',OPTIONS(10),fnew,OPTIONS(18)),
sprintf('%12.3g ',GDOLD)]);
end
XOUT=XOUT+OPTIONS(18)*SD;
if OPTIONS(5)==0
newf=GRAD'*SD+f;
GradFactor=newf'*newf;
SD=-(GRAD*GRAD'+GradFactor*(eye(nvars)))\(GRAD*f);
end
newf=GRAD'*SD+f;
XOUT=XOUT+OPTIONS(18)*SD;
EstSum=newf'*newf;
status=0;
if OPTIONS(7)==0; PCNT=1; end
else
%-------------Direction Update------------------
gdnew=gradf'*SD;
if OPTIONS(1)>0,
num=[sprintf('%5.0f %12.3g %12.3g ',OPTIONS(10),fnew,OPTIONS(18))
,sprintf('%12.3g ',gdnew)];
end
if gdnew>0 & fnew>FIRSTF
% Case 1: New function is bigger than last and gradient w.r.t. SD -ve
% ... interpolate.
how='inter';
[stepsize]=cubici1(fnew,FIRSTF,gdnew,GDOLD,OPTIONS(18));
OPTIONS(18)=0.9*stepsize;
elseif fnew<FIRSTF
% New function less than old fun. and OK for updating
% .... update and calculate new direction.
[newstep,fbest] =cubici3(fnew,FIRSTF,gdnew,GDOLD,OPTIONS(18));
if fbest>fnew,fbest=0.9*fnew; end
if gdnew<0
how='incstep';
if newstep<OPTIONS(18), newstep=(2*OPTIONS(18)+1e-4);
how=[how,'IF']; end
OPTIONS(18)=abs(newstep);
else
if OPTIONS(18)>0.9
how='int_step';
OPTIONS(18)=min([1,abs(newstep)]);
end
end
% SET DIRECTION.
% Gauss-Newton Method
temp=1;
if OPTIONS(5)==1
if OPTIONS(18)>1e-8 & cond(GRAD)<1e8
SD=GRAD'\(GRAD'*XOUT-f)-XOUT;
if SD'*gradf>eps,how='ERROR- GN not descent direction', end
temp=0;
else
if OPTIONS(1) > 0
disp('Conditioning of Gradient Poor - Switching To LM method')
end
how='CHG2LM';
OPTIONS(5)=0;
OPTIONS(18)=abs(OPTIONS(18));
end
end
if (temp)
% Levenberg_marquardt Method N.B. EstSum is the estimated sum of squares.
% GradFactor is the value of lambda.
% Estimated Residual:
if EstSum>fbest
GradFactor=GradFactor/(1+OPTIONS(18));
else
GradFactor=GradFactor+(fbest-EstSum)/(OPTIONS(18)+eps);
end
SD=-(GRAD*GRAD'+GradFactor*(eye(nvars)))\(GRAD*f);
OPTIONS(18)=1;
estf=GRAD'*SD+f;
EstSum=estf'*estf;
if OPTIONS(1)>0, num=[num,sprintf('%12.3g ',GradFactor)]; end
end
gdnew=gradf'*SD;
OLDX=XOUT;
% Save Variables
FIRSTF=fnew;
OLDG=gradf;
GDOLD=gdnew;
% If quadratic interpolation set PCNT
if OPTIONS(7)==0, PCNT=1; MATX=zeros(3,1); MATL(1)=fnew; end
else
% Halve Step-length
how='Red_Step';
if fnew==FIRSTF,
if OPTIONS(1)>0,
disp('No improvement in search direction: Terminating')
end
status=1;
else
OPTIONS(18)=OPTIONS(18)/8;
if OPTIONS(18)<1e-8
OPTIONS(18)=-OPTIONS(18);
end
end
end
XOUT=OLDX+OPTIONS(18)*SD;
if OPTIONS(1)>0,disp([num,how]),end
end %----------End of Direction Update-------------------
if OPTIONS(7)==0, PCNT=1; MATX=zeros(3,1); MATL(1)=fnew; end
% Check Termination
if max(abs(SD))< OPTIONS(2) & (gradf'*SD) < OPTIONS(3) & max(abs(gradf)) <
10*(OPTIONS(3)+OPTIONS(2))
if OPTIONS(1) > 0
disp('Optimization Terminated Successfully')
end
status=1;
elseif OPTIONS(10)>OPTIONS(14)
disp('maximum number of iterations has been exceeded');
if OPTIONS(1)>0
disp('Increase OPTIONS(14)')
end
status=1;
else
% Line search using mixed polynomial interpolation and extrapolation.
if PCNT~=0
while PCNT > 0
if etime(clock,t1) > time_out; x(:)=OLDX;; return; end
x(:) = XOUT; f(:) = eval(evalstr); OPTIONS(10)=OPTIONS(10)+1;
fnew = f'*f;
if fnew<OLDF'*OLDF, OX = XOUT; OLDF=f; end
[PCNT,MATL,MATX,steplen,fnew,how]=searchq(PCNT,fnew,OLDX,MATL,MATX,SD,GDOLD,OPTIONS(18),how);
OPTIONS(18)=steplen;
XOUT=OLDX+steplen*SD;
if fnew==FIRSTF, PCNT=0; end
end
XOUT = OX;
f=OLDF;
else
x(:)=XOUT; f(:) = eval(evalstr); OPTIONS(10)=OPTIONS(10)+1;
end
end
end
OPTIONS(8) = fnew;
XOUT=OLDX;
x(:)=XOUT;
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -