Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples

[Oscar Benjamin]

On Tue, 1 Sep 2015 18:43 Phil Hodge  wrote:

On 09/01/2015 11:14 AM, Oscar Benjamin wrote:
> Just use the next power of 2. Pure powers of 2 are the most efficient
> for FFT algorithms so it potentially works out better than finding a
> smaller but similarly composite size to pad to. Finding the next power
> of 2 is easy to code and never a bad choice.

It would be better to pad to a power of three or five, or some other
product of small primes, if the pad length would be significantly less
than padding to a power of two.  The obvious problem with padding is
that it's different from the real data, so its presence will affect the
result.  At the least, it would be good to pad with the mean value of
the original data, or to pad with zero after subtracting the mean from
the original data.



I meant performance wise it's not a bad choice. If you're concerned about
distortion then use the Bluestein algorithm as I showed.

--
Oscar
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples

[Joseph Codadeen]

Ah, looking back I see what Pierre-André did, the penny drops. Sorry I 
miss-read.I'm only interested in this part;>zeropadded_fft_A = fft(seq_A, 
n=2**(ceil(log(len(seq_A),2))+1))
>zeropadded_fft_B = fft(seq_B, n=2**(ceil(log(len(seq_B),2))+1))>You could 
>remove the "+1" above to get faster results, but then that may >lead to 
>unwanted frequencies (coming from the fact that fft assumes >periodic signals, 
>read online about zero-padding).
I did it in the following way and it seems to work. I've seen an increase in 
performance for sure.
My Python sucks by the way :)def pad_to_power_of_two(audio):
'''
Calculates next power of two and pads the audio sample array with zeros
'''# RIFF spec specifies: ckSizeA 32-bit unsigned value identifying 
the size of ckData.
max_shift = 32
shift = 1
power = len(audio) - 1
print "{0:b}".format(power)# Test if a power of two. Turn off the 
right-most bit (x & (x - 1)) and zero test
while ((power + 1) & power) and (shift <= max_shift):
power |= power >> shift
shift *= 2
print  "new power: " + "{0:b}".format(power)
power = power + 1
print  "{0:b}".format(power)

ext = np.array([0x00] * (power - len(audio)))
new_array = np.concatenate((audio, ext), axis=0)

print  "Next power of 2 greater than actual array length is " + str(power)
print  "Extending current array length (" + str(len(audio)) + ") by " + 
str(power - len(audio)) + " null bytes"
print  "New array length is " + str(len(new_array))

return new_array
I will try using the suggested methods and see if I can boost the performance 
further.Thanks.
> From: oscar.j.benja...@gmail.com
> Date: Tue, 1 Sep 2015 16:14:41 +0100
> To: numpy-discussion@scipy.org
> Subject: Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples
> 
> On 1 September 2015 at 11:38, Joseph Codadeen  wrote:
> >
> >> And while you zero-pad, you can zero-pad to a sequence that is a power of
> >> two, thus preventing awkward factorizations.
> >
> > Does numpy have an easy way to do this, i.e. for a given number, find the
> > next highest number (within a range) that could be factored into small,
> > prime numbers as Phil explained? It would help if it gave a list,
> > prioritised by number of factors.
> 
> Just use the next power of 2. Pure powers of 2 are the most efficient
> for FFT algorithms so it potentially works out better than finding a
> smaller but similarly composite size to pad to. Finding the next power
> of 2 is easy to code and never a bad choice.
> 
> To avoid the problems mentioned about zero-padding distorting the FFT
> you can use the Bluestein transform as below. This pads up to a power
> of two (greater than 2N-1) but does it in a special way that ensures
> that it is still calculating the exact (up to floating point error)
> same DFT:
> 
> from numpy import array, exp, pi, arange, concatenate
> from numpy.fft import fft, ifft
> 
> def ceilpow2(N):
> '''
> >>> ceilpow2(15)
> 16
> >>> ceilpow2(16)
> 16
> '''
> p = 1
> while p < N:
> p *= 2
> return p
> 
> def fftbs(x):
> '''
> >>> data = [1, 2, 5, 2, 5, 2, 3]
> >>> from numpy.fft import fft
> >>> from numpy import allclose
> >>> from numpy.random import randn
> >>> for n in range(1, 1000):
> ... data = randn(n)
> ... assert allclose(fft(data), fftbs(data))
> '''
> N = len(x)
> x = array(x)
> 
> n = arange(N)
> b = exp((1j*pi*n**2)/N)
> a = x * b.conjugate()
> 
> M = ceilpow2(N) * 2
> A = concatenate((a, [0] * (M - N)))
> B = concatenate((b, [0] * (M - 2*N + 1), b[:0:-1]))
> C = ifft(fft(A) * fft(B))
> c = C[:N]
> return b.conjugate() * c
> 
> if __name__ == "__main__":
> import doctest
> doctest.testmod()
> 
> 
> --
> Oscar
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

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

Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples

[Phil Hodge]

On 09/01/2015 11:14 AM, Oscar Benjamin wrote:
> Just use the next power of 2. Pure powers of 2 are the most efficient
> for FFT algorithms so it potentially works out better than finding a
> smaller but similarly composite size to pad to. Finding the next power
> of 2 is easy to code and never a bad choice.

It would be better to pad to a power of three or five, or some other 
product of small primes, if the pad length would be significantly less 
than padding to a power of two.  The obvious problem with padding is 
that it's different from the real data, so its presence will affect the 
result.  At the least, it would be good to pad with the mean value of 
the original data, or to pad with zero after subtracting the mean from 
the original data.

Phil
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] np.sign and object comparisons

[Nathaniel Smith]

On Sun, Aug 30, 2015 at 9:09 PM, Jaime Fernández del Río <
jaime.f...@gmail.com> wrote:

>
> There are three ways of fixing this that I see:
>
>1. Arbitrarily choose a value to set the return to. This is equivalent
>to choosing a default return for `cmp` for comparisons. This preserves
>behavior, but feels wrong.
>2. Similarly to how np.sign of a floating point array with nans
>returns nan for those values, return e,g, None for these cases. This
>is my preferred option.
>3. Raise an error, along the lines of the TypeError: unorderable types
>that 3.x produces for some comparisons.
>
> Having read the other replies so far -- given that no-one seems to have
any clear intuition or use cases, I guess I find option 3 somewhat
tempting... it keeps our options open until someone who actually cares
comes along with a use case to hone our intuition on, and is very safe in
the mean time.

(This was noticed in the course of routine code cleanups, right, not an
external bug report? For all we know right now, no actual user has ever
even tried to apply np.sign to an object array?)

-n

-- 
Nathaniel J. Smith -- http://vorpus.org
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] np.sign and object comparisons

[Nathaniel Smith]

On Mon, Aug 31, 2015 at 10:31 AM, Antoine Pitrou  wrote:
> On Mon, 31 Aug 2015 10:23:10 -0700
> Stephan Hoyer  wrote:
>>
>> My inclination is that return NaN would be the appropriate choice. It's
>> certainly consistent with the behavior for float dtypes -- my expectation
>> for object dtype behavior is that it works exactly like applying the
>> np.sign ufunc to each element of the array individually.
>>
>> On the other hand, I suppose there are other ways in which an object can
>> fail all those comparisons (e.g., NaT?), so I suppose we could return None.
>
> Currently:
>
 np.sign(np.timedelta64('nat'))
> numpy.timedelta64(-1)
>
> ... probably because NaT is -2**63 under the hood. But in this case
> returning NaT would sound better.

I think this is going through the np.sign timedelta64 loop, and thus
is an unrelated issue? It does look like a bug though.

-n

-- 
Nathaniel J. Smith -- http://vorpus.org
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Cythonizing some of NumPy

[Nathaniel Smith]

On Sun, Aug 30, 2015 at 2:44 PM, David Cournapeau  wrote:
> Hi there,
>
> Reading Nathaniel summary from the numpy dev meeting, it looks like there is
> a consensus on using cython in numpy for the Python-C interfaces.
>
> This has been on my radar for a long time: that was one of my rationale for
> splitting multiarray into multiple "independent" .c files half a decade ago.
> I took the opportunity of EuroScipy sprints to look back into this, but
> before looking more into it, I'd like to make sure I am not going astray:
>
> 1. The transition has to be gradual

