Non-parametric Wind Regression - speed up calculations

Share icon

Hi all,

I am working on the adaptation of a Python script I made to an Igor procedure for Non-parametric Wind Regression. First, what is NWR? it is a smoothing algorithm to better visualize the coupling of wind data (direction and speed) and pollutant concentrations. I enclosed the formulas associated with the methodology. It is basically, for each (theta,mu) couple a weighted average of concentrations where weighing coefficients are determined with two Kernel functions.

Thus, we can imagine having:
- an angle wave (theta), from 0° to 360° with a resolution of 0.5°, which makes 720 rows
- a speed wave (mu), from 0 to 20km/h with a resolution of 0.1, which makes 200 rows
The final concentration matrix will be (720,200)

So far, the procedure I wrote works pretty well, but it does the calculation cell by cell. Which takes quite a lot of time....
The ideal case would be to calculate the final matrix by only one mathemical operation (matrix multiplication), but it applies dealing with 3D matrices that are too big to handle for my computer (the concentration, ws, wd waves contains 3501 rows, leading to a 720*200*3501-element matrix). An other solution is to fragment the calculation to reduce the size of the 3D matrices.
But still, I am not very familiar with Matrix operations on Igor, and would need any kind of help!

I know my question is a pretty big thing to answer... apologies for that.

I enclose an Igor experiment in which this is:
- the NWR procedure I wrote
- some data (ws, wd, conc waves)

The main function is NWR(sigma,h). The others are here for data preparation.

Many many thanks in advance!
J-E NWR formula.png (32.48 KB) NWR_Igorexchange.pxp (1.33 MB)
Bonjour J-E,

Please post your procedure file (missing from your experiment) if you want help optimizing the code.

A.G.
WaveMetrics, Inc.

March 11, 2015 at 10:22 am - Permalink

I did not check your code or attempt to understand what you are trying to do. I note the following:

1. You need to avoid using angles in degrees. You are wasting computation time on conversion back and forth. Keep everything in radians and modify labels at the very end for the benefit of people who do not understand pi/2.

2. I do not know which data waves are important to you so I did not modify much of the code. It would be better if you wrote your code with waves as function arguments so your inputs are clear.

3. Your goal should be to convert the double loop at the end of the procedure to straight matrix operations. You should be able to do that if you create intermediate 2D waves for diff_speed and the like.

4. Replace wave assignments with MatrixOP wherever possible.

Here is a first iteration of going through your code.
Function init_angle(angle_res)
    Variable angle_res 
    Variable i
    Variable numPoints=360/angle_res
    Make/O/N=(numPoints)/D angle
    for (i=1;i<numPoints;i+=1)
        angle[i]=angle[i-1]+angle_res
    EndFor
   
    Make/O/N=(8)/D angle_ticklocator=0
    for (i=1;i<(numpnts(angle_ticklocator));i+=1)
        angle_ticklocator[i]=angle_ticklocator[i-1]+45
    EndFor
   
    angle_ticklocator=angle_ticklocator/angle_res
    Make/O/T angle_ticklabel = {"N", "NE", "E", "SE", "S", "SW", "W", "NW"}
End

Function init_speed(max_speed,speed_res)
    Variable max_speed,speed_res   
    Variable i
    Variable numPoints=max_speed/speed_res
    Make/N=(numPoints)/D speed
    for (i=1;i<numPoints;i+=1)
        speed[i]=speed[i-1]+speed_res
    EndFor
End

Function init_mx(angle,speed)
    wave angle,speed
   
    Make/N=(numpnts(angle),numpnts(speed))/D mx
    mx=nan
    //wave angle_ticklocator,angle_ticklabel
    //Display;AppendImage mx
    //ModifyGraph userticks(bottom)={angle_ticklocator,angle_ticklabel}
    //ModifyImage mx ctab= {*,*,Rainbow,1}
    //ModifyGraph axisEnab(bottom)={0,0.9}
    //ColorScale/C/N=text0/F=0/A=MC  ctab={0,100,Rainbow,1}
    //ColorScale/C/N=text0/A=RC/X=-1.65/Y=2.84
