function like accumarray or bincount

Share icon

Wed, 04/12/2017 - 10:03 pm

Does Igor have a function like MATLAB's accumarray or numpy's bincount? I essentially want to accumulate the values of 2D (complex) wave A into a 1D (complex) wave B according to indices in wave C. Therefore, wave A and wave C have the same dimensions. I have implemented this using a nested for loop, but it is somewhat slow. MatrixOp wavemap only seems to be able to sort a 1D wave into a 2D wave from the description?

function map()

    make/O/C/N=(128,128) gauss_wave = 0
    make/O/N=(128,128) map_wave = 0

    gauss_wave += cmplx(exp(-(x^2+y^2)/8192), 1-exp(-(x^2+y^2)/8192))
    map_wave = round(sqrt(x^2+y^2))
    make/O/C/N=(wavemax(map_wave)+1) accumarray

    variable i,j, index
    for(i=0; i<128; i+=1)
        for(j=0; j<128; j+=1)
            accumarray[map_wave[i][j]] = gauss_wave[i][j]
        endfor
    endfor
           

end


Is there any way to use MatrixOp or otherwise speed this up?
I guess your problem sizes are larger than the example? The example is already quick here.

Before we start optimizing the code. Does the code really do what you want?
The problem is symmetrically in i and j so I would have expected something like

accumarray[map_wave[i][j]] += gauss_wave[i][j]


as you otherwise overwrite values in the accumulation array.

April 13, 2017 at 08:04 am - Permalink

Does this not work to eliminate the loops (or am I confused by something otherwise obvious on this late Thursday after a long day)???

accumarray[map_wave[][]] = gauss_wave[p][q]

--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAH

April 13, 2017 at 04:05 pm - Permalink

Thomas, you are correct I left out the accumulation. J.J., your code would work, but I need to accumulate the values into the new array, not overwrite them as Thomas pointed out. The following code does exactly what I need:

function map()
 
    make/O/C/N=(128,128) gauss_wave = 0
    make/O/N=(128,128) map_wave = 0
 
    gauss_wave += cmplx(exp(-(x^2+y^2)/8192), 1-exp(-(x^2+y^2)/8192))
    map_wave = round(sqrt(x^2+y^2))
    make/O/C/N=(wavemax(map_wave)+1) accumarray = 0
    make/O/N=(wavemax(map_wave)+1) accumcount = 0
 
    variable i,j, index
    for(i=0; i<128; i+=1)
        for(j=0; j<128; j+=1)
            accumarray[map_wave[i][j]] += gauss_wave[i][j]
            accumcount[map_wave[i][j]] += 1
        endfor
    endfor
   
    MatrixOp/O accumarray = accumarray/accumcount
 
end


Using the MSTimer, this code takes about ~13 ms. I don't have access to my matlab setup right now but I believe the code runs in less than 1 ms. I will double check and get back.

April 14, 2017 at 02:09 pm - Permalink

Would this improve with threading?

Would this improve with some other implicit indexing trick (by example admittedly off the top of my head ...)?

make/N=128 iwave
make/N=128/D gwave
for (i ...)
       iwave = map_wave[i][p]
       gwave = gauss_wave[i][p]
       accumarray[iwave] += gwave[i]
endfor


--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAH

April 14, 2017 at 03:23 pm - Permalink

jjweimer wrote:
Would this improve with threading?



I'm not sure about threading but it should work better with MatrixOP WaveMap(). What may not be so obvious is that waveMap() works with 1D, e.g.,
make/n=10 ind2=10*(p+1)
make/n=4 ww={2,5,8,9}
matrixop/o aa=waveMap(ind2,ww)


In the case of the code above, there is no obvious reason why one can't execute:
    Duplicate/o map_wave,map1d
    Duplicate/o gauss_wave,g1d
    Redimension/n=(128*128) map1d,g1d


or possibly just redimension the original waves (faster).

April 18, 2017 at 05:44 pm - Permalink