Re: Gaussian smoother

On Fri, 7 Nov 2008, Correia, James wrote:

> All-
> I wrote some code in NCL to do a gaussian smoother but it is very SLOW. It
> does not handle missing yet.
>
> Does anyone know if there is a faster way to implement it:
>
> ;Tcc(:,:,:,:) is the array on which the smoother is operating
> ;It has dimensions of 5,224,299,49
>
> lam = 3. ; shape factor
> fac = new((/6/),float)
> do n=0,5
> ;compute the factors based on radius from central grid point
> fac(n) = exp(-1.*(abs(int2flt(n))/lam)*(abs(int2flt(n))/lam))
> end do

James,

To improve a slow script:

   - look for ways to remove calculations from a "do" loop
   - avoid nested "do" loops

You can do array arithmetic outside of a do loop, as long as the
calculations do not depend on some previous calculation in the
loop. Hence, the above code can be rewritten as:

   lam = 3. ; shape factor
   xn = ispan(0,5,1)*1.
   fac = exp(-1.*(abs(xn)/lam)*(abs(xn)/lam))

> tmp = new((/5,49/),float)
> tm1 = new((/5,49/),float)
>
> do x=5,298-5
> do y=5,223-5
> tmp(:,:) = 0.
> tm1(:,:) = 0.

Whenever you use (:,:) to reference every element of an array, you may
cause NCL to go unnecessarily into subscripting mode. If you indeed
need to do something to every element of the array, just use the array
without the (:,:)

   do x=5,298-5
   do y=5,223-5
     tmp = 0.
     tm1 = 0.

Hence lines like:

  tmp(:,:) = tmp(:,:) + tcc(:,y,x,:)*fac(rada)
  tm1(:,:) = tm1(:,:) + fac(rada)
  tcc(:,y,x,:) = tmp(:,:)/tm1(:,:)

should be rewritten as:

  tmp = tmp + tcc(:,y,x,:)*fac(rada)
  tm1 = tm1 + fac(rada)
  tcc(:,y,x,:) = tmp/tm1

> do xx=x-5,x+5
> do yy=y-5,y+5
> rad = sqrt((xx-x)*(xx-x) + (yy-y)*(yy-y))
> if(rad .lt. 6)then
> rada = floattointeger(rad)
> tmp(:,:) = tmp(:,:) + tcc(:,y,x,:)*fac(rada)
> tm1(:,:) = tm1(:,:) + fac(rada)
> end if
> end do
> end do
> tcc(:,y,x,:) = tmp(:,:)/tm1(:,:)
>end do
>end do

You can eliminate these two inner loops. The new code would look like
this:

   lam = 3. ; shape factor
   xn = ispan(0,5,1)*1.
   fac = exp(-1.*(abs(xn)/lam)*(abs(xn)/lam))

   tmp = new((/5,49/),float)
   tm1 = new((/5,49/),float)

   tmp_xy = new((/11,11/),integer)

   do x=5,298-5
     xx = ispan(x-5,x+5,1)
     xx1d = ndtooned(conform(tmp_xy,xx,0))
     do y=5,223-5
       tmp = 0.
       tm1 = 0.

       yy = ispan(y-5,y+5,1)
       yy1d = ndtooned(conform(tmp_xy,yy,1))
       rad = sqrt((xx1d-x)*(xx1d-x) + (yy1d-y)*(yy1d-y))

       ii = ind(rad.lt.6) ; Get indexes where rad < 6
       rada = floattointeger(rad(ii))
       do i=0,dimsizes(rada)-1
         tmp = tmp + tcc(:,y,x,:)*fac(rada(i))
       end do

       tm1 = sum(fac(rada))
       tcc(:,y,x,:) = tmp/tm1

       delete(ii) ; Remove before next round of loop
       delete(rada)
     end do
   end do

There are probably even more ways you can improve the above,
but hopefully you get the idea.

Using some test data, the times went from about 2.7 minutes
to 1.3 minutes. Not great, but a start.

--Mary

162.123u 0.419s 2:43.37 99.4% 0+0k 0+0io 0pf+0w

> Thanks
> jimmyc
>
>
> James Correia Jr., PhD
> Climate Physics Group
> Post. Doc.
> Pacific Northwest National Lab
>
> "Wisdom. Strength. Courage. Generosity. Each of us are born with one of
> these. It is up to us to find the other three inside of us."
> -Into the West
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sun Nov 09 2008 - 05:54:17 MST