End

constant kInvsqrt2Pi= 0.3989422804014326779399460599343

Function Kernel1(x)
    variable x
    return exp(-0.5*x*x)*kInvsqrt2Pi
End Function

Function Kernel2(x)
    variable x
    variable cx=0.75*(1-(x*x))
    if (cx<0)
        return 0
    Else
        return cx
    EndIf
End

Function NWR_init()
    wave ws, wd
   
    duplicate ws diff_angle diff_speed angle_Kernel1 speed_Kernel2, Ktot
    diff_angle=nan
    diff_speed=nan
    angle_Kernel1=nan
    speed_Kernel2=nan
    Ktot=nan
End

constant kPi180= 0.017453292519943295769236907

Function angle_vec_diff(a1,a2)
    variable a1,a2
    variable b1,b2
    b1=a1*kPi180
    b2=a2*kPi180
    variable c
    c=cos(b1)*cos(b2)+sin(b1)*sin(b2)
    return acos(c)/kPi180
End

Function NWR(sigma,h)
    variable sigma, h
   
    wave ws, wd, angle, speed, conc, mx, diff_angle, diff_speed, angle_Kernel1, speed_Kernel2, Ktot
    variable i,j
    Variable numPointsAngle=numpnts(angle)
    Variable numPointsSpeed=numpnts(speed)
    for (i=0;i<numPointsAngle;i+=1)
        diff_angle=angle_vec_diff(angle[i],wd)/sigma  // no dependence on j so factor out of loop
        angle_Kernel1=Kernel1(diff_angle)            // no dependence on j so factor out of loop
        for (j=0;j<numPointsSpeed;j+=1)
            diff_speed=(speed[j]-ws)/h
            speed_Kernel2=Kernel2(diff_speed)
            Ktot=speed_Kernel2*angle_Kernel1
            MatrixOP/O/FREE mp=(Ktot^t) x conc
            mx[i][j]=mp[0][0]/sum(Ktot)
        EndFor     
    EndFor
End


I hope this helps,

A.G.
WaveMetrics, Inc.

March 12, 2015 at 11:48 am - Permalink

Allright, thanks a lot for all this!!
As I'm starting programming on Igor, the procedure may be a little messy; but I promise, I'll keep improving!

1/ got it. I am on of these people that are not very confident with radians... :)
2/ the most important waves are the input data (ws, wd and conc), and the "angle" & "speed" waves. The other waves are intermediate for the calculation. the "mx" 2D wave is the final result.
3/ Yes, exactly, and this is the heart the problem. The first thing I could try is to remove the second loop, by doing the calculation angle-by-angle with all speed values at each time. It will lead to the handling of 2D waves which, I guess, is manageable. I'll definitly try that first. The second solution is to remove the 2 loops, but it leads to 3D-matrices that are too big... or am I missing something?
4/ I am not familiar with matrixOP and/or wave assignment... what do you mean by assign waves with MatrixOP?

Thanks again!

March 12, 2015 at 03:24 pm - Permalink

J-E Petit wrote:

3. Yes, exactly, and this is the heart the problem. The first thing I could try is to remove the second loop, by doing the calculation angle-by-angle with all speed values at each time. It will lead to the handling of 2D waves which, I guess, is manageable. I'll definitly try that first. The second solution is to remove the 2 loops, but it leads to 3D-matrices that are too big... or am I missing something?

I did not see any reason for a 3D matrix but then I may have missed something. Can you tell me what is the ultimate goal here? Are you performing some kind of interpolation? Also, your kernel is zero outside a limited domain -- you may be able to save a lot of multiplications by applying it locally.
Quote:

4. I am not familiar with matrixOP and/or wave assignment... what do you mean by assign waves with MatrixOP?

MatrixOP is a built-in operation that speeds up wave assignments. For example
Wave1=wave2+wave3
can be replaced by
MatrixOP/o wave1=wave2+wave3
Here MatrixOP creates the output wave for you (if it does not already exist) and performs the addition. The execution is significantly faster than the wave assignment above. Check out the reference documentation for MatrixOP to see more examples.

