www.gusucode.com > mbcdesign 工具箱 matlab 源码程序 > mbcdesign/@des_multimod/sumxtx.m
function [des,s, nc]=sumxtx(des,ReCalc) %SUMXTX Provide sum of x'x over entire candidate space % % S=SUMXTX(D) sums X'X over all the candidate points. Note that if the % candidate set is set to continuous then this will use some form of % continuous grid. % [S, N]=SUMXTX(D) will also return the number of points summed, N. % S will be a cell array containing the sum of X'X for each model % contained within the multimodel in D. % 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(des.store,'K') % K depends on candidate sets if (des.store.K.candstate==candstate(des) && des.store.K.modelstate==modelstate(des)) s=des.store.K.data; return end end % if design object is empty, put a dummy point in to make this work if npoints(des)==0 olddes=des; fudged=1; des=augment(des,zeros(1,nfactors(des)),'points'); else fudged=0; end % use only if no constraints if numConstraints(des)==0 s= i_RectGrid(des); else s= i_CandSet(des); end if fudged % put back olddes for storing the result des=olddes; end % store result des.store.K.data=s; des.store.K.designstate=designstate(des); des.store.K.candstate=candstate(des); des.store.K.modelstate=modelstate(des); return function [s,nc]=i_CandSet(des) % This sub-function uses the entire candidate set to calculate % s= X'*X. It can be slow nc=ncand(des); if isinf(nc) % choose suitably large number nc=300000; end % number of points to do simultaneously % 10,000 points equates to about 1.5Mb peak memory usage - not too large at all np=10000; nd= 500; niter=floor(nc./np); nover=rem(nc,np); Lsel1= 1:np; nnd= floor(np/nd); nnv= nnd*nd+1; Lsel2= 1:nd; mmdl=model(des); mdls = get(mmdl,'models'); nmdls = length(mmdl); usewait=waitbars(des) & (nmdls+niter)>1; if usewait % prevent any nested waitbars des=waitbars(des,0); % gui feedback h=xregGui.waitdlg('message','Calculating. Please Wait...'); end % calc s for each model s=cell(1,nmdls); for m=1:nmdls stmp=zeros(NumTerms(mdls{m})) ; for n=1:niter lns=indexcand(des,(Lsel1+(n-1)*np)); biglns=CalcJacob(mdls{m},lns); % The following mex function does s=s+(biglns'*biglns); % without forming transpose and taking advantage of symmetry % It seems to be faster to split up this call for i=1:nnd stmp=mx_r1update(stmp,biglns(Lsel2+(i-1)*nd,:)); end if nnv<np stmp=mx_r1update(stmp,biglns(nnv:end,:)); end if usewait h.waitbar.value= (m-1)/nmdls + n/(nmdls*(niter+1)); drawnow('expose'); end end if nover % leftover points lns=indexcand(des,((1:nover)+(niter).*np)); biglns=CalcJacob(mdls{m},lns); stmp=mx_r1update(stmp,biglns(:,:)); end if nc>0 s{m} = stmp/nc; else s{m} = stmp; end if usewait h.waitbar.value= (m/nmdls); drawnow('expose'); end end if usewait des=waitbars(des,1); delete(h); end function [s,nc]= i_RectGrid(des) % this sub function uses quadrature over [-1,1] grid nc=1; mmdl=model(des); mdls = get(mmdl,'models'); nmdls = length(mmdl); % calc s for each model s=cell(1,nmdls); fs=factorsettings(des); usewait=waitbars(des) & nmdls>1; if usewait % prevent any nested waitbars des=waitbars(des,0); % gui feedback h = xregGui.waitdlg('message', 'Calculating. Please Wait...'); end for m=1:nmdls [pev,s{m}]= MeanPredVar(InitStore(mdls{m},fs)); if usewait h.waitbar.value = m/nmdls; drawnow('expose'); end end if usewait des=waitbars(des,1); delete(h); end