Re: computing bivariate frequencies

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Tue, 8 Aug 2006 09:49:16 -0600 (MDT)

> I'm trying to compute the joint pdf of two variables C and V stored in arrays
> sized (ntime x nlat x nlon). In order to obtain the absolute frequencies, I
> need to count the number of grid cells that fall within arbitralily defined
> bins for each variable. Right now I'm attempting this using a nested do loop,
> like:
>
> do ic=0,nbin_c-1
> do iv=0,nbin_v-1
> count_c_v(ic,iv)=num(ind(c.ge.c_bin(ic).and.c.lt.c_bin(ic+1).and.
> v.ge.v_bin(iv).and.v.lt.v_bin(iv+1)))
> end do
> end do
>
> I was wondering if anybody would know how to code this more efficiently on
> ncl, in order to avoid the nested loop.
>
> ----------------------------------------------------------------

[1] I do not know why you have the "ind" function
     within in the "num" function. Not needed. Extra
     operations so slower.

[2] I can not think of an alternative to the double do loop.
     However, I did run the following: 10000 points and 2000 bins
     and it ran pretty quickly ... even with 'if' and 'print'
     within the inner most loop.

     I deliberately chose c_bin and v_bin such that there would
     be outliers. Thus, the value of 'K' will not add up to 10000.

begin
     N = 10000
     c = random_normal ( 0,5,N)
     v = random_normal ( 0,5,N)
    ;print("c="+c+" v="+v)

     nc_bin = 51 ; => 50 bins
     nv_bin = 41 ; => 40 bins
     c_bin = fspan(-10,10,nc_bin)
     v_bin = fspan(-5 ,20,nv_bin)

     count = new ( (/nc_bin,nv_bin/) , integer)
     count = 0 ; initialize to 0

     do ic=0,nc_bin-2
       do iv=0,nv_bin-2
          count(ic,iv) = num(c.ge.c_bin(ic) .and. c.lt.c_bin(ic+1) .and. \
                             v.ge.v_bin(iv) .and. v.lt.v_bin(iv+1) )
          if (count(ic,iv).gt.0) then
              print("ic="+ic+" iv="+iv+" count="+count(ic,iv) )
          end if
       end do
     end do

     fpc = (100.*count)/N ; frequency [%] (nc_bin,nv_bin)
     printVarSummary(fpc)

     K = sum( count ) ; < N
     F = sum( fpc )
     print("------")
     print("K="+K+" F[%]="+F)
     print("------")

   ; extra counting stuff

     cmin = min(c)
     cmax = max(c)
     cmin = min(v)
     cmax = max(v)

     cbmin = min(c_bin)
     cbmax = max(c_bin)
     vbmin = min(v_bin)
     vbmax = max(v_bin)
                                                 ; count corner cases
     k1 = num (c.lt.cbmin .and. v.lt.vbmin)
     k2 = num (c.lt.cbmin .and. v.ge.vbmax)
     k3 = num (c.ge.cbmax .and. v.lt.vbmin)
     k4 = num (c.ge.cbmax .and. v.ge.vbmax)
                                                 ; one outlier
     k5 = num (c.lt.cbmin .and. v.ge.vbmin .and. v.lt.vbmax)
     k6 = num (c.ge.cbmax .and. v.ge.vbmin .and. v.lt.vbmax)
     k7 = num (v.lt.vbmin .and. c.ge.cbmin .and. c.lt.cbmax)
     k8 = num (v.ge.vbmax .and. c.ge.cbmin .and. c.lt.cbmax)

     KTOT = K + k1+ k2+ k3+ k4+ k5+ k6+ k7+ k8 ; should = N

     print("k1="+k1)
     print("k2="+k2)
     print("k3="+k3)
     print("k4="+k4)
     print("k5="+k5)
     print("k6="+k6)
     print("k7="+k7)
     print("k8="+k8)

     print("------")
     print("KTOT="+KTOT)
     print("------")

     if (KTOT.eq.N) then
         print (" Lucky!!! All numbers accounted for!! :-)")
     end if
end
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Aug 08 2006 - 09:49:16 MDT

This archive was generated by hypermail 2.2.0 : Wed Aug 09 2006 - 08:25:35 MDT