March 12, 2015 at 03:48 pm - Permalink

Igor wrote:

I did not see any reason for a 3D matrix but then I may have missed something. Can you tell me what is the ultimate goal here? Are you performing some kind of interpolation? Also, your kernel is zero outside a limited domain -- you may be able to save a lot of multiplications by applying it locally.


yep, it is a sort of interpolation. The ultimate goal is to have, for any wind direction & speed, a concentration value. All this from discreet variables at the input. Does it make sense to you? I enclose an image that may be clearer than some words.
Would the Interpolate XOP be useful here?
NWR.png (165.96 KB)

March 12, 2015 at 04:35 pm - Permalink

J-E Petit wrote:

yep, it is a sort of interpolation. The ultimate goal is to have, for any wind direction & speed, a concentration value. All this from discreet variables at the input. Does it make sense to you? I enclose an image that may be clearer than some words.
Would the Interpolate XOP be useful here?

If the image on the left displays "scatter" of points were data are sampled, I'd approach the problem differently.

March 12, 2015 at 05:01 pm - Permalink

J-E Petit wrote:
Indeed, the left polar graph displays the input data. What would be your approach?


I assume you have 3 waves of data: radius, angle and some scalar data each containing N points.
1. Create a triplet wave:
Make/N=(N,3) tripletWave

2. Fill the third column with scalar data:
tripletWave[][2]=scalarWave[p]

3. compute the x and y cols from radius and angle waves:
tripletWave[][0]=radiusWave[p]*cos(angleWave[p]*pi/180)
tripletWave[][1]=radiusWave[p]*sin(angleWave[p]*pi/180)

4. determine the maximum radius value:
Variable/G rmax=WaveMax(radiusWave)

5. create a regular grid of interpolated data:
ImageInterpolate/S={-rmax,dr,rmax,-rmax,dr,rmax} voronoi tripletWave

6. display the result:
NewImage M_InterpolatedImage


Note that M_InterpolatedImage will contain NaNs outside the convex domain of the data.

HTH,

A.G.
WaveMetrics, Inc.

March 13, 2015 at 09:47 am - Permalink

I think I'm getting what you are writing. However, I get an error message when executing
ImageInterpolate/S={-rmax,dr,rmax,-rmax,dr,rmax} voronoi tripletWave
see attached image

maybe I did something wrong... but it seems like the "dr" variable is not defined. Is it related?
ImageInterpolateError.png (295.83 KB)

March 13, 2015 at 10:42 am - Permalink

(and just making sure that I understood correctly, the triplet wave is based only on the input data, right? i.e. scalarWave: concentration-containing wave; radiusWave: wind speed wave; angleWave: wind direction?)

March 13, 2015 at 10:44 am - Permalink

You need to define
Variable/G dr=...

This sets the resolution of the output and depends on your data.

March 13, 2015 at 12:41 pm - Permalink

Hum, this is weird, I set dr1 & dr2 as:
Variable/G dr1=0.1*cos(0.5*pi/180)
Variable/G dr2=0.1*sin(0.5*pi/180)


and then:
ImageInterpolate/S={-rmax,dr1,rmax,-rmax,dr2,rmax} voronoi tripletWave

but I get the same error, which is not very logical to me, because the ImageInterpolate operation doesn't need a /N flag....

March 15, 2015 at 12:26 am - Permalink

Ah, that's because the /S expects real number arguments. If you are using actual variables then they need to be passed in parenthesis as:
Variable dr=...
imageinterpolate/s={0,(dr),10,0,1,10} voronoi scatterWave

March 15, 2015 at 06:39 pm - Permalink

Okay, thanks!
Tunrs out that calculations are particularly long. Maybe I used dr values that not adapted.
Anyway, I got back to the initial approach, and tried to get rid of one loop (out of the two). It seems to work pretty well, and it is way faster than calculations through each cell (now it is row-by-row). Maybe some lines of code can be optimized?

Next step is to plot the mx wave on a polar plot :) (already in progress on the forum)
constant kPi180= 0.017453292519943295769236907
constant kInvsqrt2Pi= 0.3989422804014326779399460599343


