## How to avoid a loop that contains if statements

### Tintin Milou (view profile)

on 17 Feb 2015
Latest activity Commented on by Tintin Milou

on 20 Feb 2015

### Guillaume (view profile)

Hi all,
I've got this loop that I'm trying to get rid of, in the hope of saving on computational time.
Here is the basic idea: I have a vector x_offgrid(:,j1,j2) that contains values that are off a grid. E.g. let the grid be x=[1 3 5 7 9], then x_offgrid might have values [1.7 5.5 3 8.5]. The vector w(:,j1,j2) contains weights associated with x_offgrid. Now, I want to find w_new(:,j1,j2) that gives me weights associated with x.
For now, my idea is to follow the steps
1. Find the nearest neighbor of every value in x_offgrid in x. Call that vector x_ongrid. In my example, that is x_ongrid=[1 5 3 9]. The corresponding index is xind x_ongrid = x(xind)
2. For every value in x_offgrid(i,j1,j2), a fraction of the weight w(i,j1,j2) goes to the nearest neighbor and the other fraction goes to the opposite neighbor (xindmin and xindmax). Those fractions are given by fractiondown and fractionup.
for j1=1:N1
for j2=1:N1
z=zeros(N2);
for i = 1:N2
if x_ongrid(i,j1,j2)>x_offgrid(i,j1,j2) % nearest neighbor is above
z(i,xind(i,j1,j2)) = w(i,j1,j2)*(1-fractiondown(i,j1,j2));
z(i,xindmax(i,j1,j2)) = z(i,xindmax(i,j1,j2)) ...
+ w(i,j1,j2)*fractiondown(i,j1,j2);
else % nearest neighbor is below
z(i,xind(i,j1,j2)) = w(i,j1,j2)*(1-fractionup(i,j1,j2));
z(i,xindmin(i,j1,j2)) = z(i,xindmin(i,j1,j2)) ...
+ w(i,j1,j2)*fractionup(i,j1,j2);
end
end
w_new_tmp(:,j1,:,j2)=z;
end
end
w_new=squeeze(sum(w_new_tmp,1));
Does anybody have a general hint what I could try? Or is it not really possible to rewrite this in a more efficient way? I thought about some type of interpolation function, but right now, I don't see how that could work.
Update: Thanks to Guillaume, I've been able to get rid of the if statement, but not of the loop. Now my code looks as follows:
for j1=1:N1
for j2=1:N1
for i = 1:N2
wn(i,j1,xind(i,j1,j2)+1,j2) = w_tmp(i,j1,j2)*weightup(i,j1,j2);
wn(i,j1,xind(i,j1,j2),j2) = w_tmp(i,j1,j2)*weightdown(i,j1,j2);
end
end
end
The problem is that xind might contain the same index multiple times. I'm not sure if I can make use of accumarray because I want to keep wn is a 4-dimensional matrix.
Thanks!

Show 1 older comment
Tintin Milou

### Tintin Milou (view profile)

on 17 Feb 2015
Hi Image Analyst, you're right. I have edited the question. Let me know if it's clear now.
Guillaume

### Guillaume (view profile)

on 17 Feb 2015
A point of terminology, a vector is a 1D matrix. x_offgrid(:, j1, j2) being a 3d matrix can't be a vector.
Tintin Milou

### Tintin Milou (view profile)

on 17 Feb 2015
Ok, squeeze(x_offgrid(:,j1,j2)) is a vector.

### Guillaume (view profile)

on 17 Feb 2015
Edited by Guillaume

### Guillaume (view profile)

on 17 Feb 2015

Assuming your x is monotically increasing (if not, sort it and save the original indices to reorder the data at the end), and also assuming that x_offgrid is always within [x(1) x(end)], it think this may be what you want:
x = [1 3 5 7 9]; %also ignoring vector vs 3d array, not sure it matters
x_offgrid = [1.7 5.5 3 8.5];
w = [1 10 100 1000];
[~, xind] = histc(x_offgrid, x);
%[~, ~, xind] = histcounts(x_offgrid, x); %using new histcounts in R2014b
%x_ongrid = x(xind); %not used
weightup = w .* (x_offgrid - x(xind)) ./ (x(xind+1) - x(xind));
weightdown = w .* (x(xind+1) - x_offgrid) ./ (x(xind+1) - x(xind));
w_new = accumarray([xind xind+1]', [weightdown weightup])
Note that I've not tried to fully understand your example code (since you've not provided all the details). What I've done is find the lower point corresponding to your x_offgrid using histc. I then calculate a proportion of the weight to redistribute to the upper and lower point according to the distance of x_offgrid to these points (the weightup and weightdown). Finally, the final weight is the sum of the weightup and weightdown corresponding to each grid point. I use accumarray for that.

Tintin Milou

### Tintin Milou (view profile)

on 20 Feb 2015
Good catch. I adjusted the code to deal with out of bounds values. Now, whenever x_offgrid is either below x(1) or above x(end), I set both its upper and lower neighbor to x(1), or x(end) for that matter, so that 100% of the weight goes to that point.
Guillaume

### Guillaume (view profile)

on 20 Feb 2015
There's a few more errors in your code, which made it harder to understand. The initialisation of wn should be
N3 = size(w, 3); %and why is N2, N1 swapped?
wn = zeros(N2, N1, numel(x) + 2, N3);
%and personally, I would have the xind be the 4th dimension:
%wn = zeros([size(x_offgrid) numel(x)+2)]);
and the loop should be
for jp = 1:N3
Anyway, to obtain the same result with accumarray:
[r, c, p] = ndgrid(1:size(w, 1), 1:size(w, 2), 1:size(w, 3));
wn2 = accumarray([r(:) c(:) xind(:) p(:); r(:) c(:) xind(:)+1 p(:)], [wdown(:); wup(:)]);
Tintin Milou

### Tintin Milou (view profile)

on 20 Feb 2015
Nice. Yes, I agree the dimensions were a bit off in my example. But the accumarray function is now working, and the time savings are quite substantial! Thanks.