[Numpy-discussion] Re: Function that searches arrays for the first element that satisfies a condition

[Neal Becker]

Thanks Juan, this is really great!  I plan to make use of this right away.

On Wed, Nov 1, 2023 at 8:13 AM Juan Nunez-Iglesias  wrote:

> Have you tried timing things? Thankfully this is easy to test because the
> Python source of numba-jitted functions is available at jitted_func.py_func.
>
> In [23]: @numba.njit
> ...: def _first(arr, pred):
> ...: for i, elem in enumerate(arr):
> ...: if pred(elem):
> ...: return i
> ...:
> ...: def first(arr, pred):
> ...: _pred = numba.njit(pred)
> ...: return _first(arr, _pred)
> ...:
>
> In [24]: arr = np.random.random(100_000_000)
>
> In [25]: %timeit first(arr, lambda x: x > 5)
> 72 ms ± 1.36 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
>
> In [26]: %timeit arr + 5
> 90.3 ms ± 762 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
>
> In [27]: %timeit _first.py_func(arr, lambda x: x > 5)
> 7.8 s ± 46.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
>
> So numba gives a >100x speedup. It's still not as fast as a NumPy function
> call that doesn't have an allocation overhead:
>
> In [30]: arr2 = np.empty_like(arr, dtype=bool)
>
> In [32]: %timeit np.greater(arr, 5, out=arr2)
> 13.9 ms ± 69.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
>
> But it's certainly much better than pure Python! And it's not a huge cost
> for the flexibility.
>
> Juan.
>
> On Wed, 1 Nov 2023, at 10:42 AM, Dom Grigonis wrote:
>
> This results in a very slow code. The function calls of
>
> _pred = numba.njit(pred)
>
> are expensive and this sort of approach will be comparable to pure python
> functions.
>
> This is only recommended for sourcing functions that are not called
> frequently, but rather have a large computational content within them. In
> other words not suitable for predicates.
>
> Regards,
> DG
>
> On 1 Nov 2023, at 01:05, Juan Nunez-Iglesias  wrote:
>
> If you add a layer of indirection with Numba you can get a *very* nice API:
>
> @numba.njit
> def _first(arr, pred):
> for i, elem in enumerate(arr):
> if pred(elem):
> return i
>
> def first(arr, pred):
> _pred = numba.njit(pred)
> return _first(arr, _pred)
>
> This even works with lambdas! (TIL, thanks Numba devs!)
>
> >>> first(np.random.random(10_000_000), lambda x: x > 0.99)
> 215
>
> Since Numba has ufunc support I don't suppose it would be hard to make it
> work with an axis= argument, but I've never played with that API myself.
>
> On Tue, 31 Oct 2023, at 6:49 PM, Lev Maximov wrote:
>
> I've implemented such functions in Cython and packaged them into a library
> called numpy_illustrated 
>
> It exposes the following functions:
>
> find(a, v)  # returns the index of the first occurrence of v in a
> first_above(a, v)   # returns the index of the first element in a that is
> strictly above v
> first_nonzero(a)   # returns the index of the first nonzero element
>
> They scan the array and bail out immediately once the match is found. Have
> a significant performance gain if the element to be
> found is closer to the beginning of the array. Have roughly the same speed
> as alternative methods if the value is missing.
>
> The complete signatures of the functions look like this:
>
> find(a, v, rtol=1e-05, atol=1e-08, sorted=False, default=-1, raises=False)
> first_above(a, v, sorted=False, missing=-1, raises=False)
> first_nonzero(a, missing=-1, raises=False)
>
> This covers the most common use cases and does not accept Python callbacks
> because accepting them would nullify any speed gain
> one would expect from such a function. A Python callback can be
> implemented with Numba, but anyone who can write the callback
> in Numba has no need for a library that wraps it into a dedicated function.
>
> The library has a 100% test coverage. Code style 'black'. It should be
> easy to add functions like 'first_below' if necessary.
>
> A more detailed description of these functions can be found here
> 
> .
>
> Best regards,
>   Lev Maximov
>
> On Tue, Oct 31, 2023 at 3:50 AM Dom Grigonis 
> wrote:
>
> I juggled a bit and found pretty nice solution using numba. Which is
> probably not very robust, but proves that such thing can be optimised while
> retaining flexibility. Check if it works for your use cases and let me know
> if anything fails or if it is slow compared to what you used.
>
>
> first_true_str = """def first_true(arr, n):result = np.full((n, 
> arr.shape[1]), -1, dtype=np.int32)for j in range(arr.shape[1]):k 
> = 0for i in range(arr.shape[0]):x = arr[i:i + 1, j]   
>  if cond(x):result[k, j] = ik += 1
> if k >= n:breakreturn result"""
>
> *class* *FirstTrue*:
> CONTEXT = {'np': np}
>
> *def* __init__(self, expr):
>  

[Numpy-discussion] Re: Change definition of complex sign (and use it in copysign)

[Neal Becker]

In my opinion, with the caveat that anyone that asks for the sign of a
complex number gets what they deserve, this seems about as useful a
definition as any.

On Fri, Dec 22, 2023 at 8:23 AM  wrote:

> Hi All,
>
> A long-standing, small wart in numpy has been that the definition of sign
> for complex numbers is really useless (`np.sign(z)` gives the sign of the
> real component, unless that is zero, in which case it gives the sign of the
> imaginary component, in both cases as a complex number with zero imaginary
> part). Useless enough, in fact, that in the Array API it is suggested [1]
> that sign should return `z / |z|` (a definition consistent with those of
> reals, giving the direction in the complex plane).
>
> The question then becomes what to do. My suggestion - see
> https://github.com/numpy/numpy/pull/25441 - is to adapt the Array API
> definition for numpy 2.0, with the logic that if we don't change it in 2.0,
> when will we?
>
> Implementing it, I found no real test failures except one for
> `np.geomspace`, where it turned out that to correct the failure, the new
> definition substantially simplified the implementation.
>
> Furthermore, with the redefinition, it has become possible to extend
> ``np.copysign(x1, x2)`` to complex numbers, since it can now generally
> return ``|x1| * sign(x2)`` with the sign as defined above (with no special
> treatment for zero).
>
> Anyway, to me the main question would be whether this would break any
> workflows (though it is hard to see how it could, given that the previous
> definition was really rather useless...).
>
> Thanks,
> https://github.com/data-apis/array-api/pull/556,
> Marten
>
> [1]
> https://data-apis.org/array-api/latest/API_specification/generated/array_api.sign.html
> (and https://github.com/data-apis/array-api/pull/556, which has links to
> previous discussion)
> ___
> NumPy-Discussion mailing list -- numpy-discussion@python.org
> To unsubscribe send an email to numpy-discussion-le...@python.org
> https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
> Member address: ndbeck...@gmail.com
>


-- 
*Those who don't understand recursion are doomed to repeat it*
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] argmin of empty masked sequence not correctly handled (returns 0)

[Neal Becker]

argmin of an empty sequence throws an exception, which I think is a
good thing:
 np.argmin([])
---
ValueError

But not for a masked array:

n [24]: u = np.arange (10)

In [25]: v = np.ma.array (u, mask=np.ones(10))

In [26]: v
Out[26]:
masked_array(data=[--, --, --, --, --, --, --, --, --, --],
 mask=[ True,  True,  True,  True,  True,  True,  True,  True,
True,  True],
   fill_value=99,
dtype=int64)

In [28]: np.argmin(v)
Out[28]: 0

This is error prone and inconsistent.
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] Re: Original seed of a BitGenerator

[Neal Becker]

On Fri, Aug 2, 2024 at 9:37 AM Robert Kern  wrote:

> On Fri, Aug 2, 2024 at 1:28 AM Andrew Nelson  wrote:
>
>> When using the new `Generator`s for stochastic optimisation I sometimes
>> find myself possessing a great solution, but am wondering what path the
>> random number generation took to get to that point.
>>
>> I know that I can get the current state of the BitGenerators. However,
>> what I'd like to do is query the BitGenerator to figure out how the
>> BitGenerator was setup in the first place.
>>
>> i.e. either:
>>
>> - the seed/SeedSequence that was used to construct the BitGenerator
>>
>
> >>> rng = np.random.default_rng()
> >>> rng.bit_generator.seed_seq
> SeedSequence(
> entropy=186013007116029215180532390504704448637,
> )
>
> In some older versions of numpy, the attribute was semi-private as
> _seed_seq, if you're still using one of those.
>
> In many cases you can add your own attributes to python objects, so you
can record the seed yourself as
rng.my_seed = blah

I didn't test if the BitGenerator supports this

-- 
*Those who don't understand recursion are doomed to repeat it*
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] Set #threads from within python code

[Neal Becker]

I don't want my python code to run multi-thread.  So I can do:

MKL_NUM_THREAD=1 NUMEXPR_NUM_THREADS=1 OMP_NUM_THREADS=1 my_program...

But I don't seem to be able to achieve this effect without setting env
variables on the command line; within my_program.  Using os.environ doesn't
work.  I don't understand why.  I'd like to not have to put this on the
command line because I sometimes forget.
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Set #threads from within python code

[Neal Becker]

Ah, that probably explains it!

On Fri, Mar 24, 2017 at 7:37 AM Daπid  wrote:

> On 24 March 2017 at 12:30, Neal Becker  wrote:
> > Using os.environ doesn't
> > work.  I don't understand why.
>
> It should, I do that for other variables. Are you setting the
> variables before importing other libraries? They may only get read at
> import time.
> ___
> 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

Re: [Numpy-discussion] speed of random number generator compared to Julia

[Neal Becker]

Take a look here:
https://bashtage.github.io/ng-numpy-randomstate/doc/index.html

On Mon, Apr 3, 2017 at 9:45 AM Jaime Fernández del Río 
wrote:

