ezietsman (owner)

Revisions

gist: 226473 Download_button fork
public
Description:
Python, Cython and Fortran f2py versions of a Deeming periodogram
Public Clone URL: git://gist.github.com/226473.git
Embed All Files: show embed
build_deeming.sh #
raw 
1
2
3

#!/bin/bash
f2py -c -m deeming periodogram.f90 -lgomp
f2py -c -m deemingomp periodogram.f90 --f90flags="-fopenmp " -lgomp
cyperiodogram.pyx #
raw 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

#
# Cython version of the Deeming periodogram function. 
# See deeming.py
#
 
 
import numpy as np
cimport numpy as np
cimport cython
 
cdef extern from 'math.h':
    double sin(double theta)
    double cos(double theta)
 
 
DTYPE = np.float
ctypedef np.float_t DTYPE_t
 
 
@cython.boundscheck(False)
def periodogram(np.ndarray[DTYPE_t, ndim=1]  t not None, \
                np.ndarray[DTYPE_t, ndim=1]  m not None,\
                np.ndarray[DTYPE_t, ndim=1]  freqs not None):
 
    cdef DTYPE_t realpart = 0.0
    cdef DTYPE_t imagpart = 0.0
    cdef DTYPE_t pi = np.pi
 
    cdef int Nf = freqs.shape[0]
    cdef int Nt = t.shape[0]
 
    cdef np.ndarray[DTYPE_t, ndim=1]  amps = np.zeros([Nf], dtype=DTYPE)
    cdef int i, j
 
    for i in range(Nf):
        for j in range(Nt):
            realpart = realpart + m[j]*cos(2.0*pi*freqs[i]*t[j])
            imagpart = imagpart + m[j]*sin(2.0*pi*freqs[i]*t[j])
 
        amps[i] = 2.0*(realpart**2 + imagpart**2)**0.5/Nt
 
        realpart = 0.0
        imagpart = 0.0
 
    return amps
 
 
deeming.py #
raw 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

#!/usr/bin/env python
 
 
# Testing the efficiency of different implimentations of the Deeming
# periodogram. Its an O(N*N) algorithm for calculating a periodogram but
# it works for unevenly sampled data. Used frequently in astronomical
# time-series.
#
# cyperiodogram contains a cython implimentation
#
# periodogram.f90 contains two Fortran 90 implimentations
#   1 : Using the same array syntax as numpy
#   2 : Hand coded loops (Much faster in this case)
#
# The fortran subroutines use openmp and are wrapped using f2py. Look in
# build_deeming.sh for the compilation commands.
#
# To run this example (in Linux):
#   sh build_deeming.sh
#   python deeming.py
#
# The cython extension will be compiled when running the deeming.py file.
#
#
#
#
# Ewald Zietsman
# 5 November 2009
#
# mangled = '101,119,97,108,100,46,122,105,101,116,115,109,97,110,64,103,109,97,105,108,46,99,111,109'
# print ''.join([chr(int(c)) for c in mangled.split(',')])
#
 
import numpy as np
import deeming
import deemingomp
import time
import pyximport; pyximport.install()
import cyperiodogram as cy
 
 
# Deeming periodogram, the slow way
def periodogrampy(t,m,f):
    pi = np.pi
    amps = np.zeros(f.size, dtype='float')
    cos = np.cos
    sin = np.sin
    for i, f in enumerate(freqs):
        real = (m*cos(2*pi*f*t)).sum()
        imag = (m*sin(2*pi*f*t)).sum()
        amps[i] = real**2 + imag**2
 
    amps = 2.0*np.sqrt(amps)/t.size
    return amps
 
 
 
 
if __name__ == "__main__":
    pi = np.pi
    # Data
    t = np.linspace(0.0,0.25,5000)
    m = 1.75*np.sin(2*pi*150*t) + 0.75*np.sin(2*pi*277*t) + np.sin(2*pi*333*t)
 
    # Frequencies
    freqs = np.linspace(0.0, 1000.0, 5000)
 
    #Run the different versions
    timenow = time.time()
    amps = periodogrampy(t, m, freqs)
    print 'Python version : ', time.time() - timenow
 
 
    timenow = time.time()
    cyamps = cy.periodogram(t, m, freqs)
    print 'cython version : ', time.time() - timenow
 
    timenow = time.time()
    ampsf90 = deeming.periodogram(t, m, freqs, t.size, freqs.size, 1)
    print "Fortran 90 version : ", time.time()-timenow
 
 
    timenow = time.time()
    ampsf90_2 = deeming.periodogram2(t, m, freqs, t.size, freqs.size, 1)
    print "Fortran 90 2nd version : ", time.time()-timenow
 
    timenow = time.time()
    ampsf90omp = deemingomp.periodogram(t, m, freqs, t.size, freqs.size, 2)
    print "Fortran 90 OpenMP version : ", time.time()-timenow
 
    timenow = time.time()
    ampsf90omp_2 = deemingomp.periodogram2(t, m, freqs, t.size, freqs.size, 2)
    print "Fortran 90 OpenMP 2nd version : ", time.time()-timenow
 
    print "py == cy     ", all(amps==cyamps)
    print "py == f90    ", all(amps==ampsf90)
    print "py == f90_2  ", all(amps==ampsf90_2)
    print "py == f90omp ", all(amps==ampsf90omp)
    print "py == f90omp_2 ", all(amps==ampsf90omp_2)
 
 
 
 
 
 
 
