Re: [Numpy-discussion] Recognizing a cycle in a vector

[Oscar Benjamin]

On 3 December 2015 at 08:45, Manolo Martínez  wrote:
>> > >> Is there any way to check for cycles in this situation?
>> >
>> > > Fast fourier transform (fft)?
>> >
>> > +1 For using a discrete Fourier transform, as implemented by numpy.fft.fft.
>> > You mentioned that you sample at points which do not correspond with the
>> > period of the signal; this introduces a slight complexity in how the
>> > Fourier transform reflects information about the original signal. I attach
>> > two documents to this email with details about those (and other)
>> > complexities. There is also much information on this topic online and in
>> > signal processing books.
>>
>
> So, I thought I'd report back on what I've ended up doing. Given that
> the cycles I find in my data are usually very close to sine waves, the
> following works well enough:
>
>
> def periodic_vector(vector):
> """
> Take the FFT of a vector, and eliminate all components but the
> two main ones (i.e., the static and biggest sine amplitude) and
> compare the reconstructed wave with the original. Return true if
> close enough
> """
> rfft = np.fft.rfft(vector)
> magnitudes = np.abs(np.real(rfft))
> choice = magnitudes > sorted(magnitudes)[-3]
> newrfft = np.choose(choice, (np.zeros_like(rfft), rfft))
> newvector = np.fft.irfft(newrfft)
> return np.allclose(vector, newvector, atol=1e-2)
>
>
> This is doing the job for me at the moment, but there are, that I can
> see, a couple of things that could be improved (and surely more that I
> cannot see):
>
> 1) this func sorts the absolute value of the amplitudes to find the two
> most important  components, and this seems overkill for large vectors.
>
> 2) I'm running the inverse FFT, and comparing to the initial vector, but
> it should be possible to make a decision solely based on the size of
> terms in the FFT itself. I'm just not confident enough to design a test
> based on that.
>
> Anyway, thanks to those who pointed me in the right direction.

If what you have works out fine for you then feel free to ignore this but...

The more common way to find periodic orbits in ODEs is to pose the
question as a boundary-value problem (BVP) rather than seek orbital
patterns in the solution of an initial value problem (IVP). BVP
methods are more robust, can find unstable orbits, detect period
doubling etc. I would use the DFT method in the case that the ODEs are
of very high dimension and/or if the orbits in question are perhaps
quasi-periodic or if I can only really access the ODEs through the
output of an IVP solver. In the common case though the BVP method is
usually better. Something is written about finding periodic orbits
here:
http://www.scholarpedia.org/article/Periodic_orbit#Numerical_Methods_for_Finding_Periodic_Orbits

There are two basic ways to do this as a BVP problem: the shooting
method and the mesh. For your purposes the shooting method may suffice
so I'll briefly describe it.

Your ODEs must be of dimension at least 2 or you wouldn't have
periodic solutions. Consider x[n] to be the state vector at the nth
timestep. Suppose that x~ is part of an exact periodic orbit of the
ODE x' = f(x) with f(x) some vector field. Define P as the plane
containing x~ and normal to f(x~). The periodic orbit (if it exists)
must curve around and end up back at x~ pointing in the same
direction. For sinusoidal-ish orbits it will cross P twice, once at x~
and once somewhere else heading in something like the opposite
direction. If the orbit is more wiggly it may cross P more time but
always an even number of times before reaching x~.

Now suppose you have some guess x[0] which is close to a periodic
orbit. The true orbit should cross the plane P' generated by x[0],
f(x[0]) somewhere near x[0] pointing in approximately the same
direction. So use an ODE solver to iterate forward making x[1], x[2]
etc. until you cross the plane once and then twice coming back
somewhere near x[0]. Now you have x[n] and x[n+1] close-ish to x[0]
which lie on either side of the plane crossing in the same direction
as f(x[0]). You can now use the bisect method to take smaller and
larger timesteps from x[n] until your trajectory hits the plane
exactly. at this point your orbit is at some point x* which is on P'
near to x[0]. We now have an estimate of the period T of the orbit.

What I described in the last paragraph may be sufficient for your
purposes: if x* is sufficiently close to x[0] then you've found the
orbit and if not then it's not periodic. Usually though there is
another step: Define a function g(x, T) which takes a point x on the
plane P' and iterates the ODEs through a time T. You can put this into
a root-finder to solve g(x, T) = x for T and x. Since x is
N-dimensional we have N equations. However we have  constrained x to
lie on a plane so we have N-1 degrees of freedom in choosing x but we
also want to solve for 

