[Numpy-discussion] dtype argument description for np.array

[Kirill Balunov]

Currently in docstring the description of dtype argument for np.array says
this:

dtype : data-type, optional
> The desired data-type for the array.  If not given, then the type will
> be determined as the minimum type required to hold the objects in the
> sequence.  This argument can only be used to 'upcast' the array.  For
> downcasting, use the .astype(t) method.
>

But I found this description somewhat misleading for integer types. Is this
generally true that "the type will
be determined as the minimum type required to hold the objects in the
sequence."?

As I always thought for integer arrays first `np.int_` is assumed and if it
is not enough the bigger dtype s are used, and finally falling to `object`
dtype. This question comes from this discussion:
https://github.com/numba/numba/issues/2729.

Also somewhat (un)related (but I was always curious), why `float` type was
chosen as a default for dtype?

>>> np.array([]).dtype
dtype('float64')
>>> np.ones([]).dtype
dtype('float64')

Why np.int_ was not chosen and will this change be a good idea (in general,
without taking into account backward compatibility :-))? There is some
discussion here: https://github.com/numpy/numpy/issues/10405.

-

p.s.: For some constructors the signature looks as follows (in IPython
console):

>>> np.zeros?
Docstring:
zeros(shape, dtype=float, order='C')

Return a new array of given shape and type, filled with zeros.


>>> np.empty?
Docstring:
empty(shape, dtype=float, order='C')

Return a new array of given shape and type, without initializing entries.


But for `ones` the dtype is none instead of float and looks different:

>>> np.ones?
Signature: np.ones(shape, dtype=None, order='C')
Docstring:
Return a new array of given shape and type, filled with ones.

-

With kind regards,
-gdg
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] type promotion rules for integers

[Kirill Balunov]

 What considerations formed the basis for choosing the next type promotion
behavior in numpy:

In[2] : a = np.array([10], dtype=np.int64)
b = np.array([10], dtype=np.uint64)
(a+b).dtype
Out[2]: dtype('float64')

Why the `object` dtype was not chosen for the resulting dtype? Are
there any disadvantages in this case when choosing `object` instead of
`float64`?

With kind regards,

-gdg
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] improving arange()? introducing fma()?

[Chris Barker]

On Wed, Feb 7, 2018 at 12:09 AM, Ralf Gommers 
wrote:
>
>  It is partly a plea for some development of numerically accurate
>> functions for computing lat/lon grids from a combination of inputs: bounds,
>> counts, and resolutions.
>>
>
Can you be more specific about what problems you've run into -- I work with
lat-lon grids all the time, and have never had a problem.

float32 degrees gives you about 1 meter accuracy or better, so I can see
how losing a few digits might be an issue, though I would argue that you
maybe shouldn't use float32 if you are worried about anything close to 1m
accuracy... -- or shift to a relative coordinate system of some sort.

I have been playing around with the decimal package a bit lately,
>>
>
sigh. decimal is so often looked at a solution to a problem it isn't
designed for. lat-lon is natively Sexagesimal -- maybe we need that dtype
:-)

what you get from decimal is variable precision -- maybe a binary variable
precision lib is a better answer -- that would be a good thing to have easy
access to in numpy, but in this case, if you want better accuracy in a
computation that will end up in float32, just use float64.

and I discovered the concept of "fused multiply-add" operations for
>> improved accuracy. I have come to realize that fma operations could be used
>> to greatly improve the accuracy of linspace() and arange().
>>
>
arange() is problematic for non-integer use anyway, by its very definition
(getting the "end point" correct requires the right step, even without FP
error).

and would it really help with linspace? it's computing a delta with one
division in fp, then multiplying it by an integer (represented in fp --
why? why not keep that an integer till the multiply?).

In particular, I have been needing improved results for computing
>> latitude/longitude grids, which tend to be done in float32's to save memory
>> (at least, this is true in data I come across).
>>
>
> If you care about saving memory *and* accuracy, wouldn't it make more
> sense to do your computations in float64, and convert to float32 at the
> end?
>

that does seem to be the easy option :-)


> Now, to the crux of my problem. It is next to impossible to generate a
>> non-trivial numpy array of coordinates, even in double precision, without
>> hitting significant numerical errors.
>>
>
I'm confused, the example you posted doesn't have significant errors...


