www.gusucode.com > mbcdesign 工具箱 matlab 源码程序 > mbcdesign/@des_multimod/vaugment.m
function [des,psi] = vaugment(des,p,initpsi,DO_DESIGNTYPE) %VAUGMENT V-optimally augment design % % [D,PSI]=VAUGMENT(D,P,[INITPSI],[MODE]) augments the design D with P new % lines, in a v-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-2015 The MathWorks, Inc. and Ford Global Technologies, Inc. if nargin<3 initpsi=vcalc(des); end if nargin<4 DO_DESIGNTYPE=1; end [des,K]=sumxtx(des); mmdl=model(des); fs=factorsettings(des); mdls=get(mmdl,'models'); wts=get(mmdl,'weights'); nmdls=length(wts); % Initial Ai must be for all points. It is then iteratively % updated and hence is always correctly derived from all points. Ai = cell(1, nmdls); k = zeros(1, nmdls); for n=1:nmdls X=CalcJacob(mdls{n},fs); [~,R]= qr(X,0); ri= R\eye(size(R)); Ai{n}=ri*ri'; k(n)=NumTerms(mdls{n}); 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 for n=1:p % quick way psiupd = zeros(nc,1); for m=1:nmdls % build up div as weighted sum % quick way B= biglns{m}*Ai{m}; div=(1+sum(B.*biglns{m},2)); Bi=wts(m).*sum( B.*(B*K{m}),2 ); % check div isn't 0!! div=max(div,Bi*eps*2); % form output vector psiupd = psiupd + Bi./div; end [psiupd,i]=max(psiupd); initpsi=initpsi-psiupd; if p>1 for m=1:nmdls % Update Ai Bi=biglns{m}(i,:)*Ai{m}; divi=(1+sum(Bi.*biglns{m}(i,:),2)); Ai{m}= mx_r1update(Ai{m},Bi/sqrt(divi),1); 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'); % remove the added point from the list of options for m=1:nmdls biglns{m}(i,:)=[]; end nc=nc-1; 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+max(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 psiupd=zeros(size(lns,1),1); for l=1:nmdls biglns=CalcJacob(mdls{l},lns); B= biglns*Ai{l}; div=(1+sum(B.*biglns,2)); Bi=wts(l).*sum( B.*(B*K{l}),2 ); div=max(div,Bi*eps*2); psiupd=psiupd+Bi./div; end [localmx,locali]=max(psiupd); 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 psiupd=zeros(size(lns,1),1); for l=1:nmdls biglns=CalcJacob(mdls{l},lns); B= biglns*Ai{l}; div=(1+sum(B.*biglns,2)); Bi=wts(l).*sum( B.*(B*K{l}),2 ); div=max(div,Bi*eps*2); psiupd=psiupd+Bi./div; end [localmx,locali]=max(psiupd); if localmx>mx mx=localmx; i=inds(locali); maxln=lns(locali,:); end % update stuff initpsi=initpsi-mx; if p>1 % form appropriate B, etc for l=1:nmdls bigln=CalcJacob(mdls{l},maxln); Bi=bigln*Ai{l}; divi=wts(l).*(1+sum(Bi.*bigln,2)); % Update Ai Ai{l}= mx_r1update(Ai{l},Bi/sqrt(divi),1); 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,'v'); % update type setting end