www.gusucode.com > mbcexpr 工具箱 matlab 源码程序 > mbcexpr/@cgexpr/evaluategrid.m

    function Y = evaluategrid(obj, pGrid, func, logicalflag)
%EVALUATEGRID Evaluate expression over a grid
%
%  Y = EVALUATEGRID(OBJ, GRIDPTRS, QUANTITY) evaluates the expression OBJ
%  over a grid.  GRIDPTRS is a list of pointers to inputs to OBJ that you
%  want to create a grid over.  The output Y will be returned as an N-D
%  array with the dimensions in the same order as GRIDPTRS.   QUANTITY is
%  an optional parameter to specify which aspect of the expression to
%  evaluate and may be one of 'value', 'pev' or 'constraint' or a function
%  handle to a function with the syntax Y = FUNC(OBJ).
%
%  If QUANTITY is not specified then 'value' is assumed.
%
%  Y = EVALUATEGRID(OBJ, GRIDPTRS, QUANTITY, LOGICAL) indicates whether the
%  function evaluation will return logical values (LOGICAL=true).
%
%  All other pointers in the inports list must not contain vectors unless
%  they are dependent on one of the pointers in GRIDPTRS.  All of the
%  GRIDPTRS must be independent of each other.
%
%  See also CGEXPR/EVALUATE, CGEXPR/EVALUATESEQUENTIAL,
%           CGEXPR/EVALUATEAVERAGE.

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


if nargin<2
    error(message('mbc:cgexpr:InvalidArgument'));
end

if nargin<4
    logicalflag = false;
end
if nargin<3
    func = @i_eval;
end

if length(pGrid)<2
    % transfer to standard evaluate mechanism as that is faster
    Y = evaluate(obj, func, logicalflag);
    
    % Check that output is the correct length and scalar expand if required
    if length(pGrid)>0
        outLen = length(pGrid.getvalue);
    else
        outLen = 1;
    end
    if length(Y)~=outLen
        Y(1:outLen,1) = Y;
    end
    return
end

% Check evaluation function
if ischar(func)
    nQuant = strcmpi(func, {'value', 'pev', 'constraint'});
    if ~any(nQuant)
        % Convert the string to a function handle
        func = str2func(func);
    elseif nQuant(1)
        func = @i_eval;
    elseif nQuant(2)
        func = @peveval;
    else
        func = @ceval;
    end
end

% Check that the inputs are all independent of each other
if ~cgisindependentvars(pGrid)
    error(message('mbc:cgexpr:InvalidArgument2'));
end

% Find the pointers that are not gridding variables or dependent on the
% gridding variables (this includes the grid pointers themselves).  All of
% these remaining pointers must be scalars.
pInputs = getinports(obj);
pInputs = pInputs(cgisindependentvars(pInputs, pGrid));bScalar = pveceval(pInputs, @isscalar);
if ~all([bScalar{:}])
    error(message('mbc:cgexpr:InvalidArgument3'));   
end


dValues = pveceval(pGrid, @getvalue);
dims = cellfun('length', dValues);
nElements  = prod(dims);

% Check that the ndgrid won't consume too much memory.  If it looks big
% then the calculation is split into chunks which takes longer but will
% avoid failing or causing disk thrashing.  Note that the actual memory
% consumed is in part determined by the complexity of the expression.

% LIMIT is the limit on the number of input elements
LIMIT = 2^16;
BLOCKSIZE = 2^15;
if (nElements*length(pGrid)) < LIMIT
    dGridValues = cell(size(dValues));
    if length(dGridValues)>1
        [dGridValues{1:length(dGridValues)}] = ndgrid(dValues{:});
        for n = 1:length(dGridValues)
            dGridValues{n} = dGridValues{n}(:);
        end
    else
        dGridValues{1} = dValues{1};
    end
    passign(pGrid, parrayeval(pGrid, @setvalue, {dGridValues}));
    % Pre-allocating the output automatically corrects the rare cases where
    % the output of i_eval is scalar even when the grid inputs are vectors.
    if logicalflag
        Y = false(nElements, 1);
    else
        Y = zeros(nElements, 1);
    end
    Y(:) = func(obj);
else
    dGridValues = cell(size(dValues));
    cs = cset_grid(dValues);
    nBlocks = floor(nElements./BLOCKSIZE);
    nLeftover = mod(nElements, BLOCKSIZE);
    if logicalflag
        Y = false(nElements, 1);
    else
        Y = zeros(nElements, 1);
    end
    nIndex = uint32(1):uint32(BLOCKSIZE);
    for n = 1:nBlocks
        dPoints = partialset(cs, nIndex);
        for m = 1:length(dGridValues)
            dGridValues{m} = dPoints(:,m);            
        end
        passign(pGrid, parrayeval(pGrid, @setvalue, {dGridValues}));
        Y(nIndex) = func(obj);
        nIndex = nIndex+BLOCKSIZE;
    end
    if nLeftover
        % Final block is smaller
        dPoints = partialset(cs, nIndex(1:nLeftover));
        for m = 1:length(dGridValues)
            dGridValues{m} = dPoints(:,m);            
        end
        passign(pGrid, parrayeval(pGrid, @setvalue, {dGridValues}));
        Y(nIndex(1:nLeftover)) = func(obj);
    end
end

passign(pGrid, parrayeval(pGrid, @setvalue, {dValues}));
if length(dims)>1
    Y = reshape(Y, dims);
end