> Which has lead me down the path of using the decimal package (which
>> doesn't play very nicely with numpy because of the lack of casting rules
>> for it). Consider the following:
>> ```
>> $ cat test_fma.py
>> from __future__ import print_function
>> import numpy as np
>> res = np.float32(0.01)
>> cnt = 7001
>> x0 = np.float32(-115.0)
>> x1 = res * cnt + x0
>> print("res * cnt + x0 = %.16f" % x1)
>> x = np.arange(-115.0, -44.99 + (res / 2), 0.01, dtype='float32')
>> print("len(arange()): %d  arange()[-1]: %16f" % (len(x), x[-1]))
>> x = np.linspace(-115.0, -44.99, cnt, dtype='float32')
>> print("linspace()[-1]: %.16f" % x[-1])
>>
>> $ python test_fma.py
>> res * cnt + x0 = -44.9900015648454428
>> len(arange()): 7002  arange()[-1]:   -44.975044
>> linspace()[-1]: -44.9900016784667969
>> ```
>> arange just produces silly results (puts out an extra element... adding
>> half of the resolution is typically mentioned as a solution on mailing
>> lists to get around arange()'s limitations -- I personally don't do this).
>>
>
The real solution is "don't do that" arange is not the right tool for the
job.

Then there is this:

res * cnt + x0 = -44.9900015648454428
linspace()[-1]: -44.9900016784667969

that's as good as you are ever going to get with 32 bit floats...

Though I just noticed something about your numbers -- there should be a
nice even base ten delta if you have 7001 gaps -- but linspace produces N
points, not N gaps -- so maybe you want:


In [*17*]: l = np.linspace(-115.0, -44.99, 7002)


In [*18*]: l[:5]

Out[*18*]: array([-115.  , -114.99, -114.98, -114.97, -114.96])


In [*19*]: l[-5:]

Out[*19*]: array([-45.03, -45.02, -45.01, -45.  , -44.99])


or, in float32 -- not as pretty:


In [*20*]: l = np.linspace(-115.0, -44.99, 7002, dtype=np.float32)


In [*21*]: l[:5]

Out[*21*]:

array([-115., -114.98999786, -114.98000336, -114.97000122,

   -114.9508], dtype=float32)


In [*22*]: l[-5:]

Out[*22*]: array([-45.02999878, -45.0246, -45.00999832, -45.,
-44.99000168], dtype=float32)


but still as good as you get with float32, and exactly the same result as
computing in float64 and converting:



In [*25*]: l = np.linspace(-115.0, -44.99, 7002).astype(np.float32)


In [*26*]: l[:5]

Out[*26*]:

array([-115., -114.98999786, -114.98000336, -114.97000122,

   -114.9508], dtype=float32)


In [*27*]: l[-5:]

Out[*27*]: array([-45.02999878, -45.0246, -45.00999832, -45.,
-44.99000168], dtype=float32)



>> So, does it make any sense to improve arange by utilizing fma() under the
>> hood?
>>
>

Re: [Numpy-discussion] improving arange()? introducing fma()?

[Benjamin Root]

On Fri, Feb 9, 2018 at 12:19 PM, Chris Barker  wrote:

> On Wed, Feb 7, 2018 at 12:09 AM, Ralf Gommers 
> wrote:
>>
>>  It is partly a plea for some development of numerically accurate
>>> functions for computing lat/lon grids from a combination of inputs: bounds,
>>> counts, and resolutions.
>>>
>>
> Can you be more specific about what problems you've run into -- I work
> with lat-lon grids all the time, and have never had a problem.
>
> float32 degrees gives you about 1 meter accuracy or better, so I can see
> how losing a few digits might be an issue, though I would argue that you
> maybe shouldn't use float32 if you are worried about anything close to 1m
> accuracy... -- or shift to a relative coordinate system of some sort.
>