> On Mon, Apr 3, 2017 at 3:20 PM, Pierre Haessig 
> wrote:
>
> Hello,
> Le 30/03/2017 à 13:31, Pierre Haessig a écrit :
>
> []
>
> But how come Julia is 4-5x faster since Numpy uses C implementation for
> the entire process ? (Mersenne Twister -> uniform double -> Box-Muller
> transform to get a Gaussian
> https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/randomkit.c).
> Also I noticed that Julia uses a different algorithm (Ziggurat Method
> from Marsaglia and Tsang ,
> https://github.com/JuliaLang/julia/blob/master/base/random.jl#L700) but
> this doesn't explain the difference for uniform rng.
>
> Any ideas?
>
>
> This
> 
>  says
> that Julia uses this library
> , which is
> different from the home brewed version of the Mersenne twister in NumPy.
> The second link I posted claims their speed comes from generating double
> precision numbers directly, rather than generating random bytes that have
> to be converted to doubles, as is the case of NumPy through this magical
> incantation
> .
> They also throw the SIMD acronym around, which likely means their random
> number generation is parallelized.
>
> My guess is that most of the speed-up comes from the SIMD parallelization:
> the Mersenne algorithm does a lot of work
> 
>  to
> produce 32 random bits, so that likely dominates over a couple of
> arithmetic operations, even if divisions are involved.
>
> Jaime
>
> Do you think Stackoverflow would be a better place for my question?
>
> best,
>
> Pierre
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion
>
>
>
>
> --
> (\__/)
> ( O.o)
> ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes
> de dominación mundial.
> ___
> 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

Re: [Numpy-discussion] speed of random number generator compared to Julia

[Neal Becker]

I think the intention is that this is the next gen of numpy randomstate,
and will eventually be merged in.

On Mon, Apr 3, 2017 at 11:47 AM Pierre Haessig 
wrote:

>
> Le 03/04/2017 à 15:52, Neal Becker a écrit :
> > Take a look here:
> > https://bashtage.github.io/ng-numpy-randomstate/doc/index.html
> Thanks for the pointer. A very feature-full random generator package.
>
> So it is indeed possible to have in Python/Numpy both the "advanced"
> Mersenne Twister (dSFMT) at the lower level and the Ziggurat algorithm
> for Gaussian transform on top. Perfect!
>
> In an ideal world, this would be implemented by default in Numpy, but I
> understand that this would break the reproducibility of existing codes.
>
> best,
> Pierre
> ___
> 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

Re: [Numpy-discussion] proposal: smaller representation of string arrays

[Neal Becker]

I'm no unicode expert, but can't we truncate unicode strings so that only
valid characters are included?

On Thu, Apr 20, 2017 at 1:32 PM Chris Barker  wrote:

> On Thu, Apr 20, 2017 at 9:47 AM, Anne Archibald  > wrote:
>
>> Is there any reason not to support all Unicode encodings that python
>> does, with the same names and semantics? This would surely be the simplest
>> to understand.
>>
>
> I think it should support all fixed-length encodings, but not the
> non-fixed length ones -- they just don't fit well into the numpy data model.
>
>
>> Also, if latin1 is to going to be the only practical 8-bit encoding,
>> maybe check with some non-Western users to make sure it's not going to
>> wreck their lives? I'd have selected ASCII as an encoding to treat
>> specially, if any, because Unicode already does that and the consequences
>> are familiar. (I'm used to writing and reading French without accents
>> because it's passed through ASCII, for example.)
>>
>
> latin-1 (or latin-9) only makes things better than ASCII -- it buys most
> of the accented characters for the European language and some symbols that
> are nice to have (I use the degree symbol a lot...). And it is ASCII
> compatible -- so there is NO reason to choose ASCII over Latin-*
>
> Which does no good for non-latin languages -- so we need to hear from the
> community -- is there a substantial demand for a non-latin one-byte per
> character encoding?
>
>
>> Variable-length encodings, of which UTF-8 is obviously the one that makes
>> good handling essential, are indeed more complicated. But is it strictly
>> necessary that string arrays hold fixed-length *strings*, or can the
>> encoding length be fixed instead? That is, currently if you try to assign a
>> longer string than will fit, the string is truncated to the number of
>> characters in the data type.
>>
>
> we could do that, yes, but an improperly truncated "string" becomes
> invalid -- just seems like a recipe for bugs that won't be found in testing.
>
> memory is cheap, compressing is fast -- we really shouldn't get hung up on
> this!
>
> Note: if you are storing a LOT of text (which I have no idea why you would
> use numpy anyway), then the memory size might matter, but then
> semi-arbitrary truncation would probably matter, too.
>
> I expect most text storage in numpy arrays is things like names of
> datasets, ids, etc, etc -- not massive amounts of text -- so storage space
> really isn't critical. but having an id or something unexpectedly truncated
> could be bad.
>
> I think practical experience has shown us that people do not handle
> "mostly fixed length but once in awhile not" text well -- see the nightmare
> of UTF-16 on Windows. Granted, utf-8 is multi-byte far more often, so
> errors are far more likely to be found in tests (why would you use utf-8 is
> all your data are in ascii???). but still -- why invite hard to test for
> errors?
>
> Final point -- as Julian suggests, one reason to support utf-8 is for
> interoperability with other systems -- but that makes errors more of an
> issue -- if it doesn't pass through the numpy truncation machinery, invalid
> data could easily get put in a numpy array.
>
> -CHB
>
>  it would allow UTF-8 to be used just the way it usually is - as an
>> encoding that's almost 8-bit.
>>
>
> ouch! that perception is the route to way too many errors! it is by no
> means almost 8-bit, unless your data are almost ascii -- in which case, use
> latin-1 for pity's sake!
>
> This highlights my point though -- if we support UTF-8, people WILL use
> it, and only test it with mostly-ascii text, and not find the bugs that
> will crop up later.
>
> All this said, it seems to me that the important use cases for string
>> arrays involve interaction with existing binary formats, so people who have
>> to deal with such data should have the final say. (My own closest approach
>> to this is the FITS format, which is restricted by the standard to ASCII.)
>>
>
> yup -- not sure we'll get much guidance here though -- netdf does not
> solve this problem well, either.
>
> But if you are pulling, say, a utf-8 encoded string out of a netcdf file
> -- it's probably better to pull it out as bytes and pass it through the
> python decoding/encoding machinery than pasting the bytes straight to a
> numpy array and hope that the encoding and truncation are correct.
>
> -CHB
>
>
> --
>
> Christopher Barker, Ph.D.
> Oceanographer
>
> Emergency Response Division
> NOAA/NOS/OR&R(206) 526-6959   voice
> 7600 Sand Point Way NE   (206) 526-6329   fax
> Seattle, WA  98115   (206) 526-6317   main reception
>
> chris.bar...@noaa.gov
> ___
> 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

Re: [Numpy-discussion] proposal: smaller representation of string arrays

[Neal Becker]

So while compression+ucs-4 might be OK for out-of-core representation, what
about in-core?  blosc+ucs-4?  I don't think that works for mmap, does it?

On Thu, Apr 27, 2017 at 7:11 AM Francesc Alted  wrote:

> 2017-04-27 3:34 GMT+02:00 Stephan Hoyer :
>
>> On Wed, Apr 26, 2017 at 4:49 PM, Nathaniel Smith  wrote:
>>
>>> It's worthwhile enough that both major HDF5 bindings don't support
>>> Unicode arrays, despite user requests for years. The sticking point seems
>>> to be the difference between HDF5's view of a Unicode string array (defined
>>> in size by the bytes of UTF-8 data) and numpy's current view of a Unicode
>>> string array (because of UCS-4, defined by the number of
>>> characters/codepoints/whatever). So there are HDF5 files out there that
>>> none of our HDF5 bindings can read, and it is impossible to write certain
>>> data efficiently.
>>>
>>>
>>> I would really like to hear more from the authors of these libraries
>>> about what exactly it is they feel they're missing. Is it that they want
>>> numpy to enforce the length limit early, to catch errors when the array is
>>> modified instead of when they go to write it to the file? Is it that they
>>> really want an O(1) way to look at a array and know the maximum number of
>>> bytes needed to represent it in utf-8? Is it that utf8<->utf-32 conversion
>>> is really annoying and files that need it are rare so they haven't had the
>>> motivation to implement it? My impression is similar to Julian's: you
>>> *could* implement HDF5 fixed-length utf-8 <-> numpy U arrays with a few
>>> dozen lines of code, which is nothing compared to all the other hoops these
>>> libraries are already jumping through, so if this is really the roadblock
>>> then I must be missing something.
>>>
>>
>> I actually agree with you. I think it's mostly a matter of convenience
>> that h5py matched up HDF5 dtypes with numpy dtypes:
>> fixed width ASCII -> np.string_/bytes
>> variable length ASCII -> object arrays of np.string_/bytes
>> variable length UTF-8 -> object arrays of unicode
>>
>> This was tenable in a Python 2 world, but on Python 3 it's broken and
>> there's not an easy fix.
>>
>> We absolutely could fix h5py by mapping everything to object arrays of
>> Python unicode strings, as has been discussed (
>> https://github.com/h5py/h5py/pull/871). For fixed width UTF-8, this
>> would be a fine but non-ideal solution, since there is currently no fixed
>> width UTF-8 support.
>>
>> For fixed width ASCII arrays, this would mean increased convenience for
>> Python 3 users, at the price of decreased convenience for Python 2 users
>> (arrays now contain boxed Python objects), unless we made the h5py behavior
>> dependent on the version of Python. Hence, we're back here, waiting for
>> better dtypes for encoded strings.
>>
>> So for HDF5, I see good use cases for ASCII-with-surrogateescape (for
>> handling ASCII arrays as strings) and UTF-8 with length equal to the number
>> of bytes.
>>
>
> Well, I'll say upfront that I have not read this discussion in the fully,
> but apparently some opinions from developers of HDF5 Python packages would
> be welcome here, so here I go :) ​
>
> As a long-time developer of one of the Python HDF5 packages (PyTables), I
> have always been of the opinion that plain ASCII (for byte strings) and
> UCS-4 (for Unicode) encoding would be the appropriate dtypes​ for storing
> large amounts of data, most specially for disk storage (but also using
> compressed in-memory containers).  My rational is that, although UCS-4 may
> require way too much space, compression would reduce that to basically the
> space that is required by compressed UTF-8 (I won't go into detail, but
> basically this is possible by using the shuffle filter).
>
> I remember advocating for UCS-4 adoption in the HDF5 library many years
> ago (2007?), but I had no success and UTF-8 was decided to be the best
> candidate.  So, the boat with HDF5 using UTF-8 sailed many years ago, and I
> don't think there is a go back (not even adding UCS-4 support on it,
> although I continue to think it would be a good idea).  So, I suppose that
> if HDF5 is found to be an important format for NumPy users (and I think
> this is the case), a solution for representing Unicode characters by using
> UTF-8 in NumPy would be desirable (at the risk of making the implementation
> more complex).
>
> ​Francesc
> ​
>
>>
>> ___
>> NumPy-Discussion mailing list
>> NumPy-Discussion@python.org
>> https://mail.python.org/mailman/listinfo/numpy-discussion
>>
>>
>
>
> --
> Francesc Alted
> ___
> 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 numpy arrays to Matlab?

[Neal Becker]

I've searched but haven't found any decent answer.  I need to call Matlab
from python.  Matlab has a python module for this purpose, but it doesn't
understand numpy AFAICT.  What solutions are there for efficiently
interfacing numpy arrays to Matlab?

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

Re: [Numpy-discussion] Interface numpy arrays to Matlab?

[Neal Becker]

Transplant sounds interesting, I think I could use this.  I don't
understand though why nobody has used a more direct approach?  Matlab has
their python API
https://www.mathworks.com/help/matlab/matlab-engine-for-python.html.  This
will pass Matlab arrays to/from python as some kind of opaque blob.  I
would guess that inside every Matlab array is a numpy array crying to be
freed - in both cases an array is a block of memory together with shape and
stride information.  So I would hope a direct conversion could be done, at
least via C API if not directly with python numpy API.  But it seems nobody
has done this, so maybe it's not that simple?


On Mon, Aug 28, 2017 at 5:32 PM Gregory Lee  wrote:

> I have not used Transplant, but it sounds fairly similar to
> Python-matlab-bridge.  We currently optionally call Matlab via
> Python-matlab-bridge in some of the the tests for the PyWavelets package.
>
> https://arokem.github.io/python-matlab-bridge/
> https://github.com/arokem/python-matlab-bridge
>
> I would be interested in hearing about the benefits/drawbacks relative to
> Transplant if there is anyone who has used both.
>
>
> On Mon, Aug 28, 2017 at 4:29 PM, CJ Carey 
> wrote:
>
>> Looks like Transplant can handle this use-case.
>>
>> Blog post: http://bastibe.de/2015-11-03-matlab-engine-performance.html
>> GitHub link: https://github.com/bastibe/transplant
>>
>> I haven't given it a try myself, but it looks promising.
>>
>> On Mon, Aug 28, 2017 at 4:21 PM, Stephan Hoyer  wrote:
>>
>>> If you can use Octave instead of Matlab, I've had a very good experience
>>> with Oct2Py:
>>> https://github.com/blink1073/oct2py
>>>
>>> On Mon, Aug 28, 2017 at 12:20 PM, Neal Becker 
>>> wrote:
>>>
>>>> I've searched but haven't found any decent answer.  I need to call
>>>> Matlab from python.  Matlab has a python module for this purpose, but it
>>>> doesn't understand numpy AFAICT.  What solutions are there for efficiently
>>>> interfacing numpy arrays to Matlab?
>>>>
>>>> Thanks,
>>>> Neal
>>>>
>>>> ___
>>>> 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 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 mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Moving NumPy's PRNG Forward

[Neal Becker]

On Fri, Jan 19, 2018 at 6:13 PM Nathaniel Smith  wrote:

> ...



> I agree that relaxing our policy would be better than the status quo.
> Before making any decisions, though, I'd like to make sure we
> understand the alternatives and their trade-offs. Specifically, I
> think the main alternative would be the following approach to
> versioning:
>
> 1) make RandomState's state be a tuple (underlying RNG algorithm,
> underlying RNG state, distribution version)
> 2) zero-argument initialization/seeding, like RandomState() or
> rstate.seed(), sets the state to: (our recommended RNG algorithm,
> os.urandom(...), version=LATEST_VERSION)
> 3) for backcompat, single-argument seeding like RandomState(123) or
> rstate.seed(123), sets the state to: (mersenne twister,
> expand_mt_seed(123), version=0)
> 4) also allow seeding to explicitly control all the parameters, like
> RandomState(PCG_XSL_RR(123), version=12) or whatever
> 5) the distribution functions are implemented like:
>
> def normal(*args, **kwargs):
> if self.version < 3:
> return self._normal_box_muller(*args, **kwargs)
> elif self.version < 8:
> return self._normal_ziggurat_v1(*args, **kwargs)
> else:  # version >= 8
> return self._normal_ziggurat_v2(*args, **kwargs)
>

I like this suggestion, but I suggest to modify it so that a zero-argument
initialization or 1-argument seeding will initialize to a global value,
which would  default to backcompat, but could be changed.  Then my old code
would by default produce the same old results, but adding 1 line at the top
switches to faster code if I want.
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Prototype CorePRNG implementation feedback

[Neal Becker]

What is the syntax to construct an initialized generator?
RandomGenerator(Xoroshiro128(my_seed))?

On Thu, Feb 22, 2018 at 6:06 AM Kevin Sheppard 
wrote:

> I have implemented a working prototype of a pluggable RandomState.  The
> repo is located at https://github.com/bashtage/core-prng.
>
> The main design has a small PRNG (currently only SplitMix64 and
> Xoroshiro128 are implemented) that does not have any public functions
> available to generate random numbers.  These objects only expose methods to
> get and set state and to jump the PRNG.
>
> The random number generator is exposed via an opaque structure wrapped in
> a PyCapsule.  This structure has two void pointers, one to the state and
> one to the function that generates the next uint 64.  I think in the
> long-run it might make sense to expose a few additional interfaces, such as
> the next random double (makes sense for dSFMT) or the next random uint32
> (makes sense for MT19937).  For now, I have deliberately kept it simple.
>
> The user-facing object is called `RandomGenerator` and it can be
> initialized using the syntax `RandomGenerator` or
> `RandomGenerator(Xoroshiro128())`.  The object exposes user-facing methods
> (currently only `random_integer` and `random_double`).
>
> A basic demo:
>
> from core_prng.generator import RandomGenerator
> # Splitmix 64 for now
> RandomGenerator().random_integer()
>
> from core_prng.xoroshiro128 import Xoroshiro128
> RandomGenerator(Xoroshiro128()).random_integer()
>
> A few questions that have come up:
>
>1. Should a passed core PRNG be initialized or not. The syntax 
> RandomGenerator(Xoroshiro128)
>looks nicer than RandomGenerator(Xoroshiro128())
>2. Which low-level methods should be part of the core PRNG?  These
>would all be scalar generators that would be consumed by RandomGenerator
>and exposed through void *. I could imagine some subset of the
>following:
>   - uint64
>   - uint32
>   - double
>   - float
>   - uint53 (applicable to PRNGs that use doubles as default type)
>3. Since RandomState is taken, a new name is needed for the
>user-facing object.  RandomGenerator is a placeholder but ideally,
>something meaningful could be chosen.
>4. I can't see how locking can be usefully implemented in the core
>PRNGs without a large penalty. This means that these are not threadsafe.
>This isn't a massive problem but they do access the state (get or set)
>although with GIL.
>5. Should state be a property -- this isn't Java...
>
>
> Other feedback is of course also welcome.
>
> Kevin
>
>
> ___
> 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] update to numpy-1.5.0 gives new warnings from scipy

[Neal Becker]

After update to numpy-1.5.0, I'm getting warnings from scipy.
These probably come from my code using convolve.  Does scipy need updating?

/home/nbecker/.local/lib/python3.6/site-packages/scipy/fftpack/basic.py:160:
FutureWarning: Using a non-tuple sequence for multidimensional indexing is
deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this
will be interpreted as an array index, `arr[np.array(seq)]`, which will
result either in an error or a different result.
  z[index] = x
/home/nbecker/.local/lib/python3.6/site-packages/scipy/signal/signaltools.py:491:
FutureWarning: Using a non-tuple sequence for multidimensional indexing is
deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this
will be interpreted as an array index, `arr[np.array(seq)]`, which will
result either in an error or a different result.
  return x[reverse].conj()
/home/nbecker/.local/lib/python3.6/site-packages/scipy/signal/signaltools.py:251:
FutureWarning: Using a non-tuple sequence for multidimensional indexing is
deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this
will be interpreted as an array index, `arr[np.array(seq)]`, which will
result either in an error or a different result.
  in1zpadded[sc] = in1.copy()
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Add pybind11 to docs about writing binding code

[Neal Becker]

The only way I've seen xtensor used is *with* pybind11

On Fri, Aug 17, 2018 at 5:45 PM Robert Kern  wrote:

> On Fri, Aug 17, 2018 at 2:00 AM Hans Dembinski 
> wrote:
>
>> Hi Sylvain,
>>
>> On 16. Aug 2018, at 13:29, Sylvain Corlay 
>> wrote:
>>
>> Actually, xtensor-python does a lot more in terms of numpy bindings, as
>> it uses the C APIs of numpy directly for a number of things.
>>
>> Plus, the integration into the xtensor expression system lets you do
>> things such as view / broadcasting / newaxis / ufuncs directly from the C++
>> side (and all that is in the cheat sheets).
>>
>> ok, good, but my point was different. The page in question is about
>> Python as a glue language. The other solutions on that site are general
>> purpose binding solutions for any kind of C++ code, while xtensor-python is
>> xtensor-specific. xtensor in turn is a library that mimics the numpy API in
>> C++.
>>
>
> Even if you don't use the numpy-mimicking parts of the xtensor API,
> xtensor-python is a probably a net improvement over pybind11 for
> communicating arrays back and forth across the C++/Python boundary. Even if
> the rest of your C++ code doesn't use xtensor, you could profitably use
> xtensor-python at the interface. Also, though the article is generally
> framed as using Python as a glue language (i.e. communicating with existing
> C/C++/Fortran code), it is also relevant for the use case where you are
> writing the C/C++/Fortran code from scratch (perhaps just accelerating
> small kernels or whatever). Talking about the available options for that
> use case is perfectly on-topic for that article.
>
> You don't have to be the one that writes it, though, if you just want to
> cover pybind11.
>
> --
> 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

Re: [Numpy-discussion] Add pybind11 to docs about writing binding code

[Neal Becker]

I'm confused, do you have a link or example showing how to use
xtensor-python without pybind11?
Thanks,
Neal

On Mon, Aug 20, 2018 at 5:50 PM Hans Dembinski 
wrote:

> Dear Robert,
>
> > On 17. Aug 2018, at 23:44, Robert Kern  wrote:
> >
> > Even if you don't use the numpy-mimicking parts of the xtensor API,
> xtensor-python is a probably a net improvement over pybind11 for
> communicating arrays back and forth across the C++/Python boundary. Even if
> the rest of your C++ code doesn't use xtensor, you could profitably use
> xtensor-python at the interface. Also, though the article is generally
> framed as using Python as a glue language (i.e. communicating with existing
> C/C++/Fortran code), it is also relevant for the use case where you are
> writing the C/C++/Fortran code from scratch (perhaps just accelerating
> small kernels or whatever). Talking about the available options for that
> use case is perfectly on-topic for that article.
>
> no objections here, xtensor should be highlighted in the pybind11 part for
> these reasons. I just think it should not be a separate section.
>
> Best regards,
> Hans
> ___
> 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

Re: [Numpy-discussion] Ia `-fno-strict-aliasing` still required?

[Neal Becker]

Define "doing fine", does it pass all tests?

On Fri, Jan 4, 2019 at 12:36 PM Charles R Harris 
wrote:

> Hi All,
>
> Just asking if the `-fno-strict-aliasing` flag is still required for gcc.
> Supposedly `-fstrict-aliasing` is enabled by default with optimization
> levels >= `-O2` and that used to result in errors. However, I have noticed
> that numpy wheels are being built without the `-fno-strict-aliasing` flag
> and doing fine. Same on my local machine. So maybe the compiler has gotten
> better at knowing when pointers may be aliased? I cannot find informative
> documentation on how `-fstrict-aliasing` is actually implemented, so
> thought I'd ask here if someone knows more about it.
>
> Chuck
> ___
> 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

Re: [Numpy-discussion] Add guaranteed no-copy to array creation and reshape?

[Neal Becker]

constants are easier to support for autocompletion than strings.  My
current env (emacs) will (usually) autocomplete the former, but not the
latter.

On Thu, Jan 10, 2019 at 2:21 PM Feng Yu  wrote:

> Hi Todd,
>
> I agree a flag is more suitable than classes.
>
> I would add another bonus of a flag than a function argument is to avoid
> massive contamination of function signatures for a global variation of
> behavior that affects many functions.
>
> Yu
>
> On Wed, Jan 9, 2019 at 11:34 PM Todd  wrote:
>
>> On Mon, Jan 7, 2019, 14:22 Feng Yu >
>>> Hi,
>>>
>>> Was it ever brought up the possibility of a new array class (ndrefonly,
>>> ndview) that is strictly no copy?
>>>
>>> All operations on ndrefonly will return ndrefonly and if the operation
>>> cannot be completed without making a copy, it shall throw an error.
>>>
>>> On the implementation there are two choices if we use subclasses:
>>>
>>> - ndrefonly can be a subclass of ndarray. The pattern would be subclass
>>> limiting functionality of super, but ndrefonly is a ndarray.
>>> - ndarray as a subclass of ndarray. Subclass supplements functionality
>>> of super. : ndarray will not throw an error when a copy is necessary.
>>> However ndarray is not a ndarray.
>>>
>>> If we want to be wild they do not even need to be subclasses of each
>>> other, or maybe they shall both be subclasses of something more
>>> fundamental.
>>>
>>> - Yu
>>>
>>
>> I would prefer a flag for this.  Someone can make an array read-only by
>> setting `arr.flags.writable=False`.  So along those lines, we could have a
>> `arr.flags.copyable` flag that if set to `False` would result in an error
>> of any operation tried to copy the data.
>>
>>> ___
>> 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 mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] overhauling numpy.random and randomgen Message-ID:

