-
Notifications
You must be signed in to change notification settings - Fork 0
/
hpd.py
executable file
·138 lines (110 loc) · 4.61 KB
/
hpd.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
#! /usr/bin/env python
# -*- coding: utf-8 -*-
# filename: hpd.py
# Copyright 2008-2009 Stefano Costa <[email protected]>
# Copyright 2008 David Laban <[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/>.
from collections import namedtuple
from copy import copy
from numpy import asarray
def findsorted(n, array):
'''Return sorted array and index of n inside array.'''
a = asarray(array)
a.sort()
i = a.searchsorted(n)
return a, i
def prev(n, array):
'''Find interval between n and its previous, inside array.'''
a,i = findsorted(n, array)
if i-1 < 0:
prev = None
else:
prev = a[i-1]
return prev
def next(n, array):
'''Find interval between n and its next, inside array.'''
a,i = findsorted(n, array)
try:
next = a[i+1]
except IndexError:
next = None
return next
def alsuren_hpd(calibrated_age, alpha):
'''Return year spans that have the required Highest Probability Density.'''
hpd_curve = calibrated_age.copy()
# sort rows by second column in inverse order
hpd_sorted = hpd_curve[hpd_curve[:,1].argsort(),][::-1]
hpd_cumsum = hpd_sorted[:,1].cumsum()
# normalised values
hpd_cumsum /= hpd_cumsum[-1]
threshold_index = hpd_cumsum.searchsorted(1 - alpha)
threshold_p = hpd_sorted[threshold_index][1]
threshold_index = calibrated_age[:,1] > threshold_p
hpd = list(hpd_curve[threshold_index,0])
confidence_intervals = list()
for i in hpd:
# ^ is the XOR operator
if (prev(i,hpd_curve[:,0]) not in hpd) ^ (next(i,hpd_curve[:,0]) not in hpd):
confidence_intervals.append(i)
return asarray(confidence_intervals).reshape(len(confidence_intervals)//2,2)
def confidence_percent(years, calibrated_age):
'''Return HPD as percent value for a given span of years.'''
percent_curve = calibrated_age.copy()
percent_curve[:,1] /= percent_curve[:,1].sum()
percent_sorted = percent_curve[percent_curve[:,0].argsort(),]
year1_index = percent_sorted[:,0].searchsorted([years[0]])
year2_index = percent_sorted[:,0].searchsorted([years[1]])
indices = [ percent_sorted[:,0].searchsorted([year]) for year in years ]
indices.sort()
min_year = indices[0][0]
max_year = indices[1][0]
confidence_interval = percent_sorted[min_year:max_year+1,1]
percent_result = confidence_interval.sum()
return float(percent_result)
class ConfIntv(namedtuple('ConfIntv', ['from_year', 'to_year', 'conf_perc'])):
__slots__ = ()
def __format__(self, fmt='bp'):
def ad_bc_prefix(year, prefixes='ad'):
'''Return a string with BC/AD prefix and the given year.'''
if prefixes == 'ad':
neg, pos = ('BC', 'AD')
elif prefixes == 'ce':
neg, pos = ('BCE', 'CE')
if year < 0:
yearf = '{0} {1:.0f}'.format(neg, abs(year))
else:
yearf = '{0} {1:.0f}'.format(pos, year)
return yearf
fmt = fmt if len(fmt) > 0 else 'bp' # default fmt is ''
if fmt == 'bp':
f = '{from_year:.0f} CalBP - {to_year:.0f} CalBP ({conf_perc:2.1%})'.format(**self._asdict())
elif fmt == 'ad' or 'ce':
ad = {'from_year': ad_bc_prefix(1950 - self.from_year, fmt),
'to_year': ad_bc_prefix(1950 - self.to_year, fmt),
'conf_perc': self.conf_perc}
f = '{from_year} - {to_year} ({conf_perc:2.1%})'.format(**ad)
return f
__str__ = __format__
class ConfIntvList(list):
def __format__(self, fmt):
return '\n'.join(ci.__format__(fmt) for ci in self)
__str__ = __format__
def hpd_interval(calibrated_age, alpha):
'''Wrapper around other functions, returns a single object.'''
res = ConfIntvList()
intervals = alsuren_hpd(calibrated_age, alpha)
for interval in intervals:
percent = confidence_percent(interval, calibrated_age)
res.append(ConfIntv(interval[0], interval[1], percent))
return res