-
Notifications
You must be signed in to change notification settings - Fork 0
/
spd.py
executable file
·332 lines (235 loc) · 10.8 KB
/
spd.py
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
# -*- coding: utf-8 -*-
# filename: spd.py
# Copyright 2017 Mario Gutiérrez-Roig <[email protected]>
# Copyright 2020 Stefano Costa <[email protected]>
#
# This file is part of IOSACal, the IOSA Radiocarbon Calibration Library.
# IOSACal is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# IOSACal is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with IOSACal. If not, see <http://www.gnu.org/licenses/>.
import numpy as np
from thanados.models.iosacal.core import CalibrationCurve, RadiocarbonDetermination
def spdsum(spdlist, norm=True):
""" Sums several SPDs stored in a list.
Arguments:
norm -- If 'True' final SPD will be normalized.
"""
spd = spdlist[0].copy()
maxi = np.max([np.max(np.array(x).T[0]) for x in spdlist])
mini = np.min([np.min(np.array(x).T[0]) for x in spdlist])
# Creates a void vector where perform the sum
xrange = np.arange(mini, maxi + 1, 1)
yrange = np.zeros((int)(maxi) + 1 - (int)(mini))
for d in spdlist:
# Reshapes every array by adding zeros at the head and tail
y = np.lib.pad(
d.T[1],
(np.int(d.T[0][0] - mini), np.int(maxi - d.T[0][-1])),
"constant",
constant_values=0,
)
# Summing over all histogram
yrange += y
# Normalizating the SPD in calendar scale
if norm == True:
yrange = yrange / np.sum(yrange)
spd.resize(len(xrange), 2, refcheck=False)
spd.T[0] = xrange
spd.T[1] = yrange
spd.ndates = np.sum([x.ndates for x in spdlist])
return spd
class SPD(np.ndarray):
""" A Sum of Probability Distributions of calibrated dates """
def __new__(cls, caldates, norm=True):
""" Initializes the SPD class from calibrated dates
Arguments:
caldates -- List with all calibrated dates (after R.calibrate()) in a list
norm -- If 'True' SPD will be normalized
"""
# Finds maximum and minimum value of calibrated dates
maxi = np.max([np.max(np.array(x).T[0]) for x in caldates])
mini = np.min([np.min(np.array(x).T[0]) for x in caldates])
# Creates a void vector where perform the sum
xrange = np.arange(mini, maxi + 1, 1)
yrange = np.zeros((int)(maxi) + 1 - (int)(mini))
for d in caldates:
# Disentangles and reverses the arrays
x = np.array(d).T[0][::-1]
y = np.array(d).T[1][::-1]
# Reshapes every array by adding zeros at the head and tail
yy = np.lib.pad(
y,
(np.int(x[0] - mini), np.int(maxi - x[-1])),
"constant",
constant_values=0,
)
# Summing over all histogram
yrange += yy
# Normalizating the SPD in calendar scale
if norm == True:
yrange = yrange / np.sum(yrange)
obj = np.asarray(np.column_stack((xrange, yrange))).view(cls)
obj.ndates = len(caldates)
return obj
def __array_finalize__(self, obj):
# see InfoArray.__array_finalize__ for comments
if obj is None:
return
self.ndates = getattr(obj, "ndates", None)
def normalize(self):
""" Normalizes the SPD """
y = np.array(self.T[1])
self.T[1] /= np.sum(y)
def taphcorr(self, function="Surovell", factor=1):
""" Performs a taphonomic correction over the SPD
Taphonomic Correction can be made according two different curves
indicated by the variable 'function': Surovell's or Williams' correction:
Surovell's correction:
n(t) = 5.726442 * 10^6 * (t + 2176.4)^-1.3925309
Surovell, T. A., Finley, J. B., Smith, G. M., Brantingham, P. J.,
& Kelly, R. (2009). Correcting temporal frequency distributions for
taphonomic bias. Journal of Archaeological Science, 36(8), 1715-1724.
Williams' correction:
n(t) = 2.107 * 10^7 * (t + 2754)^-1.526
Williams, A. N. (2012). The use of summed radiocarbon probability
distributions in archaeology: a review of methods.
Journal of Archaeological Science, 39(3), 578-589.
The taphonomic correction is not applied homogeneously, but taking into
account the percentage of open sites, rock shelters and caves samples
in each age. This percentage correction is given by the vector
"factor". If factor=None, then correction is applied homogeneously.
Arguments:
factor -- Percentage of correction between 0 and 1.
"""
corr = self.copy()
if function == "Surovell":
taph = 5726442 * np.power((self.T[0] + 2176.4), -1.3925309)
elif function == "Williams":
taph = 21070000 * np.power((self.T[0] + 2754), -1.526)
else:
print("ERROR: Invalid argument for 'function' variable")
sys.exit()
corr_y = [c / t if t > 0 else c for c, t in zip(np.array(self.T[1]), taph)]
corr_y = np.array(corr_y) * np.sum(self.T[1]) / np.sum(corr_y)
corr.T[1] = corr_y * factor + (1. - factor) * self.T[1]
return corr
def rollmean(self, m):
""" Performs a Rolling Mean of a 2*m+1 window (-m to +m)
Note:
- Returns a SPD of smaller size.
- Variable m should always be greater than 0
"""
a = self.copy()
n = 2 * m + 1
ret = np.cumsum(a.T[1], dtype=float)
ret[n:] = ret[n:] - ret[:-n]
a.T[1][m:-m] = ret[n - 1:] / n
return a[m:-m]
def simdates(self, N, calcurve_name, c14std, seed=0):
""" Simulated dates drawn from SPD curve
Performs a Bootstrap reseampling on the original SPD and generates a
new list of simulated dates after backcalibrating. The errors are
assigned at random from the original distribution.
Arguments:
N -- Number of dates to simulate
seed -- Seed for random generator (if any)
calcurve_name -- Calibration curve name for backcalibrating
c14std -- List with all original errors for generating the simulated
std.
"""
np.random.seed = seed
SPDcum = np.cumsum(self.T[1])
calcurve = CalibrationCurve(calcurve_name)
dates = []
for i in range(N):
ran = np.random.random()
idx = (np.abs(SPDcum - ran)).argmin()
date = self[idx][0]
idx2 = (np.abs(calcurve.T[0] - date)).argmin()
c14date = calcurve[idx2][1]
std = np.random.choice(c14std)
idname = "sim_" + str(i)
dates.append([idname, c14date, std])
return dates
def simSPD(self, calcurve_name, c14std, seed=0):
""" Simulated SPD curve generated from the original SPD
Performs a Bootstrap reseampling on the original SPD and generates a
new list of simulated dates that are backcalibrated in the 14C axis,
then those dates are calibrated again for constructing the SPD. The
number of simulated dates is the same than the original curve.
Arguments:
calcurve_name -- Calibration curve name for backcalibrating
c14std -- List with all original errors for generating the simulated
std.
seed -- Seed for random generator (if any)
"""
simdateslist = self.simdates(self.ndates, calcurve_name, c14std, seed)
simdates = [RadiocarbonDetermination(x[1], x[2], x[0]) for x in simdateslist]
# Calibrated date-by-date in order to save memory
caldates = []
for x in simdates:
cal_r = x.calibrate(calcurve_name, norm=True)
del cal_r.calibration_curve
del cal_r.radiocarbon_sample
del cal_r.intervals
del cal_r.median
caldates.append(cal_r)
simulatedSPD = SPD(caldates)
return simulatedSPD
class FreqHist(np.ndarray):
""" A Frequency Histogram of calibrated dates """
def __new__(cls, caldates, bins):
""" Initializes the FreqHist class from calibrated dates
Arguments:
caldates -- List with all calibrated dates (after R.calibrate()) in a list
binwidth -- Array with the bins
"""
meds = [x.median for x in caldates]
freqs, bins = np.histogram(meds, bins)
avbins = 0.5 * (bins[1:] + bins[:-1])
obj = np.asarray(np.column_stack((avbins, freqs))).view(cls)
return obj
def taphcorr(self, function="Surovell", factor=None):
""" Performs a taphonomic correction over the Frequency Histogram
Taphonomic Correction can be made according two different curves
indicated by the variable 'function': Surovell's or Williams' correction:
Surovell's correction:
n(t) = 5.726442 * 10^6 * (t + 2176.4)^-1.3925309
Surovell, T. A., Finley, J. B., Smith, G. M., Brantingham, P. J.,
& Kelly, R. (2009). Correcting temporal frequency distributions for
taphonomic bias. Journal of Archaeological Science, 36(8), 1715-1724.
Williams' correction:
n(t) = 2.107 * 10^7 * (t + 2754)^-1.526
Williams, A. N. (2012). The use of summed radiocarbon probability
distributions in archaeology: a review of methods.
Journal of Archaeological Science, 39(3), 578-589.
The taphonomic correction is not applied homogeneously, but taking into
account the percentage of open sites, rock shelters and caves samples
in each age. This percentage correction is given by the vector
"factor". If factor=None, then correction is applied homogeneously.
Arguments:
factor -- Percentage of correction in each bin.
"""
corr = self.copy()
if function == "Surovell":
taph = 5726442 * np.power((self.T[0] + 2176.4), -1.3925309)
elif function == "Williams":
taph = 21070000 * np.power((self.T[0] + 2754), -1.526)
else:
print("ERROR: Invalid argument for 'function' variable")
sys.exit()
if factor is None:
taph = np.array(taph)
else:
taph = np.array(taph) * factor
corr_y = [c / t if t > 0 else c for c, t in zip(np.array(self.T[1]), taph)]
corr_y = np.array(corr_y) * np.sum(self.T[1]) / np.sum(corr_y)
corr.T[1] = corr_y
return corr