Re: [Numpy-discussion] Recognizing a cycle in a vector

[Manolo Martínez]

Dear Oscar, 
 >
> > This is doing the job for me at the moment, but there are, that I can
> > see, a couple of things that could be improved (and surely more that I
> > cannot see):
 
> If what you have works out fine for you then feel free to ignore this but...
> [snip]

Talk about things I cannot see :) Thanks a lot for that very detailed
explanation. I will certainly look into settting my problem up as a BVP.

Incidentally, is there any modern textbook on numerical solving of ODEs
that you could recommend? 

Thanks again,
Manolo
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Recognizing a cycle in a vector

[Oscar Benjamin]

On 3 December 2015 at 11:58, Manolo Martínez  wrote:
>> > This is doing the job for me at the moment, but there are, that I can
>> > see, a couple of things that could be improved (and surely more that I
>> > cannot see):
>
>> If what you have works out fine for you then feel free to ignore this but...
>> [snip]
>
> Talk about things I cannot see :) Thanks a lot for that very detailed
> explanation. I will certainly look into settting my problem up as a BVP.
>
> Incidentally, is there any modern textbook on numerical solving of ODEs
> that you could recommend?

Not particularly. The shooting and mesh methods are described in
chapter 17 of the Numerical Recipes in C book which is relatively
accessible:
http://apps.nrbook.com/c/index.html

In terms of out of the box software I can recommend auto and xpp. Each
is esoteric and comes with a clunky interface. XPP has a strange GUI
and auto is controlled through Python bindings using IPython as
frontend.

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

Re: [Numpy-discussion] Recognizing a cycle in a vector

[Manolo Martínez]

> > >> Is there any way to check for cycles in this situation?
> > 
> > > Fast fourier transform (fft)?
> > 
> > +1 For using a discrete Fourier transform, as implemented by numpy.fft.fft.
> > You mentioned that you sample at points which do not correspond with the
> > period of the signal; this introduces a slight complexity in how the
> > Fourier transform reflects information about the original signal. I attach
> > two documents to this email with details about those (and other)
> > complexities. There is also much information on this topic online and in
> > signal processing books.
> 

So, I thought I'd report back on what I've ended up doing. Given that
the cycles I find in my data are usually very close to sine waves, the
following works well enough:


def periodic_vector(vector):
"""
Take the FFT of a vector, and eliminate all components but the
two main ones (i.e., the static and biggest sine amplitude) and
compare the reconstructed wave with the original. Return true if
close enough 
"""
rfft = np.fft.rfft(vector)
magnitudes = np.abs(np.real(rfft))
choice = magnitudes > sorted(magnitudes)[-3]
newrfft = np.choose(choice, (np.zeros_like(rfft), rfft))
newvector = np.fft.irfft(newrfft)
return np.allclose(vector, newvector, atol=1e-2)


This is doing the job for me at the moment, but there are, that I can
see, a couple of things that could be improved (and surely more that I
cannot see): 

1) this func sorts the absolute value of the amplitudes to find the two
most important  components, and this seems overkill for large vectors. 

2) I'm running the inverse FFT, and comparing to the initial vector, but
it should be possible to make a decision solely based on the size of
terms in the FFT itself. I'm just not confident enough to design a test
based on that.

Anyway, thanks to those who pointed me in the right direction.

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

Re: [Numpy-discussion] Recognizing a cycle in a vector

[Manolo Martínez]

On 12/03/15 at 05:39am, Eric Firing wrote:
> On 2015/12/02 10:45 PM, Manolo Martínez wrote:
> >1) this func sorts the absolute value of the amplitudes to find the two
> >most important  components, and this seems overkill for large vectors.
> 
> Try
> 
> inds = np.argpartition(-np.abs(ft), 2)[:2]
> 
> Now inds holds the indices of the two largest components.
> 

That is better than sorted indeed. Thanks,

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

Re: [Numpy-discussion] Recognizing a cycle in a vector

[Eric Firing]

On 2015/12/02 10:45 PM, Manolo Martínez wrote:

1) this func sorts the absolute value of the amplitudes to find the two
most important  components, and this seems overkill for large vectors.


Try

inds = np.argpartition(-np.abs(ft), 2)[:2]

