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

    function [des,psi] = daugment(des,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.


mmdl=model(des);

fs=factorsettings(des);
mdls=get(mmdl,'models');
wts=get(mmdl,'weights');
nmdls=length(wts);

ri = cell(1,nmdls);
k = zeros(1, nmdls);
for n=nmdls:-1:1
    X=CalcJacob(mdls{n},fs);
    [Q, R]=qr(X,0);
    ri{n}=R\eye(size(R));
    ri{n}= chol(ri{n}*ri{n}');
    k(n)=NumTerms(mdls{n});
end

if nargin<3
    initpsi=dcalc(des);
end
if nargin<4
    DO_DESIGNTYPE=1;
end

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

usewait=waitbars(des);
if usewait
    % prevent any nested waitbars
    des=waitbars(des,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(des);
end

% 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)

if nc*(nf+sum(k))<2.5e6;
    % get the whole candidate list in one go, don't need to keep doing x2fx
    if allowreps(des)
        lns=candidates(des,'constrained');
    else
        lns=candidates(des,'constrained','noreplacement');
    end
    % multiple biglns copies, unfortunately
    biglns=cell(1,nmdls);
    for n=1:nmdls
        biglns{n}=CalcJacob(mdls{n},lns)';
    end
    % get rid of lns - reclaim valuable memory
    clear lns

    for n=1:p
        if allowreps(des)
            coef = ones(nc,1);
        else
            coef = ones(nc+1-n,1);
        end
        for m=1:nmdls
            d= ri{m}*biglns{m};
            coef = coef.*(1+sum((d).^2)).^(wts(m)./k(m))';
        end

        [delta,i]=max(coef);

        % update stuff
        initpsi=initpsi+log(delta);

        if p>1
            for m=1:nmdls
                % recalc the necessary bits
                di=ri{m}*biglns{m}(:,i);
                ri{m}=cholupdate(ri{m},(ri{m}'*di)/sqrt(1+sum(di.^2)),'-');
                if ~allowreps(des)
                    % remove point from candidate list
                    biglns{m}(:,i)=[];
                end
            end
        end

        if ~allowreps(des)
            % 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
            des=augment(des,i,'indexnorep');
        else
            des=augment(des,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+sum(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(des);

    reps=allowreps(des);

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

        for m=1:nchunks
            if reps
                [lns,inds]=indexcand(des,((m-1).*chunksz)+ind);
            else
                [lns,inds]=indexcand(des,setxor(((m-1).*chunksz)+ind,intersect(((m-1).*chunksz)+ind,des_ind)));
            end

            coef = ones(size(lns,1),1);
            for l=1:nmdls
                biglns=CalcJacob(mdls{l},lns);
                d= ri{l}*biglns';
                coef = coef.*((1+sum((d).^2)).^(wts(l)./k(l)))';
            end
            [localmx,locali]=max(coef);
            if localmx>mx
                mx=localmx;
                i=inds(locali);
                maxln=lns(locali,:);
            end
            if usewait
                h.waitbar.value= n - 1 + m/(nchunks+1);
            end
        end
        % do last chunk
        if reps
            [lns,inds]=indexcand(des,((nc-lastchunksz+1):nc));
        else
            [lns,inds]=indexcand(des,setxor(((nc-lastchunksz+1):nc),intersect(((nc-lastchunksz+1):nc),des_ind)));
        end

        coef=ones(size(lns,1),1);
        for l=1:nmdls
            biglns=CalcJacob(mdls{l},lns);
            d= ri{l}*biglns';
            coef = coef.*(1+sum((d).^2)).^(wts(l)./k(l))';
        end
        [localmx,locali]=max(coef);
        if localmx>mx
            mx=localmx;
            i=inds(locali);
            maxln=lns(locali,:);
        end

        % update stuff
        initpsi=initpsi+log(mx);
        if p>1
            for l=1:nmdls
                bigmaxln=CalcJacob(mdls{l},maxln);
                % form correct d for updating
                di= ri{l}*bigmaxln';
                % update Ai for next iteration
                ri{l}= cholupdate(ri{l},(ri{l}'*di)/sqrt(1+sum(di.^2)),'-');
            end
        end
        des = augment(des,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
    des=waitbars(des,1);
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