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)]);