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

    function [smod,psi]=daugment(smod,p,initpsi,DO_DESIGNTYPE)
%DAUGMENT D-optimally augment design
%
%  [D,PSI]=DAUGMENT(D,P,[INITPSI],[MODE]) augments the design D with P new
%  lines, in a d-optimal manner.  The optional INITPSI provides the
%  starting value of psi and saves it being calculated.  The optional MODE
%  can be set to 0 which will prevent design type information updates.  D
%  must pass a rankcheck else this routine will fail.

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


mdl=model(smod);

X=CalcJacob(mdl,factorsettings(smod));
[Q, R]=qr(X,0);
ri=R\eye(size(R));
ri= chol(ri*ri');

if nargin<4
    DO_DESIGNTYPE=1;
end
if nargin<3
    initpsi=(2*sum(log(abs(diag(R)))))./(size(X,2));
end

% need to do different things depending on whether we are
% replacing or not
if allowreps(smod)
    nc=ncand(smod);
else
    nc=ncandleft(smod);
end
k=NumTerms(mdl);
nf=nfactors(smod);

% test to see how big the matrices are going to be. If the final
% memory footprint is > ~20 Meg then we'll split into chunks (slower)

usewait=waitbars(smod);
if usewait
    % prevent any nested waitbars
    smod=waitbars(smod,0);
    % gui feedback
    h=xregGui.waitdlg('title','MBC Toolbox','message',['Augmenting design: 0/' sprintf('%d',p) ' points added.']);
    h.waitbar.max=p;
end

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

if nc*(nf+k)<2.5e6;
    % get the whole candidate list in one go, don't need to keep doing x2fx
    if allowreps(smod)
        lns=candidates(smod,'constrained');
    else
        lns=candidates(smod,'constrained','noreplacement');
    end
    biglns=CalcJacob(mdl,lns);

    for n=1:p
        d= ri*biglns';

        [mx,i]=max(sum((d).^2));

        % update stuff
        initpsi=initpsi+log(1+mx)./k;
        if p>1
            % update Ai for next iteration
            ri=cholupdate(ri,(ri'*d(:,i))/sqrt(1+mx),'-');
        end
        if ~allowreps(smod)
            % update design object
            % i is the index into lns, not into the full candidate set...
            % lns is the (constrained) candidate set, minus the design points
            % therefore need to do a 'noreplacement' indexcand inside design/augment
            smod=augment(smod,i,'indexnorep');
            % remove the added point from the list of options
            biglns(i,:)=[];
        else
            smod=augment(smod,i,'index');
        end
        if usewait
            h.waitbar.value=n;
            h.message= sprintf('Augmenting design: %d/%d points added.',n,p);
        end
    end
    psi=initpsi;

else
    % chunks of 0.5e6./(nf+k) :-(
    chunksz=fix(0.5e6./(nf+k));
    nchunks=fix(nc./chunksz);
    % calculate leftover (last) chunk size
    lastchunksz=nc-(nchunks.*chunksz);
    % if the last chunk is too small then it may cause problems, so add it to the next-to-last
    if lastchunksz<p
        nchunks=nchunks-1;
        lastchunksz=chunksz+lastchunksz;
    end
    % design indices used if ~reps
    des_ind=designindex(smod);

    reps=allowreps(smod);

    for n=1:p
        mx=0;
        i=[];
        maxln=[];
        % basic index vector
        ind=1:chunksz;

        for m=1:nchunks
            if reps
                [lns,inds]=indexcand(smod,((m-1).*chunksz)+ind);
            else
                [lns,inds]=indexcand(smod,setxor(((m-1).*chunksz)+ind,intersect(((m-1).*chunksz)+ind,des_ind)));
            end
            biglns=CalcJacob(mdl,lns);
            d= ri*biglns';
            [localmx,locali]=max(sum((d).^2));
            if localmx>mx
                mx=localmx;
                i=inds(locali);
                maxln=biglns(locali,:);
            end
            if usewait
                h.waitbar.value=(n-1)+m/(nchunks+1);
            end
        end
        % do last chunk
        if reps
            [lns,inds]=indexcand(smod,((nc-lastchunksz+1):nc));
        else
            [lns,inds]=indexcand(smod,setxor(((nc-lastchunksz+1):nc),intersect(((nc-lastchunksz+1):nc),des_ind)));
        end
        biglns=CalcJacob(mdl,lns);
        d= ri*biglns';
        [localmx,locali]=max(sum((d).^2));
        if localmx>mx
            mx=localmx;
            i=inds(locali);
            maxln=biglns(locali,:);
        end

        % update stuff
        initpsi=initpsi+log(1+mx)./k;

        if p>1
            % form correct d for updating
            d= ri*maxln';
            % update Ai for next iteration
            ri= cholupdate(ri,(ri'*d)/sqrt(1+mx),'-');
        end

        smod = augment(smod,i,'absindex');

        if ~reps
            % reduce nc and the last chunk size
            nc = nc-1;
            lastchunksz = lastchunksz-1;
        end
        if usewait
            h.waitbar.value=n;
            h.message= sprintf('Augmenting design: %d/%d points added.',n,p);
        end
    end
    psi = initpsi;
end

if usewait
    delete(h);
    % turn waitbars back on
    smod=waitbars(smod,1);
end
if DO_DESIGNTYPE
    % update design type
    smod=DesignType(smod,TP,INFO);         % reset object to initial setting
    smod=UpdateDesignType(smod,'d');       % update type setting
end