Re: [Numpy-discussion] is_triangular, is_diagonal, is_symmetric et al. in NumPy or SciPy linalg

[Evgeni Burovski]

Hi Ilhan,

Overall I think something like this would be great. However, I wonder
if you considered having a specialized container with a structure tag
instead of trying to discover the structure. If it's a container, it
can neatly wrap various lapack storage schemes and dispatch to an
appropriate lapack functionality. Possibly even sparse storage
schemes. And it seems a bit more robust than trying to discover the
structure (e.g. what about off-band elements of  \sim 1e-16 etc).

The next question is of course if this should live in scipy/numpy
.linalg or as a separate repo, at least for some time (maybe in the
scipy organization?). So that it can iterate faster, among other
things.
(I'd be interested in contributing FWIW)

Cheers,

Evgeni


On Wed, Jun 30, 2021 at 1:22 AM Ilhan Polat  wrote:
>
> Dear all,
>
> I'm writing some helper Cythpm functions for scipy.linalg which is kinda 
> performant and usable. And there is still quite some wiggle room for more.
>
> In many linalg routines there is a lot of performance benefit if the 
> structure can be discovered in a cheap and reliable way at the outset. For 
> example if symmetric then eig can delegate to eigh or if triangular then 
> triangular solvers can be used in linalg.solve and lstsq so forth
>
> Here is the Cythonized version for Jupyter notebook to paste to discover the 
> lower/upper bandwidth of square array A that competes well with A != 0 just 
> to use some low level function (note the latter returns an array hence more 
> cost is involved) There is a higher level supervisor function that checks 
> C-contiguousness otherwise specializes to different versions of it
>
> Initial cell
>
> %load_ext Cython
> %load_ext line_profiler
> import cython
> import line_profiler
>
> Then another cell
>
> %%cython
> # cython: language_level=3
> # cython: linetrace=True
> # cython: binding = True
> # distutils: define_macros=CYTHON_TRACE=1
> # distutils: define_macros=CYTHON_TRACE_NOGIL=1
>
> cimport cython
> cimport numpy as cnp
> import numpy as np
> import line_profiler
> ctypedef fused np_numeric_t:
> cnp.int8_t
> cnp.int16_t
> cnp.int32_t
> cnp.int64_t
> cnp.uint8_t
> cnp.uint16_t
> cnp.uint32_t
> cnp.uint64_t
> cnp.float32_t
> cnp.float64_t
> cnp.complex64_t
> cnp.complex128_t
> cnp.int_t
> cnp.long_t
> cnp.longlong_t
> cnp.uint_t
> cnp.ulong_t
> cnp.ulonglong_t
> cnp.intp_t
> cnp.uintp_t
> cnp.float_t
> cnp.double_t
> cnp.longdouble_t
>
>
> @cython.linetrace(True)
> @cython.initializedcheck(False)
> @cython.boundscheck(False)
> @cython.wraparound(False)
> cpdef inline (int, int) band_check_internal(np_numeric_t[:, ::1]A):
> cdef Py_ssize_t n = A.shape[0], lower_band = 0, upper_band = 0, r, c
> cdef np_numeric_t zero = 0
>
> for r in xrange(n):
> # Only bother if outside the existing band:
> for c in xrange(r-lower_band):
> if A[r, c] != zero:
> lower_band = r - c
> break
>
> for c in xrange(n - 1, r + upper_band, -1):
> if A[r, c] != zero:
> upper_band = c - r
> break
>
> return lower_band, upper_band
>
> Final cell for use-case ---
>
> # Make arbitrary lower-banded array
> n = 50 # array size
> k = 3 # k'th subdiagonal
> R = np.zeros([n, n], dtype=np.float32)
> R[[x for x in range(n)], [x for x in range(n)]] = 1
> R[[x for x in range(n-1)], [x for x in range(1,n)]] = 1
> R[[x for x in range(1,n)], [x for x in range(n-1)]] = 1
> R[[x for x in range(k,n)], [x for x in range(n-k)]] = 2
>
> Some very haphazardly put together metrics
>
> %timeit band_check_internal(R)
> 2.59 µs ± 84.7 ns per loop (mean ± std. dev. of 7 runs, 10 loops each)
>
> %timeit np.linalg.solve(R, zzz)
> 824 µs ± 6.24 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
>
> %timeit R != 0.
> 1.65 µs ± 43.1 ns per loop (mean ± std. dev. of 7 runs, 100 loops each)
>
> So the worst case cost is negligible in general (note that the given code is 
> slower as it uses the fused type however if I go with tempita standalone 
> version is faster)
>
> Two questions:
>
> 1) This is missing np.half/float16 functionality since any arithmetic with 
> float16 is might not be reliable including nonzero check. IS it safe to view 
> it as np.uint16 and use that specialization? I'm not sure about the sign bit 
> hence the question. I can leave this out since almost all linalg suite 
> rejects this datatype due to well-known lack of supprt.
>
> 2) Should this be in NumPy or SciPy linalg? It is quite relevant to be on 
> SciPy but then again this stuff is purely about array structures. But if the 
> opinion is for NumPy then I would need a volunteer because NumPy codebase 
> flies way above my head.
>
>
> All feedback welcome
>
> Best
> ilhan
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@python.org
> https://

