MatrixOp wave averaging and getting it really fast

Share icon

Hi,

I want to average some columns from some waves. These waves have many points.
I'm currently using fWaveAverage but need to get something quicker (IP7).

I've come up with the following test procedures (number of points are realistic, number of waves are larger):

Function DoStuff()

    variable timer

    variable numPoints = 2e5
    variable numCols   = 12
    variable numWaves  = 3

    Make/N=(numPoints, numCols)/O input1, input2, input3

    input1[][] = enoise(p*q)
    input2[][] = enoise(p*q)
    input3[][] = enoise(p*q)

    /// #1
    timer = startmstimer

    MatrixOp/O input1Single = col(input1, 1)
    MatrixOp/O input2Single = col(input2, 5)
    MatrixOp/O input3Single = col(input3, 10)

    fWaveAverage("input1Single;input2Single;input3Single;", "", 0, 0, "root:average1", "")
    Wave average1

    print stopmstimer(timer)/1e6

    Make/O/N=(numWaves, numPoints) matrix

    /// #2
    timer = startmstimer

    Multithread matrix[0][] = input1[q][1]
    Multithread matrix[1][] = input2[q][5]
    Multithread matrix[2][] = input3[q][10]

    MatrixOp/O average2 = averageCols(matrix)^t
   
    printf "%g, %g\r", stopmstimer(timer)/1e6, EqualWaves(average1, average2, 1)

    /// #3
    timer = startmstimer

    MatrixOp/FREE data = col(input1, 1)
    MultiThread matrix[0][] = data[q]

    MatrixOp/FREE data = col(input2, 5)
    MultiThread matrix[1][] = data[q]

    MatrixOp/FREE data = col(input3, 10)
    MultiThread matrix[2][] = data[q]

    MatrixOp/O average3 = averageCols(matrix)^t

    printf "%g, %g\r", stopmstimer(timer)/1e6, EqualWaves(average1, average3, 1)

    /// #4
    timer = startmstimer

    MatrixOP/FREE matrix = setRow(matrix, 0, col(input1, 1))
    MatrixOP/FREE matrix = setRow(matrix, 1, col(input2, 5))
    MatrixOP/FREE matrix = setRow(matrix, 2, col(input3, 10))

    MatrixOp/O average4 = averageCols(matrix)^t

    printf "%g, %g\r", stopmstimer(timer)/1e6, EqualWaves(average1, average4, 1)

        /// #5
    timer = startmstimer

    MatrixOp/O average5 = averageCols(setRow(setRow(setRow(matrix, 1, col(input2, 5)), 0, col(input1, 1)), 2, col(input3, 10)))^t

    printf "%g, %g\r", stopmstimer(timer)/1e6, EqualWaves(average1, average5, 1)
End


with the results

•dostuff() 0.023102 0.00567959, 1 0.0080306, 1 0.00736661, 1 0.0038697, 1

So 5 looks best, but the problem is that neither the column from the input waves nor the number of input waves are known beforehand.
Is there some way of construction the MatrixOP line from example 5 at runtime? Using something like Execute "MatrixOP" is a bit messy for my taste.

Thanks a lot for reading that far ;)
Thomas
Hello Thomas,

I don't like version (5) because it is difficult to read. I am also confused as to why you chose to use setRow() with ^t when setCol() is an available option and is more efficient for large number of points.

I also don't know that it is more efficient to combine the 3 columns as opposed to executing something like:
MatrixOP/O aa=averagecols(col(inwave,colnum))


A.G.

May 13, 2016 at 01:39 pm - Permalink

Igor wrote:
I don't like version (5) because it is difficult to read.


I'm not fond of that either, very much write only code.

Igor wrote:

I am also confused as to why you chose to use setRow() with ^t when setCol() is an available option and is more efficient for large number of points.


Well in the end I want to have a 1D wave with the averages and averageCols returns a 1xm wave, so I thought transposing it is the right thing to do.
Also if I use setCol I have to transpose matrix again as I want to average not each input wave but the same points from all input waves.

I've added a two more tests:

#include <Waves Average>

