Index x values of one wave using y values of another

Share icon

I have one wave of 1000 rows = output.
I would like to to set all of the y values in the rows of output whose nearest x values correspond to the y values of a second wave (of 50 rows; raster).

I would like to create a function and to avoid using a for loop because I will use this later in Optimize:
    variable i
    for(i=0; i<numpnts(raster)-1; i+=1)
        variable outputPnt = x2pnt(output, raster[i])
        output[outputPnt] = 1
    endfor
Maybe someone else can come up with a clever way to do that, but as it is I can't see a way to avoid the loop. A small optimization would be to take the call to numpnts() out of the for loop and assign it to a variable:
    variable i
    variable npnts = numpnts(raster)
    for(i=0; i<npnts-1; i+=1)
        ...

Igor's compiler isn't very smart- it will execute the numpnts() function every time through the loop.

Did you intend to terminate your loop one point shy of the end?

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

February 15, 2017 at 10:29 am - Permalink

How about something like this:

Function Test1()
    Make/O/N=2000 Output=p^2
    Make/O/N=200 Raster=100+p*2
    Raster[]=Test2(Output, Raster[p])
    Display Output
end
 
Function Test2(Output, XValue)
Wave Output
Variable XValue
    Output[x2pnt(Output, XValue)]=1
    Return XValue
end


I used the raster wave to make one calculation per data point in the wave, instead of the for loop. You could even MultiThread it if the Raster wave has enough data points for it to be worthwhile. I'm not sure what would happen if you try to MultiThread it and two points in Raster point to the same x-value in Output.

February 16, 2017 at 01:37 am - Permalink

olelytken wrote:
I used the raster wave to make one calculation per data point in the wave, instead of the for loop. You could even MultiThread it if the Raster wave has enough data points for it to be worthwhile. I'm not sure what would happen if you try to MultiThread it and two points in Raster point to the same x-value in Output.


Clever solution, but testing on my machine showed that it was actually almost twice as slow as the original implementation. Is that due to the wave lookup every time Test2 is called? I didn't try multithreading.

Another way I found to slow it down is to pre-fill a wave with the "x2pnt" calculation outside of the loop, then look up those values inside the loop. Curious. Here's my code:
Function wavePrep()
    Make/O/N=5000000 Output=p^2
    SetScale/P x 0,1.1,"", Output
    Make/O/N=20000 Raster=100+sqrt(p*50+1)*p/4
    Make/O/N=20000 RasterPnts=x2pnt(output,raster[p])
end
 
Function testA()
    variable i, timerrefnum, microsec
    wave Output,Raster
    variable pnts=numpnts(raster)
    timerrefnum=startMStimer
    for(i=0; i<pnts-1; i+=1)
        variable outputPnt = x2pnt(output, raster[i])
        output[outputPnt] = 1
    endfor
    microsec=Stopmstimer(timerrefnum)
    Print microSec, "microseconds"
end
 
Function testB()
    variable i, timerrefnum, microsec
    wave Output,RasterPnts
    variable pnts=numpnts(rasterPnts)
    timerrefnum=startMStimer
    for(i=0; i<pnts-1; i+=1)
        output[RasterPnts[i]] = 1
    endfor
    microsec=Stopmstimer(timerrefnum)
    Print microSec, "microseconds"
end
 
Function Test1()
    variable timerrefnum, microsec
    wave Output,Raster
    timerrefnum=startMSTimer
    Raster[]=Test2(Output, Raster[p])
    microsec=Stopmstimer(timerrefnum)
    Print microSec, "microseconds"
end
 
Function Test2(Output, XValue)
    Wave Output
    Variable XValue
    Output[x2pnt(Output, XValue)]=1
    Return Value
end

Timing results:
•waveprep(); testA() //Original method
4706.18 microseconds
•waveprep(); testB() //Original with pre-filled x2pnt
7494.86 microseconds
•waveprep(); test1() //Olelytken method
8225.75 microseconds

I feel that there should be a way to have the pre-filled x2pnt wave (RasterPnts) set the indexed points in Output with one command, perhaps with MatrixOP waveIndexSet, but I couldn't figure out how.

February 16, 2017 at 02:30 pm - Permalink

Here are faster versions of TestA() (original speed 3480 us, new speed 2330 us) and TestB() (original speed 4270 us new speed 1426 us). The pnts-1 calculation in the loop slows it down, so I moved that outside the loop. Also first declaring a variable with the value x2pnt(output, raster[i]) is for some reason much faster than the single command output[x2pnt(output, raster[i])] = 1. I guess a conversion between different number types is taking place?

Function testA()    //  3480 us, 3400 us, 2330 us
    variable i, timerrefnum, microsec, outputPnt
    wave Output,Raster
    //variable pnts=numpnts(raster)
    variable pnts=numpnts(raster)-1
    timerrefnum=startMStimer
    for(i=0; i<pnts; i+=1)
        //variable outputPnt = x2pnt(output, raster[i])
        outputPnt = x2pnt(output, raster[i])
        output[outputPnt] = 1
        //output[x2pnt(output, raster[i])] = 1      //  5000 us
    endfor
    microsec=Stopmstimer(timerrefnum)
    Print microSec, "microseconds"
end
 
 
Function testB() // 4270 us, 1426 us
    variable i, timerrefnum, microsec, outputPnt
    wave Output,RasterPnts
    variable pnts=numpnts(rasterPnts)-1
    timerrefnum=startMStimer
    for(i=0; i<pnts; i+=1)
        //output[RasterPnts[i]] = 1
        outputPnt=RasterPnts[i]
        output[outputPnt] = 1
    endfor
    microsec=Stopmstimer(timerrefnum)
    Print microSec, "microseconds"
end


I usually find that wave indexing is faster than a for loop and if I do a simple test like the one below that is also the case. The clean for loop executes in 1313 us and the wave indexing in 262 us. I'm not sure what the significant difference is between that and the above examples.

Function Test()
 
    Print(" ")
    Print(" ")
 
    Variable i=0
    Variable t=StartMSTimer
    for (i=0; i<20000; i+=1)
    endfor
    Print(StopMSTimer(t))
 
    Make/O/N=2000 TestWave
    t=StartMSTimer
    TestWave[]=TestFunc(TestWave, 1)
    Print(StopMSTimer(t))
end
 
Function TestFunc(TestWave, Value)
    Wave TestWave
    Variable Value
    Return 0
end

February 17, 2017 at 01:19 am - Permalink

olelytken wrote: 

Function Test2(Output, XValue)
        Wave Output
        Variable XValue
    Output[x2pnt(Output, XValue)]=1
    Return XValue
end


Wow. This is very clever. Thanks.

February 21, 2017 at 07:34 am - Permalink