Now inds holds the indices of the two largest components.

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

Re: [Numpy-discussion] When to stop supporting Python 2.6?

[Eric Firing]

On 2015/12/03 12:47 PM, Charles R Harris wrote:

Hi All,

Thought I would raise the topic apropos this post

. There is not a great advantage to dropping 2.6, OTOH, 2.7 has more
features (memoryview) and we could clean up the code a bit.

Along the same lines, dropping support for Python 3.2 would allow more
cleanup. In fact, I'd like to get to 3.4 as soon as possible, but don't
know what would be a reasonable schedule. The Python 3 series might be
easier to move forward on, as I think that Python 3 is just now starting
to become the dominant version in some areas.

Chuck



Chuck,

I would support dropping the old versions now.  As a related data point, 
matplotlib is testing master on 2.7, 3.4, and 3.5--no more 2.6 and 3.3.


Eric

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

[Numpy-discussion] When to stop supporting Python 2.6?

[Charles R Harris]

Hi All,

Thought I would raise the topic apropos this post

. There is not a great advantage to dropping 2.6, OTOH, 2.7 has more
features (memoryview) and we could clean up the code a bit.

Along the same lines, dropping support for Python 3.2 would allow more
cleanup. In fact, I'd like to get to 3.4 as soon as possible, but don't
know what would be a reasonable schedule. The Python 3 series might be
easier to move forward on, as I think that Python 3 is just now starting to
become the dominant version in some areas.

Chuck

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

Re: [Numpy-discussion] When to stop supporting Python 2.6?

[Bryan Van de Ven]

> On Dec 3, 2015, at 4:59 PM, Eric Firing  wrote:
> 
> Chuck,
> 
> I would support dropping the old versions now.  As a related data point, 
> matplotlib is testing master on 2.7, 3.4, and 3.5--no more 2.6 and 3.3.

Ditto for Bokeh. 
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] When to stop supporting Python 2.6?

[Jeff Reback]

pandas is going to drop
2.6 and 3.3 next release at end of Jan 

(3.2 dropped in 0.17, in October)



I can be reached on my cell 917-971-6387
> On Dec 3, 2015, at 6:00 PM, Bryan Van de Ven  wrote:
> 
> 
>> On Dec 3, 2015, at 4:59 PM, Eric Firing  wrote:
>> 
>> Chuck,
>> 
>> I would support dropping the old versions now.  As a related data point, 
>> matplotlib is testing master on 2.7, 3.4, and 3.5--no more 2.6 and 3.3.
> 
> Ditto for Bokeh. 
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] future of f2py and Fortran90+

[David Verelst]

f90wrap [1] extends the functionality of f2py, and can automatically
generate sensible wrappers for certain cases.
[1] https://github.com/jameskermode/f90wrap

On 15 July 2015 at 03:45, Sturla Molden  wrote:

> Eric Firing  wrote:
>
> > I'm curious: has anyone been looking into what it would take to enable
> > f2py to handle modern Fortran in general?  And into prospects for
> > getting such an effort funded?
>
> No need. Use Cython and Fortran 2003 ISO C bindings. That is the only
> portable way to interop between Fortran and C (including CPython) anyway.
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] f2py, numpy.distutils and multiple Fortran source files

[David Verelst]

Hi,

For the wafo [1] package we are trying to include the extension compilation
process in setup.py [2] by using setuptools and numpy.distutils [3]. Some
of the extensions have one Fortran interface source file, but it depends on
several other Fortran sources (modules). The manual compilation process
would go as follows:

gfortran -fPIC -c source_01.f
gfortran -fPIC -c source_02.f
f2py -m module_name -c source_01.o source_02.o source_interface.f

Can this workflow be incorporated into setuptools/numpy.distutils?
Something along the lines as:

from numpy.distutils.core import setup, Extension
ext = Extension('module.name',
   depends=['source_01.f', 'source_02.f'],
   sources=['source_interface.f'])

(note that the above does not work)

[1] https://github.com/wafo-project/pywafo
[2] https://github.com/wafo-project/pywafo/blob/pipinstall/setup.py
[3] http://docs.scipy.org/doc/numpy/reference/distutils.html

Regards,
David
​
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] future of f2py and Fortran90+

[Yuxiang Wang]

Too add to Sturla - I think this is what he mentioned but in more details:

http://www.fortran90.org/src/best-practices.html#interfacing-with-python

