Divide a 3D wave by a 1D wave layer by layer

Share icon

Hello everyone,
I am trying to improve the speed of a function which divides each layer of a 3D wave by a scaler read from a 1D wave with the same number of pnts as the number of layers in the 3D wave. I cant figure out the correct MatrixOp expression or if it is even possible to do this. I can do the manipulation with the p,q,r indexers just fine:


Function testfunc()
	Make/O/N=(256,256,600) images = 0
        Make/O/N=600 wav2=p+1
        Duplicate/O images result
	multithread result=images[p][q][r]/wav2[r]
End

any suggestions?
Why couldn't you make 
wav2
3d and just make the division?

Make/O/N=(256,256,600) w_images = r // each layer has the value from 0 to 600, i.e. r0 = 0, r1 = 1, etc
Make/O/N=(256,256,600) w_wav2 = r/2 //
matrixop/o w_res = w_images / w_wav2 // w_res now contains 2 everywhere except first layer


idk which is faster.

edit:

function test_time()
	variable v_t0 
	variable i=0
	
	
	make/o/n=(256,256,600) w_images = r
	make/o/n=600 w_beam = p/2
	make/o/n=(256,256,600) w_res = 0
	make/o/n=(256,256,600) w_div = r/2
	
	v_t0 = ticks
	for (i=0; i<10; i+=1)
		// multithread w_res = w_images[p][q][r] / w_beam[r] // 1
		// w_res = w_images[p][q][r] / w_beam[r] // 2
		// matrixop/o w_res = w_images / w_div // 3
	endfor
	print (ticks-v_t0)/60/i
end




•test_time()
  1.22833
•test_time()
  4.17167
•test_time()
  0.216667


best,
_sk

March 9, 2018 at 01:34 am - Permalink
I think this does the job

Function Test()
//	Test Function

	Make/O/N=(256,256,600) Images
	
	Make/O/FREE/N=(256) Wave1
	Make/O/FREE/N=(1,1,600) Wave3
	
	FastOP Wave1=1
	MultiThread Wave3=r
	FastOP Wave3=Wave3+1
      
	MatrixOP/O Result=Images / (Wave1 x Wave3 x Wave1^t)
end

March 9, 2018 at 01:52 am - Permalink
thanks _sk and olelytken for your help! turns out the clever matrix math is the fastest even though it requires more lines of code. just changes the time check a little to account for memory allocation time. images and I_total have been previously made. eventually this code will be used for data sets around 10 gb so speed is important:

edit: removed the "= 0" from the make w_res command and everything sped up but option 4 is still the fastest by ~4x


function test()
	Make/O/N=(256,256,600) images = 100
	Make/O/N=(600) I_total = p/100+1
	
	variable v_t0 = ticks
	
	////Option 1:
	//Make/O/N=(256,256,600) w_res = 0
	//Make/FREE/O/N=(256,256,600) wave1 = I_total[r]
	
	//MatrixOp/O w_res = images/wave1
	////Time: 2.97 sec
	////===================================
	
	////Option 2:
//	Make/O/N=(256,256,600) w_res = 0
//	multithread w_res = images[p][q][r]/i_total[r]
	////Time: 1.54 sec
	////===================================
	
	////Option 3:
	//Make/O/N=(256,256,600) w_res = 0
	//w_res = images[p][q][r]/i_total[r]
	////Time: 3.54 sec
	////===================================
	
	////Option 4:
	Make/FREE/O/N=(256) w1
	Make/FREE/O/N=(1,1,600) w2
	
	FastOp w1 = 1
	Multithread w2 = I_total[r]
	
	MatrixOp/O w_res = Images / (w1 x w2 x w1^t)
	////Time: 0.09975 sec
	////===================================
	
	print (ticks-v_t0)/60.15
	
end

March 9, 2018 at 08:43 am - Permalink
Unless I am missing something, it should be pretty simple:


// create sample data:
make/n=(10,20,30) ddd=1+z
make/n=(1,1,30)/o eee=1+z
// Remember that matrixop works on a layer by layer basis:
matrixop/o aa=ddd*rec(eee)

March 9, 2018 at 12:13 pm - Permalink
Igor wrote: Unless I am missing something, it should be pretty simple:


// create sample data:
make/n=(10,20,30) ddd=1+z
make/n=(1,1,30)/o eee=1+z
// Remember that matrixop works on a layer by layer basis:
matrixop/o aa=ddd*rec(eee)


I get an error when I run that "While executing MatrixOP the following error occurred: Bad MatrixOPs token."

PS: I use Igor 6.37

March 19, 2018 at 02:55 am - Permalink
olelytken wrote:
I get an error when I run that "While executing MatrixOP the following error occurred: Bad MatrixOPs token."
PS: I use Igor 6.37


Igor 6.x is a bit old and does not support the full range of options that are available in IP7 or IP8. I really recommend upgrading to IP7. If you upgrade now you will get a free upgrade to IP8.

A.G.

March 19, 2018 at 10:31 am - Permalink