The issue isn't so much the accuracy of the coordinates themselves. I am
only worried about 1km resolution (which is approximately 0.01 degrees at
mid-latitudes). My concern is with consistent *construction* of a
coordinate grid with even spacing. As it stands right now. If I provide a
corner coordinate, a resolution, and the number of pixels, the result is
not terrible (indeed, this is the approach used by gdal/rasterio). If I
have start/end coordinates and the number of pixels, the result is not bad,
either (use linspace). But, if I have start/end coordinates and a
resolution, then determining the number of pixels from that is actually
tricky to get right in the general case, especially with float32 and large
grids, and especially if the bounding box specified isn't exactly divisible
by the resolution.


>
> I have been playing around with the decimal package a bit lately,
>>>
>>
> sigh. decimal is so often looked at a solution to a problem it isn't
> designed for. lat-lon is natively Sexagesimal -- maybe we need that dtype
> :-)
>
> what you get from decimal is variable precision -- maybe a binary variable
> precision lib is a better answer -- that would be a good thing to have easy
> access to in numpy, but in this case, if you want better accuracy in a
> computation that will end up in float32, just use float64.
>

I am not concerned about computing distances or anything like that, I am
trying to properly construct my grid. I need consistent results regardless
of which way the grid is specified (start/end/count, start/res/count,
start/end/res). I have found that loading up the grid specs (using in a
config file or command-line) using the Decimal class allows me to exactly
and consistently represent the grid specification, and gets me most of the
way there. But the problems with arange() is frustrating, and I have to
have extra logic to go around that and over to linspace() instead.


>
> and I discovered the concept of "fused multiply-add" operations for
>>> improved accuracy. I have come to realize that fma operations could be used
>>> to greatly improve the accuracy of linspace() and arange().
>>>
>>
> arange() is problematic for non-integer use anyway, by its very definition
> (getting the "end point" correct requires the right step, even without FP
> error).
>
> and would it really help with linspace? it's computing a delta with one
> division in fp, then multiplying it by an integer (represented in fp --
> why? why not keep that an integer till the multiply?).
>

Sorry, that was a left-over from a previous draft of my email after I
discovered that linspace's accuracy was on par with fma(). And while
arange() has inherent problems, it can still be made better than it is now.
In fact, I haven't investigated this, but I did recently discover some unit
tests of mine started to fail after a numpy upgrade, and traced it back to
a reduction in the accuracy of a usage of arange() with float32s. So,
something got worse at some point, which means we could still get accuracy
back if we can figure out what changed.


>
> In particular, I have been needing improved results for computing
>>> latitude/longitude grids, which tend to be done in float32's to save memory
>>> (at least, this is true in data I come across).
>>>
>>
>> If you care about saving memory *and* accuracy, wouldn't it make more
>> sense to do your computations in float64, and convert to float32 at the
>> end?
>>
>
> that does seem to be the easy option :-)
>

Kinda missing the point, isn't it? Isn't that like saying "convert all your
data to float64s prior to calling np.mean()"? That's ridiculous. Instead,
we made np.mean() upcast the inner-loop operation, and even allow an option
to specify what the dtype that should be used for the aggregator.


>
>
>> Now, to the crux of my problem. It is next to impossible to generate a
>>> non-trivial numpy array of coordinates, even in double precision, without
>>> hitting significant numerical errors.
>>>
>>
> I'm confused, the example you posted doesn't have significant errors...
>

Hmm, "errors" was the wrong word. "Differences between methods" might be
more along the lines of what I was thinking. Remember, I am looking for
consistency.


>
>
>> Which has lead me down the 

Re: [Numpy-discussion] improving arange()? introducing fma()?

[Matthew Harrigan]

I apologize if I'm missing something basic, but why are floats being
accumulated in the first place?  Can't arange and linspace operations with
floats be done internally similar to `start + np.arange(num_steps) *
step_size`?  I.e. always accumulate (really increment) integers to limit
errors.

On Fri, Feb 9, 2018 at 3:43 PM, Benjamin Root  wrote:

>
>
> On Fri, Feb 9, 2018 at 12:19 PM, Chris Barker 
> wrote:
>
>> On Wed, Feb 7, 2018 at 12:09 AM, Ralf Gommers 
>> wrote:
>>>
>>>  It is partly a plea for some development of numerically accurate
 functions for computing lat/lon grids from a combination of inputs: bounds,
 counts, and resolutions.