Shawn

On Tue, Jul 14, 2015 at 9:45 PM, Sturla Molden  wrote:
> Eric Firing  wrote:
>
>> I'm curious: has anyone been looking into what it would take to enable
>> f2py to handle modern Fortran in general?  And into prospects for
>> getting such an effort funded?
>
> No need. Use Cython and Fortran 2003 ISO C bindings. That is the only
> portable way to interop between Fortran and C (including CPython) anyway.
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion



-- 
Yuxiang "Shawn" Wang
Gerling Haptics Lab
University of Virginia
yw...@virginia.edu
+1 (434) 284-0836
https://sites.google.com/a/virginia.edu/yw5aj/
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] future of f2py and Fortran90+

[Eric Firing]

On 2015/12/03 11:08 AM, Yuxiang Wang wrote:

Too add to Sturla - I think this is what he mentioned but in more details:

http://www.fortran90.org/src/best-practices.html#interfacing-with-python


Right, but for each function that requires writing two wrappers, one in 
Fortran and a second one in cython.  Even though they are very simple, 
this would be cumbersome for a library with more than a few functions. 
Therefore I think there is still a place for f2py and f90wrap, and I am 
happy to see development continuing at least on the latter.


Eric



Shawn

On Tue, Jul 14, 2015 at 9:45 PM, Sturla Molden  wrote:

Eric Firing  wrote:


I'm curious: has anyone been looking into what it would take to enable
f2py to handle modern Fortran in general?  And into prospects for
getting such an effort funded?


No need. Use Cython and Fortran 2003 ISO C bindings. That is the only
portable way to interop between Fortran and C (including CPython) anyway.

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






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

Re: [Numpy-discussion] future of f2py and Fortran90+

[Tim Cera]

On Tue, Jul 14, 2015 at 10:13 PM Sturla Molden 
wrote:

> Eric Firing  wrote:
>
> > I'm curious: has anyone been looking into what it would take to enable
> > f2py to handle modern Fortran in general?  And into prospects for
> > getting such an effort funded?
>
> No need. Use Cython and Fortran 2003 ISO C bindings. That is the only
> portable way to interop between Fortran and C (including CPython) anyway.
>

For my wdmtoolbox I have a f2py wrapped Fortran 77 library.  Works great on
Linux, but because of the 'C' wrapper that f2py creates you run into the
same problem as Cython on Windows.

https://github.com/cython/cython/wiki/CythonExtensionsOnWindows

I guess if you first find a Windows machine, get it all setup, you may not
have to futz with it too much, but I just don't want to do it for such a
niche package.  I probably have a handful of users.

So I've been toying with the idea to use ctypes + iso_c_binding.  I could
then use MinGW on Windows to compile the code for all versions of Python
that have ctypes.  I've tested this approach on a few functions and it
works, but far from done.

My $0.02.

Kindest regards,
Tim
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] f2py, numpy.distutils and multiple Fortran source files

[Tim Cera]

On Thu, Dec 3, 2015 at 4:07 PM David Verelst 
wrote:

> Hi,
>
> For the wafo [1] package we are trying to include the extension
> compilation process in setup.py [2] by using setuptools and
> numpy.distutils [3]. Some of the extensions have one Fortran interface
> source file, but it depends on several other Fortran sources (modules). The
> manual compilation process would go as follows:
>
> gfortran -fPIC -c source_01.f
> gfortran -fPIC -c source_02.f
> f2py -m module_name -c source_01.o source_02.o source_interface.f
>
> Can this workflow be incorporated into setuptools/numpy.distutils?
> Something along the lines as:
>
> from numpy.distutils.core import setup, Extension
> ext = Extension('module.name',
>depends=['source_01.f', 'source_02.f'],
>sources=['source_interface.f'])
>
> (note that the above does not work)
>
> [1] https://github.com/wafo-project/pywafo
> [2] https://github.com/wafo-project/pywafo/blob/pipinstall/setup.py
> [3] http://docs.scipy.org/doc/numpy/reference/distutils.html
>
> Regards,
> David
> ​
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion


This might be helpful:
https://github.com/timcera/wdmtoolbox/blob/master/setup.py

Looks like I created the *.pyf file and include that in sources.  I think I
only used f2py to create the pyf file and not directly part of the
compilation process.  If memory serves the Extension function knows what to
do with the pyf file.

Kindest regards,
Tim
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion