www.gusucode.com > mbcexpr 工具箱 matlab 源码程序 > mbcexpr/+cgsimfill/@Feature/Smooth.m
function [Z,Cost,msg]= Smooth(F,ModelValues,fOutFcn) %SMOOTH feature fill with first order regulariser % % [Z,Cost]= Smooth(F,ModelValues,fOutFcn) % Copyright 2005-2015 The MathWorks, Inc. % use medium-scale algorithm with jacobian calculated with numjac [V0,JPattern,Bnds,A,C]= getOptimValues(F); n= size(JPattern,1); % make matrix of second derivatives As= cell(1,length(F.Tables)); for i=1:length(As) T= F.Tables(i); As{i}= del2(T.Object,T.Indices); end Aw= F.SmoothingFactor*blkdiag(As{:}); F.JacWorkFac= []; F.JacWorkG= []; msg = ''; if isempty(A) fopts= optimset(optimset('lsqnonlin'),... 'Algorithm','levenberg-marquardt',... 'OutputFcn',fOutFcn,... 'Display','none',... 'Jacobian','on'); if isempty(fOutFcn) fopts= optimset(fopts,'Display','iter'); end if any(isfinite(Bnds(:))) % have to use large scale algorithm for bounded problems fopts= optimset(fopts,'Algorithm','trust-region-reflective'); if length(ModelValues)<numel(V0) msg = 'Insufficient data points for bounded fill problems'; return end end % smoothing [Z,~,es]= lsqnonlin(@iCostSimFillSmooth,V0(:),Bnds(:,1),Bnds(:,2),fopts,... F,JPattern,ModelValues,Aw); else % Hard gradient constraints + smoothing fopts= optimset(optimset('fmincon'),... 'Algorithm','interior-point',... 'OutputFcn',fOutFcn,... 'Display','none',... 'DerivativeCheck','off',... 'GradConstr','on',... 'GradObj','on'); if isempty(fOutFcn) fopts= optimset(fopts,'Display','iter'); end qopt= optimset(optimset('quadprog'),... 'Algorithm','interior-point-convex',... 'display','none'); % find a feasible starting point close to initial value n= length(V0); dV =quadprog(eye(n),zeros(1,n),A,C-A*V0,[],[],... Bnds(:,1)-V0,Bnds(:,2)-V0,zeros(n,1),qopt); V0= V0+dV; fopts= optimset(fopts,'TolFun',F.Tolerance); Z= fmincon(@iCostSimFillConstraints,V0(:),A,C,[],[],Bnds(:,1),Bnds(:,2),'',fopts,... F,JPattern,ModelValues,Aw); es= iCostSimFillSmooth(Z,F,JPattern,ModelValues,Aw); end FTY= es(1:n); Cost= FTY'*FTY; Av= Aw*Z; solnmsg= sprintf('%8.5g %8.5g\n',[sqrt((FTY'*FTY)/n) sqrt((Av'*Av)/length(Z))]); if isempty(fOutFcn) fprintf('%s\n',solnmsg); else fOutFcn([],[],'message',solnmsg); end function [es,Js]= iCostSimFillSmooth(V,F,JPattern,ModelValues,Aw) FTY= CostSimFill(F,V,ModelValues); es= [FTY; Aw*V]; if nargout>1 % numjac uses scoped variables F.JacWorkFac and F.JacWorkG % as 'reference' V = V(:); Threshold = max(abs(V)*1e-8,1e-6); % set typical value to 1 where V is 0 otherwise numjac is poor TypicalY = V; TypicalY(V==0) = 1; [J,F.JacWorkFac,F.JacWorkG]= numjac(@CostSimFill,F,V(:),FTY(:),... {Threshold,TypicalY},F.JacWorkFac,0,JPattern,F.JacWorkG,... ModelValues); % medium scale algorithm needs a full matrix Js= full([J; Aw]); end function [c,g]= iCostSimFillConstraints(V,F,JPattern,ModelValues,Aw) if nargout==1 FTY = iCostSimFillSmooth(V,F,JPattern,ModelValues,Aw); else % numjac uses scoped variables JacWorkFac and JacG % as 'reference' [FTY,J]= iCostSimFillSmooth(V,F,JPattern,ModelValues,Aw); g= full(2*FTY'*J); end c= FTY'*FTY;