fgridmin.m

function [i,i_] = fgridmin(y)
% Search an N-dimensional array y for points that may be proximate to a global minimum
% of the interpolated array; return the proximate points' N-D subscripts (i) in y. Also
% optionally return the subscripts (i_) for discrete local minima in y.
%
% syntax:
% i = fgridmin(y);
% [i,i_] = fgridmin(y);
%
% input arg:
%
% y: N-dimensional array, real double (See validation conditions below.)
%
% output args:
%
% i: size-[M, N] matrix of integer
%
% i_: size-[M_,N] matrix of integer (optional output)
%
% N is defined above such that size(y,N) > 1 and size(y,N+m) == 1 for all m > 0. If y is
% a scalar (case N = 0), or if size(y,m) < 3 for any m <= N, then i will be empty.
% Otherwise, i will be a size-[M,N] matrix whose rows represent N-dimensional array
% subscripts of points in y that may be proximate to a global minimum of the
% interpolated y array. (M will be zero if none is found.) The subscripts will only
% address interior points in y, i.e., 1 < i(m,j) < size(y,j) for all m, j. Also, if y
% contains any NaNs, the NaN points and points adjacent to NaNs will be excluded from
% the search. (Use NaNs to represent "missing data".)
%
% If the i_ result is requested, it will similarly contain subscripts for discrete local
% minima in y. The local minima are determined by comparing adjcent y points, including
% diagonal adjacencies. (If two adjacent points have equal y values, only one will be
% counted as a local minimum.) The i result will comprise a subset of the i_ rows.
%
% Version 04/09/2006
% Author: Kenneth C. Johnson
% software.kjinnovation.com
%
% See ZipInterp.pdf, Appendix ("Interpolator application example - global
% optimization"), on the KJ Innovation website for an application example illustrating
% the use of fgridmin.m.
%
% Get y size; define N.
size_y = size(y);
N = length(size_y);
if size_y(N)==1
N = N-1; % column vector
if size_y(N)==1
N = N-1; % scalar
end
end

% If y has no interior points return an empty subscript list.
if N==0 || any(size_y(1:N) < 3)
i = zeros(0,N);
i_ = i;
return
end

% Construct a list of index offsets for adjacent point comparisons. ("Adjacent" includes
% diagonal adjacencies.) Each index offset corresponds to a size-[1,N] subscript offset
% with elements 0, +1, or -1 (not all zero), and with the last non-zero subscript equal
% to +1. (This only includes half of the comparison offsets; the others are obtained by
% sign-inverting the offset list.)
%
% Also construct indices (n) for interior points in y.
%
offset = [];
n = 0;
stride = 1;
for j = 1:N
offset = [offset;[-offset;0;offset]+stride];
n_ = (1:size_y(j)-2)*stride;
n = repmat(n,size(n_))+repmat(n_,size(n));
n = n(:);
stride = stride*size_y(j);
end
n = n+1;

% Select indices of discrete local minima in y; restrict n to selected indices.
for m = 1:length(offset)
% Note: The first comparison below is >=, and second comparison is >, so that if two
% adjacent points have equal y values only one will be counted as a local minimum.
% Negative logic is used (~(...<...). ~(...>=...)) to filter out NaNs and points
% adjacent to NaNs. (Comparisons to NaN always yield false.)
n(~(y(n)<y(n-offset(m)))) = [];
n(~(y(n)<=y(n+offset(m)))) = [];
if isempty(n)
% no local minima
i = zeros(0,N);
i_ = i;
return
end
end

if nargout >= 2
% Convert flat indices n to N-D subscripts i_.
[i_{1:N}] = ind2sub(size_y,n);
i_ = cell2mat(i_);
end

% Determine an estimated upper bound (d) on the variation of the interpolated y from
% each y minimum within a half-step subscript range from each minimum. The d estimate
% is half the maximum difference between the y minimum and any adjacent y value, and the
% interpolated y values proximate to the minima will be assumed to be between y(n)-d and
% y(n)+d.
y_n = y(n);
d = y_n;
for m = 1:length(offset)
d = max(d,y(n-offset(m)));
d = max(d,y(n+offset(m)));
end
d = 0.5*(d-y_n);

% A conservative upper bound on the interpolated global minimum is min(y(n)+d); select
% local minima whose proximate lower bounds (y(n)-d) are within this limit.
n(y_n-d > min(y_n+d)) = [];
clear y_n d

% Convert flat indices n to N-D subscripts i.
[i{1:N}] = ind2sub(size_y,n);
i = cell2mat(i);