Function init_angle(angle_res)
    Variable angle_res 
    Variable i
    Make/N=(360/angle_res)/D angle
    for (i=1;i<(numpnts(angle));i+=1)
        angle[i]=angle[i-1]+angle_res
    EndFor
    Make/N=(8)/D angle_ticklocator
    for (i=1;i<(numpnts(angle_ticklocator));i+=1)
        angle_ticklocator[i]=angle_ticklocator[i-1]+45
    EndFor
    angle_ticklocator=angle_ticklocator/angle_res
    Make/O/T angle_ticklabel = {"N", "NE", "E", "SE", "S", "SW", "W", "NW"}
End Function

Function init_speed(max_speed,speed_res)
    Variable max_speed,speed_res   
    Variable i
    Make/N=(max_speed/speed_res)/D speed
    for (i=1;i<(numpnts(speed));i+=1)
        speed[i]=speed[i-1]+speed_res
    EndFor
End Function

Function init_mx(angle,speed)
    wave angle,speed
    Make/N=(numpnts(angle),numpnts(speed))/D mx, mx_jp
    mx=nan
    mx_jp=nan
    wave angle_ticklocator,angle_ticklabel
    Display/N=NWR_graph as "Wind apportioned concentrations";AppendImage mx
    ModifyGraph/W=NWR_graph userticks(bottom)={angle_ticklocator,angle_ticklabel}
    ModifyImage/W=NWR_graph mx ctab= {*,*,Rainbow,1}
    ModifyGraph/W=NWR_graph axisEnab(bottom)={0,0.9}
    ModifyGraph/W=NWR_graph fStyle(bottom)=1,fSize(bottom)=14,fStyle(left)=1,fSize(left)=14
    ColorScale/W=NWR_graph/C/N=cb_NWR/F=0/A=RC/X=-9.61/Y=0.51 heightPct=110,image=mx,fsize=14
    ColorScale/W=NWR_graph/C/N=cb_NWR "concentration"
    SetScale/P y 0,0.5,"Wind speed (km/h)", mx
   
    Display/N=JP_graph as "Joint probability";AppendImage mx_jp
    ModifyGraph/W=JP_graph userticks(bottom)={angle_ticklocator,angle_ticklabel}
    ModifyImage/W=JP_graph mx_jp ctab= {*,*,Rainbow,1}
    ModifyGraph/W=JP_graph axisEnab(bottom)={0,0.9}
    ModifyGraph/W=JP_graph fStyle(bottom)=1,fSize(bottom)=14,fStyle(left)=1,fSize(left)=14
    SetScale/P y 0,0.5,"Wind speed (km/h)", mx_jp
   
End Function

Function Kernel1(x)
    variable x
    return exp(-0.5*x*x)*kInvsqrt2Pi
End Function

Function Kernel2(x)
    variable x
    variable cx=0.75*(1-(x*x))
    if (cx<0)
        return 0
    Else
        return cx
    EndIf
End

Function NWR_init()
    wave ws, wd
    duplicate ws diff_angle diff_speed angle_Kernel1 speed_Kernel2, Ktot
    diff_angle=nan
    diff_speed=nan
    angle_Kernel1=nan
    speed_Kernel2=nan
    Ktot=nan
End Function

Function angle_vec_diff(a1,a2)
    variable a1,a2
    variable b1,b2
    b1=a1*kPi180
    b2=a2*kPi180
    variable c
    c=cos(b1)*cos(b2)+sin(b1)*sin(b2)
    return acos(c)/kPi180
End Function

Function meshspeed(x)
    variable x
    wave speed, ws
    Make/N=(numpnts(ws),numpnts(speed)) meshspeed_mx
    variable i
    for(i=0;i<numpnts(speed);i+=1)
        meshspeed_mx[][i]=Kernel2((speed[i]-ws[p])/x)
    endfor
    duplicate meshspeed_mx,Ktot1
End Function