[Neal Becker]

The boost_random c++ library uses the terms 'generators' and
'distributions'.  Distributions are applied to generators.

On Fri, Apr 19, 2019 at 7:54 AM Kevin Sheppard
 wrote:
>
> >  Rather than "base RNG", what about calling these classes a "random source"
> or "random stream"? In particular, I would suggest defining two Python
> classes:
> > - np.random.Generator as a less redundant name for what is currently called
> RandomGenerator
> > - np.random.Source or np.random.Stream as an abstract base class for what
> are currently called "base RNGs"
>
> Naming is definitely hard.  Simple RNGs are currently called basic RNGs which 
> was inspired by mkl-random. 'source' sounds OK to me, but sort of hides the 
> fact that these are the actual Psuedo RNGs. `stream` has a technical meaning 
> (a single PRNG make produce multiple independent streams) and IMO should be 
> avoided since this might lead to confusion.  Perhaps source_rng (or in docs 
> Source RNG)?
>
> RandomGenerator is actually RandomTransformer, but I didn't like the latter 
> name.
>
> > There are also a couple of convenience attributes in the user-facing API
> that I would suggest refining:
> >   - The "brng" attribute of RandomGenerator is not a very descriptive name. 
> > I
> would prefer "stream" or "source", or the more explicit "base_rng" if we
> stick with that term.
> >   - I don't think we need the "generator" property on base RNG objects. It 
> > is
> fine to require writing np.random.Generator(base) instead. Looking at the
> implementation, .generator caches the RandomGenerator objects it creates on
> the base RNG, which creates a reference cycle. Yes, Python can garbage
> collect reference cycles, but this is still a muddled data model.
>
> The attribute name should match the final (descriptive) name, whatever it is. 
>  In RandomGen I am using the `basic_rng` attribute name, but this could be 
> `source`.  I also use a property so that the attribute can have a docstring 
> attached for use in IPython. I think this is more user-friendly.
>
> I think dropping the `generator` property on the basic RNGs is reasonable.  
> It was a convenience but is awkward, and I always understood that it creates 
> a cycle.
>
> > Finally, why do we expose the np.random.gen object? I thought part of the
> idea with the new API was to avoid global mutable state.
>
> Module level functions are essential for quick experiments and should be 
> provided.  The only difference here is that the singleton `seed`  and `state` 
> are no longer exposed so that it isn't possible (using the exposed API) to 
> set the seed.
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion



-- 
Those who don't understand recursion are doomed to repeat it
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Moving forward with value based casting

[Neal Becker]

On Wed, Jun 5, 2019 at 4:42 PM Sebastian Berg
 wrote:

I think the best approach is that if the user gave unambiguous types
as inputs to operators then the output should be the same dtype, or
type corresponding to the common promotion type of the inputs.

If the input type is not specified, I agree with the suggestion here:

> 3. Just associate python float with float64 and python integers with
>long/int64 and force users to always type them explicitly if they
>need to.

Explicit is better than implicit

-- 
Those who don't understand recursion are doomed to repeat it
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] float16 but no complex32?

[Neal Becker]

I see a float16 dtype but not complex32.  Is this an oversight?

Thanks,
Neal

-- 
Those who don't understand recursion are doomed to repeat it
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Low-level API for Random

[Neal Becker]

I have used C-api in the past, and would like to see a convenient and
stable way to do this.  Currently I'm using randomgen, but calling
(from c++)
to the python api.  The inefficiency is amortized by generating and
caching batches of results.

I thought randomgen was supposed to be the future of numpy random, so
I've based on that.

On Fri, Sep 20, 2019 at 6:08 AM Ralf Gommers  wrote:
>
>
>
> On Fri, Sep 20, 2019 at 5:29 AM Robert Kern  wrote:
>>
>> On Thu, Sep 19, 2019 at 11:04 PM Ralf Gommers  wrote:
>>>
>>>
>>>
>>> On Thu, Sep 19, 2019 at 4:53 PM Robert Kern  wrote:

 On Thu, Sep 19, 2019 at 5:24 AM Ralf Gommers  
 wrote:
>
>
> On Thu, Sep 19, 2019 at 10:28 AM Kevin Sheppard 
>  wrote:
>>
>> There are some users of the NumPy C code in randomkit.  This was never 
>> officially supported.  There has been a long open issue to provide this 
>> officially.
>>
>> When I wrote randomgen I supplied .pdx files that make it simpler to 
>> write Cython code that uses the components.  The lower-level API has not 
>> had much scrutiny and is in need of a clean-up.   I thought this would 
>> also encourage users to extend the random machinery themselves as part 
>> of their project or code so as to minimize the requests for new (exotic) 
>> distributions to be included in Generator.
>>
>> Most of the generator functions follow a pattern random_DISTRIBUTION.  
>> Some have a bit more name mangling which can easily be cleaned up (like 
>> ranomd_gauss_zig, which should become PREFIX_standard_normal).
>>
>> Ralf Gommers suggested unprefixed names.
>
>
> I suggested that the names should match the Python API, which I think 
> isn't quite the same. The Python API doesn't contain things like "gamma", 
> "t" or "f".


 As the implementations evolve, they aren't going to match one-to-one 100%. 
 The implementations are shared by the legacy RandomState. When we update 
 an algorithm, we'll need to make a new function with the better algorithm 
 for Generator to use, then we'll have two C functions roughly 
 corresponding to the same method name (albeit on different classes). C 
 doesn't give us as many namespace options as Python. We could rely on 
 conventional prefixes to distinguish between the two classes of function 
 (e.g. legacy_normal vs random_normal).
>>>
>>>
>>> That seems simple and clear
>>>
 There are times when it would be nice to be more descriptive about the 
 algorithm difference (e.g. random_normal_polar vs random_normal_ziggurat),
>>>
>>>
>>> We decided against versioning algorithms in NEP 19, so an update to an 
>>> algorithm would mean we'd want to get rid of the older version (unless it's 
>>> still in use by legacy). So AFAICT we'd never have both random_normal_polar 
>>> and random_normal_ziggurat present at the same time?
>>
>>
>> Well, we must because one's used by the legacy RandomState and one's used by 
>> Generator. :-)
>>
>>>
>>> I may be missing your point here, but if we have in Python 
>>> `Generator.normal` and can switch its implementation from polar to ziggurat 
>>> or vice versa without any deprecation, then why would we want to switch 
>>> names in the C API?
>>
>>
>> I didn't mean to suggest that we'd have an unbounded number of functions as 
>> we improve the algorithms, just that we might have 2 once we decide to 
>> change something about the algorithm. We need 2 to support both the improved 
>> algorithm in Generator and the legacy algorithm in RandomState. The current 
>> implementation of the C function would be copied to a new name (`legacy_foo` 
>> or whatever), then we'd make RandomState use that frozen copy, then we make 
>> the desired modifications to the main function that Generator is referencing 
>> (`random_foo`).
>>
>> Or we could just make those legacy copies now so that people get to use them 
>> explicitly under the legacy names, whatever they are, and we can feel more 
>> free to modify the main implementations. I suggested this earlier, but 
>> convinced myself that it wasn't strictly necessary. But then I admit I was 
>> more focused on the Python API stability than any promises about the 
>> C/Cython API.
>>
>> We might end up with more than 2 implementations if we need to change 
>> something about the function signature, for whatever reason, and we want to 
>> retain C/Cython API compatibility with older code. The C functions aren't 
>> necessarily going to be one-to-one to the Generator methods. They're just 
>> part of the implementation. So for example, if we wanted to, say, precompute 
>> some intermediate values from the given scalar parameters so we don't have 
>> to recompute them for each element of the `size`-large requested output, we 
>> might do that in one C function and pass those intermediate values as 
>> arguments to the C function that does the actual sampling. So we'd have two 
>

Re: [Numpy-discussion] Low-level API for Random

[Neal Becker]

I'm using the low-level generator.  In this example I need to generate
small random integers of defined bit widths (e.g., 2 bit).  So
I get 64-bit uniform random uintegers, and cache the values, returning
them n-bits (e.g. 2 bits) at a time to the caller.

On Fri, Sep 20, 2019 at 9:03 AM Matti Picus  wrote:
>
> On 20/9/19 2:18 pm, Neal Becker wrote:
> > I have used C-api in the past, and would like to see a convenient and
> > stable way to do this.  Currently I'm using randomgen, but calling
> > (from c++)
> > to the python api.  The inefficiency is amortized by generating and
> > caching batches of results.
> >
> > I thought randomgen was supposed to be the future of numpy random, so
> > I've based on that.
> >
>
> It would be good to have actual users tell us what APIs they need.
>
> Are you using the BitGenerators or only the higher level Generator
> functions?
>
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion



-- 
Those who don't understand recursion are doomed to repeat it
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] argmax() indexes to value

[Neal Becker]

max(axis=1)?

On Wed, Oct 30, 2019, 7:33 PM Daniele Nicolodi  wrote:

> Hello,
>
> this is a very basic question, but I cannot find a satisfying answer.
> Assume a is a 2D array and that I get the index of the maximum value
> along the second dimension:
>
> i = a.argmax(axis=1)
>
> Is there a better way to get the value of the maximum array entries
> along the second axis other than:
>
> v = a[np.arange(len(a)), i]
>
> ??
>
> Thank you.
>
> Cheers,
> Daniele
> ___
> 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

Re: [Numpy-discussion] manylinux upgrade for numpy wheels

[Neal Becker]

Slightly off topic perhaps, it is recommended to perform custom compilation
for best performance, yet is there an
easy way to do this?  I don't think a simple pip will do.

On Wed, Feb 5, 2020 at 4:07 AM Matthew Brett 
wrote:

> Hi,
>
> On Tue, Feb 4, 2020 at 10:38 PM Nathaniel Smith  wrote:
> >
> > Pretty sure the 2010 and 2014 images both have much newer compilers than
> that.
> >
> > There are still a lot of users on CentOS 6, so I'd still stick to 2010
> for now on x86_64 at least. We could potentially start adding 2014 wheels
> for the other platforms where we currently don't ship wheels – gotta be
> better than nothing, right?
> >
> > There probably still is some tail of end users whose pip is too old to
> know about 2010 wheels. I don't know how big that tail is. If we wanted to
> be really careful, we could ship both manylinux1 and manylinux2010 wheels
> for a bit – pip will automatically pick the latest one it recognizes – and
> see what the download numbers look like.
>
> That all sounds right to me too.
>
> Cheers,
>
> Matthew
>
> > On Tue, Feb 4, 2020, 13:18 Charles R Harris 
> wrote:
> >>
> >> Hi All,
> >>
> >> Thought now would be a good time to decide on upgrading manylinux for
> the 1.19 release so that we can make sure that everything works as
> expected. The choices are
> >>
> >> manylinux1 -- CentOS 5, currently used, gcc 4.2 (in practice 4.5), only
> supports i686, x86_64.
> >> manylinux2010 -- CentOS 6, gcc 4.5, only supports i686, x86_64.
> >> manylinux2014 -- CentOS 7, gcc 4.8, supports many more architectures.
> >>
> >> The main advantage of manylinux2014 is that it supports many new
> architectures, some of which we are already testing against. The main
> disadvantage is that it requires pip >= 19.x, which may not be much of a
> problem 4 months from now but will undoubtedly cause some installation
> problems. Unfortunately, the compiler remains archaic, but folks interested
> in performance should be using a performance oriented distribution or
> compiling for their native architecture.
> >>
> >> Chuck
> >>
> >> ___
> >> 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 mailing list
> NumPy-Discussion@python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion
>


