www.gusucode.com > mbcdesign 工具箱 matlab 源码程序 > mbcdesign/@des_linearmod/aaugment.m
function [smod,psi]=aaugment(smod, p, initpsi,DO_DESIGNTYPE) %AAUGMENT A-optimally augment design % % [D,PSI]=AAUGMENT(D,P,[INITPSI],[MODE]) augments the design D with P new % lines, in an a-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. if nargin<4 DO_DESIGNTYPE=1; end if nargin<3 initpsi=acalc(smod); end mdl=model(smod); % Initial Ai must be for all points. It is then iteratively % updated and hence is always correctly derived from all points. X=CalcJacob(mdl,factorsettings(smod)); [Q,R]= qr(X,0); ri= R\eye(size(R)); Ai=ri*ri'; % X used for augmentation search is taken from the candidate set % 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 % quick way B= biglns*Ai; div=(1+sum(B.*biglns,2)); Bi=sum( B.*B,2 ); div=max(div,Bi*eps*2); [delpsi,i]=max(Bi./div); initpsi=initpsi-delpsi; if p>1 % Update Ai Ai= mx_r1update(Ai,B(i,:)/sqrt(div(i)),1); 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=[]; maxdiv=[]; % 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); B= biglns*Ai; div=(1+sum(B.*biglns,2)); Bi=sum( B.*B,2 ); div=max(div,Bi*eps*2); [localmx,locali]=max(Bi./div); if localmx>mx mx=localmx; i = inds(locali); maxln=biglns(locali,:); maxdiv=div(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); B= biglns*Ai; div=(1+sum(B.*biglns,2)); Bi=sum( B.*B,2 ); div=max(div,Bi*eps*2); [localmx,locali]=max(Bi./div); if localmx>mx mx=localmx; i=inds(locali); maxln=biglns(locali,:); maxdiv=div(locali); end % update stuff initpsi=initpsi-mx; if p>1 % form appropriate B, etc B=maxln*Ai; % Update Ai Ai= mx_r1update(Ai,B/sqrt(maxdiv),1); 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,'a'); % update type setting end