# Re: Strange behavior in algebric computations

From: Mateus Teixeira <mateus.teixeira_at_nyahnyahspammersnyahnyah>
Date: Thu Oct 01 2009 - 16:36:35 MDT

Hi Mary,

Thanks for the help.

Just to inform you:
I have tested this expression putting 0.5 at its end and the same
error arises. At first time I thought that was the division by 2 the
source of error, but it seems that its position defines the result.

Putting 1/2 or 0.5 before dh, eliminates the problem.

Best regards,

Mateus

2009/10/1 Mary Haley <haley@ucar.edu>:
> 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
>>  _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
>>  _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
>>>> 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
>> 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