-
Notifications
You must be signed in to change notification settings - Fork 4
Expand file tree
/
Copy pathhistogram2d.pyx
More file actions
97 lines (86 loc) · 3.33 KB
/
histogram2d.pyx
File metadata and controls
97 lines (86 loc) · 3.33 KB
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
# -*- coding: utf-8 -*-
#
# Licensed under the terms of the BSD 3-Clause
# (see plotpy/__init__.py for details)
"""2D-Histogram algorithm"""
cimport cython
cimport numpy as cnp
from libc.math cimport log
cnp.import_array()
cdef inline double double_max(double a, double b): return a if a >= b else b
cdef inline double double_min(double a, double b): return a if a <= b else b
@cython.profile(False)
@cython.boundscheck(False)
def histogram2d(cnp.ndarray[double, ndim=1] X, cnp.ndarray[double, ndim=1] Y,
double i0, double i1, double j0, double j1,
cnp.ndarray[double, ndim=2] data, logscale):
"""Compute 2-D Histogram from data X, Y"""
cdef double cx, cy, nmax, ix, iy
cdef unsigned int i
cdef Py_ssize_t n = X.shape[0]
cdef Py_ssize_t nx = data.shape[1]
cdef Py_ssize_t ny = data.shape[0]
cx = nx/(i1-i0)
cy = ny/(j1-j0)
for i in range(n):
# Centered bins => - .5
ix = (X[i] - i0) * cx - .5
iy = (Y[i] - j0) * cy - .5
if ix >= 0 and ix <= nx-1 and iy >= 0 and iy <= ny-1:
data[<int> iy, <int> ix] += 1
nmax = 0.
if logscale:
for j in range(ny):
for i in range(nx):
data[j, i] = log(1+data[j, i])
nmax = double_max(nmax, data[j, i])
else:
for j in range(ny):
for i in range(nx):
nmax = double_max(nmax, data[j, i])
return nmax
@cython.profile(False)
@cython.boundscheck(False)
def histogram2d_func(cnp.ndarray[double, ndim=1] X,
cnp.ndarray[double, ndim=1] Y,
cnp.ndarray[double, ndim=1] Z,
double i0, double i1, double j0, double j1,
cnp.ndarray[double, ndim=2] data_tmp,
cnp.ndarray[double, ndim=2] data, int computation):
"""Compute 2-D Histogram from data X, Y"""
cdef double cx, cy, nmax, ix, iy
cdef unsigned int i, u, v
cdef Py_ssize_t n = X.shape[0]
cdef Py_ssize_t nx = data.shape[1]
cdef Py_ssize_t ny = data.shape[0]
cx = nx/(i1-i0)
cy = ny/(j1-j0)
for i in range(n):
# Centered bins => - .5
ix = (X[i] - i0) * cx - .5
iy = (Y[i] - j0) * cy - .5
if ix >= 0 and ix <= nx-1 and iy >= 0 and iy <= ny-1:
u, v = <int> iy, <int> ix
if computation == 0: # max
data_tmp[u, v] += 1
data[u, v] = double_max(data[u, v], Z[i])
elif computation == 1: # min
data_tmp[u, v] += 1
data[u, v] = double_min(data[u, v], Z[i])
elif computation == 2: # sum
data_tmp[u, v] += 1
data[u, v] += Z[i]
elif computation == 3: # prod
data_tmp[u, v] += 1
data[u, v] *= Z[i]
elif computation == 4: # avg
data_tmp[u, v] += 1
data[u, v] += (Z[i]-data[u, v])/data_tmp[u, v]
elif computation == 5: # argmin
if data[u, v] > Z[i]:
data_tmp[u, v] = i
data[u, v] = Z[i]
elif computation == 6: # argmax
if data[u, v] < Z[i]:
data_tmp[u, v] = i
data[u, v] = Z[i]