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;