periodogram.f90 #
raw 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

!
!  See deeming.py for info. This file contains 2 Fortran implimentations
!  of the deeming periodogram algorithm
subroutine periodogram(time, value, freqs, nt, nf , numprocs, amps)
 
! periodogram calculates the Deeming periodogram of the datapoints contained
! in the arrays time, value at the frequencies contained in freqs
! This subroutine is meant for f2py wrapping as well as a testbed of openMP
! within a python environment
 
    use omp_lib
    implicit none
 
 
    ! sizes of time and frequency arrays
    integer :: nt, nf, id
    ! arrays
    double precision time(nt), value(nt)
    double precision freqs(nf), amps(nf)
    integer :: numprocs, f
    double precision :: pi, realpart, imagpart
 
    pi = 3.1415926535897931D0
    realpart = 0.0D0
    imagpart = 0.0D0
 
! Add compiler directives so f2py knows what to do
 
!f2py intent(in) time
!f2py intent(in) value
 
!f2py intent(in) nt
!f2py depend(nt) time
!f2py depend(nt) value
 
!f2py intent(in) freqs
!f2py intent(in) nf
!f2py depend(nf) freqs
 
!f2py intent(in) numprocs
!f2py intent(out) amps
 
    call OMP_SET_NUM_THREADS(numprocs)
 
    !$OMP PARALLEL DO &
    !$OMP DEFAULT(SHARED) PRIVATE(f,realpart, imagpart)
    do f = 1, nf
        ! calculate real and imaginary parts
        realpart = sum((value*dcos(2.0D0*pi*freqs(f)*time)))
        imagpart = sum((value*dsin(2.0D0*pi*freqs(f)*time)))
 
        ! calculate the amplitude
        amps(f) = 2.0D0*dsqrt(realpart**2 + imagpart**2)/nt
    end do
    !$OMP END PARALLEL DO
 
 
 
end subroutine
 
 
subroutine periodogram2(time, value, freqs, nt, nf , numprocs, amps)
 
! periodogram calculates the Deeming periodogram of the datapoints contained
! in the arrays time, value at the frequencies contained in freqs
! This subroutine is meant for f2py wrapping as well as a testbed of openMP
! within a python environment
 
    use omp_lib
    implicit none
 
 
    ! sizes of time and frequency arrays
    integer :: nt, nf, id
    ! arrays
    double precision time(nt), value(nt)
    double precision freqs(nf), amps(nf)
    integer :: numprocs, f, d
    double precision :: pi, realpart, imagpart
 
    pi = 3.1415926535897931D0
    realpart = 0.0D0
    imagpart = 0.0D0
 
! Add compiler directives so f2py knows what to do
 
!f2py threadsafe
!f2py intent(in) time
!f2py intent(in) value
 
!f2py intent(in) nt
!f2py depend(nt) time
!f2py depend(nt) value
 
!f2py intent(in) freqs
!f2py intent(in) nf
!f2py depend(nf) freqs
 
!f2py intent(in) numprocs
!f2py intent(out) amps
 
    call OMP_SET_NUM_THREADS(numprocs)
 
    !$OMP PARALLEL DO &
    !$OMP DEFAULT(SHARED) PRIVATE(f,realpart, imagpart)
    do f = 1, nf
        ! calculate real and imaginary parts
        do d = 1, nt
            realpart = realpart + value(d)*dcos(2.0D0*pi*freqs(f)*time(d))
            imagpart = imagpart + value(d)*dsin(2.0D0*pi*freqs(f)*time(d))
        end do
 
        ! calculate the amplitude
        amps(f) = 2.0D0*dsqrt(realpart*realpart + imagpart*imagpart)/nt
 
        realpart = 0.0D0
        imagpart = 0.0D0
    end do
    !$OMP END PARALLEL DO
 
end subroutine