-- 
*Those who don't understand recursion are doomed to repeat it*
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Tensor Developer Summit

[Neal Becker]

Stefan van der Walt wrote:

> Hi all,
> 
> This has been mentioned on the community calls, but not on the mailing
> list, so a reminder about the Tensor Developer Summit happening at March
> in Berkeley:
> 
> https://xd-con.org/tensor-2020/
> 
> We would love to have developers and advanced users of NumPy (or other
> array libraries with Python interfaces) attend.  Registration closes 20
> February.
> 
> Best regards,
> Stéfan
Sounds like an exciting group!  Will it be streamed?

Thanks,
Neal

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

[Numpy-discussion] reseed random generator (1.19)

[Neal Becker]

Consider the following:

from numpy.random import default_rng
rs = default_rng()

Now how do I re-seed the generator?
I thought perhaps rs.bit_generator.seed(), but there is no such attribute.

Thanks,
Neal

-- 
*Those who don't understand recursion are doomed to repeat it*
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Numpy FFT normalization options issue (addition of new option)

[Neal Becker]

Honestly, I don't find "forward" very informative.  There isn't any real
convention on whether FFT of IFFT have any normalization.
To the best of my experience, either forward or inverse could be normalized
by 1/N, or each normalized by 1/sqrt(N), or neither
could be normalized.  I will say my expertise is in signal processing and
communications.

Perhaps
norm = {full, half, none} would be clearest to me.

Thanks,
Neal

On Sat, Jun 27, 2020 at 10:40 AM Sebastian Berg 
wrote:

> On Fri, 2020-06-26 at 21:53 -0700, leofang wrote:
> > Hi all,
> >
> >
> > Since I brought this issue from CuPy to Numpy, I'd like to see a
> > decision
> > made sooner than later so that downstream libraries like SciPy and
> > CuPy can
> > act accordingly. I think norm='forward' is fine. If there're still
> > people
> > unhappy with it after my reply, I'd suggest norm='reverse'. It has
> > the same
> > meaning, but is less confusing (than 'inverse' or other choices on
> > the
> > table) to me.
> >
>
> I expect "forward" is good (if I misread something please correct me),
> and I think we can go ahead with it, sorry for the delay.  However, I
> have send an email to scipy-dev, since we should give them at least a
> heads-up, and if you do not mind, I would wait a few days to actually
> merge (although we can also simply reverse, as long as CuPy does not
> have a release with it).
>
> It might be nice to expand the kwarg docs slightly with a sentence for
> each normalization mode?  Refering to `np.fft` docs is good, but if we
> can squeeze in a short refresher and refer there for details/formula it
> would be nicer.
> I feel "forward" is very intuitive, but only after pointing out that it
> is related to whether the fft or ifft has the normalization factor.
>
> Cheers,
>
> Sebastian
>
>
> >
> > Best,
> > Leo
> >
> >
> >
> > --
> > Sent from: http://numpy-discussion.10968.n7.nabble.com/
> > ___
> > 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
>


-- 
*Those who don't understand recursion are doomed to repeat it*
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] reseed random generator (1.19)

[Neal Becker]

I was using this to reset the generator, in order to repeat the same
sequence again for testing purposes.

On Wed, Jun 24, 2020 at 6:40 PM Robert Kern  wrote:

> On Wed, Jun 24, 2020 at 3:31 PM Neal Becker  wrote:
>
>> Consider the following:
>>
>> from numpy.random import default_rng
>> rs = default_rng()
>>
>> Now how do I re-seed the generator?
>> I thought perhaps rs.bit_generator.seed(), but there is no such attribute.
>>
>
> In general, reseeding an existing generator instance is not a good
> practice. What effect are you trying to accomplish? I assume that you are
> asking this because you are currently using `RandomState.seed()`. In what
> circumstances?
>
> The raw `bit_generator.state` property *can* be assigned to, in order to
> support some advanced use cases (mostly involving de/serialization and
> similar kinds of meta-programming tasks). It's also been helpful for me to
> construct worst-case scenarios for testing parallel streams. But it quite
> deliberately bypasses the notion of deriving the state from a
> human-friendly seed number.
>
> --
> Robert Kern
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion
>


-- 
*Those who don't understand recursion are doomed to repeat it*
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] reseed random generator (1.19)

[Neal Becker]

On Sat, Jul 4, 2020 at 1:56 PM Robert Kern  wrote:


>
> 3. Is there a way of telling the number of draws a generator did?
>>
>> The use case is to checkpoint the number of draws and `.advance` the
>> bit generator when resuming from the checkpoint. (The runs are longer
>> then the batch queue limits).
>>
>
> There are computations you can do on the internal state of PCG64 and
> Philox to get this information, but not in general, no. I do recommend
> serializing the Generator or BitGenerator (or at least the BitGenerator's
> .state property, which is a nice JSONable dict for PCG64) for checkpointing
> purposes. Among other things, there is a cached uint32 for when odd numbers
> of uint32s are drawn that you might need to handle. The state of the
> default PCG64 is much smaller than MT19937. It's less work and more
> reliable than computing that distance and storing the original seed and the
> distance.
>
> --
> Robert Kern
>

Sorry, you lost me here.  If I want to save, restore the state of a
generator, can I use pickle/unpickle?


-- 
*Those who don't understand recursion are doomed to repeat it*
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] C-coded dot 1000x faster than numpy?

[Neal Becker]

I have code that performs dot product of a 2D matrix of size (on the
order of) [1000,16] with a vector of size [1000].  The matrix is
float64 and the vector is complex128.  I was using numpy.dot but it
turned out to be a bottleneck.

So I coded dot2x1 in c++ (using xtensor-python just for the
interface).  No fancy simd was used, unless g++ did it on it's own.

On a simple benchmark using timeit I find my hand-coded routine is on
the order of 1000x faster than numpy?  Here is the test code:
My custom c++ code is dot2x1.  I'm not copying it here because it has
some dependencies.  Any idea what is going on?

import numpy as np

from dot2x1 import dot2x1

a = np.ones ((1000,16))
b = np.array([ 0.80311816+0.80311816j,  0.80311816-0.80311816j,
   -0.80311816+0.80311816j, -0.80311816-0.80311816j,
1.09707981+0.29396165j,  1.09707981-0.29396165j,
   -1.09707981+0.29396165j, -1.09707981-0.29396165j,
0.29396165+1.09707981j,  0.29396165-1.09707981j,
   -0.29396165+1.09707981j, -0.29396165-1.09707981j,
0.25495815+0.25495815j,  0.25495815-0.25495815j,
   -0.25495815+0.25495815j, -0.25495815-0.25495815j])

def F1():
d = dot2x1 (a, b)

def F2():
d = np.dot (a, b)

from timeit import timeit
print (timeit ('F1()', globals=globals(), number=1000))
print (timeit ('F2()', globals=globals(), number=1000))

In [13]: 0.013910860987380147 << 1st timeit
28.608758996007964  << 2nd timeit
-- 
Those who don't understand recursion are doomed to repeat it
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] C-coded dot 1000x faster than numpy?

[Neal Becker]

One suspect is that it seems the numpy version was multi-threading.
This isn't useful here, because I'm running parallel monte-carlo
simulations using all cores.  Perhaps this is perversely slowing
things down?  I don't know how to account for 1000x  slowdown though.

On Tue, Feb 23, 2021 at 1:40 PM Roman Yurchak  wrote:
>
> For the first benchmark apparently A.dot(B) with A real and B complex is
> a known issue performance wise https://github.com/numpy/numpy/issues/10468
>
> In general, it might be worth trying different BLAS backends. For
> instance, if you install numpy from conda-forge you should be able to
> switch between OpenBLAS, MKL and BLIS:
> https://conda-forge.org/docs/maintainer/knowledge_base.html#switching-blas-implementation
>
> Roman
>
> On 23/02/2021 19:32, Andrea Gavana wrote:
> > Hi,
> >
> > On Tue, 23 Feb 2021 at 19.11, Neal Becker  > <mailto:ndbeck...@gmail.com>> wrote:
> >
> > I have code that performs dot product of a 2D matrix of size (on the
> > order of) [1000,16] with a vector of size [1000].  The matrix is
> > float64 and the vector is complex128.  I was using numpy.dot but it
> > turned out to be a bottleneck.
> >
> > So I coded dot2x1 in c++ (using xtensor-python just for the
> > interface).  No fancy simd was used, unless g++ did it on it's own.
> >
> > On a simple benchmark using timeit I find my hand-coded routine is on
> > the order of 1000x faster than numpy?  Here is the test code:
> > My custom c++ code is dot2x1.  I'm not copying it here because it has
> > some dependencies.  Any idea what is going on?
> >
> >
> >
> > I had a similar experience - albeit with an older numpy and Python 2.7,
> > so my comments are easily outdated and irrelevant. This was on Windows
> > 10 64 bit, way more than plenty RAM.
> >
> > It took me forever to find out that numpy.dot was the culprit, and I
> > ended up using fortran + f2py. Even with the overhead of having to go
> > through f2py bridge, the fortran dot_product was several times faster.
> >
> > Sorry if It doesn’t help much.
> >
> > Andrea.
> >
> >
> >
> >
> > import numpy as np
> >
> > from dot2x1 import dot2x1
> >
> > a = np.ones ((1000,16))
> > b = np.array([ 0.80311816+0.80311816j,  0.80311816-0.80311816j,
> > -0.80311816+0.80311816j, -0.80311816-0.80311816j,
> >  1.09707981+0.29396165j,  1.09707981-0.29396165j,
> > -1.09707981+0.29396165j, -1.09707981-0.29396165j,
> >  0.29396165+1.09707981j,  0.29396165-1.09707981j,
> > -0.29396165+1.09707981j, -0.29396165-1.09707981j,
> >  0.25495815+0.25495815j,  0.25495815-0.25495815j,
> > -0.25495815+0.25495815j, -0.25495815-0.25495815j])
> >
> > def F1():
> >  d = dot2x1 (a, b)
> >
> > def F2():
> >  d = np.dot (a, b)
> >
> > from timeit import timeit
> > print (timeit ('F1()', globals=globals(), number=1000))
> > print (timeit ('F2()', globals=globals(), number=1000))
> >
> > In [13]: 0.013910860987380147 << 1st timeit
> > 28.608758996007964  << 2nd timeit
> > --
> > Those who don't understand recursion are doomed to repeat it
> > ___
> > NumPy-Discussion mailing list
> > NumPy-Discussion@python.org <mailto:NumPy-Discussion@python.org>
> > https://mail.python.org/mailman/listinfo/numpy-discussion
> > <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 mailing list
> NumPy-Discussion@python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion



-- 
Those who don't understand recursion are doomed to repeat it
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] C-coded dot 1000x faster than numpy?

[Neal Becker]

