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