Re: [Numpy-discussion] Interface for wrapping random distribution functions (__array_function__, __ufunc__, ?)

[Robert Kern]

On Wed, Jun 30, 2021 at 3:34 PM Yoann Mocquin  wrote:

>
>
> Hello there,
>
> here is a feature request for the possibility to wrap numpy random
> distributions to be wrapped by a mechanism like __array_function__ or
> __array_ufunc__. You can find the GH issue at :
> https://github.com/numpy/numpy/issues/19382.
>
>
>
> This post follows [this one](https://github.com/numpy/numpy/issues/18902).
> I figured I should open another issue for `np.random.normal`, but this
> applies for all distributions I guess.
>
> ## Feature
>
> Basicaly, I would like that `numpy.random.*` distributions could trigger
> an interface when passed instances from custom classes, "à la"
> `__array_ufunc__` or `__array_function__`. Here is an example concept :
>

 The main problem is that these are not functions but methods on a hidden
global `RandomState` instance. I haven't kept up fully with all of the
efforts in the `__array_*__` proposals, but I do see that methods are
explicitly excluded from the `__array_function__` mechanism:

  https://numpy.org/neps/nep-0018-array-function-protocol.html#non-goals

The `RandomState` functionality (and thus all of these aliases in
`numpy.random`) are now frozen in functionality (per NEP 19), so we will
not be adding this functionality to them. If the `__array_function__`
developments eventually work out a good way to wrap methods, then we can
think about using that on `Generator`, but I suspect that will not be
straightforward.

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

Re: [Numpy-discussion] Interface for wrapping random distribution functions (__array_function__, __ufunc__, ?)

[Kevin Sheppard]

There is an issue that suggests reimplementing many of the Generator
methods as ufuncs so they they would inherit all of the free features of
ufuncs. It seems like a doable but large project for some useful but
nonessential features.

Kevin


On Wed, Jun 30, 2021, 20:49 Robert Kern  wrote:

> On Wed, Jun 30, 2021 at 3:34 PM Yoann Mocquin  wrote:
>
>>
>>
>> Hello there,
>>
>> here is a feature request for the possibility to wrap numpy random
>> distributions to be wrapped by a mechanism like __array_function__ or
>> __array_ufunc__. You can find the GH issue at :
>> https://github.com/numpy/numpy/issues/19382.
>>
>>
>>
>> This post follows [this one](https://github.com/numpy/numpy/issues/18902).
>> I figured I should open another issue for `np.random.normal`, but this
>> applies for all distributions I guess.
>>
>> ## Feature
>>
>> Basicaly, I would like that `numpy.random.*` distributions could trigger
>> an interface when passed instances from custom classes, "à la"
>> `__array_ufunc__` or `__array_function__`. Here is an example concept :
>>
>
>  The main problem is that these are not functions but methods on a hidden
> global `RandomState` instance. I haven't kept up fully with all of the
> efforts in the `__array_*__` proposals, but I do see that methods are
> explicitly excluded from the `__array_function__` mechanism:
>
>   https://numpy.org/neps/nep-0018-array-function-protocol.html#non-goals
>
> The `RandomState` functionality (and thus all of these aliases in
> `numpy.random`) are now frozen in functionality (per NEP 19), so we will
> not be adding this functionality to them. If the `__array_function__`
> developments eventually work out a good way to wrap methods, then we can
> think about using that on `Generator`, but I suspect that will not be
> straightforward.
>
> --
> Robert Kern
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion
>
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] Interface for wrapping random distribution functions (__array_function__, __ufunc__, ?)

[Yoann Mocquin]

Hello there,

here is a feature request for the possibility to wrap numpy random 
distributions to be wrapped by a mechanism like __array_function__ or 
__array_ufunc__. You can find the GH issue at : 
https://github.com/numpy/numpy/issues/19382 
.



This post follows [this one](https://github.com/numpy/numpy/issues/18902 
). I figured I should open another 
issue for `np.random.normal`, but this applies for all distributions I guess.

## Feature

Basicaly, I would like that `numpy.random.*` distributions could trigger an 
interface when passed instances from custom classes, "à la" `__array_ufunc__` 
or `__array_function__`. Here is an example concept : 

```python
arr = np.arange(10)

# Reference result 
print(np.mean(arr))
print(np.random.normal(arr))

custom_obj = MyArrayLike(arr)
print(np.mean(custom_obj))   # OK : np.mean will trigger 
__array_function__ interface
print(np.random.normal(custom_obj))  # KO : np.random.normal will "only" try to 
cast the object to float
```

And here is a MWE : 
```python
import numpy as np
np.random.seed(1234)

HANDLED_FUNCTIONS = {}

class NumericalLabeled():
def __init__(self, value, label=""):
self.value = value
self.label = label

def __repr__(self):
return "NumericalLabelled<"+str(self.value) + "," + self.label+">"

def __array_function__(self, func, types, args, kwargs):
if func not in HANDLED_FUNCTIONS:
return NotImplemented
return HANDLED_FUNCTIONS[func](*args, **kwargs)


def make_numericallabelled(x, label=""):
"""
Helper function to cast anything into a NumericalLabelled object.
"""
if isinstance(x, NumericalLabeled):
return x
else:
return NumericalLabeled(x, label=label)

# Numpy functions
# Override functions - used with __array_function__
def implements(np_function):
def decorator(func):
HANDLED_FUNCTIONS[np_function] = func
return func
return decorator


@implements(np.random.normal)
def np_random_normal(loc=0.0, scale=1.0, **kwargs):
# cast both loc and scale into Numericallabelled
loc = make_numericallabelled(loc)
scale = make_numericallabelled(scale)
# check their label is "compatible"
if not loc.label == scale.label:
raise ValueError
return NumericalLabeled(np.random.rand(loc=loc.value,
   scale=scale.value, **kwargs), 
loc.label+scale.label)

@implements(np.mean)
def np_mean(a, *args, **kwargs):
return NumericalLabeled(np.mean(a.value, *args, **kwargs),
a.label)



def main():
# reference result for standard array
arr = np.arange(10)
print(np.mean(arr))
print(np.random.normal(arr))

# array-like object
num_labeled = NumericalLabeled(arr, "toto")
print(np.mean(num_labeled))
try:
print(np.random.normal(num_labeled))
except Exception as e:
print(e)

main()
```
which results in 
```
4.5
[ 0.47143516 -0.19097569  3.43270697  2.6873481   3.27941127  5.88716294
  6.85958841  6.3634765   8.01569637  6.75731505]
NumericalLabelled<4.5,toto> 
float() argument must be a string or a number, not 'NumericalLabeled'
```

Since the distribution functions accept array as input, I would expect them to 
be wrappable like many other array functions.

### Versions
```
Python : 3.8.5 (default, Sep  4 2020, 02:22:02) 
[Clang 10.0.0 ]
Numpy : 1.21.0
```


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

[Numpy-discussion] `keepdims=True` for argmin/argmx and C-API `PyArray_ArgMaxWithKeepdims`

[Sebastian Berg]

Hi all,

The PR https://github.com/numpy/numpy/pull/19211 proposes to extend
argmin and argmax with a `keepdims=False` keyword-only argument.

This is a standard argument in NumPy, so it is a small API addition.

The PR  also proposes to add:

* `PyArray_ArgMinWithKeepdims`
* `PyArray_ArgMaxWithKeepdims`

in the C-API.  We have barely extended the C-API in a very long time,
so if anyone has concerns, we could pull that out again [1].


Otherwise, this should go in soon, and we will have `keepdims` for both
of those functions in the next release :).

Cheers,

Sebastian



[1] I do not see this is much of a maintenance concern, since the
original function is just a one-line wrapper of the new one.

The API is fairly large and it probably is not used much. So it doesn't
feel important to add to me.  Overally, I just don't have a preference.

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