Code covered by the BSD License

# find_idx

### Bruno Luong (view profile)

20 Feb 2009 (Updated )

Fractional binning

idx=find_idx(xi, xgrid, options, InputCheck)
```function idx=find_idx(xi, xgrid, options, InputCheck)
% function idx=find_idx(xi, xgrid)
% function idx=find_idx(xi, xgrid, 'extrap')
%
% Purpose: bins of array 'xi' in the container 'xgrid'
%
% INPUT
%   xgrid: is a sorted array (ascending or descending)
%   xi: array of points
%   CALL idx = find_idx(..., InputCheck)
%      InputCheck [false]: optional boolean flag, force check whereas
%                          xgrid is sorted
% OUTPUT
%   idx: fractional indices of xi in xgrid, same size as xi
%       - The first element starts from 1
%       - The fractional part is calculated as linear interpolation
%         between two brakets
% Algorithm: dichotomy, complexity of m.log(n), where m is number of data
%   points (xi) and n is number of bins (xgrid)
%
% Note 1: This function is equivalent to
%   idx = interp1(xgrid,(1:length(xgrid)),xi,'linear',options)
%   except that overflow values are clipped (interp1 returns NaN).
%   Speed improvement is about 5 times
% Note 2: round(idx) is equivalent to interp1(..., 'nearest')
% Note 3: floor(idx) is equivalent to [trash idx] = histc(xi, xgrid)
%        except overflow values
%
% Call mex function findidxmex.c,
%    compile command: >> mex -O -v findidxmex.c
%
% Example: bining 1 billions data points into 1000 intervals
%   xgrid=cumsum(3*rand(1,1000));
%   xi=rand(1,1e6)*2000;
%   idx=find_idx(xi, xgrid,'extrap');
%
% Author: Bruno Luong
% Original: 19/Feb/2009
%           20/Feb/2009, spare called interp1 when compiled
%              mex engine does not exist

persistent mexexist

% Check if mex file exists
if isempty(mexexist)
mexexist = false;
mexloc = which('findidxmex');
if ~isempty(mexloc)
[pname fname ext] = fileparts(mexloc);
mexexist = ~strcmpi(ext,'.m');
end
if ~mexexist
warning('find_idx:MexNotCompiled', ...
'find_idx: mex engine is not compiled, interp1 invokes');
fprintf('\tRecommendation action: compile then reset by\n');
fprintf('\t>> mex -O -v findidxmex.c\n');
fprintf('\t>> clear find_idx');

end
end

if nargin<4 || isempty(InputCheck)
InputCheck=false; % Not checking by default
end

if InputCheck
if ~isvector(xgrid) || (~issorted(xgrid) && ~issorted(-xgrid))
error('find_idx: xgrid must be sorted vector');
elseif any(diff(xgrid)==0)
error('find_idx: xgrid elements must be distinct');
end
end

nx=length(xgrid);

descending = (nx>0) && xgrid(nx)<xgrid(1);
if descending % decending-order
xgrid = xgrid(end:-1:1); % flip
end

% Cast if necessary
if ~strcmp(class(xgrid),'double')
xgrid = double(xgrid);
end
if ~strcmp(class(xi),'double')
xi = double(xi);
end

% Call the mex engine (dichotomy search)
if mexexist
idx=findidxmex(xgrid(:),xi(:));
else % spare solution
idx=interp1(xgrid,(1:length(xgrid)),xi,'linear','extrap');
end
if isempty(idx) && ~isempty(xi) % problem
error('find_idx: out of memory');
end

if descending % decending-order
idx = nx+1-idx;
end

% Hack the code here if you want different behaviour for overflowed data
if nargin<3 || ~strcmpi(options,'extrap')
idx=max(min(idx,nx),1); % clipping
end

% Reshape the output
idx=reshape(idx,size(xi));

```