Mercurial > hg > python-allan

annotate _allan.pyx @ 0:3dcd4cfd8f82

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import
author Daniele Nicolodi <daniele.nicolodi@obspm.fr>
date Fri, 18 Apr 2014 18:44:54 +0200
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
1 import numpy as np
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
2 cimport numpy as np
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
3 cimport cython
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
4
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
5 @cython.boundscheck(False)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
6 @cython.wraparound(False)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
7 @cython.cdivision(True)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
8 def _xmvar(np.ndarray[np.double_t, ndim=1] x, int m):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
9
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
10 cdef double c
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
11 cdef double v
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
12 cdef int n = x.shape[0]
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
13 cdef unsigned int j
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
14 cdef unsigned int i
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
15
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
16 c = 0
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
17 for j in range(n - 3 * m + 1):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
18 v = 0
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
19 for i in range(j, j + m):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
20 v += (x[i + 2 * m] - 2 * x[i + m] + x[i]) / m
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
21 c += v * v
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
22 return c / (2.0 * m*m * (n - 3 * m + 1))
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
23
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
24
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
25 @cython.boundscheck(False)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
26 @cython.wraparound(False)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
27 @cython.cdivision(True)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
28 def _xmvar2(np.ndarray[np.double_t, ndim=1] x, np.ndarray[np.int_t, ndim=1] mm):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
29
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
30 cdef double z
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
31 cdef double zi
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
32 cdef unsigned int n = x.shape[0]
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
33 cdef unsigned int i
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
34 cdef unsigned int k
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
35 cdef unsigned int m
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
36 cdef np.ndarray[np.double_t, ndim=1] w = np.empty([n + 1, ], np.double)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
37 cdef np.ndarray[np.double_t, ndim=1] mvar = np.empty([mm.shape[0], ], np.double)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
38
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
39 mvar.fill(np.nan)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
40
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
41 w[0] = 0.0
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
42 for i in range(n):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
43 w[i + 1] = w[i] + x[i]
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
44
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
45 for k, m in enumerate(mm):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
46 if 3 * m > n:
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
47 continue
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
48 z = 0
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
49 for i in range(n - 3 * m + 1):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
50 zi = (w[i + 3 * m] - 3 * w[i + 2 * m] + 3 * w[i + m] - w[i]) / m
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
51 z += zi * zi
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
52 mvar[k] = z / (2.0 * m*m * (n - 3 * m + 1))
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
53
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
54 return mvar
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
55
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
56
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
57 @cython.boundscheck(False)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
58 @cython.wraparound(False)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
59 @cython.cdivision(True)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
60 def _xavar(double[:] x, int[:] mm):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
61
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
62 cdef double z
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
63 cdef double zi
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
64 cdef unsigned int n = x.shape[0]
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
65 cdef unsigned int i
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
66 cdef unsigned int k
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
67 cdef unsigned int m
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
68 cdef int N
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
69
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
70 r = np.empty([mm.shape[0], ], np.double)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
71 r.fill(np.nan)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
72 cdef double[::1] avar = r
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
73
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
74 for k, m in enumerate(mm):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
75 if m < 1:
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
76 continue
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
77 N = (n // m) - 2
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
78 if N < 1:
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
79 continue
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
80 z = 0
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
81 for i in range(N):
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
82 zi = x[i * m] - 2 * x[(i + 1) * m] + x[(i + 2) * m]
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
83 z += (zi * zi)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
84 avar[k] = z / (2.0 * m*m * N)
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
85
Daniele Nicolodi <daniele.nicolodi@obspm.fr>
parents:
diff changeset
86 return r