-
Notifications
You must be signed in to change notification settings - Fork 70
/
estimator.pyx
374 lines (309 loc) · 11.3 KB
/
estimator.pyx
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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
# distutils: language = c++
"""
This is the variogram estimater, implemented in cython.
"""
import numpy as np
from cython.parallel import parallel, prange
try:
cimport openmp
except:
...
cimport numpy as np
from libc.math cimport M_PI, acos, atan2, cos, fabs, isnan, pow, sin, sqrt
def set_num_threads(num_threads):
cdef int num_threads_c = 1
if num_threads is None:
try:
num_threads_c = openmp.omp_get_num_procs()
except:
...
else:
num_threads_c = num_threads
return num_threads_c
cdef inline double dist_euclid(
const int dim,
const double[:, :] pos,
const int i,
const int j,
) nogil:
cdef int d
cdef double dist_squared = 0.0
for d in range(dim):
dist_squared += ((pos[d, i] - pos[d, j]) * (pos[d, i] - pos[d, j]))
return sqrt(dist_squared)
cdef inline double dist_haversine(
const int dim,
const double[:, :] pos,
const int i,
const int j,
) nogil:
# pos holds lat-lon in deg
cdef double deg_2_rad = M_PI / 180.0
cdef double diff_lat = (pos[0, j] - pos[0, i]) * deg_2_rad
cdef double diff_lon = (pos[1, j] - pos[1, i]) * deg_2_rad
cdef double arg = (
pow(sin(diff_lat/2.0), 2) +
cos(pos[0, i]*deg_2_rad) *
cos(pos[0, j]*deg_2_rad) *
pow(sin(diff_lon/2.0), 2)
)
return 2.0 * atan2(sqrt(arg), sqrt(1.0-arg))
ctypedef double (*_dist_func)(
const int,
const double[:, :],
const int,
const int,
) nogil
cdef inline bint dir_test(
const int dim,
const double[:, :] pos,
const double dist,
const double[:, :] direction,
const double angles_tol,
const double bandwidth,
const int i,
const int j,
const int d,
) nogil:
cdef double s_prod = 0.0 # scalar product
cdef double b_dist = 0.0 # band-distance
cdef double tmp # temporary variable
cdef int k
cdef bint in_band = True
cdef bint in_angle = True
# scalar-product calculation for bandwidth projection and angle calculation
for k in range(dim):
s_prod += (pos[k, i] - pos[k, j]) * direction[d, k]
# calculate band-distance by projection of point-pair-vec to direction line
if bandwidth > 0.0:
for k in range(dim):
tmp = (pos[k, i] - pos[k, j]) - s_prod * direction[d, k]
b_dist += tmp * tmp
in_band = sqrt(b_dist) < bandwidth
# allow repeating points (dist = 0)
if dist > 0.0:
# use smallest angle by taking absolute value for arccos angle formula
tmp = fabs(s_prod) / dist
if tmp < 1.0: # else same direction (prevent numerical errors)
in_angle = acos(tmp) < angles_tol
return in_band and in_angle
cdef inline double estimator_matheron(const double f_diff) nogil:
return f_diff * f_diff
cdef inline double estimator_cressie(const double f_diff) nogil:
return sqrt(fabs(f_diff))
ctypedef double (*_estimator_func)(const double) nogil
cdef inline void normalization_matheron(
double[:] variogram,
long[:] counts,
):
cdef int i
for i in range(variogram.shape[0]):
# avoid division by zero
variogram[i] /= (2. * max(counts[i], 1))
cdef inline void normalization_cressie(
double[:] variogram,
long[:] counts,
):
cdef int i
cdef long cnt
for i in range(variogram.shape[0]):
# avoid division by zero
cnt = max(counts[i], 1)
variogram[i] = (
0.5 * (1./cnt * variogram[i])**4 /
(0.457 + 0.494 / cnt + 0.045 / cnt**2)
)
ctypedef void (*_normalization_func)(
double[:],
long[:],
)
cdef inline void normalization_matheron_vec(
double[:, :] variogram,
long[:, :] counts,
):
cdef int d
for d in range(variogram.shape[0]):
normalization_matheron(variogram[d, :], counts[d, :])
cdef inline void normalization_cressie_vec(
double[:, :] variogram,
long[:, :] counts,
):
cdef int d
for d in range(variogram.shape[0]):
normalization_cressie(variogram[d, :], counts[d, :])
ctypedef void (*_normalization_func_vec)(
double[:, :],
long[:, :],
)
cdef _estimator_func choose_estimator_func(str estimator_type):
cdef _estimator_func estimator_func
if estimator_type == 'm':
estimator_func = estimator_matheron
else: # estimator_type == 'c'
estimator_func = estimator_cressie
return estimator_func
cdef _normalization_func choose_estimator_normalization(str estimator_type):
cdef _normalization_func normalization_func
if estimator_type == 'm':
normalization_func = normalization_matheron
else: # estimator_type == 'c'
normalization_func = normalization_cressie
return normalization_func
cdef _normalization_func_vec choose_estimator_normalization_vec(str estimator_type):
cdef _normalization_func_vec normalization_func_vec
if estimator_type == 'm':
normalization_func_vec = normalization_matheron_vec
else: # estimator_type == 'c'
normalization_func_vec = normalization_cressie_vec
return normalization_func_vec
def directional(
const double[:, :] f,
const double[:] bin_edges,
const double[:, :] pos,
const double[:, :] direction, # should be normed
const double angles_tol=M_PI/8.0,
const double bandwidth=-1.0, # negative values to turn of bandwidth search
const bint separate_dirs=False, # whether the direction bands don't overlap
str estimator_type='m',
const int num_threads=1,
):
if pos.shape[1] != f.shape[1]:
raise ValueError(f'len(pos) = {pos.shape[1]} != len(f) = {f.shape[1])}')
if bin_edges.shape[0] < 2:
raise ValueError('len(bin_edges) too small')
if angles_tol <= 0:
raise ValueError('tolerance for angle search masks must be > 0')
cdef _estimator_func estimator_func = choose_estimator_func(estimator_type)
cdef _normalization_func_vec normalization_func_vec = (
choose_estimator_normalization_vec(estimator_type)
)
cdef int dim = pos.shape[0]
cdef int d_max = direction.shape[0]
cdef int i_max = bin_edges.shape[0] - 1
cdef int j_max = pos.shape[1] - 1
cdef int k_max = pos.shape[1]
cdef int f_max = f.shape[0]
cdef double[:, :] variogram = np.zeros((d_max, len(bin_edges)-1))
cdef long[:, :] counts = np.zeros((d_max, len(bin_edges)-1), dtype=long)
cdef int i, j, k, m, d
cdef double dist
cdef int num_threads_c = set_num_threads(num_threads)
for i in prange(i_max, nogil=True, num_threads=num_threads_c):
for j in range(j_max):
for k in range(j+1, k_max):
dist = dist_euclid(dim, pos, j, k)
if dist < bin_edges[i] or dist >= bin_edges[i+1]:
continue # skip if not in current bin
for d in range(d_max):
if not dir_test(
dim, pos, dist, direction, angles_tol, bandwidth, k, j, d
):
continue # skip if not in current direction
for m in range(f_max):
# skip no data values
if not (isnan(f[m, k]) or isnan(f[m, j])):
counts[d, i] += 1
variogram[d, i] += estimator_func(f[m, k] - f[m, j])
# once we found a fitting direction
# break the search if directions are separated
if separate_dirs:
break
normalization_func_vec(variogram, counts)
return np.asarray(variogram), np.asarray(counts)
def unstructured(
const double[:, :] f,
const double[:] bin_edges,
const double[:, :] pos,
str estimator_type='m',
str distance_type='e',
const int num_threads=1,
):
cdef int dim = pos.shape[0]
cdef _dist_func distance
if distance_type == 'e':
distance = dist_euclid
else:
distance = dist_haversine
if dim != 2:
raise ValueError(f'Haversine: dim = {dim} != 2')
if pos.shape[1] != f.shape[1]:
raise ValueError(f'len(pos) = {pos.shape[1]} != len(f) = {f.shape[1])}')
if bin_edges.shape[0] < 2:
raise ValueError('len(bin_edges) too small')
cdef _estimator_func estimator_func = choose_estimator_func(estimator_type)
cdef _normalization_func normalization_func = (
choose_estimator_normalization(estimator_type)
)
cdef int i_max = bin_edges.shape[0] - 1
cdef int j_max = pos.shape[1] - 1
cdef int k_max = pos.shape[1]
cdef int f_max = f.shape[0]
cdef double[:] variogram = np.zeros(len(bin_edges)-1)
cdef long[:] counts = np.zeros(len(bin_edges)-1, dtype=long)
cdef int i, j, k, m
cdef double dist
cdef int num_threads_c = set_num_threads(num_threads)
for i in prange(i_max, nogil=True, num_threads=num_threads_c):
for j in range(j_max):
for k in range(j+1, k_max):
dist = distance(dim, pos, j, k)
if dist < bin_edges[i] or dist >= bin_edges[i+1]:
continue # skip if not in current bin
for m in range(f_max):
# skip no data values
if not (isnan(f[m, k]) or isnan(f[m, j])):
counts[i] += 1
variogram[i] += estimator_func(f[m, k] - f[m, j])
normalization_func(variogram, counts)
return np.asarray(variogram), np.asarray(counts)
def structured(
const double[:, :] f,
str estimator_type='m',
const int num_threads=1,
):
cdef _estimator_func estimator_func = choose_estimator_func(estimator_type)
cdef _normalization_func normalization_func = (
choose_estimator_normalization(estimator_type)
)
cdef int i_max = f.shape[0] - 1
cdef int j_max = f.shape[1]
cdef int k_max = i_max + 1
cdef double[:] variogram = np.zeros(k_max)
cdef long[:] counts = np.zeros(k_max, dtype=long)
cdef int i, j, k
cdef int num_threads_c = set_num_threads(num_threads)
with nogil, parallel(num_threads=num_threads_c):
for i in range(i_max):
for j in range(j_max):
for k in prange(1, k_max-i):
counts[k] += 1
variogram[k] += estimator_func(f[i, j] - f[i+k, j])
normalization_func(variogram, counts)
return np.asarray(variogram)
def ma_structured(
const double[:, :] f,
const bint[:, :] mask,
str estimator_type='m',
const int num_threads=1,
):
cdef _estimator_func estimator_func = choose_estimator_func(estimator_type)
cdef _normalization_func normalization_func = (
choose_estimator_normalization(estimator_type)
)
cdef int i_max = f.shape[0] - 1
cdef int j_max = f.shape[1]
cdef int k_max = i_max + 1
cdef double[:] variogram = np.zeros(k_max)
cdef long[:] counts = np.zeros(k_max, dtype=long)
cdef int i, j, k
cdef int num_threads_c = set_num_threads(num_threads)
with nogil, parallel(num_threads=num_threads_c):
for i in range(i_max):
for j in range(j_max):
for k in prange(1, k_max-i):
if not mask[i, j] and not mask[i+k, j]:
counts[k] += 1
variogram[k] += estimator_func(f[i, j] - f[i+k, j])
normalization_func(variogram, counts)
return np.asarray(variogram)