Re: Strange behavior in algebric computations

Hi Mateus,

I think I have to disagree with Mary on this. I do not think it is a
bug, but a natural consequence of NCL's rules for the propagation of
missing values. The only way i think it can be considered a bug is in
that the initial choices of default missing values for NCL values were
a bit myopic. The default missing values should have been outside the
range of "normal" calculations. Likewise we should be encouraging
users to set their own missing values outside the range they expect
their computations to encompass, as Mary does here when she recommends
you use 1e20 rather than -32767 as the missing value. If you look at
this section of the NCL reference manual:

http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclExpressions.shtml#MissingValues

you will see that in the evaluation of any expression if any variable
has a _FillValue attribute it is propagated to the result of the
expression with elements on the left side of each operation prevailing
if there are different fill values. NCL has no way to know that in
this case an intermediate result is just accidently equal to the
_FillValue of the left hand variable in the calculation. And of course
once you have a missing value it is "viral" -- all subsequent
calculations are short-circuited and the missing value is propagated
as the result.

The bottom line is that this behavior results from a fundamental
property of NCL's data model that is built into all the arithmetic
operations. I don't think it could be changed without a radical
reworking of much significant code.

  -dave

On Oct 1, 2009, at 3:17 PM, Mary Haley wrote:

> Mateus,
>
> I've tracked down the source of the problem, and I think it's a bug
> on our end, which I've reported.
>
> The issue is that when dh=2114, and the two v's equal -10.7 and -4.8,
> then this part of the calcuation is equal to exactly -32767, which is
> the missing value for "dh":
>
> dh * ( v(c,t,k-1,j,i) + v(c,t,k,j,i) )
>
> For some reason this makes the rest of the calculation missing.
>
> When you changed it so that it was multiplying by 0.5 first, then 0.5
> got multiplied to dh before the rest of the calculation was done, so
> you didn't end up with -32767.
>
> You've already found one work-around. Another one is to change the
> missing value of "dh" to something other than -32767, like 1e20.
>
> I'll let you know if I hear about a fix, or whether this is a not a
> bug.
>
> --Mary
>
> On Oct 1, 2009, at 10:18 AM, Mateus Teixeira wrote:
>
>> Hi Mary,
>>
>> Here are the results of your codes:
>>
>> ncl 0> su = 0.
>> ncl 1> dh = 1271
>> ncl 2> v1 = 0.699997
>> ncl 3> v2 = -4.8
>> ncl 4> su1 = su + dh * (v1 + v2) / 2
>> ncl 5> su2 = su + 0.5 * dh * (v1 + v2)
>> ncl 6> print("su1 = " + su1)
>> (0) su1 = -2605.55
>> ncl 7> print("su2 = " + su2)
>> (0) su2 = -2605.55
>>
>> ncl 10> su = -2605.55
>> ncl 11> dh = 2114
>> ncl 12> v1 = -4.8
>> ncl 13> v2 = -10.7
>> ncl 14> su1 = su + dh * (v1 + v2) / 2
>> ncl 15> su2 = su + 0.5 * dh * (v1 + v2)
>> ncl 16> print("su1 = " + su1)
>> (0) su1 = -18989.1
>> ncl 17> print("su2 = " + su2)
>> (0) su2 = -18989.1
>>
>> As I mentioned in the previous email, with both formulae the result
>> is
>> the same, when computed by hand, but in the loop it goes to a missing
>> value.
>>
>> And here are the outputs from printVarSummary for all variables used
>> in the formula:
>>
>>
>> Variable: dh
>> Type: float
>> Total Size: 4 bytes
>> 1 values
>> Number of Dimensions: 1
>> Dimensions and sizes: [1]
>> Coordinates:
>> Number Of Attributes: 1
>> _FillValue : -32767
>>
>> Variable: sdh
>> Type: float
>> Total Size: 4 bytes
>> 1 values
>> Number of Dimensions: 1
>> Dimensions and sizes: [1]
>> Coordinates:
>> Number Of Attributes: 1
>> _FillValue : -32767
>>
>>
>> Variable: sv
>> Type: float
>> Total Size: 4 bytes
>> 1 values
>> Number of Dimensions: 1
>> Dimensions and sizes: [1]
>> Coordinates:
>> Number Of Attributes: 1
>> _FillValue : -32767
>>
>>
>> Variable: v (subsection)
>> Type: float
>> Total Size: 4 bytes
>> 1 values
>> Number of Dimensions: 1
>> Dimensions and sizes: [1]
>> Coordinates:
>> Number Of Attributes: 8
>> lon : 302.5
>> lat : -35
>> nivel : 250
>> tempo : 1
>> padrao : 1
>> _FillValue : -32767
>> long_name : vento meridional
>> units : m/s
>>
>>
>> Variable: hgt (subsection)
>> Type: float
>> Total Size: 4 bytes
>> 1 values
>> Number of Dimensions: 1
>> Dimensions and sizes: [1]
>> Coordinates:
>> Number Of Attributes: 8
>> lon : 302.5
>> lat : -35
>> nivel : 300
>> tempo : 1
>> padrao : 1
>> _FillValue : -32767
>> long_name : altura geopotencial
>> units : m
>>
>>
>> And here is the loop
>>
>> sdh = 0.
>> su = 0.
>> sv = 0.
>> dh = 0.
>> do k = k10000, k3000
>> dh = hgt(c,t,k-1,j,i) - hgt(c,t,k,j,i)
>> sdh = sdh + dh
>> su = su + 0.5 * dh * ( u(c,t,k-1,j,i) + u(c,t,k,j,i) )
>> sv = sv + 0.5 * dh * ( v(c,t,k-1,j,i) + v(c,t,k,j,i) )
>> end do
>>
>> And all these lines are inside other loops for c, t, i, and j. I
>> searched for missing values in the data (hgt, u, and v), founding
>> nothing. It is very strange. I thought it was due a error
>> accumulations, but the variables involved are equaled to zero
>> before k
>> loop.
>>
>> Best regards,
>>
>> Mateus
>>
>>
>>
>> 2009/10/1 Mary Haley <haley@ucar.edu>:
>>> Hi Mateus,
>>>
>>> I was unable to reproduce your problem. Can you try this code and
>>> let me
>>> know what you get:
>>>
>>> su = 0.
>>> dh = 1271
>>> v1 = 0.699997
>>> v2 = -4.8
>>> su1 = su + dh * (v1 + v2) / 2
>>> su2 = su + 0.5 * dh * (v1 + v2)
>>> print("su1 = " + su1)
>>> print("su2 = " + su2)
>>>
>>> su = -2605.55
>>> dh = 2114
>>> v1 = -4.8
>>> v2 = -10.7
>>> su1 = su + dh * (v1 + v2) / 2
>>> su2 = su + 0.5 * dh * (v1 + v2)
>>> print("su1 = " + su1)
>>> print("su2 = " + su2)
>>>
>>> My first thought is that this has something to do with the types
>>> of your
>>> variables, and the order of which the arithmetic is being done
>>> on them.
>>>
>>> Can you do a "printVarSummary" on all your variables and send me
>>> the output
>>> so I have the correct types?
>>>
>>> --Mary
>>>
>>>
>>> On Oct 1, 2009, at 8:13 AM, Mateus Teixeira wrote:
>>>
>>>> Dear NCL users,
>>>>
>>>> I'm doing some algebric computations that are having strange
>>>> behavior
>>>> depending on way the formula is written. When written this way
>>>>
>>>> su = su + dh * ( u(c,t,k-1,j,i) + u(c,t,k,j,i) ) / 2
>>>>
>>>> I'm getting a missing value (the formula is inside a do loop), as
>>>> shown
>>>> below:
>>>>
>>>> (0) sv(before)=0; dh=1271; v(k-1)=0.699997; v(k)=-4.8
>>>> (0) sv(after)=-2605.55
>>>> (0) sv(before)=-2605.55; dh=2114; v(k-1)=-4.8; v(k)=-10.7
>>>> (0) sv(after)=-32767 <= MISSING VALUE
>>>>
>>>> But when the formula is written as
>>>>
>>>> su = su + 0.5 * dh * ( u(c,t,k-1,j,i) + u(c,t,k,j,i) )
>>>>
>>>> I haven't a missingl value, as shown below for the same stop
>>>> point showed
>>>> above
>>>>
>>>> (0) sv(before)=0; dh=1271; v(k-1)=0.699997; v(k)=-4.8
>>>> (0) sv(after)=-2605.55
>>>> (0) sv(before)=-2605.55; dh=2114; v(k-1)=-4.8; v(k)=-10.7
>>>> (0) sv(after)=-18989.1 <=CORRECT VALUE
>>>>
>>>> I made these same computations by hand in NCL with both formulae
>>>> and
>>>> get the correct value given by the second way of the formula.
>>>>
>>>> Can anyone imagine what is causing this? I thought that the float
>>>> representation of the 1/2 could be the source of the error, but the
>>>> difference between missing value and the correct value is not
>>>> little.
>>>>
>>>> Best regards,
>>>>
>>>> --
>>>> Mateus da Silva Teixeira
>>>> Meteorologista
>>>> Instituto de Pesquisas Meteorológicas - IPMet
>>>> Universidade Estadual Paulista - UNESP
>>>> Av. Luis Edmundo Carrijo Coube, 14-01 - Cx. Postal 281 - CEP
>>>> 17033-360
>>>> Bauru - SP - Brasil
>>>> Fone: +55 14 3103-6030 / fax: 3203-3649
>>>>
>>>> Registered Linux User #466740 (http://counter.li.org/)
>>>> _______________________________________________
>>>> ncl-talk mailing list
>>>> List instructions, subscriber options, unsubscribe:
>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>>
>>
>>
>>
>> --
>> Mateus da Silva Teixeira
>> Meteorologista
>> Instituto de Pesquisas Meteorológicas - IPMet
>> Universidade Estadual Paulista - UNESP
>> Av. Luis Edmundo Carrijo Coube, 14-01 - Cx. Postal 281 - CEP
>> 17033-360
>> Bauru - SP - Brasil
>> Fone: +55 14 3103-6030 / fax: 3203-3649
>>
>> Registered Linux User #466740 (http://counter.li.org/)
>> _______________________________________________
>> 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

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Oct 1 18:29:26 2009