Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Comparing changes

Choose two branches to see what's changed or to start a new pull request. If you need to, you can also compare across forks.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also compare across forks.
base fork: petigura/FFA
base: a3595b3193
...
head fork: petigura/FFA
compare: eb34a5c5af
Checking mergeability… Don't worry, you can still create the pull request.
  • 10 commits
  • 2 files changed
  • 0 commit comments
  • 1 contributor
Showing with 67 additions and 69 deletions.
  1. +43 −52 FFABench_cy.py
  2. +24 −17 FFA_cy.pyx
View
95 FFABench_cy.py
@@ -1,57 +1,48 @@
import numpy as np
from numpy import ma
import FFA_cy as FFA
-from matplotlib.pylab import *
-def seg(X,P0):
- XW = FFA.XWrap2(X,P0,pow2=True)
- M = XW.shape[0] # number of rows
-
- idCol = np.arange(P0,dtype=int) # id of each column
- idRow = np.arange(M,dtype=int) # id of each row
- P = P0 + idRow.astype(float) / (M - 1)
+def FFABench():
+ def seg(X,P0):
+ XW = FFA.XWrap2(X,P0,pow2=True)
+ M = XW.shape[0] # number of rows
+
+ idCol = np.arange(P0,dtype=int) # id of each column
+ idRow = np.arange(M,dtype=int) # id of each row
+ P = P0 + idRow.astype(float) / (M - 1)
+
+ XW.fill_value=0
+ data = XW.filled()
+ mask = (~XW.mask).astype(int)
+
+ sumF = FFA.FFA(data) # Sum of array elements folded on P0, P0 + i/(1-M)
+ countF = FFA.FFA(mask) # Number of valid data points
+ meanF = sumF/countF
+
+ names = ['mean','count','s2n','P']
+ dtype = zip(names,[float]*len(names) )
+ rep = np.empty(M,dtype=dtype)
+
+ # Take maximum epoch
+ idColMa = meanF.argmax(axis=1)
+ rep['mean'] = meanF[idRow,idColMa]
+ rep['count'] = countF[idRow,idColMa]
+ rep['s2n'] = rep['mean']*np.sqrt(rep['count'])
+ rep['P'] = P
+
+ return rep
+
+ X = np.load('pulse_train_data.npy')
+ Xmask = np.load('pulse_train_mask.npy')
+ X = ma.masked_array(X,Xmask,fill_value=0)
+
+ X = X[:10000] # Modify this to change execution time.
+
+ Pmin,Pmax = 250,2500
+ PGrid = np.arange(Pmin,Pmax)
+
+ func = lambda P0: seg(X,P0)
+ rep = map(func,PGrid)
+ rep = np.hstack(rep)
+ return rep
- XW.fill_value=0
- data = XW.filled()
- mask = (~XW.mask).astype(int)
-
- sumF = FFA.FFA(data) # Sum of array elements folded on P0, P0 + i/(1-M)
- countF = FFA.FFA(mask) # Number of valid data points
- meanF = sumF/countF
-
- names = ['mean','count','s2n','P']
- dtype = zip(names,[float]*len(names) )
- rep = np.empty(M,dtype=dtype)
-
- # Take maximum epoch
- idColMa = meanF.argmax(axis=1)
- rep['mean'] = meanF[idRow,idColMa]
- rep['count'] = countF[idRow,idColMa]
- rep['s2n'] = rep['mean']*np.sqrt(rep['count'])
- rep['P'] = P
-
- return rep
-
-X = np.load('pulse_train_data.npy')
-Xmask = np.load('pulse_train_mask.npy')
-X = ma.masked_array(X,Xmask,fill_value=0)
-
-X = X[:10000] # Modify this to change execution time.
-
-Pmin,Pmax = 250,2500
-PGrid = np.arange(Pmin,Pmax)
-
-func = lambda P0: seg(X,P0)
-rep = map(func,PGrid)
-rep = np.hstack(rep)
-
-fig,axL = subplots(nrows=2)
-axL[0].plot(X,label='Time Series')
-axL[1].plot(rep['P'],rep['s2n'],label='Periodogram')
-sca(axL[1])
-xlabel('Period')
-ylabel('S/N')
-legend()
-
-draw()
-show()
View
41 FFA_cy.pyx
@@ -5,6 +5,7 @@ cimport cython
cimport numpy as cnp
# cython: cdivision = True
+@cython.profile(True)
def FFA(XW):
"""
Fast Folding Algorithm
@@ -53,9 +54,10 @@ def FFA(XW):
XWFS = XW.copy()
for stage in range(1,nStage+1):
- XWFS = FFAShiftAdd(XWFS,stage)
+ XWFS = FFAShiftAdd(XWFS.astype(float),stage)
return XWFS
+@cython.profile(True)
def FFAButterfly(stage):
"""
FFA Butterfly
@@ -87,7 +89,9 @@ def FFAButterfly(stage):
@cython.wraparound(False)
@cython.boundscheck(False)
-def FFAGroupShiftAdd(group00,
+@cython.profile(True)
+@cython.cdivision(True)
+def FFAGroupShiftAdd(cnp.ndarray[cnp.float64_t, ndim=2,mode='c'] group0,
cnp.ndarray[cnp.int64_t, ndim=1] Arow,
cnp.ndarray[cnp.int64_t, ndim=1] Brow,
cnp.ndarray[cnp.int64_t, ndim=1] Bshft,
@@ -106,15 +110,11 @@ def FFAGroupShiftAdd(group00,
"""
maxShft = max(Bshft)
- group00 = np.hstack( [group00 , group00[:,: maxShft]] )
- cdef cnp.ndarray[cnp.float64_t, ndim=2] group0 = group00
-
-
cdef int iRow,iCol,iA,iB,Bs
- cdef int nRowGroup = group00.shape[0]
- cdef int nColGroup = group00.shape[1]
-
- cdef cnp.ndarray[cnp.float64_t, ndim=2] group
+ cdef int nRowGroup = group0.shape[0]
+ cdef int nColGroup = group0.shape[1]
+
+ cdef cnp.ndarray[cnp.float64_t, ndim=2] group = np.zeros((nRowGroup,nColGroup))
# Grow group by the maximum shift value
# Loop over rows in group
@@ -124,11 +124,14 @@ def FFAGroupShiftAdd(group00,
Bs = Bshft[iRow]
# Loop over the columns in the group
for iCol in range(nColGroup):
- group[iRow,iCol] = group0[iA,iCol] + group0[iB,Bs+iCol]
-
+ iBCol = (iCol + Bs + nColGroup) % nColGroup
+ group[iRow,iCol] += group0[iA,iCol]
+ group[iRow,iCol] += group0[iB,iBCol]
return group
-def FFAShiftAdd(XW0,stage):
+@cython.profile(True)
+def FFAShiftAdd(cnp.ndarray[cnp.float64_t, ndim=2] XW0,
+ int stage):
"""
FFA Shift and Add
@@ -155,16 +158,20 @@ def FFAShiftAdd(XW0,stage):
[ 0., 0., 1., 1.],
[ 0., 0., 2., 0.]])
"""
- nRow = XW0.shape[0]
+ cdef int nRow,nCol,nGroup,nRowGroup,iGroup,start,stop
+ cdef cnp.ndarray[cnp.float64_t, ndim=2] XW
+
+ nRow = XW0.shape[0]
+ nCol = XW0.shape[1]
nRowGroup = 2**stage
nGroup = nRow/nRowGroup
- XW = np.empty(XW0.shape)
+ XW = np.zeros((nRow,nCol),dtype=float)
+
Arow,Brow,Bshft = FFAButterfly(stage)
for iGroup in range(nGroup):
start = iGroup*nRowGroup
stop = (iGroup+1)*nRowGroup
- sG = slice(start,stop)
- XW[sG] = FFAGroupShiftAdd(XW0[sG],Arow,Brow,Bshft)
+ XW[start:stop] = FFAGroupShiftAdd(XW0[start:stop],Arow,Brow,Bshft)
return XW

No commit comments for this range

Something went wrong with that request. Please try again.