Faster finding of 1D linear interpolation nodes and weights for each element in ND matrix
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
In a problem I'm working on now, I compute some values in a matrix x and I then for each element in x need to find the index of the closest element below in a monotonically increasing vector X as well as the relative proximity of the x elements to the first elements on their either side. (This is essentially linear interpolation without doing the actual interpolation.) I'm doing this maaaany times so I really super extra interested in it being as fast as possible.
I have written a function locate that I can call with some example data:
X = linspace(5, 300, 40);
x = randi(310, 5, 6, 7);
[ii, weights] = locate(x, X);
I have written two versions of locate. The first is for exposition and the second is my best attempt at speeding up the computations. Do you have any suggestions or alternative approaches for how I could accelerate performance further?
1. Exposition
function [ii, weights] = locate(x, X)
% LOCATE Locate first node on grid below a given value.
%
% [ii, weights] = locate(x, X) returns the first node in X that is below
% each element in x and the relative proximities to the two closest nodes.
%
% X must be a monotonically increasing vector. x is a matrix (of any
% order).
% Preallocate
ii = ones(size(x)); % Indices of first node below (or 1 if no nodes below)
weights = zeros([2, size(x)]); % Relative proximity of the two closest nodes
% Find indices and compute weights
for ix = 1:numel(x)
if x(ix) <= X(1)
ii(ix) = 1;
weights(:, ix) = [1; 0];
elseif x(ix) >= X(end)
ii(ix) = length(X) - 1;
weights(:, ix) = [0; 1];
else
ii(ix) = find(X <= x(ix), 1, 'last');
weights(:, ix) = ...
[X(ii(ix) + 1) - x(ix); x(ix) - X(ii(ix))] / (X(ii(ix) + 1) - X(ii(ix)));
end
end
end
2. Best attempt
function [ii, weights] = locate(x, X)
% LOCATE Locate first node on grid below a given value.
%
% [ii, weights] = locate(x, X) returns the first node in X that is below
% each element in x and the relative proximities to the two closest nodes.
%
% X must be a monotonically increasing vector. x is a matrix (of any
% order).
% Preallocate
ii = ones(size(x)); % Indices of first node below (or 1 if no nodes below)
weights = zeros([2, size(x)]); % Relative proximity of the two closest nodes
% Find indices
for iX = 1:length(X) - 1
ii(X(iX) <= x) = iX;
end
% Find weights
below = x <= X(1);
weights(1, below) = 1; % All mass on the first node
weights(2, below) = 0;
above = x >= X(end);
weights(1, above) = 0;
weights(2, above) = 1; % All mass on the last node
interior = ~below & ~above;
xInterior = x(interior)';
iiInterior = ii(interior);
XBelow = X(iiInterior)';
XAbove = X(iiInterior + 1)';
weights(:, interior) = ...
[XAbove - xInterior; xInterior - XBelow] ./ (XAbove - XBelow);
end
4 comentarios
Respuestas (0)
Ver también
Categorías
Más información sobre Numerical Integration and Differentiation en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!