# MatrixOp wave averaging and getting it really fast

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

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)

MatrixOp/FREE data = col(input2, 5)

MatrixOp/FREE data = col(input3, 10)

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.
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

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)

MatrixOp/FREE data = col(input2, 5)

MatrixOp/FREE data = col(input3, 10)

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 ```
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

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
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.

@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.