How to avoid a loop that contains if statements

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!

[Image Analyst]

Explain, in words, what this does. What is TMz, xval, yp, w, etadown, wn, etaup, etc. Right now it just looks like alphabet soup - a jumble of random variables and I have trouble following it because there are no comments.

[Tintin Milou]

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

[Guillaume]

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]

Ok, squeeze(x_offgrid(:,j1,j2)) is a vector.

Answer 1

[Guillaume]

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]

That looks great, Guillaume. Thanks! You even guessed my weights right. There are two issues:

1. x_offgrid can potentially have values that are out of bounds (i.e. either below min(x) or above max(x). But that's a minor adjustment.
2. Eventually, x is part of a 3 dimensional matrix X. Can accumarray deal with 3 dimensional matrices? Or do I still have to loop over j1 and j2 (that means, using your code within two for loops and setting x=X(:,j1,j2) and w=W(:,j1,j2))

[Guillaume]

For 1., you can extend the x (and w) array each way, for example:

x = [-Inf x Inf];
w = [0 w 0];

You just then need to deal with the edges in the weight calculation. I'm not sure what weight you want to assign to these values outside x. Possibly:

 weightdown(xind == 1) = 0;
 weightup(xind == 1) = w(1);
 weightup(xind == numel(w)) = 0;
 weightdown(xind == numel(w)) = w(end);

For 2., I don't see where the dimensionality of the variables matter in the processing. I would just flatten the inputs into vectors, and reshape at the end:

w_new = reshape(w_new, size(x_offgrid));

Otherwise, yes, you can create 3d matrices straight out of accumarray, the columns of the first argument correspond to each dimension. I've no idea what you want to do though.

[Tintin Milou]

Thanks, Guillaume. So yes, I solved the first issue. But the accumarray issue is still a mystery to me. Even if I flatten x_offgrid, I have an issue with w. Since

[~, xind] = histc(x_offgrid(:),[-inf; x; inf]);
xind = min(N2-1,max(1,xind-1)); 
max(xind) = N2,

but

size(w(:)) = N2*N1*N1,

the indicators in xind will never reach anything behind the N2th position in w.

Plus, I want the new w to be a 4D matrix. So here's the code that I want to use accumarray on:

    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

[Guillaume]

I'm afraid I'm not understanding the business with the dimensions. What do they represent?

What I understood is that you have some matrices x_offgrid of a given size and shape, and w of the same size and shape. You also have a vector (any other shape wouldn't make much sense?) of x that you want to move the offgrid data to. In the process, the weight w would be redistributed to the ongrid points. With these premises, the shape of the original data does not matter and at the end, you just do:

x_ongrid = reshape(x_ongrid, size(x_offgrid));
new_w = reshape(new_w, size(x_offgrid));

I don't see where an additional dimension would come from. So can you explain what is going on. A concrete example with small matrices of the right shape would be very helpful.

[Tintin Milou]

Here's a minimal example. I'm not sure how to get from wn2 to wn.

    % Set up variables
    x = [1 3 5 7 9]'; 
    x_offgrid = cat(3,[1.7 5.5; 3.2 7.9; 4.3 9.2],[1.7 5.5; 3.2 7.9; 4.3 9.2]-.5);  

    N1=size(x_offgrid,2); 
    N2=size(x_offgrid,1); 

    w_tmp = [.1 .2; .2 .2; .1 .2]; TM = [.5 .5; .5 .5];
    w = bsxfun(@times,w_tmp,reshape(TM,1,N1,N1));
    wn = zeros(N2,N1,N2,N1); 

    % Find lower and upper neighbor
    [~, xind] = histc(x_offgrid,[-inf; x; inf]); 
    xln=x(max(1,xind-1)); xun=x(min(numel(x),xind)); 
    % for out of bounds values, xln=xun=boundary

    % Calculate weights
    weightup = max(0,min(1,(x_offgrid - xln) ./ (xun - xln)));
    wup =  w.*weightup;
    wdown = w.*(1-weightup);

    % Find new weights
    for j=1:N1
        for jp=1:N1 
            for i = 1:N2    
                wn(i,j,xind(i,j,jp)+1,jp) = wup(i,j,jp);
                wn(i,j,xind(i,j,jp),jp) = wdown(i,j,jp);    
            end     
        end
    end     

    % Alternative approach using accumarray
    x_offgrid = x_offgrid(:);
    w = w(:);

    [~, xind] = histc(x_offgrid,[-inf; x; inf]); 
    xind = min(N2-1,max(1,xind-1)); % xind is the lower neighbor
    xln=x(xind); xun=x(xind+1);

    weightup = max(0,min(1,(x_offgrid - xln) ./ (xun - xln)));
    wup =  w.*weightup;
    wdown = w.*(1-weightup);

    wn2 = accumarray([xind; xind+1], [wdown; wup])

    % Check equality
    min(abs(wn2-reshape(sum(sum(sum(wn,1),2),4),size(wn2))))

    % What I need is wn, not wn2
    reshape(wn2,N2,N1,N2,N1) % ???

[Guillaume]

I'm afraid I really don't understand the purpose of this line:

xind = min(N2-1,max(1,xind-1)); % xind is the lower neighbor

And as a result, I don't get the code that follows. What does xind represents after the line?

[Tintin Milou]

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]

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]

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.