MatrixOp wave averaging and getting it really fast

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

Igor

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

thomas_braun

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

olelytken

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

Igor

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

thomas_braun

@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