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