www.gusucode.com > mbcdesign 工具箱 matlab 源码程序 > mbcdesign/@des_multimod/ddelete.m

    function [des, psi]=ddelete(des,initpsi,p,DO_DESIGNTYPE)
%DDELETE D-optimal deletion
%
%  [D,PSI]=DDELETE(D,INITPSI,P) deletes P lines from the design D using
%  d-optimality.  A new design object and the new d-opimality criteria,
%  PSI, are returned.

%  Copyright 2000-2004 The MathWorks, Inc. and Ford Global Technologies, Inc.


if nargin<4
    DO_DESIGNTYPE=1;
end

mmdl=model(des);
fs= factorsettings(des);
mdls= get(mmdl,'models');
wts= get(mmdl,'weights');
fp=freepoints(des);
nmdls=length(wts);
ri = cell(1, nmdls);
Xsmall = cell(1, nmdls);
Nt = zeros(1, nmdls);
for n=1:nmdls
    X=CalcJacob(mdls{n},fs);
    [Q,R]=qr(X,0);
    ri{n}=R\eye(size(R));
    ri{n}= chol(ri{n}*ri{n}');

    Xsmall{n}= X(fp,:)';
    Nt(n)= size(ri{n},1);
end

if DO_DESIGNTYPE
    % need to remember start point
    [TP,INFO]= DesignType(des);
end

Np=size(fs,1);
psinew=initpsi;
for n=1:p
    % get new optimality criteria for if each line were removed
    % and test for maximum

    % alternative implementation - don't save whole d matrices for later on
    coef = ones(Np+1-n,1);
    for m=1:nmdls
        d= ri{m}*Xsmall{m};
        coef = coef.*((1-sum(d.^2, 1)).^(wts(m)./Nt(m)))';
    end

    [delta,i]=max(coef);

    psinew=psinew+log(delta);

    if p>1
        for m=1:nmdls
            % apply changes to each stored matrix
            % recalc the necessary bits
            di=ri{m}*Xsmall{m}(:,i);
            ri{m}=cholupdate(ri{m},(ri{m}'*di)/sqrt(1-sum(di.^2)));

            Xsmall{m}(:,i)=[];
        end
    end
    des=delete(des,'indexed',i,'changeable');
end

if DO_DESIGNTYPE
    % update design type
    des=DesignType(des,TP,INFO);         % reset object to initial setting
    des=UpdateDesignType(des,'d');       % update type setting
end

psi=psinew;