>>>
>> Can you be more specific about what problems you've run into -- I work
>> with lat-lon grids all the time, and have never had a problem.
>>
>> float32 degrees gives you about 1 meter accuracy or better, so I can see
>> how losing a few digits might be an issue, though I would argue that you
>> maybe shouldn't use float32 if you are worried about anything close to 1m
>> accuracy... -- or shift to a relative coordinate system of some sort.
>>
>
> The issue isn't so much the accuracy of the coordinates themselves. I am
> only worried about 1km resolution (which is approximately 0.01 degrees at
> mid-latitudes). My concern is with consistent *construction* of a
> coordinate grid with even spacing. As it stands right now. If I provide a
> corner coordinate, a resolution, and the number of pixels, the result is
> not terrible (indeed, this is the approach used by gdal/rasterio). If I
> have start/end coordinates and the number of pixels, the result is not bad,
> either (use linspace). But, if I have start/end coordinates and a
> resolution, then determining the number of pixels from that is actually
> tricky to get right in the general case, especially with float32 and large
> grids, and especially if the bounding box specified isn't exactly divisible
> by the resolution.
>
>
>>
>> I have been playing around with the decimal package a bit lately,

>>>
>> sigh. decimal is so often looked at a solution to a problem it isn't
>> designed for. lat-lon is natively Sexagesimal -- maybe we need that dtype
>> :-)
>>
>> what you get from decimal is variable precision -- maybe a binary
>> variable precision lib is a better answer -- that would be a good thing to
>> have easy access to in numpy, but in this case, if you want better accuracy
>> in a computation that will end up in float32, just use float64.
>>
>
> I am not concerned about computing distances or anything like that, I am
> trying to properly construct my grid. I need consistent results regardless
> of which way the grid is specified (start/end/count, start/res/count,
> start/end/res). I have found that loading up the grid specs (using in a
> config file or command-line) using the Decimal class allows me to exactly
> and consistently represent the grid specification, and gets me most of the
> way there. But the problems with arange() is frustrating, and I have to
> have extra logic to go around that and over to linspace() instead.
>
>
>>
>> and I discovered the concept of "fused multiply-add" operations for
 improved accuracy. I have come to realize that fma operations could be used
 to greatly improve the accuracy of linspace() and arange().

>>>
>> arange() is problematic for non-integer use anyway, by its very
>> definition (getting the "end point" correct requires the right step, even
>> without FP error).
>>
>> and would it really help with linspace? it's computing a delta with one
>> division in fp, then multiplying it by an integer (represented in fp --
>> why? why not keep that an integer till the multiply?).
>>
>
> Sorry, that was a left-over from a previous draft of my email after I
> discovered that linspace's accuracy was on par with fma(). And while
> arange() has inherent problems, it can still be made better than it is now.
> In fact, I haven't investigated this, but I did recently discover some unit
> tests of mine started to fail after a numpy upgrade, and traced it back to
> a reduction in the accuracy of a usage of arange() with float32s. So,
> something got worse at some point, which means we could still get accuracy
> back if we can figure out what changed.
>
>
>>
>> In particular, I have been needing improved results for computing
 latitude/longitude grids, which tend to be done in float32's to save memory
 (at least, this is true in data I come across).

>>>
>>> If you care about saving memory *and* accuracy, wouldn't it make more
>>> sense to do your computations in float64, and convert to float32 at the
>>> end?
>>>
>>
>> that does seem to be the easy option :-)
>>
>
> Kinda missing the point, isn't it? Isn't that like saying "convert all
> your data to float64s prior to calling np.mean()"? That's ridiculous.
> Instead, we made np.mean() upcast the inner-loop operation, and even allow
> an option to specify what the dtype that should be used for the a

Re: [Numpy-discussion] improving arange()? introducing fma()?

[Eric Wieser]

Can’t arange and linspace operations with floats be done internally

Yes, and they probably should be - they’re done this way as a hack because
the api exposed for custom dtypes is here
,
(example implementation here
)
- essentially, you give it the first two elements of the array, and ask it
to fill in the rest.
​

On Fri, 9 Feb 2018 at 13:17 Matthew Harrigan 
wrote:

> I apologize if I'm missing something basic, but why are floats being
> accumulated in the first place?  Can't arange and linspace operations with
> floats be done internally similar to `start + np.arange(num_steps) *
> step_size`?  I.e. always accumulate (really increment) integers to limit
> errors.
>
> On Fri, Feb 9, 2018 at 3:43 PM, Benjamin Root 
> wrote:
>
>>
>>
>> On Fri, Feb 9, 2018 at 12:19 PM, Chris Barker 
>> wrote:
>>
>>> On Wed, Feb 7, 2018 at 12:09 AM, Ralf Gommers 
>>> wrote:

  It is partly a plea for some development of numerically accurate
> functions for computing lat/lon grids from a combination of inputs: 
> bounds,
> counts, and resolutions.
>

>>> Can you be more specific about what problems you've run into -- I work
>>> with lat-lon grids all the time, and have never had a problem.
>>>
>>> float32 degrees gives you about 1 meter accuracy or better, so I can see
>>> how losing a few digits might be an issue, though I would argue that you
>>> maybe shouldn't use float32 if you are worried about anything close to 1m
>>> accuracy... -- or shift to a relative coordinate system of some sort.
>>>
>>
>> The issue isn't so much the accuracy of the coordinates themselves. I am
>> only worried about 1km resolution (which is approximately 0.01 degrees at
>> mid-latitudes). My concern is with consistent *construction* of a
>> coordinate grid with even spacing. As it stands right now. If I provide a
>> corner coordinate, a resolution, and the number of pixels, the result is
>> not terrible (indeed, this is the approach used by gdal/rasterio). If I
>> have start/end coordinates and the number of pixels, the result is not bad,
>> either (use linspace). But, if I have start/end coordinates and a
>> resolution, then determining the number of pixels from that is actually
>> tricky to get right in the general case, especially with float32 and large
>> grids, and especially if the bounding box specified isn't exactly divisible
>> by the resolution.
>>
>>
>>>
>>> I have been playing around with the decimal package a bit lately,
>

>>> sigh. decimal is so often looked at a solution to a problem it isn't
>>> designed for. lat-lon is natively Sexagesimal -- maybe we need that dtype
>>> :-)
>>>
>>> what you get from decimal is variable precision -- maybe a binary
>>> variable precision lib is a better answer -- that would be a good thing to
>>> have easy access to in numpy, but in this case, if you want better accuracy
>>> in a computation that will end up in float32, just use float64.
>>>
>>
>> I am not concerned about computing distances or anything like that, I am
>> trying to properly construct my grid. I need consistent results regardless
>> of which way the grid is specified (start/end/count, start/res/count,
>> start/end/res). I have found that loading up the grid specs (using in a
>> config file or command-line) using the Decimal class allows me to exactly
>> and consistently represent the grid specification, and gets me most of the
>> way there. But the problems with arange() is frustrating, and I have to
>> have extra logic to go around that and over to linspace() instead.
>>
>>
>>>
>>> and I discovered the concept of "fused multiply-add" operations for
> improved accuracy. I have come to realize that fma operations could be 
> used
> to greatly improve the accuracy of linspace() and arange().
>

>>> arange() is problematic for non-integer use anyway, by its very
>>> definition (getting the "end point" correct requires the right step, even
>>> without FP error).
>>>
>>> and would it really help with linspace? it's computing a delta with one
>>> division in fp, then multiplying it by an integer (represented in fp --
>>> why? why not keep that an integer till the multiply?).
>>>
>>
>> Sorry, that was a left-over from a previous draft of my email after I
>> discovered that linspace's accuracy was on par with fma(). And while
>> arange() has inherent problems, it can still be made better than it is now.
>> In fact, I haven't investigated this, but I did recently discover some unit
>> tests of mine started to fail after a numpy upgrade, and traced it back to
>> a reduction in the accuracy of a usage of arange() with float32s. So,
>> something got worse at some point, which means we could still get accuracy
>> back if we can fi

Re: [Numpy-discussion] improving arange()? introducing fma()?

[Benjamin Root]

Interesting...

```
static void
@NAME@_fill(@type@ *buffer, npy_intp length, void *NPY_UNUSED(ignored))
{
npy_intp i;
@type@ start = buffer[0];
@type@ delta = buffer[1];
delta -= start;
for (i = 2; i < length; ++i) {
buffer[i] = start + i*delta;
}
}
```

