Coverage for pesummary/utils/array.py: 87.1%
140 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-12-09 22:34 +0000
1# Licensed under an MIT style license -- see LICENSE.md
3import numpy as np
4from pesummary.utils.decorators import deprecation
6__author__ = ["Charlie Hoy <charlie.hoy@ligo.org>"]
9def _2DArray(input_array, likelihood=None, prior=None, weights=None):
10 """Helper function for initialising multiple Array objects.
12 Parameters
13 ----------
14 input_array: np.ndarray, list
15 input list/array
16 likelihood: list, optional
17 log likelihood samples to use when calculating maximum likelihood
18 prior: list, optional
19 log prior samples to use when calculating maximum posterior
20 weights: list, optional
21 weights to use for the samples
23 Returns
24 -------
25 array: list
26 list of Array objects of size len(input_array)
27 """
28 obj = np.atleast_2d(input_array)
29 if obj.ndim > 2:
30 raise ValueError("Input must be one- or two-dimensional")
31 if not obj.shape[-1]:
32 standard_deviation, minimum, maximum = [None], [None], [None]
33 else:
34 standard_deviation = np.std(obj, axis=1)
35 minimum = np.min(obj, axis=1)
36 maximum = np.max(obj, axis=1)
37 try:
38 maxL = Array._maxL(obj.T, likelihood=likelihood)
39 except Exception:
40 maxL = None
41 try:
42 maxP = Array._maxP(obj.T, log_likelihood=likelihood, log_prior=prior)
43 except Exception:
44 maxP = None
45 return [
46 Array(
47 _array, minimum=minimum[num], maximum=maximum[num],
48 standard_deviation=standard_deviation[num],
49 maxL=maxL[num] if maxL is not None else None,
50 maxP=maxP[num] if maxP is not None else None, weights=weights
51 ) for num, _array in enumerate(obj)
52 ]
55class Array(np.ndarray):
56 """Class to add extra functions and methods to np.ndarray
58 Parameters
59 ----------
60 input_aray: list/array
61 input list/array
63 Attributes
64 ----------
65 median: float
66 median of the input array
67 mean: float
68 mean of the input array
69 key_data: dict
70 dictionary containing the key data associated with the array
71 """
72 __slots__ = [
73 "standard_deviation", "minimum", "maximum", "maxL", "maxP", "weights"
74 ]
76 def __new__(
77 cls, input_array, standard_deviation=None, minimum=None, maximum=None,
78 maxL=None, maxP=None, likelihood=None, prior=None,
79 weights=None
80 ):
81 obj = np.asarray(input_array).view(cls)
82 obj.weights = weights
83 mapping = {
84 "standard_deviation": [standard_deviation, np.std, {}],
85 "minimum": [minimum, np.min, {}],
86 "maximum": [maximum, np.max, {}],
87 "maxL": [maxL, cls._maxL, {"likelihood": likelihood}],
88 "maxP": [
89 maxP, cls._maxP,
90 {"log_likelihood": likelihood, "log_prior": prior}
91 ],
92 }
93 for attr, item in mapping.items():
94 if item[0] is None:
95 try:
96 setattr(obj, attr, item[1](obj, **item[2]))
97 except Exception:
98 setattr(obj, attr, None)
99 else:
100 setattr(obj, attr, item[0])
101 return obj
103 def __reduce__(self):
104 pickled_state = super(Array, self).__reduce__()
105 new_state = pickled_state[2] + tuple(
106 [getattr(self, i) for i in self.__slots__]
107 )
108 return (pickled_state[0], pickled_state[1], new_state)
110 def __setstate__(self, state):
111 self.standard_deviation = state[-6]
112 self.minimum = state[-5]
113 self.maximum = state[-4]
114 self.maxL = state[-3]
115 self.maxP = state[-2]
116 self.weights = state[-1]
117 super(Array, self).__setstate__(state[0:-6])
119 def average(self, type="mean"):
120 """Return the average of the array
122 Parameters
123 ----------
124 type: str
125 the method to average the array
126 """
127 if type == "mean":
128 return self._mean(self, weights=self.weights)
129 elif type == "median":
130 return self._median(self, weights=self.weights)
131 else:
132 return None
134 @staticmethod
135 def _mean(array, weights=None):
136 """Compute the mean from a set of weighted samples
138 Parameters
139 ----------
140 array: np.ndarray
141 input array
142 weights: np.ndarray, optional
143 list of weights associated with each sample
144 """
145 if weights is None:
146 return np.mean(array)
147 weights = np.array(weights).flatten() / float(sum(weights))
148 return float(np.dot(np.array(array), weights))
150 @staticmethod
151 def _median(array, weights=None):
152 """Compute the median from a set of weighted samples
154 Parameters
155 ----------
156 array: np.ndarray
157 input array
158 weights: np.ndarray, optional
159 list of weights associated with each sample
160 """
161 from pesummary.utils.credible_interval import credible_interval
162 if weights is None:
163 return np.median(array)
164 return credible_interval(array, 50, weights=weights)
166 @staticmethod
167 def _maxL(array, likelihood=None):
168 """Return the maximum likelihood value of the array
170 Parameters
171 ----------
172 array: np.ndarray
173 input array
174 likelihood: np.ndarray, optional
175 likelihoods associated with each sample
176 """
177 if likelihood is not None:
178 likelihood = list(likelihood)
179 ind = likelihood.index(np.max(likelihood))
180 return array[ind]
181 return None
183 @staticmethod
184 def _maxP(array, log_likelihood=None, log_prior=None):
185 """Return the maximum posterior value of the array
187 Parameters
188 ----------
189 array: np.ndarray
190 input array
191 log_likelihood: np.ndarray, optional
192 log likelihoods associated with each sample
193 log_prior: np.ndarray, optional
194 log prior associated with each sample
195 """
196 if any(param is None for param in [log_likelihood, log_prior]):
197 return None
198 likelihood = np.array(log_likelihood)
199 prior = np.array(log_prior)
200 posterior = likelihood + prior
201 ind = np.argmax(posterior)
202 return array[ind]
204 def to_dtype(self, _dtype):
205 return _dtype(self)
207 @property
208 def key_data(self):
209 return self._key_data(self)
211 @staticmethod
212 def _key_data(
213 array, header=[
214 "mean", "median", "std", "maxL", "maxP", "5th percentile",
215 "95th percentile"
216 ]
217 ):
218 """Return a dictionary containing the key data associated with the
219 array
221 Parameters
222 ----------
223 array: np.ndarray
224 input array
225 header: list
226 list of properties you wish to return
227 """
228 def safe_dtype_change(array, _dtype):
229 if array is not None:
230 if isinstance(array, Array):
231 return array.to_dtype(_dtype)
232 else:
233 return _dtype(array)
234 return array
236 mydict = {}
237 for key in header:
238 if not hasattr(np.ndarray, key):
239 try:
240 _value = safe_dtype_change(getattr(array, key), float)
241 except AttributeError:
242 if key == "5th percentile":
243 _value = safe_dtype_change(
244 array.credible_interval(percentile=5), float
245 )
246 elif key == "95th percentile":
247 _value = safe_dtype_change(
248 array.credible_interval(percentile=95), float
249 )
250 else:
251 _value = safe_dtype_change(
252 array.average(type=key), float
253 )
254 else:
255 if key == "std":
256 _value = safe_dtype_change(array.standard_deviation, float)
257 else:
258 _value = safe_dtype_change(array.average(type=key), float)
259 mydict[key] = _value
260 return mydict
262 @staticmethod
263 @deprecation(
264 "The Array.percentile method will be deprecated in future releases. Please "
265 "use pesummary.utils.credible_interval.credible_interval instead"
266 )
267 def percentile(array, weights=None, percentile=None):
268 from pesummary.utils.credible_interval import credible_interval
269 return credible_interval(array, percentile, weights=weights)
271 @deprecation(
272 "The Array.confidence_interval method will be deprecated in future releases. "
273 "Please use Array.credible_interval instead"
274 )
275 def confidence_interval(self, percentile=None):
276 return self.credible_interval(percentile=percentile)
278 def credible_interval(self, percentile=None):
279 """Return the credible interval of the array
281 Parameters
282 ----------
283 percentile: int/list, optional
284 Percentile or sequence of percentiles to compute, which must be
285 between 0 and 100 inclusive
286 """
287 from pesummary.utils.credible_interval import credible_interval
288 if percentile is not None:
289 return credible_interval(self, percentile, weights=self.weights)
290 return self.credible_interval(percentile=[5, 95])
292 def two_sided_credible_interval(self, percentile=90):
293 """Return the two-sided credible interval of the array
295 Parameters
296 ----------
297 percentile: float, optional
298 Percentile to compute. Must be between 0 and 100 inclusive. Default 90
299 """
300 from pesummary.utils.credible_interval import two_sided_credible_interval
301 return two_sided_credible_interval(self, percentile, weights=self.weights)
303 def __array_finalize__(self, obj):
304 if obj is None:
305 return
306 self.standard_deviation = getattr(obj, 'standard_deviation', None)
307 self.minimum = getattr(obj, 'minimum', None)
308 self.maximum = getattr(obj, 'maximum', None)
309 self.maxL = getattr(obj, 'maxL', None)
310 self.maxP = getattr(obj, 'maxP', None)
311 self.weights = getattr(obj, 'weights', None)
313 def to_numpy(self):
314 """Convert Array object to a numpy.ndarray
316 Returns
317 -------
318 data: np.ndarray/tuple
319 return stored data as a np.ndarray. If weights are stored in the
320 Array object, return a tuple containing the stored data as a
321 np.ndarray and the weights as a np.ndarray
322 """
323 _array = np.asarray(self)
324 if self.weights is None:
325 return _array
326 return _array, np.asarray(self.weights)