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

    function [psi,smod,X]=gcalc(smod,ReCalc,Method,NumPts)
% DES_LINEARMOD/GCALC  G-optimal value
%   PSI=GCALC(D) returns the g-optimality value for the
%   design object D.
%   See also: DCALC,VCAlC

%  Copyright 2000-2015 The MathWorks, Inc. and Ford Global Technologies, Inc.



% Created 8/11/99

[OK,smod]=rankcheck(smod);
if ~OK
    psi=[];
   return
end

if nargin < 2
   ReCalc= 0;
end

s= store(smod);

% search store for a valid copy
if ~ReCalc && isfield(s,'gpsi')
   % psi depends on design and candidate sets
   if (s.gpsi.designstate==designstate(smod)) && (s.gpsi.candstate==candstate(smod)) && ...
         (s.gpsi.modelstate==modelstate(smod))
      psi= s.gpsi.data;
      X=[];
      return
   end
end

if nargin<3
   Method='Sample';
end
switch lower(Method)
case 'sample'
	% sample from candidate set and optimise
   if nargin<4
      NumPts=[4 10];
   elseif length(NumPts)<2
      % number of groups
      NumPts= [NumPts(1) 10];
   end
   if nargout==3
      [psi,X]= i_samp(smod,NumPts(1),NumPts(2));
   else
      psi= i_samp(smod,NumPts(1),NumPts(2));
   end
case 'full'
	% ful candidate set
   [psi,X]= i_fullg(smod,1);
end   

% store result
s.gpsi.data=psi;
s.gpsi.designstate=designstate(smod);
s.gpsi.candstate=candstate(smod);
s.gpsi.modelstate=modelstate(smod);

smod= store(smod,s);

nm=inputname(1);
if nargout<2 && ~isempty(nm)
   assignin('caller',nm,smod);
end



function [psi,Xmax]= i_fullg(smod,NumPts)
% find G over full candidate set 

smod = InitStore(smod);

nc=ncand(smod);
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
np=10000;

niter=floor(nc./np);
nover=rem(nc,np);
if NumPts==1
   maxg=0;
   gvect=zeros(np,1);
else
   maxg=[];
   Xmax=[];
end
for n=1:niter
   lns=indexcand(smod,((1:np)+(n-1).*np));
	gvect= evalpev(smod,lns);
	
   if NumPts==1
      if max(gvect)>maxg
         [maxg,ind]=max(gvect);
         Xmax= lns(ind,:);
      end  
   else
      X= [lns;Xmax];
      P= [gvect;maxg];
      [dum,ind]= sort(P);
      maxg= P(ind(end:-1:end-NumPts+1));
      Xmax= X(ind(end:-1:end-NumPts+1),:);
   end
end
if nover
   % last load of points (<10000)
   lns=indexcand(smod,((1:nover)+(niter).*np));
	gvect= evalpev(smod,lns);
   if NumPts==1
      if max(gvect)>maxg
         [maxg,ind]=max(gvect);
         Xmax= lns(ind,:);
      end
   else
      X= [lns;Xmax];
      P= [gvect;maxg];
      [dum,ind]= sort(P);
      maxg= P(ind(end:-1:end-NumPts+1));
      Xmax= X(ind(end:-1:end-NumPts+1),:);
   end   
end
psi=maxg;

return




function [psi,G]= i_samp(smod,NumPts,nGroups)

m= model(smod);
smod = InitStore(smod);

nc= ncand(smod);
if nc> nGroups*1e4
   % make sure design points are used
   X=factorsettings(smod);
   P= evalpev(smod,X);
   [dum,ind]= sort(P);
   Y= P(ind(end:-1:end-NumPts+1));
   % try nGroups of 10000
   for i= 1:nGroups
      % take random sample from candidate space
      j=unidrnd(fix(nc/1e4));
      sind= j:fix(nc/1e4):nc;
      x=indexcand(smod,sind);
      X= [x; X];
      P= [evalpev(smod,x); Y];
      [dum,ind]= sort(P);
      Y= P(ind(end:-1:end-NumPts+1));
      X= X(ind(end:-1:end-NumPts+1),:);
   end   
elseif nc==0
    psi = 0;
    G = zeros(1, nfactors(smod)+1);
    return
else
   % do a full search if small candidate set
   [Y,X]=i_fullg(smod,NumPts);
end

% oprimise on best 'NumPts' values
fopts= optimset(optimset('fmincon'),...
    'largescale','off',...
    'Algorithm','active-set',...
    'MaxFunEvals',500*size(X,2),...
    'display','none');

% Constraints
Tgt = gettarget( m );
[conFcn, A, b] = getOptimConstraints( smod ); 

% Objective function
objFcn = @(x) gopt( x, smod );

for i=1:size(X,1)
	% perform nonlinear optimisations to search for 
   x0= X(i,:);
   [xopt, f, exitflag,output] = fmincon( objFcn, x0, A, b, [], [], Tgt(:,1), Tgt(:,2), conFcn, fopts );
   if exitflag>0
       % only update if optimization has been successful
       X(i,:)=xopt;
       Y(i)=evalpev(smod,xopt);
   end
end
psi= max(Y);
if nargout==2
	% determine bigist NumPts pevs
   [X,ind]=unique(X,'rows');
   Y= Y(ind);
   [Y,ind]= sort(Y);
   Y= Y(end:-1:1);
   X= X(ind(end:-1:1),:);
   j=[];
   for i=2:length(Y);
      if norm(X(i-1,:)-X(i,:))<1e-4
         j=[j i];
      end
   end
   if ~isempty(j)
      Y(j)=[];
      X(j,:)=[];
   end
   G= [Y(:) X];
end