Function NWR3(angle_res,speed_max,speed_res,sigma,h,conc)
    variable angle_res,speed_max,speed_res,sigma, h
    wave conc
    wave ws, wd
    init_angle(angle_res)
    init_speed(speed_max,speed_res)
    wave angle, speed
    init_mx(angle,speed)
    wave mx, mx_jp
    NWR_init()
    wave diff_angle, diff_speed, angle_Kernel1, speed_Kernel2, Ktot
    meshspeed(h)
    wave Ktot1,meshspeed_mx
    variable i
    Variable numPointsAngle=numpnts(angle)
    Variable numPointsSpeed=numpnts(speed)
    Variable kNsigh=numpnts(conc)*sigma*h
    Make/N=(1,numPointsSpeed)/D temp
    for (i=0;i<numPointsAngle;i+=1)
        Ktot1[][]=meshspeed_mx[p][q]
        diff_angle=angle_vec_diff(angle[i],wd)/sigma
        angle_Kernel1=Kernel1(diff_angle)
        Ktot1[][]*=angle_Kernel1[p]
        MatrixOP/O mp=(Ktot1^t) x conc
        MatrixOP/O colSums = SumCols(Ktot1)
        MatrixTranspose mp
        temp=(mp[0][q])/colSums[0][q]
        mx[i][]=temp[0][q]
        mx_jp[i][]=colSums[0][q]/kNsigh
    EndFor 
End

March 17, 2015 at 09:49 am - Permalink

J-E Petit wrote:
Okay, thanks!
Tunrs out that calculations are particularly long. Maybe I used dr values that not adapted.


There are two components for that calculation: triangulation and interpolation. The value of (dr) will only affect the interpolation stage and its time-dependence is linear for each dimension. The triangulation is O(N^2) where N is the number of rows in your triplet wave.

Looking at the last function in your code, here are a few more suggestions you might want to try:
Function NWR3(angle_res,speed_max,speed_res,sigma,h,conc)
    variable angle_res,speed_max,speed_res,sigma, h
    wave conc
    wave ws, wd
    init_angle(angle_res)
    init_speed(speed_max,speed_res)
    wave angle, speed
    init_mx(angle,speed)
    wave mx, mx_jp
    NWR_init()
    wave diff_angle, diff_speed, angle_Kernel1, speed_Kernel2, Ktot
    meshspeed(h)
    wave Ktot1,meshspeed_mx
    variable i
    Variable numPointsAngle=numpnts(angle)
    Variable numPointsSpeed=numpnts(speed)
    Variable kNsigh=numpnts(conc)*sigma*h
    Make/N=(1,numPointsSpeed)/D temp
    // convert to radians only once!
    MatrixOP/O wds=wd*kPi180
    MatrixOP/O angles=angle*kPi180
    // create a single variable that eliminates one multiplication.
    Variable invSigma=1/(sigma*kPi180)
    Variable angles_i
    for (i=0;i<numPointsAngle;i+=1)
        // Duplicate is more efficient though in IP7 (see comment below)
        // this instruction could be skipped completely.
        // Ktot1[][]=meshspeed_mx[p][q]
        Duplicate/O meshspeed_mx,Ktot1

        // create a variable to pass MatrixOP.
        angles_i=angles[i]

        // eliminate function call:
        // diff_angle=angle_vec_diff(angles[i],wds)/sigma
        MatrixOP/O diff_angle=invSigma*acos(cos(angles_i)*cos(wds)+sin(angles_i)*sin(wds))

        // eliminate function call:
        // angle_Kernel1=Kernel1(diff_angle)
        MatrixOP/O angle_Kernel1=exp(-0.5*diff_angle*diff_angle)*kInvsqrt2Pi
        Ktot1[][]*=angle_Kernel1[p] // In IP7 this would be MatrixOP/o Ktot1=scaleRows(Ktot1,angle_Kernel1)
        MatrixOP/O mp=((Ktot1^t) x conc)^t
        MatrixOP/O colSums = SumCols(Ktot1)
        // MatrixTranspose mp
        temp=(mp[0][q])/colSums[0][q]
        mx[i][]=temp[0][q]
        mx_jp[i][]=colSums[0][q]/kNsigh
    EndFor 
End

March 17, 2015 at 10:34 am - Permalink