www.gusucode.com > mbcdesign 工具箱 matlab 源码程序 > mbcdesign/@des_linearmod/sumxtx.m
function [smod,s, nc]=sumxtx(smod,ReCalc) %SUMXTX Provide sum of x'x over entire candidate space % % [D,s]=SUMXTX(D) sums X'X over all the candidate points Note that if th % candidate set is set to continuous then this will use some form of % continuous grid. [D,s, nc]=SUMXTX(D) will also return the number of points % summed, N. % Copyright 2000-2015 The MathWorks, Inc. and Ford Global Technologies, Inc. if nargin==1 ReCalc=0; end % nc is not used any more nc=1; % search store for a valid copy if ~ReCalc && isfield(smod.store,'K') % K depends on candidate sets if (smod.store.K.candstate==candstate(smod) && smod.store.K.modelstate==modelstate(smod)) s=smod.store.K.data; return end end % if design object is empty, put a dummy point in to make this work if npoints(smod)==0 oldsmod=smod; fudged=1; smod=augment(smod,zeros(1,nfactors(smod)),'points'); else fudged=0; end % use only if no constraints if numConstraints(smod)==0 s= i_RectGrid(smod); else s= i_CandSet(smod); end if fudged % put back oldsmod for storing the result smod=oldsmod; end % store result smod.store.K.data=s; smod.store.K.designstate=designstate(smod); smod.store.K.candstate=candstate(smod); smod.store.K.modelstate=modelstate(smod); return function [s,nc]=i_CandSet(smod) % This sub-function uses the entire candidate set to calculate % s= X'*X. It can be slow nc=ncand(smod); if isinf(nc) % choose suitably large number nc=300000; end usewait=(waitbars(smod) && (nc*nfactors(smod))>50000); if usewait % prevent any nested waitbars smod=waitbars(smod,0); % gui feedback h=xregGui.waitdlg('title','MBC Toolbox','message','Calculating. Please wait...'); end mdl=model(smod); % number of points to do simultaneously % 10,000 points equates to about 1.5Mb peak memory usage np=10000; nd= 500; niter=floor(nc./np); nover=rem(nc,np); s=zeros(NumTerms(mdl)) ; if nc==0 % Can return now with s = 0 return end Lsel1= 1:np; nnd= floor(np/nd); nnv= nnd*nd+1; Lsel2= 1:nd; if usewait h.Waitbar.Max=niter*(nnd+1)+1; end for n=1:niter lns=indexcand(smod,(Lsel1+(n-1)*np)); biglns=CalcJacob(mdl,lns); % The following mex function does s=s+(biglns'*biglns); % without forming transpose and taking advantage of symmetry %s=mx_r1update(s,biglns(:,t)); % It seems to be faster to split up this call for i=1:nnd s=mx_r1update(s,biglns(Lsel2+(i-1)*nd,:)); if usewait h.waitbar.value= h.waitbar.value+1; drawnow('expose'); end end if nnv<np s=mx_r1update(s,biglns(nnv:end,:)); end if usewait h.waitbar.value= h.waitbar.value+1; drawnow('expose'); end end % leftover points if nover lns=indexcand(smod,((1:nover)+(niter).*np)); biglns=CalcJacob(mdl,lns); s=mx_r1update(s,biglns(:,:)); end if usewait h.waitbar.value= h.waitbar.value+1; drawnow('expose'); end s= s/nc; if usewait delete(h); % turn waitbars back on smod=waitbars(smod,1); end function [s,nc]= i_RectGrid(smod) % this sub function uses quadrature over [-1,1] grid mdl=model(smod); nc=1; ord= get(mdl,'order'); if any(ord>3) || size(mdl,1)*prod(ord+1) > 10e6 [s,nc]=i_CandSet(smod); else mdl= InitStore(mdl,factorsettings(smod)); [pev,s]= MeanPredVar(mdl); end