I'm using fedora 33 standard numpy.
ldd says:

/usr/lib64/python3.9/site-packages/numpy/core/_multiarray_umath.cpython-39-x86_64-linux-gnu.so:
linux-vdso.so.1 (0x7ffdd1487000)
libflexiblas.so.3 => /lib64/libflexiblas.so.3 (0x7f0512787000)

So whatever flexiblas is doing controls blas.

On Tue, Feb 23, 2021 at 1:51 PM Neal Becker  wrote:
>
> One suspect is that it seems the numpy version was multi-threading.
> This isn't useful here, because I'm running parallel monte-carlo
> simulations using all cores.  Perhaps this is perversely slowing
> things down?  I don't know how to account for 1000x  slowdown though.
>
> On Tue, Feb 23, 2021 at 1:40 PM Roman Yurchak  wrote:
> >
> > For the first benchmark apparently A.dot(B) with A real and B complex is
> > a known issue performance wise https://github.com/numpy/numpy/issues/10468
> >
> > In general, it might be worth trying different BLAS backends. For
> > instance, if you install numpy from conda-forge you should be able to
> > switch between OpenBLAS, MKL and BLIS:
> > https://conda-forge.org/docs/maintainer/knowledge_base.html#switching-blas-implementation
> >
> > Roman
> >
> > On 23/02/2021 19:32, Andrea Gavana wrote:
> > > Hi,
> > >
> > > On Tue, 23 Feb 2021 at 19.11, Neal Becker  > > <mailto:ndbeck...@gmail.com>> wrote:
> > >
> > > I have code that performs dot product of a 2D matrix of size (on the
> > > order of) [1000,16] with a vector of size [1000].  The matrix is
> > > float64 and the vector is complex128.  I was using numpy.dot but it
> > > turned out to be a bottleneck.
> > >
> > > So I coded dot2x1 in c++ (using xtensor-python just for the
> > > interface).  No fancy simd was used, unless g++ did it on it's own.
> > >
> > > On a simple benchmark using timeit I find my hand-coded routine is on
> > > the order of 1000x faster than numpy?  Here is the test code:
> > > My custom c++ code is dot2x1.  I'm not copying it here because it has
> > > some dependencies.  Any idea what is going on?
> > >
> > >
> > >
> > > I had a similar experience - albeit with an older numpy and Python 2.7,
> > > so my comments are easily outdated and irrelevant. This was on Windows
> > > 10 64 bit, way more than plenty RAM.
> > >
> > > It took me forever to find out that numpy.dot was the culprit, and I
> > > ended up using fortran + f2py. Even with the overhead of having to go
> > > through f2py bridge, the fortran dot_product was several times faster.
> > >
> > > Sorry if It doesn’t help much.
> > >
> > > Andrea.
> > >
> > >
> > >
> > >
> > > import numpy as np
> > >
> > > from dot2x1 import dot2x1
> > >
> > > a = np.ones ((1000,16))
> > > b = np.array([ 0.80311816+0.80311816j,  0.80311816-0.80311816j,
> > > -0.80311816+0.80311816j, -0.80311816-0.80311816j,
> > >  1.09707981+0.29396165j,  1.09707981-0.29396165j,
> > > -1.09707981+0.29396165j, -1.09707981-0.29396165j,
> > >  0.29396165+1.09707981j,  0.29396165-1.09707981j,
> > > -0.29396165+1.09707981j, -0.29396165-1.09707981j,
> > >  0.25495815+0.25495815j,  0.25495815-0.25495815j,
> > > -0.25495815+0.25495815j, -0.25495815-0.25495815j])
> > >
> > > def F1():
> > >  d = dot2x1 (a, b)
> > >
> > > def F2():
> > >  d = np.dot (a, b)
> > >
> > > from timeit import timeit
> > > print (timeit ('F1()', globals=globals(), number=1000))
> > > print (timeit ('F2()', globals=globals(), number=1000))
> > >
> > > In [13]: 0.013910860987380147 << 1st timeit
> > > 28.608758996007964  << 2nd timeit
> > > --
> > > Those who don't understand recursion are doomed to repeat it
> > > ___
> > > NumPy-Discussion mailing list
> > > NumPy-Discussion@python.org <mailto:NumPy-Discussion@python.org>
> > > https://mail.python.org/mailman/listinfo/numpy-discussion
> > > <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 mailing list
> > NumPy-Discussion@python.org
> > https://mail.python.org/mailman/listinfo/numpy-discussion
>
>
>
> --
> Those who don't understand recursion are doomed to repeat it



-- 
Those who don't understand recursion are doomed to repeat it
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] C-coded dot 1000x faster than numpy?

[Neal Becker]

See my earlier email - this is fedora 33, python3.9.

I'm using fedora 33 standard numpy.
ldd says:

/usr/lib64/python3.9/site-packages/numpy/core/_multiarray_umath.cpython-39-x86_64-linux-gnu.so:
linux-vdso.so.1 (0x7ffdd1487000)
libflexiblas.so.3 => /lib64/libflexiblas.so.3 (0x7f0512787000)

So whatever flexiblas is doing controls blas.
flexiblas print
FlexiBLAS, version 3.0.4
Copyright (C) 2014, 2015, 2016, 2017, 2018, 2019, 2020 Martin Koehler
and others.
This is free software; see the source code for copying conditions.
There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.


Configured BLAS libraries:
System-wide (/etc/flexiblasrc):

System-wide from config directory (/etc/flexiblasrc.d/)
 OPENBLAS-OPENMP
   library = libflexiblas_openblas-openmp.so
   comment =
 NETLIB
   library = libflexiblas_netlib.so
   comment =
 ATLAS
   library = libflexiblas_atlas.so
   comment =

User config (/home/nbecker/.flexiblasrc):

Host config (/home/nbecker/.flexiblasrc.nbecker8):

Available hooks:

Backend and hook search paths:
  /usr/lib64/flexiblas/

Default BLAS:
System:   OPENBLAS-OPENMP
User: (none)
Host: (none)
Active Default: OPENBLAS-OPENMP (System)
Runtime properties:
   verbose = 0 (System)

So it looks  to me it is using openblas-openmp.


On Tue, Feb 23, 2021 at 8:00 PM Charles R Harris
 wrote:
>
>
>
> On Tue, Feb 23, 2021 at 5:47 PM Charles R Harris  
> wrote:
>>
>>
>>
>> On Tue, Feb 23, 2021 at 11:10 AM Neal Becker  wrote:
>>>
>>> I have code that performs dot product of a 2D matrix of size (on the
>>> order of) [1000,16] with a vector of size [1000].  The matrix is
>>> float64 and the vector is complex128.  I was using numpy.dot but it
>>> turned out to be a bottleneck.
>>>
>>> So I coded dot2x1 in c++ (using xtensor-python just for the
>>> interface).  No fancy simd was used, unless g++ did it on it's own.
>>>
>>> On a simple benchmark using timeit I find my hand-coded routine is on
>>> the order of 1000x faster than numpy?  Here is the test code:
>>> My custom c++ code is dot2x1.  I'm not copying it here because it has
>>> some dependencies.  Any idea what is going on?
>>>
>>> import numpy as np
>>>
>>> from dot2x1 import dot2x1
>>>
>>> a = np.ones ((1000,16))
>>> b = np.array([ 0.80311816+0.80311816j,  0.80311816-0.80311816j,
>>>-0.80311816+0.80311816j, -0.80311816-0.80311816j,
>>> 1.09707981+0.29396165j,  1.09707981-0.29396165j,
>>>-1.09707981+0.29396165j, -1.09707981-0.29396165j,
>>> 0.29396165+1.09707981j,  0.29396165-1.09707981j,
>>>-0.29396165+1.09707981j, -0.29396165-1.09707981j,
>>> 0.25495815+0.25495815j,  0.25495815-0.25495815j,
>>>-0.25495815+0.25495815j, -0.25495815-0.25495815j])
>>>
>>> def F1():
>>> d = dot2x1 (a, b)
>>>
>>> def F2():
>>> d = np.dot (a, b)
>>>
>>> from timeit import timeit
>>> print (timeit ('F1()', globals=globals(), number=1000))
>>> print (timeit ('F2()', globals=globals(), number=1000))
>>>
>>> In [13]: 0.013910860987380147 << 1st timeit
>>> 28.608758996007964  << 2nd timeit
>>
>>
>> I'm going to guess threading, although huge pages can also be a problem on a 
>> machine under heavy load running other processes. Call overhead may also 
>> matter for such small matrices.
>>
>
> What BLAS library are you using. I get much better results using an 8 year 
> old i5 and ATLAS.
>
> Chuck
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion



--
Those who don't understand recursion are doomed to repeat it
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] C-coded dot 1000x faster than numpy?

[Neal Becker]

Supposedly can control through env variables but I didn't see any effect

On Wed, Feb 24, 2021, 10:12 AM Charles R Harris 
wrote:

>
>
> On Wed, Feb 24, 2021 at 8:02 AM Charles R Harris <
> charlesr.har...@gmail.com> wrote:
>
>>
>>
>> On Wed, Feb 24, 2021 at 5:36 AM Neal Becker  wrote:
>>
>>> See my earlier email - this is fedora 33, python3.9.
>>>
>>> I'm using fedora 33 standard numpy.
>>> ldd says:
>>>
>>> /usr/lib64/python3.9/site-packages/numpy/core/_
>>> multiarray_umath.cpython-39-x86_64-linux-gnu.so:
>>> linux-vdso.so.1 (0x7ffdd1487000)
>>> libflexiblas.so.3 => /lib64/libflexiblas.so.3 (0x7f0512787000)
>>>
>>> So whatever flexiblas is doing controls blas.
>>> flexiblas print
>>> FlexiBLAS, version 3.0.4
>>> Copyright (C) 2014, 2015, 2016, 2017, 2018, 2019, 2020 Martin Koehler
>>> and others.
>>> This is free software; see the source code for copying conditions.
>>> There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
>>> FITNESS FOR A PARTICULAR PURPOSE.
>>>
>>>
>>> Configured BLAS libraries:
>>> System-wide (/etc/flexiblasrc):
>>>
>>> System-wide from config directory (/etc/flexiblasrc.d/)
>>>  OPENBLAS-OPENMP
>>>library = libflexiblas_openblas-openmp.so
>>>comment =
>>>  NETLIB
>>>library = libflexiblas_netlib.so
>>>comment =
>>>  ATLAS
>>>library = libflexiblas_atlas.so
>>>comment =
>>>
>>> User config (/home/nbecker/.flexiblasrc):
>>>
>>> Host config (/home/nbecker/.flexiblasrc.nbecker8):
>>>
>>> Available hooks:
>>>
>>> Backend and hook search paths:
>>>   /usr/lib64/flexiblas/
>>>
>>> Default BLAS:
>>> System:   OPENBLAS-OPENMP
>>> User: (none)
>>> Host: (none)
>>> Active Default: OPENBLAS-OPENMP (System)
>>> Runtime properties:
>>>verbose = 0 (System)
>>>
>>> So it looks  to me it is using openblas-openmp.
>>>
>>>
>> ISTR that there have been problems with openmp. There are a ton of
>> OpenBLAS versions available in fedora 33. Just available via flexiblas
>>
>>
>>1. flexiblas-openblas-openmp.x86_64 : FlexiBLAS wrappers for OpenBLAS
>>2. flexiblas-openblas-openmp.i686 : FlexiBLAS wrappers for OpenBLAS
>>3. flexiblas-openblas-openmp64.x86_64 : FlexiBLAS wrappers for
>>OpenBLAS (64-bit)
>>4. flexiblas-openblas-serial.x86_64 : FlexiBLAS wrappers for OpenBLAS
>>5. flexiblas-openblas-serial64.x86_64 : FlexiBLAS wrappers for
>>OpenBLAS (64-bit)
>>6. flexiblas-openblas-threads.x86_64 : FlexiBLAS wrappers for OpenBLAS
>>7. flexiblas-openblas-threads64.x86_64 : FlexiBLAS wrappers for
>>OpenBLAS (64-bit)
>>
>> I am not sure how to make use of flexiblas, but would explore that. We
>> might need to do something with distutils to interoperate with it or maybe
>> you can control it though site.cfg. There are 12 versions available in
>> total. I would suggest trying serial or pthreads.
>>
>>
> Seems to be controlled in the /etc directory:
>
> /etc/flexiblas64rc.d/openblas-openmp64.conf
>
> On my machine it looks like openmp64 is the system default.
>
> Chuck
> ___
> 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] Pi day easter egg

