www.gusucode.com > mbcdesign 工具箱 matlab 源码程序 > mbcdesign/@cset_lhs/private/pr_selectlhs.m
function inds=pr_selectlhs(obj,N,nf,tpflg,meth,optimmeth,guiflag,nperms,strat,strat_lvls,doSymmetry) %PR_SELECTLHS Select a Latin Hypercube % % INDS=PR_SELECTLHS(N,nf,tpflg,costmeth,optimmeth,GUIFLAG,Nperms,stratify) % creates nf sets of permutations for the latin hypercube sampling % candidate set. tpflg is 0/1/2/3 corresponding to UINT8/16/32/double for % the indices. optimmeth parameter is currently ignored (only 'random' % allowed) % Copyright 2000-2011 The MathWorks, Inc. and Ford Global Technologies, Inc. % permutations are created here, then discrepancy checks are done % in mex function. The one with lowest variance is chosen. The random % seeds used for each permutation are saved and then used to recreate % the permutation indices. if nargin<8 strat=zeros(1,nf); end if nargin<7 nperms=500; end if nargin<6 guiflag=0; end % initialise inds cls = {'uint8', 'uint16', 'uint32', 'double'}; inds = zeros(N, nf, cls{tpflg+1}); minrandst=[]; % decide which subfunc to use for generating points if any(strat) genfunc=@i_geninds_strat; elseif doSymmetry genfunc=@i_geninds_sym; else genfunc=@i_geninds; end % Decide whether to show a waitbar or not. We only do this if the % computation is expected to take longer than about 0.5 seconds on a slow % machine if guiflag % Use same rules that the preview check uses. guiflag = ~checkPreview(obj); end h = []; Nperbar_counter = 0; Nperbar = 1; if guiflag % create a waitbar h=xregGui.waitdlg('message','Selecting Latin Hypercube. Please wait...'); h.waitbar.max = nperms; drawnow; % work out divisor to make sure we don't do too many waitbar calls Nperbar = floor(nperms/100); end switch lower(meth) case 'discrepancy' ndiscrep_boxes = 100+20*nf; discrep_frac=0.2; % width of box discrep_width = floor(discrep_frac*N); minvar=inf; for n=1:nperms % initialise random state [inds,randst]=feval(genfunc,inds,N,nf,strat); vars=qdiscrep(obj,inds,ndiscrep_boxes,discrep_width); if vars<minvar minvar=vars; minrandst=randst; end Nperbar_counter = i_incrementwait(h, guiflag, Nperbar, Nperbar_counter, n); end case 'minimax' minvar=inf; for n=1:nperms % initialise random state [inds,randst]=feval(genfunc,inds,N,nf,strat,strat_lvls); vars= qmaxdist(obj,inds); % mx_distance in max distance mode if vars<minvar minvar=vars; minrandst=randst; end Nperbar_counter = i_incrementwait(h, guiflag, Nperbar, Nperbar_counter, n); end case 'maximin' minvar=-1; for n=1:nperms % initialise random state [inds,randst]=feval(genfunc,inds,N,nf,strat,strat_lvls); vars= qmindist(obj,inds); % mx_distance in min distance mode if vars>minvar minvar=vars; minrandst=randst; end Nperbar_counter = i_incrementwait(h, guiflag, Nperbar, Nperbar_counter, n); end case 'cdfvariance' minvar=inf; for n=1:nperms % initialise random state [inds,randst]=feval(genfunc,inds,N,nf,strat,strat_lvls); vars= qcdfvar(obj,inds); if vars<minvar minvar=vars; minrandst=randst; end Nperbar_counter = i_incrementwait(h, guiflag, Nperbar, Nperbar_counter, n); end case 'cdfmaximum' minvar=inf; for n=1:nperms % initialise random state [inds,randst]=feval(genfunc,inds,N,nf,strat,strat_lvls); vars= qmaxcdfvar(obj,inds); if vars<minvar minvar=vars; minrandst=randst; end Nperbar_counter = i_incrementwait(h, guiflag, Nperbar, Nperbar_counter, n); end otherwise error(message('mbc:cset_lhs:InvalidState')); end if ~isempty(minrandst) % reconstruct best indices. although this resets the state of the default % stream, we have not reset it to a fixed state, only to a state visited % by this function in one of the above loops. dflt = RandStream.getGlobalStream; dflt.State = minrandst; inds=feval(genfunc,inds,N,nf,strat,strat_lvls); else warning(message('mbc:cset_lhs:InvalidState1')); end if guiflag delete(h); end function Counter = i_incrementwait(hWait, guiflag, Nperbar, Counter, n) if guiflag if Counter>=Nperbar hWait.waitbar.value = n; drawnow; Counter=0; else Counter = Counter+1; end end function [inds,rnd]=i_geninds(inds,N,nf,strat,strat_lvls) dflt=RandStream.getGlobalStream; rnd=dflt.State; % initialise random state for m=1:nf inds(:,m)=randperm(N)'; end function [inds,rnd]=i_geninds_strat(inds,N,nf,strat,strat_lvls) dflt=RandStream.getGlobalStream; rnd=dflt.State; % initialise random state for m=1:nf tmp=randperm(N)'; if strat(m) if strat(m)==-1 % use the specified levels Nstrat = length(strat_lvls{m}); tmp = (ceil(((Nstrat).*tmp./N))); % Nstrat stratified levels in 1--> Nstrat region tmp = strat_lvls{m}(tmp); % index into the 1--->N domain tmp = tmp(:); else tmp=((N-1).*(ceil(((strat(m)).*tmp./N))-1)./(strat(m)-1))+1; end end inds(:,m)=tmp; end function [inds,rnd]=i_geninds_sym(inds,N,nf,strat,strat_lvls) dflt=RandStream.getGlobalStream; rnd=dflt.State; halfN=floor(N/2); if (halfN)==(N/2) for m=1:nf inds(:,m)=i_symrandperm(N)'; end else for m=1:nf tmp=i_symrandperm(N)'; inds(1:halfN,m)=tmp(1:halfN); inds((halfN+2):end,m)=tmp(halfN+1:end); inds(halfN+1,m)=halfN+1; end end function vect=i_symrandperm(N) % generate a random permutation with forced symmetry, for i_geninds_sym % ONLY WORKS FOR N=EVEN halfN=floor(N/2); vect=rand(1,halfN); [nul,vect]=sort([vect 1-vect(end:-1:1)]);