Function DoStuff()
 
    variable timer
 
    variable numPoints = 2e5
    variable numCols   = 12
    variable numWaves  = 3
 
    Make/N=(numPoints, numCols)/O input1, input2, input3
 
    input1[][] = enoise(p*q)
    input2[][] = enoise(p*q)
    input3[][] = enoise(p*q)
 
    /// #1
    timer = startmstimer
 
    MatrixOp/O input1Single = col(input1, 1)
    MatrixOp/O input2Single = col(input2, 5)
    MatrixOp/O input3Single = col(input3, 10)
 
    fWaveAverage("input1Single;input2Single;input3Single;", "", 0, 0, "root:average1", "")
    Wave average1
 
    printf "1: %g, %g\r", stopmstimer(timer)/1e6, 1
 
    Make/O/N=(numWaves, numPoints) matrix
 
    /// #2
    timer = startmstimer
 
    Multithread matrix[0][] = input1[q][1]
    Multithread matrix[1][] = input2[q][5]
    Multithread matrix[2][] = input3[q][10]
 
    MatrixOp/O average2 = averageCols(matrix)^t
 
    printf "2: %g, %g\r", stopmstimer(timer)/1e6, EqualWaves(average1, average2, 1)
 
    /// #3
    timer = startmstimer
 
    MatrixOp/FREE data = col(input1, 1)
    MultiThread matrix[0][] = data[q]
 
    MatrixOp/FREE data = col(input2, 5)
    MultiThread matrix[1][] = data[q]
 
    MatrixOp/FREE data = col(input3, 10)
    MultiThread matrix[2][] = data[q]
 
    MatrixOp/O average3 = averageCols(matrix)^t
 
    printf "3: %g, %g\r", stopmstimer(timer)/1e6, EqualWaves(average1, average3, 1)
 
    /// #4
    timer = startmstimer
 
    MatrixOP/FREE matrix = setRow(matrix, 0, col(input1, 1))
    MatrixOP/FREE matrix = setRow(matrix, 1, col(input2, 5))
    MatrixOP/FREE matrix = setRow(matrix, 2, col(input3, 10))
 
    MatrixOp/O average4 = averageCols(matrix)^t
 
    printf "4: %g, %g\r", stopmstimer(timer)/1e6, EqualWaves(average1, average4, 1)
 
    /// #5
    timer = startmstimer
 
    MatrixOp/O average5 = averageCols(setRow(setRow(setRow(matrix, 1, col(input2, 5)), 0, col(input1, 1)), 2, col(input3, 10)))^t
 
    printf "5: %g, %g\r", stopmstimer(timer)/1e6, EqualWaves(average1, average5, 1)

    Make/O/N=(numPoints, numWaves) matrix

    /// #6
    timer = startmstimer

    MatrixOP/FREE matrix = setCol(matrix, 0, col(input1, 1))
    MatrixOP/FREE matrix = setCol(matrix, 1, col(input2, 5))
    MatrixOP/FREE matrix = setCol(matrix, 2, col(input3, 10))

    MatrixOp/O average6 = averageCols(matrix^t)^t
 
    printf "6: %g, %g\r", stopmstimer(timer)/1e6, EqualWaves(average1, average6, 1)

    /// #7
    timer = startmstimer

    MatrixOp/O average7 = averageCols(setCol(setCol(setCol(matrix, 1, col(input2, 5)), 0, col(input1, 1)), 2, col(input3, 10))^t)^t

    printf "7: %g, %g\r", stopmstimer(timer)/1e6, EqualWaves(average1, average7, 1)
End


and these now run with

•dostuff() 1: 0.0247083, 1 2: 0.00551995, 1 3: 0.00792857, 1 4: 0.00731161, 1 5: 0.00384668, 1 6: 0.00865464, 1 7: 0.00467676, 1

May 16, 2016 at 12:07 pm - Permalink

This is how I would do it, but I am using Igor 6.37 not Igor 7.

It was a little difficult to understand what you wanted; I hope I got it right.

Generally I noticed that FastOP is faster than MatrixOP for some of the operations. I wanted to MultiThread the heaviest calculation, but apparently that has almost no effect on MatrixOP calculations. I also noticed that it was faster to do one operation at a time with FastOP and MatrixOP instead of multiple operations, e.g. it's faster to do a=a+b and then a=a+c, than it is to do them in one step a=a+b+c.

If you want to use EqualWaves, I suggest doing it outside the timed sequence, since EqualWaves took a not insignificant amount of time to calculate when I tested it.

Function DoStuff()
 
    variable numRows=2e5
    variable numCols=12
    variable numWaves=30
 
    //  Creates the input waves
    Make/O/FREE/WAVE/N=(numWaves) InputWaves
    MultiThread InputWaves[]=DoStuffCreateInputWaves(p, numRows, numCols)       //  2.6 seconds
   
    //  Starts the timer
    Variable t=StartMSTimer

    //  Sums the rows of each input wave, one at a time. I meant to use multithreading, but it showed no improvement in speed
    Duplicate/O/FREE/WAVE InputWaves, SumRowWaves       //  0.02 ms
    SumRowWaves[]=DoStuffSumRows(InputWaves[p])         //  121 ms

    //  Adds all the summed rows together
    Wave SumRowWave=SumRowWaves[0]
    Duplicate/O SumRowWave, AverageWave //  0.5 ms
    Variable i=0
    for (i=1; i<numWaves; i+=1)
        Wave SumRowWave=SumRowWaves[i]
        FastOP AverageWave=AverageWave+SumRowWave
    endfor  //  6.5 ms
   
    //  Normalizes with the total number of columns. If the number of columns are not the same for all waves, you will have to calculate the sum instead
    Variable a=1/(numWaves*numCols)
    FastOP AverageWave=(a)*AverageWave      //  0.14 ms
   
    //  Displays the elapsed time in miliseconds. 127 ms overall for 30 waves
    Print(StopMSTimer(t)/1000)
    Print(" ")
end

ThreadSafe Function/WAVE DoStuffCreateInputWaves(WaveNum, numRows, numCols)
Variable WaveNum, numRows, numCols
    Make/O/FREE/N=(numRows, numCols) InputWave
    InputWave[][] = enoise(p*q)
    Return InputWave
end

Function/WAVE DoStuffSumRows(InputWave)
Wave InputWave
    MatrixOP/O/FREE SumRowWave=sumRows(InputWave)
    Return SumRowWave
end

May 17, 2016 at 04:49 am - Permalink

thomas_braun wrote:

Well in the end I want to have a 1D wave with the averages and averageCols returns a 1xm wave, so I thought transposing it is the right thing to do.
Also if I use setCol I have to transpose matrix again as I want to average not each input wave but the same points from all input waves.


FWIW, I just added one further optimization in MatrixOP for the special case of transposing 1D arrays.

A.G.

May 17, 2016 at 11:03 am - Permalink

@olelytken: Thanks for your solution, but I want to average along the number of input waves and not along the number of points in each wave. Regarding EqualWaves, I *think* it is executed after stopmstimer.

@AG: Thanks, I'll give it a try.

May 20, 2016 at 03:55 pm - Permalink