Yes, definitely.

> 2. The obvious way I can think of allowing cython in multiarray is modifying
> multiarray such as cython "owns" the PyMODINIT_FUNC and the module
> PyModuleDef table.

The seems like a plausible place to start.

In the longer run, I think we'll need to figure out a strategy to have
source code divided over multiple .pyx files (for the same reason we
want multiple .c files -- it'll just be impossible to work with
otherwise). And this will be difficult for annoying technical reasons,
since we definitely do *not* want to increase the API surface exposed
by multiarray.so, so we will need to compile these multiple .pyx and
.c files into a single module, and have them talk to each other via
internal interfaces. But Cython is currently very insistent that every
.pyx file should be its own extension module, and the interface
between different files should be via public APIs.

I spent some time poking at this, and I think it's possible but will
take a few kluges at least initially. IIRC the tricky points I noticed
are:

- For everything except the top-level .pyx file, we'd need to call the
generated module initialization functions "by hand", and have a bit of
utility code to let us access the symbol tables for the resulting
modules

- We'd need some preprocessor hack (or something?) to prevent the
non-main module initialization functions from being exposed at the .so
level (like 'cdef extern from "foo.h"', 'foo.h' re#defines
PyMODINIT_FUNC to remove the visibility declaration)

- By default 'cdef' functions are name-mangled, which is annoying if
you want to be able to do direct C calls between different .pyx and .c
files. You can fix this by adding a 'public' declaration to your cdef
function. But 'public' also adds dllexport stuff which would need to
be hacked out as per above.

I think the best strategy for this is to do whatever horrible things
are necessary to get an initial version working (on a branch, of
course), and then once that's done assess what changes we want to ask
the cython folks for to let us eliminate the gross parts.

(Insisting on compiling everything into the same .so will probably
also help at some point in avoiding Cython-Related Binary Size Blowup
Syndrome (CRBSBS), because the masses of boilerplate could in
principle be shared between the different files. I think some modern
linkers are even clever enough to eliminate this kind of duplicate
code automatically, since C++ suffers from a similar problem.)

> 3. We start using cython for the parts that are mostly menial refcount work.
> Things like functions in calculation.c are obvious candidates.
>
> Step 2 should not be disruptive, and does not look like a lot of work: there
> are < 60 methods in the table, and most of them should be fairly
> straightforward to cythonize. At worse, we could just keep them as is
> outside cython and just "export" them in cython.
>
> Does that sound like an acceptable plan ?
>
> If so, I will start working on a PR to work on 2.

Makes sense to me!

-n

-- 
Nathaniel J. Smith -- http://vorpus.org
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] UTC-based datetime64

[Chris Barker]

> Googling for a way to print UTC out of the box, the best thing I could
> find is:
>
> In [40]: [str(i.item()) for i in np.array([t], dtype="datetime64[s]")]
> Out[40]: ['2015-08-26 11:52:10']
>
> Now, is there a better way to specify that I want the datetimes printed
> always in UTC?
>

maybe, but it's a kludge no matter how you do it :-(

It's been in the list for a while to fix this, but hasn't happened yet.
Partly due to me not finishing writing out the notes from SciPy 2014 where
we discussed a way forward on this.

Sorry,

-Chris




> Thanks,
> --
> Francesc Alted
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> 
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>