[Neal Becker]

There's a little pi day easter egg for all math fans.  Google for pi to
find it.
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] Indexing question

[Neal Becker]

This seems like something that can be done with indexing, but I
haven't found the solution.

out is a 2D array is initialized to zeros.  x is a 1D array whose
values correspond to the columns of out.  For each row in out, set
out[row,x[row]] = 1.  Here is working code:
def orthogonal_mod (x, nbits):
out = np.zeros ((len(x), 1

Re: [Numpy-discussion] Indexing question

[Neal Becker]

Thanks!

On Thu, May 20, 2021 at 9:53 AM Robert Kern  wrote:
>
> On Thu, May 20, 2021 at 9:47 AM Neal Becker  wrote:
>>
>> This seems like something that can be done with indexing, but I
>> haven't found the solution.
>>
>> out is a 2D array is initialized to zeros.  x is a 1D array whose
>> values correspond to the columns of out.  For each row in out, set
>> out[row,x[row]] = 1.  Here is working code:
>> def orthogonal_mod (x, nbits):
>> out = np.zeros ((len(x), 1<> for e in range (len (x)):
>> out[e,x[e]] = 1
>> return out
>>
>> Any idea to do this without an explicit python loop?
>
>
>
> i = np.arange(len(x))
> j = x[i]
> out[i, j] = 1
>
> --
> Robert Kern
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion



-- 
Those who don't understand recursion are doomed to repeat it
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] EHN: Discusions about 'add numpy.topk'

[Neal Becker]

Topk is a bad choice imo.  I initially parsed it as to_pk, and had no idea
what that was, although sounded a lot like a scipy signal function.
Nlargest would be very obvious.

On Sun, May 30, 2021, 7:50 AM Alan G. Isaac  wrote:

> Mathematica and Julia both seem relevant here.
> Mma has TakeLargest (and Wolfram tends to think hard about names).
> https://reference.wolfram.com/language/ref/TakeLargest.html
> Julia's closest comparable is perhaps partialsortperm:
> https://docs.julialang.org/en/v1/base/sort/#Base.Sort.partialsortperm
> Alan Isaac
>
>
>
> On 5/30/2021 4:40 AM, kang...@mail.ustc.edu.cn wrote:
> > Hi, Thanks for reply, I present some details below:
> ___
> 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] Add count (and dtype) to packbits

[Neal Becker]

In my application I need to pack bits of a specified group size into
integral values.
Currently np.packbits only packs into full bytes.
For example, I might have a string of bits encoded as a np.uint8
vector with each uint8 item specifying a single bit 1/0.  I want to
encode them 4 bits at a time into a np.uint32 vector.

python code to implement this:

---
def pack_bits (inp, bits_per_word, dir=1, dtype=np.int32):
assert bits_per_word <= np.dtype(dtype).itemsize * 8
assert len(inp) % bits_per_word == 0
out = np.empty (len (inp)//bits_per_word, dtype=dtype)
i = 0
o = 0
while i < len(inp):
ret = 0
for b in range (bits_per_word):
if dir > 0:
ret |= inp[i] << b
else:
ret |= inp[i] << (bits_per_word - b - 1)
i += 1
out[o] = ret
o += 1
return out
---

It looks like unpackbits has a "count" parameter but packbits does not.
Also would be good to be able to specify an output dtype.
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Add count (and dtype) to packbits

[Neal Becker]

Well that's just the point, I wanted to consider group size > 8.

On Wed, Jul 21, 2021 at 8:53 AM Andras Deak  wrote:
>
> On Wed, Jul 21, 2021 at 2:40 PM Neal Becker  wrote:
>>
>> In my application I need to pack bits of a specified group size into
>> integral values.
>> Currently np.packbits only packs into full bytes.
>> For example, I might have a string of bits encoded as a np.uint8
>> vector with each uint8 item specifying a single bit 1/0.  I want to
>> encode them 4 bits at a time into a np.uint32 vector.
>>
>> python code to implement this:
>>
>> ---
>> def pack_bits (inp, bits_per_word, dir=1, dtype=np.int32):
>> assert bits_per_word <= np.dtype(dtype).itemsize * 8
>> assert len(inp) % bits_per_word == 0
>> out = np.empty (len (inp)//bits_per_word, dtype=dtype)
>> i = 0
>> o = 0
>> while i < len(inp):
>> ret = 0
>> for b in range (bits_per_word):
>> if dir > 0:
>> ret |= inp[i] << b
>> else:
>> ret |= inp[i] << (bits_per_word - b - 1)
>> i += 1
>> out[o] = ret
>> o += 1
>> return out
>> ---
>
>
> Can't you just `packbits` into a uint8 array and then convert that to uint32? 
> If I change `dtype` in your code from `np.int32` to `np.uint32` (as you 
> mentioned in your email) I can do this:
>
> rng = np.random.default_rng()
> arr = (rng.uniform(size=32) < 0.5).astype(np.uint8)
> group_size = 4
> original = pack_bits(arr, group_size, dtype=np.uint32)
> new = np.packbits(arr.reshape(-1, group_size), axis=-1, 
> bitorder='little').ravel().astype(np.uint32)
> print(np.array_equal(new, original))
> # True
>
> There could be edge cases where the result dtype is too small, but I haven't 
> thought about that part of the problem. I assume this would work as long as 
> `group_size <= 8`.
>
> András
>
>>
>> It looks like unpackbits has a "count" parameter but packbits does not.
>> Also would be good to be able to specify an output dtype.
>> ___
>> 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



-- 
Those who don't understand recursion are doomed to repeat it
___
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] Re: A bite of C++

[Neal Becker]

A couple of questions from a quick casual reading

1. radixsort (void *start...)
Do we really need void*?  We don't know the type of start at compile time?

2. reinterpret_cast start (related to #1).

3. reinterpret_cast(malloc(num * sizeof(T)));
A C-ism.  Would it work to use  new T[num]?

On Tue, Sep 28, 2021 at 10:03 PM Sebastian Berg
 wrote:
>
> On Wed, 2021-08-25 at 17:48 +0200, Serge Guelton wrote:
> > Hi folks,
> >
> > https://github.com/numpy/numpy/pull/19713 showcases what *could* be a
> > first step
> > toward getting rid of generated C code within numpy, in favor of some
> > C++ code,
> > coupled with a single macro trick.
>
>
> It seems time to pick up this discussion again.  I think we were
> cautiously in favor of trying this out in the limited proposed form.
>
> Aside from ensuring we do not use C++ exceptions/run-time environment.
> I.e. compile with `-fno-exception` and `-nostdlib`, the PR is probably
> largely ready.
>
> Are there any other bigger issues that need to be discussed?
>
>
> The goal for now seems to allow this type of function templates (C++ as
> a glorified templating language).  We may be able to enforce this with
> clang-tidy (not sure if that is necessary right away though).
> There are parts of NumPy (mainly `fft`) where more C++ features could
> be helpful but that should probably be discussed separately.
>
>
> The current solution ends up using no X-Macro anymore, which means
> listing the C symbols like:
>
> NPY_NO_EXPORT int radixsort_byte(void* vec, npy_intp cnt, void 
> *NPY_UNUSED(null)) { return radixsort(vec, cnt); }
> NPY_NO_EXPORT int radixsort_ubyte(void* vec, npy_intp cnt, void 
> *NPY_UNUSED(null)) { return radixsort(vec, cnt); }
> NPY_NO_EXPORT int radixsort_short(void* vec, npy_intp cnt, void 
> *NPY_UNUSED(null)) { return radixsort(vec, cnt); }
> NPY_NO_EXPORT int radixsort_ushort(void* vec, npy_intp cnt, void 
> *NPY_UNUSED(null)) { return radixsort(vec, cnt); }
> NPY_NO_EXPORT int radixsort_int(void* vec, npy_intp cnt, void 
> *NPY_UNUSED(null)) { return radixsort(vec, cnt); }
> NPY_NO_EXPORT int radixsort_uint(void* vec, npy_intp cnt, void 
> *NPY_UNUSED(null)) { return radixsort(vec, cnt); }
> NPY_NO_EXPORT int radixsort_long(void* vec, npy_intp cnt, void 
> *NPY_UNUSED(null)) { return radixsort(vec, cnt); }
> NPY_NO_EXPORT int radixsort_ulong(void* vec, npy_intp cnt, void 
> *NPY_UNUSED(null)) { return radixsort(vec, cnt); }
>
> I guess we could try to write some macro's to make this shorter,
> although I am fine with delaying this, since the main point is probably
> getting a brief example and the build changes in.
>
> It also defines the following type tags:
> https://github.com/numpy/numpy/pull/19713/files#diff-4de4beaecd7b0a1f4d887221191843a54b0c62f13eea18b8716c54bec9657f20
>
> Any comments to the approach (especially from anyone more familiar with
> C++) would be much appreciated!
>
> Cheers,
>
> Sebastian
>
>
>
> >
> > Basically, templated code is an easy and robust way to replace
> > generated code
> > (the C++ compiler becomes the code generator when instantiating
> > code), and a
> > single X-macro takes care of the glue with the C world.
> >
> > Some changes in distutils were needed to cope with C++-specific
> > flags, and
> > extensions that consist in mixed C and C++ code.
> >
> > I've kept the change as minimal as possible to ease the (potential)
> > transition
> > and keep the C++ code close to the C code. This led to less idiomatic
> > C++ code,
> > but I value a "correct first" approach. There's an on-going effort by
> > seiko2plus
> > to remove that C layer, I acknowledge this would bring even more C++
> > code, but
> > that looks orthogonal to me (and a very good second step!)
> >
> > All lights are green for the PR, let's assume it's a solid ground for
> > discussion :-)
> > So, well, what do you think? Should we go forward?
> >
> > Potential follow-ups :
> >
> > - do we want to use -nostdlib, to be sure we don't bring any C++
> > runtime dep?
> > - what about -fno-exception, -fno-rtti?
> > - coding style?
> > - (I'm-not-a-farseer-I-don-t-know-all-topics)
> >
> > ___
> > NumPy-Discussion mailing list
> > NumPy-Discussion@python.org
> > https://mail.python.org/mailman/listinfo/numpy-discussion
> >
>
> ___
> NumPy-Discussion mailing list -- numpy-discussion@python.org
> To unsubscribe send an email to numpy-discussion-le...@python.org
> https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
> Member address: ndbeck...@gmail.com



-- 
Those who don't understand recursion are doomed to repeat it
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] Re: Code formatters

[Neal Becker]

I've also used uncrustify for c/c++.  Of course this plays hell with
revision control.

On Sun, Nov 14, 2021 at 6:29 PM Juan Nunez-Iglesias  wrote:
>
>
>
> On 15 Nov 2021, at 8:23 am, Stefan van der Walt  wrote:
>
> On Sun, Nov 14, 2021, at 09:13, Charles R Harris wrote:
>
> The black formatter is much improved in its latest version and I think good 
> enough to start using. The main drawbacks that I see are:
>
> all operators, including '*' and '/',  get spaces around them,
> very long strings are not broken into multiple lines,
> lists, tuples, and function signatures are either on one line, or broken into 
> multiple lines of one element/argument each,
> the formatting of extended logical expressions could be improved to emphasize 
> the priority of 'and' over 'or' operators
>
>
> We've also been having a conversation around mathematical formatting here: 
> https://discuss.scientific-python.org/t/how-to-format-mathematical-expressions/62/8
>
> I tried yapf recently, and was pleased with the output.  One concern about 
> yapf was that it has many configuration options: but the only important thing 
> is that you fix the knobs, then you simply have a different version of black.
>
>
> In my experience, while none of these tools are perfect, not having to have 
> discussions around formatting is completely worth it!
>
>
> +1 on everything Stéfan said. I never liked black’s formatting, but I have 
> *absolutely* appreciated having zero discussions/push commits/code 
> suggestions to deal with formatting in napari. I have since added yapf to my 
> own repos with a config I like *and* added yapf auto-formatting-on-save to my 
> VSCode, and I don’t even have to have formatting discussions with *myself* 
> anymore. 😂 It’s very liberating!
>
> For reference, here’s my yapf config:
>
> https://github.com/jni/skan/blob/74507344b4cd4453cc43b4dbd0b5742fc08eb5a0/.style.yapf
>
> As Stéfan said, fix the knobs (yours might be different), then forget about 
> it!
>
> Oh, and yes, yapf does allow formatting only the diff. I agree that 
> reformatting the entire code base is problematic.
>
> Juan.
> ___
> NumPy-Discussion mailing list -- numpy-discussion@python.org
> To unsubscribe send an email to numpy-discussion-le...@python.org
> https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
> Member address: ndbeck...@gmail.com



-- 
Those who don't understand recursion are doomed to repeat it
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] Re: Exporting numpy arrays to binary JSON (BJData) for better portability

[Neal Becker]

>
>
> the loading time (from an nvme drive, Ubuntu 18.04, python 3.6.9, numpy
> 1.19.5) for each file is listed below:
>
> 0.179s  eye1e4.npy (mmap_mode=None)
> 0.001s  eye1e4.npy (mmap_mode=r)
> 0.718s  eye1e4_bjd_raw_ndsyntax.jdb
> 1.474s  eye1e4_bjd_zlib.jdb
> 0.635s  eye1e4_bjd_lzma.jdb
>
>
> clearly, mmapped loading is the fastest option without a surprise; it is
> true that the raw bjdata file is about 5x slower than npy loading, but
> given the main chunk of the data are stored identically (as contiguous
> buffer), I suppose with some optimization of the decoder, the gap between
> the two can be substantially shortened. The longer loading time of
> zlib/lzma (and similarly saving times) reflects a trade-off between smaller
> file sizes and time for compression/decompression/disk-IO.
>
> I think the load time for mmap may be deceptive, it isn't actually loading
> anything, just mapping to memory.  Maybe a better benchmark is to actually
> process the data, e.g., find the mean which would require reading the
> values.
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] broadcasting question

[Neal Becker]

I have a function F
def F(a, b):
c = a * b

Initially, a is a scalar, b[240,3000].  No problem.
Later I want to use F, where a[240] is a vector.  I want to allow both the
scalar and vector cases.  So I write:

def F(a,b):
  a = np.atleast_1d(a)
  c = a[:,None] * b

This now works for scalar a or vector a.  But this solutions seems
inelegant, and somewhat fragile.  Suppose later we want to allow
a[240,3000], a 2d array matching b.

Certainly don't want to write code like:
if a.ndim == 0:...

Is there a more elegant/robust approach?

Thanks,
Neal
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] Re: broadcasting question