So, the second value is computed using the delta arange was given, but then
tries to get the delta back, which incurs errors:
```
>>> a = np.float32(-115)
>>> delta = np.float32(0.01)
>>> b = a + delta
>>> new_delta = b - a
>>> "%.16f" % delta
'0.009997764826'
>>> "%.16f" % new_delta
'0.0100021362304688'
```

Also, right there is a good example of where the use of fma() could be of
value.

Cheers!
Ben Root


On Fri, Feb 9, 2018 at 4:56 PM, Eric Wieser 
wrote:

> Can’t arange and linspace operations with floats be done internally
>
> Yes, and they probably should be - they’re done this way as a hack because
> the api exposed for custom dtypes is here
> ,
> (example implementation here
> )
> - essentially, you give it the first two elements of the array, and ask it
> to fill in the rest.
> ​
>
> On Fri, 9 Feb 2018 at 13:17 Matthew Harrigan 
> wrote:
>
>> I apologize if I'm missing something basic, but why are floats being
>> accumulated in the first place?  Can't arange and linspace operations with
>> floats be done internally similar to `start + np.arange(num_steps) *
>> step_size`?  I.e. always accumulate (really increment) integers to limit
>> errors.
>>
>> On Fri, Feb 9, 2018 at 3:43 PM, Benjamin Root 
>> wrote:
>>
>>>
>>>
>>> On Fri, Feb 9, 2018 at 12:19 PM, Chris Barker 
>>> wrote:
>>>
 On Wed, Feb 7, 2018 at 12:09 AM, Ralf Gommers 
 wrote:
>
>  It is partly a plea for some development of numerically accurate
>> functions for computing lat/lon grids from a combination of inputs: 
>> bounds,
>> counts, and resolutions.
>>
>
 Can you be more specific about what problems you've run into -- I work
 with lat-lon grids all the time, and have never had a problem.

 float32 degrees gives you about 1 meter accuracy or better, so I can
 see how losing a few digits might be an issue, though I would argue that
 you maybe shouldn't use float32 if you are worried about anything close to
 1m accuracy... -- or shift to a relative coordinate system of some sort.

>>>
>>> The issue isn't so much the accuracy of the coordinates themselves. I am
>>> only worried about 1km resolution (which is approximately 0.01 degrees at
>>> mid-latitudes). My concern is with consistent *construction* of a
>>> coordinate grid with even spacing. As it stands right now. If I provide a
>>> corner coordinate, a resolution, and the number of pixels, the result is
>>> not terrible (indeed, this is the approach used by gdal/rasterio). If I
>>> have start/end coordinates and the number of pixels, the result is not bad,
>>> either (use linspace). But, if I have start/end coordinates and a
>>> resolution, then determining the number of pixels from that is actually
>>> tricky to get right in the general case, especially with float32 and large
>>> grids, and especially if the bounding box specified isn't exactly divisible
>>> by the resolution.
>>>
>>>

 I have been playing around with the decimal package a bit lately,
>>
>
 sigh. decimal is so often looked at a solution to a problem it isn't
 designed for. lat-lon is natively Sexagesimal -- maybe we need that dtype
 :-)

 what you get from decimal is variable precision -- maybe a binary
 variable precision lib is a better answer -- that would be a good thing to
 have easy access to in numpy, but in this case, if you want better accuracy
 in a computation that will end up in float32, just use float64.

>>>
>>> I am not concerned about computing distances or anything like that, I am
>>> trying to properly construct my grid. I need consistent results regardless
>>> of which way the grid is specified (start/end/count, start/res/count,
>>> start/end/res). I have found that loading up the grid specs (using in a
>>> config file or command-line) using the Decimal class allows me to exactly
>>> and consistently represent the grid specification, and gets me most of the
>>> way there. But the problems with arange() is frustrating, and I have to
>>> have extra logic to go around that and over to linspace() instead.
>>>
>>>

 and I discovered the concept of "fused multiply-add" operations for
>> improved accuracy. I have come to realize that fma operations could be 
>> used
>> to greatly improve the accuracy of linspace() and arange().
>>
>
 arange() is problematic for non-integer use anyway, by its very
 definition (getting the "end point" correct