Re: missing values in the gamma function

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Tue Feb 05 2013 - 14:25:58 MST

Bummer .... Looks like a bug in the interface to the
'under-the-hood' SLATEC function. Here is
a work around
---------------------------------------------------------------

undef ("gamma_fix")
function gamma_fix(z:numeric)
local x, dimx, rankx, i, x1d
begin
   x = z

   if (isatt(x,"_FillValue")) then
       if (any(ismissing(x))) then
           dimx = dimsizes(x)
           rankx = dimsizes(dimx)

           if (rankx.eq.1) then
               i = ind(.not.ismissing(x))
               x(i) = gamma(x(i))
               return(x)
           else
               x1d = ndtooned(x)
               i = ind(.not.ismissing(x1d))
               x1d(i) = gamma(x1d(i))
               return( onedtond(x1d, dimx) )
           end if
       else
            x = gamma(x)
            return(x)
       end if
   else
       x = gamma(x)
       return(x)
   end if
end

;========= MAIN =======================================

    x = (/ 0.5, 0.33333333, 0.25, 0.2, 2, 3, 4 /)
    gx_a = gamma_fix(x)
    print(gx_a)
    print("================")

    x(1) = 1e20
    x(5) = 1e20
    print("================")
    x@_FillValue = 1e20
    gx_b = gamma_fix(x)
    print(gx_b)

    print("================")
    x@_FillValue = -999.
    gx_c = gamma_fix(x)
    print(gx_c)

;-------------------------------------------------
    nx = dimsizes(x)

    q = conform_dims((/2,3,nx/), x, 2)
    gq = gamma_fix(q)
    print(gq)

On 02/05/2013 08:55 AM, Nick Guy wrote:
> Hello all,
>
> I am passing an array of values to the gamma function as below. Note the -999 causes an error:
>
> ncl 5> a=(/5,2,-999/)
> ncl 6> print(gamma(a))
> ==>SLATEC/NCL: DGAMMA: X SO SMALL GAMMA UNDERFLOWS: NERR= 2 : LEVEL= 1
> (0) 24
> (1) 1
> (2) 0
>
> Now I set the -999 as a missing values and try again, but the error is still returned:
> ncl 7> a@_FillValue=-999
> ncl 8> print(gamma(a))
> ==>SLATEC/NCL: DGAMMA: X SO SMALL GAMMA UNDERFLOWS: NERR= 2 : LEVEL= 1
> (0) 24
> (1) 1
> (2) 0
>
> So I also try using the where command to only process the "good" values:
> ncl 10> b=where(ismissing(a),a@_FillValue,gamma(a))
> ==>SLATEC/NCL: DGAMMA: X SO SMALL GAMMA UNDERFLOWS: NERR= 2 : LEVEL= 1
> ncl 11> print(b)
>
>
> Variable: b
> Type: float
> Total Size: 12 bytes
> 3 values
> Number of Dimensions: 1
> Dimensions and sizes: [3]
> Coordinates:
> (0) 24
> (1) 1
> (2) -999
>
> So the missing value is always processed (if the default missing float value is used, a big gamma underflow occurs). I assume this may be part of the DCDFLIB library and that all array numbers must be passed? In the end NCL does what the command tells it to do, but the resulting output is annoying more than anything.
>
> Is there any way to cut the error by not processing the missing values?
> Thanks,
>
> Nick Guy, Ph.D.
> NOAA National Severe Storms Laboratory
> NOAA/NSSL/WRDD
> 120 David L. Boren Blvd.
> Norman, OK 73072
> (405)325-6287
> nick.guy@noaa.gov
>
>
>
>
>
> _______________________________________________
> 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 Tue Feb 5 14:26:10 2013

This archive was generated by hypermail 2.1.8 : Wed Feb 06 2013 - 16:37:00 MST