[Neal Becker]

On Thu, Mar 23, 2023 at 5:21 AM Sebastian Berg 
wrote:

> On Wed, 2023-03-22 at 12:00 -0400, Robert Kern wrote:
> > On Wed, Mar 22, 2023 at 9:34 AM Neal Becker 
> > wrote:
> >
> > > I have a function F
> > > def F(a, b):
> > > c = a * b
> > >
> > > Initially, a is a scalar, b[240,3000].  No problem.
> > > Later I want to use F, where a[240] is a vector.  I want to allow
> > > both the
> > > scalar and vector cases.  So I write:
> > >
> > > def F(a,b):
> > >   a = np.atleast_1d(a)
> > >   c = a[:,None] * b
> > >
> > > This now works for scalar a or vector a.  But this solutions seems
> > > inelegant, and somewhat fragile.  Suppose later we want to allow
> > > a[240,3000], a 2d array matching b.
> > >
> > > Certainly don't want to write code like:
> > > if a.ndim == 0:...
> > >
> > > Is there a more elegant/robust approach?
> > >
> >
> > I would leave it as `c = a * b` and simply record in the docstring
> > that `a`
> > and `b` should be broadcastable. Yes, that means that the user will
> > have to
> > write `F(a[:, np.newaxis], b)` for that one case, and that looks a
> > little
> > ugly, but overall it's less cognitive load on the user to just reuse
> > the
> > common convention of broadcasting than to record the special case.
>
>
> I will note that it is not hard to insert the new axes.
> `np.expand_dims` may be convenient.  many functions (ufuncs) also have
> the `outer` version which does this: `np.add.outer()`, etc.
>
> However, I agree.  Unless the use-case exceedingly clear about
> requiring "outer" behavior.  "outer" behavior is uncommon for functions
> in the NumPy world and broadcasting is what users will generally expect
> (and that includes future self).
>
> - Sebastian
>
> Thanks for the advice!  On reflection, I agree.
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] Re: mean_std function returning both mean and std

[Neal Becker]

On a somewhat related note, I usually find I need to compute stats
incrementally.  To do this, a stat object is created so batches of samples
can be fed to it sequentially.

I used to use an implementation based on boost::accumulator for this.  More
recently I'm using my own c++ code based on xtensor, exposed to python with
xtensor-python and pybind11.

The basic technique to find 2nd order stats is to keep 2 running sums,
sum(x) and sum(x**2).

It would be useful to have functionality for incremental stats like this in
numpy, as well as other incremental operations (e.g., histogram).  I
frequently find I need to process large amounts of data in small batches at
a time, generated by iterative monte-carlo simulations, for example.
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] np.ndenumerate doesn't obey mask?

[Neal Becker via NumPy-Discussion]

I was using ndenuerate with a masked array, and it seems that the mask is
ignored.  Is this true?  If so, isn't that a bug?
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] Re: NumPy 2.2.0 Released

[Neal Becker via NumPy-Discussion]

Where can I find more information on improvements to stringdtype?

On Sun, Dec 8, 2024, 11:25 AM Charles R Harris via NumPy-Discussion <
numpy-discussion@python.org> wrote:

> Hi All,
>
> On behalf of the NumPy team, I'm pleased to announce the release of NumPy
> 2.2.0. The NumPy 2.2.0 release is a short release that brings us back
> into sync with the usual twice yearly release cycle. There have been a
> number of small cleanups, as well as work bringing the new StringDType to
> completion and improving support for free threaded Python. Highlights are:
>
>- New functions `matvec` and `vecmat`, see below.
>- Many improved annotations.
>- Improved support for the new StringDType.
>- Improved support for free threaded Python
>- Fixes for f2py
>
> This release supports Python 3.10-3.13. Wheels can be downloaded from PyPI
> ; source archives, release notes,
> and wheel hashes are available on Github
> .
>
> Cheers,
>
> Charles Harris
> ___
> NumPy-Discussion mailing list -- numpy-discussion@python.org
> To unsubscribe send an email to numpy-discussion-le...@python.org
> https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
> Member address: ndbeck...@gmail.com
>
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com

[Numpy-discussion] Re: ENH: Add saturating arithmetic functions

[Neal Becker via NumPy-Discussion]

On Mon, Apr 28, 2025 at 4:52 AM Ralf Gommers via NumPy-Discussion <
numpy-discussion@python.org> wrote:

>
>
> On Sun, Apr 27, 2025 at 7:31 AM Carlos Martin 
> wrote:
>
>> Saturating arithmetic (
>> https://en.wikipedia.org/wiki/Saturation_arithmetic) is important in
>> digital signal processing and other areas.
>>
>> Feature request: Add saturating arithmetic functions for the following
>> basic operations:
>>
>> - addition (C++ counterpart:
>> https://en.cppreference.com/w/cpp/numeric/add_sat)
>> - subtraction (C++ counterpart:
>> https://en.cppreference.com/w/cpp/numeric/sub_sat)
>> - multiplication (C++ counterpart:
>> https://en.cppreference.com/w/cpp/numeric/mul_sat)
>> - division (C++ counterpart:
>> https://en.cppreference.com/w/cpp/numeric/div_sat)
>> - casting (C++ counterpart:
>> https://en.cppreference.com/w/cpp/numeric/saturate_cast)
>> - negation
>>
>> I've implemented these for JAX at
>> https://gist.github.com/carlosgmartin/b32fa6fed3aa82f83dfbaac4b6345672.
>>
>> Corresponding issue: https://github.com/jax-ml/jax/issues/26566.
>>
>
> Thanks for the proposal Carlos. On
> https://github.com/numpy/numpy/issues/28837 Matti suggested that it may
> be possible to implement this as a casting mode. If so, then it may be
> feasible to have this functionality inside NumPy. If it would require new
> API beyond that, it seems a bit niche for NumPy, but if the amount of extra
> code/maintenance is reasonable and it's only one extra casting mode, then
> it does seem reasonable to me.
>
> Cheers,
> Ralf
>

I have worked in this area (DSP hardware) and have developed and used tools
for such algorithms.   I think this is part of a larger topic of fixed
point arithmetic.  There are 2 aspects to fixed point arithmetic that we
can consider.  One is the interpretation of a binary point.  The other
aspect is the handling of arithmetic overflow for a finite field of bits.
It is this latter aspect that I think is currently being discussed.  When
studying overflow effects in DSP algorithms it is useful to be able to
choose between at least 3 different behaviors: 1) wraparound 2) saturation
3) throw error.  In hardware design saturation logic has a cost (in both
space and time) and so should only be added when needed.  In my own work I
have implemented (in c++) data types that allow switching overflow behavior
as a parameter of the class, allowing the designer to study the effect of
the different choices.  The "throw error" exists to find whether the
potential exists for overflow and so indicate whether there is a need for
any additional circuitry.

Perhaps this should be considered as a part of a more comprehensive effort
to implement fixed-point arithmetic?

Thanks,
Neal
___
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com