-- 

Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR(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



-- 

Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR(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@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples

[Charles R Harris]

On Tue, Sep 1, 2015 at 12:06 PM, Oscar Benjamin 
wrote:

>
> On Tue, 1 Sep 2015 18:43 Phil Hodge  wrote:
>
> On 09/01/2015 11:14 AM, Oscar Benjamin wrote:
> > Just use the next power of 2. Pure powers of 2 are the most efficient
> > for FFT algorithms so it potentially works out better than finding a
> > smaller but similarly composite size to pad to. Finding the next power
> > of 2 is easy to code and never a bad choice.
>
> It would be better to pad to a power of three or five, or some other
> product of small primes, if the pad length would be significantly less
> than padding to a power of two.  The obvious problem with padding is
> that it's different from the real data, so its presence will affect the
> result.  At the least, it would be good to pad with the mean value of
> the original data, or to pad with zero after subtracting the mean from
> the original data.
>
>
>
> I meant performance wise it's not a bad choice. If you're concerned about
> distortion then use the Bluestein algorithm as I showed.
>

I don't see the problem with padding. After all, the spectrum is an
estimate and the effect of the padding, as well as the use of windowing is
well understood. However, I would recommend subtracting the average before
padding, as otherwise the DC frequency is likely to be of large amplitude
and its side lobes will mask the lower frequencies.


> --
> Oscar
>
>
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] UTC-based datetime64

[Jim Benson]

On Tue, Sep 01, 2015 at 05:28:56PM -0700, Chris Barker wrote:
> > Googling for a way to print UTC out of the box, the best thing I could
> > find is:
> >
> > In [40]: [str(i.item()) for i in np.array([t], dtype="datetime64[s]")]
> > Out[40]: ['2015-08-26 11:52:10']
> >
> > Now, is there a better way to specify that I want the datetimes printed
> > always in UTC?
> >
> 
> maybe, but it's a kludge no matter how you do it :-(
> 
> It's been in the list for a while to fix this, but hasn't happened yet.
> Partly due to me not finishing writing out the notes from SciPy 2014 where
> we discussed a way forward on this.
> 
> Sorry,
> 
> -Chris
> 

Plans for getting times in UTC1 at the ms accuracy?

Sorry, a joke, i couldn't resist.

Cheers,

Jim

> 
> 
> 
> > Thanks,
> > --
> > Francesc Alted
> >
> > ___
> > NumPy-Discussion mailing list
> > NumPy-Discussion@scipy.org
> > 
> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >
> >
> 
> 
> -- 
> 
> Christopher Barker, Ph.D.
> Oceanographer
> 
> Emergency Response Division
> NOAA/NOS/OR(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
> 
> 
> 
> -- 
> 
> Christopher Barker, Ph.D.
> Oceanographer
> 
> Emergency Response Division
> NOAA/NOS/OR(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@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Notes from the numpy dev meeting at scipy 2015

[Nathaniel Smith]

On Sun, Aug 30, 2015 at 9:12 PM, Marten van Kerkwijk
 wrote:
> Hi Nathaniel, others,
>
> I read the discussion of plans with interest. One item that struck me is
> that while there are great plans to have a proper extensible and presumably
> subclassable dtype, it is discouraged to subclass ndarray itself (rather, it
> is encouraged to use a broader array interface). From my experience with
> astropy in both Quantity (an ndarray subclass), Time (a separate class
> containing high precision times using two ndarray float64), and Table
> (initially holding structured arrays, but now sets of Columns, which
> themselves are ndarray subclasses), I'm not convinced the broader, new
> containers approach is that much preferable. Rather, it leads to a lot of
> boiler-plate code to reimplement things ndarray does already (since one is
> effectively just calling the methods on the underlying arrays).
>
> I also think the idea that a dtype becomes something that also contains a
> unit is a bit odd. Shouldn't dtype just be about how data is stored? Why
> include meta-data such as units?
>
> Instead, I think a quantity is most logically seen as numbers with a unit,
> just like masked arrays are numbers with masks, and variables numbers with
> uncertainties. Each of these cases adds extra information in different
> forms, and all are quite easily thought of as subclasses of ndarray where
> all operations do the normal operation, plus some extra work to keep the
> extra information up to date.

The intuition behind the array/dtype split is that an array is just a
container: it knows how to shuffle bytes around, be reshaped, indexed,
etc., but it knows nothing about the meaning of the items it holds --
as far as it's concerned, each entry is just an opaque binary blobs.
If it wants to actually do anything with these blobs, it has to ask
the dtype for help.

The dtype, OTOH, knows how to interpret these blobs, and (in
cooperation with ufuncs) to perform operations on them, but it doesn't
need to know how they're stored, or about slicing or anything like
that -- all that's the container's job.

Think about it this way: does it make sense to have a sparse array of
numbers-with-units? how about a blosc-style compressed array of
numbers-with-units? If yes, then numbers-with-units are a special kind
of dtype, not a special kind of array.

Another way of getting this intuition: if I have 8 bytes, that could
be an int64, or it could be a float64. Which one it is doesn't affect
how it's stored at all -- either way it's stored as a chunk of 8
arbitrary bytes. What it affects is how we *interpret* these bytes --
e.g. there is one function called "int64 addition" which takes two 8
byte chunks and returns a new 8 byte chunk as the result, and a second
function called "float64 addition" which takes those same two 8 byte
chunks and returns a different one. The dtype tells you which of these
operations should be used for a particular array. What's special about
a float64-with-units? Well, it's 8 bytes, but the addition operation
is different from regular float64 addition: it has to do some extra
checks and possibly unit conversions. This is exactly what the ufunc
dtype dispatch and casting system is there for.

This also solves your problem with having to write lots of boilerplate
code, b/c if this is a dtype then it means you can just use the actual
ndarray class directly without subclassing or anything :-).

> Anyway, my suggestion would be to *encourage* rather than discourage ndarray
> subclassing, and help this by making ndarray (even) better.

So, we very much need robust support for
objects-that-quack-like-an-array that are *not* ndarrays, because
ndarray subclasses are forced to use ndarray-style strided in-memory
storage, and there's huge demand for objects that expose an array-like
interface but that use a different storage strategy underneath: sparse
arrays, compressed arrays (like blosc), out-of-core arrays,
computed-on-demand arrays (like dask), distributed arrays, etc. etc.

And once we have solid support for duck-arrays and for user-defined
dtypes (as discussed above), then those two things remove a huge
amount of the motivation for subclassing ndarray.

At the same time, ndarray subclassing is... nearly unmaintainable,
AFAICT. The problem with subclassing is that you're basically taking
some interface, making a copy of it, and then monkeypatching the copy.
As you would expect, this is intrinsically very fragile, because it
breaks abstraction barriers. Suddenly things that used to be
implementation details -- like which methods are implemented in terms
of which other methods -- become part of the public API. And there's
never been any coherent, documentable theory of how ndarray
subclassing is *supposed* to work, so in practice it's just a bunch of
ad hoc hooks designed around the needs of np.matrix and np.ma. We get
a regular stream of bug reports asking us to tweak things one way or
another, and 

Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples

[Joseph Codadeen]

Hi,
I cannot see how the following would work  when it is np.fft.fft() that takes a 
long time based on the length of data. In my case my data is non-periodic.> 
from numpy.fft import fft> from numpy.random import rand> from math import log, 
ceil> seq_A = rand(2649674)> seq_B = rand(2646070)> fft_A = fft(seq_A) #Long> 
fft_B = fft(seq_B)>zeropadded_fft_A = fft(seq_A, 
n=2**(ceil(log(len(seq_A),2))+1))>zeropadded_fft_B = fft(seq_B, 
n=2**(ceil(log(len(seq_B),2))+1))Ideally I need to calculate a favourable 
length to pad, prior to calling np.fft.fft() as Stéfan pointed out by 
calculating the factors. 
In [6]: from sympy import factorint 
In [7]: max(factorint(2646070)) Out[7]: 367 
In [8]: max(factorint(2649674)) Out[8]: 1324837 I will try adding some code to 
calculate the next power of two above my array length as you suggest;> And 
while you zero-pad, you can zero-pad to a sequence that is a power of two, thus 
preventing awkward factorizations.Does numpy have an easy way to do this, i.e. 
for a given number, find the next highest number (within a range) that could be 
factored into small, prime numbers as Phil explained? It would help if it gave 
a list, prioritised by number of factors.Thanks,Joseph
Date: Fri, 28 Aug 2015 17:26:32 -0700
From: stef...@berkeley.edu
To: numpy-discussion@scipy.org
Subject: Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples

On Aug 28, 2015 5:17 PM, "Pierre-Andre Noel"  
wrote:

>

> I had in mind the use of FFT to do convolutions (

> https://en.wikipedia.org/wiki/Convolution_theorem ). If you do not

> zero-pad properly, then the end of the signal may "bleed" on the

> beginning, and vice versa.
Ah, gotcha! All these things should also be handled nicely in 
scipy.signal.fftconvolve.
Stéfan 

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

Re: [Numpy-discussion] Cythonizing some of NumPy

[David Cournapeau]

On Tue, Sep 1, 2015 at 8:16 AM, Nathaniel Smith  wrote:

> On Sun, Aug 30, 2015 at 2:44 PM, David Cournapeau 
> wrote:
> > Hi there,
> >
> > Reading Nathaniel summary from the numpy dev meeting, it looks like
> there is
> > a consensus on using cython in numpy for the Python-C interfaces.
> >
> > This has been on my radar for a long time: that was one of my rationale
> for
> > splitting multiarray into multiple "independent" .c files half a decade
> ago.
> > I took the opportunity of EuroScipy sprints to look back into this, but
> > before looking more into it, I'd like to make sure I am not going astray:
> >
> > 1. The transition has to be gradual
>
> Yes, definitely.
>
> > 2. The obvious way I can think of allowing cython in multiarray is
> modifying
> > multiarray such as cython "owns" the PyMODINIT_FUNC and the module
> > PyModuleDef table.
>
> The seems like a plausible place to start.
>
> In the longer run, I think we'll need to figure out a strategy to have
> source code divided over multiple .pyx files (for the same reason we
> want multiple .c files -- it'll just be impossible to work with
> otherwise). And this will be difficult for annoying technical reasons,
> since we definitely do *not* want to increase the API surface exposed
> by multiarray.so, so we will need to compile these multiple .pyx and
> .c files into a single module, and have them talk to each other via
> internal interfaces. But Cython is currently very insistent that every
> .pyx file should be its own extension module, and the interface
> between different files should be via public APIs.
>
> I spent some time poking at this, and I think it's possible but will
> take a few kluges at least initially. IIRC the tricky points I noticed
> are:
>
> - For everything except the top-level .pyx file, we'd need to call the
> generated module initialization functions "by hand", and have a bit of
> utility code to let us access the symbol tables for the resulting
> modules
>
> - We'd need some preprocessor hack (or something?) to prevent the
> non-main module initialization functions from being exposed at the .so
> level (like 'cdef extern from "foo.h"', 'foo.h' re#defines
> PyMODINIT_FUNC to remove the visibility declaration)
>
> - By default 'cdef' functions are name-mangled, which is annoying if
> you want to be able to do direct C calls between different .pyx and .c
> files. You can fix this by adding a 'public' declaration to your cdef
> function. But 'public' also adds dllexport stuff which would need to
> be hacked out as per above.
>
> I think the best strategy for this is to do whatever horrible things
> are necessary to get an initial version working (on a branch, of
> course), and then once that's done assess what changes we want to ask
> the cython folks for to let us eliminate the gross parts.
>

Agreed.

Regarding multiple cython .pyx and symbol pollution, I think it would be
fine to have an internal API with the required prefix (say `_npy_cpy_`) in
a core library, and control the exported symbols at the .so level. This is
how many large libraries work in practice (e.g. MKL), and is a model well
understood by library users.

I will start the cythonize process without caring about any of that though:
one large .pyx file, and everything build together by putting everything in
one .so. That will avoid having to fight both cython and distutils at the
same time :)

David

>
> (Insisting on compiling everything into the same .so will probably
> also help at some point in avoiding Cython-Related Binary Size Blowup
> Syndrome (CRBSBS), because the masses of boilerplate could in
> principle be shared between the different files. I think some modern
> linkers are even clever enough to eliminate this kind of duplicate
> code automatically, since C++ suffers from a similar problem.)
>
> > 3. We start using cython for the parts that are mostly menial refcount
> work.
> > Things like functions in calculation.c are obvious candidates.
> >
> > Step 2 should not be disruptive, and does not look like a lot of work:
> there
> > are < 60 methods in the table, and most of them should be fairly
> > straightforward to cythonize. At worse, we could just keep them as is
> > outside cython and just "export" them in cython.
> >
> > Does that sound like an acceptable plan ?
> >
> > If so, I will start working on a PR to work on 2.
>
> Makes sense to me!
>
> -n
>
> --
> Nathaniel J. Smith -- http://vorpus.org
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] np.sign and object comparisons

[Jaime Fernández del Río]

On Mon, Aug 31, 2015 at 11:49 PM, Nathaniel Smith  wrote:

> On Sun, Aug 30, 2015 at 9:09 PM, Jaime Fernández del Río <
> jaime.f...@gmail.com> wrote:
>
>>
>> There are three ways of fixing this that I see:
>>
>>1. Arbitrarily choose a value to set the return to. This is
>>equivalent to choosing a default return for `cmp` for comparisons. This
>>preserves behavior, but feels wrong.
>>2. Similarly to how np.sign of a floating point array with nans
>>returns nan for those values, return e,g, None for these cases. This
>>is my preferred option.
>>3. Raise an error, along the lines of the TypeError: unorderable types
>>that 3.x produces for some comparisons.
>>
>> Having read the other replies so far -- given that no-one seems to have
> any clear intuition or use cases, I guess I find option 3 somewhat
> tempting... it keeps our options open until someone who actually cares
> comes along with a use case to hone our intuition on, and is very safe in
> the mean time.
>
> (This was noticed in the course of routine code cleanups, right, not an
> external bug report? For all we know right now, no actual user has ever
> even tried to apply np.sign to an object array?)
>

We do have a user that tried np.sign on an object array, and discovered
that our Py3K object comparison was crap:
https://github.com/numpy/numpy/issues/6229

No report of anyone trying np.sign on anything other than numbers that we
know of, though.

I'm starting to think that, given the lack of agreement, I thinking I am
going to agree with you that raising an error may be the better option,
because it's the least likely to break people's code if we later find we
need to change it.

Jaime

-- 
(\__/)
( 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@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples

[Oscar Benjamin]

On 1 September 2015 at 11:38, Joseph Codadeen  wrote:
>
>> And while you zero-pad, you can zero-pad to a sequence that is a power of
>> two, thus preventing awkward factorizations.
>
> Does numpy have an easy way to do this, i.e. for a given number, find the
> next highest number (within a range) that could be factored into small,
> prime numbers as Phil explained? It would help if it gave a list,
> prioritised by number of factors.

Just use the next power of 2. Pure powers of 2 are the most efficient
for FFT algorithms so it potentially works out better than finding a
smaller but similarly composite size to pad to. Finding the next power
of 2 is easy to code and never a bad choice.

To avoid the problems mentioned about zero-padding distorting the FFT
you can use the Bluestein transform as below. This pads up to a power
of two (greater than 2N-1) but does it in a special way that ensures
that it is still calculating the exact (up to floating point error)
same DFT:

from numpy import array, exp, pi, arange, concatenate
from numpy.fft import fft, ifft

def ceilpow2(N):
'''
>>> ceilpow2(15)
16
>>> ceilpow2(16)
16
'''
p = 1
while p < N:
p *= 2
return p

def fftbs(x):
'''
>>> data = [1, 2, 5, 2, 5, 2, 3]
>>> from numpy.fft import fft
>>> from numpy import allclose
>>> from numpy.random import randn
>>> for n in range(1, 1000):
... data = randn(n)
... assert allclose(fft(data), fftbs(data))
'''
N = len(x)
x = array(x)

n = arange(N)
b = exp((1j*pi*n**2)/N)
a = x * b.conjugate()

M = ceilpow2(N) * 2
A = concatenate((a, [0] * (M - N)))
B = concatenate((b, [0] * (M - 2*N + 1), b[:0:-1]))
C = ifft(fft(A) * fft(B))
c = C[:N]
return b.conjugate() * c

if __name__ == "__main__":
import doctest
doctest.testmod()


--
Oscar
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion