Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/cardinal_pythonlib/maths_numpy.py : 27%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1#!/usr/bin/env python
2# cardinal_pythonlib/maths_numpy.py
4"""
5===============================================================================
7 Original code copyright (C) 2009-2021 Rudolf Cardinal (rudolf@pobox.com).
9 This file is part of cardinal_pythonlib.
11 Licensed under the Apache License, Version 2.0 (the "License");
12 you may not use this file except in compliance with the License.
13 You may obtain a copy of the License at
15 https://www.apache.org/licenses/LICENSE-2.0
17 Unless required by applicable law or agreed to in writing, software
18 distributed under the License is distributed on an "AS IS" BASIS,
19 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
20 See the License for the specific language governing permissions and
21 limitations under the License.
23===============================================================================
25**Miscellaneous mathematical functions that use Numpy** (which can be slow to
26load).
28"""
30# =============================================================================
31# Imports
32# =============================================================================
34from collections import Counter
35import random
36import sys
37from typing import List, Optional, Union
39import numpy as np # pip install numpy
41from cardinal_pythonlib.logs import get_brace_style_log_with_null_handler
43log = get_brace_style_log_with_null_handler(__name__)
46# =============================================================================
47# Constants
48# =============================================================================
50# sys.float_info.max_10_exp:
51# largest integer x such that 10 ** x is representable as a float
52# sys.float_info.max_exp:
53# largest integer x such that float(sys.float_info.radix) ** x is
54# representable as a float... and that's 2.0 ** x
55# But what we want is the largest integer x such that e ** x = math.exp(x)
56# is representable as a float, and that is:
58MAX_E_EXPONENT = int(np.log(sys.float_info.max)) # typically, 709
61# =============================================================================
62# Softmax
63# =============================================================================
65def softmax(x: np.ndarray,
66 b: float = 1.0) -> np.ndarray:
67 r"""
68 Standard softmax function:
70 .. math::
72 P_i = \frac {e ^ {\beta \cdot x_i}} { \sum_{i}{\beta \cdot x_i} }
74 Args:
75 x: vector (``numpy.array``) of values
76 b: exploration parameter :math:`\beta`, or inverse temperature
77 [Daw2009], or :math:`1/t`; see below
79 Returns:
80 vector of probabilities corresponding to the input values
82 where:
84 - :math:`t` is temperature (towards infinity: all actions equally likely;
85 towards zero: probability of action with highest value tends to 1)
86 - Temperature is not used directly as optimizers may take it to zero,
87 giving an infinity; use inverse temperature instead.
88 - [Daw2009] Daw ND, "Trial-by-trial data analysis using computational
89 methods", 2009/2011; in "Decision Making, Affect, and Learning: Attention
90 and Performance XXIII"; Delgado MR, Phelps EA, Robbins TW (eds),
91 Oxford University Press.
93 """
94 n = len(x)
95 if b == 0.0:
96 # e^0 / sum(a bunch of e^0) = 1 / n
97 return np.full(n, 1 / n)
99 constant = np.mean(x)
100 products = x * b - constant
101 # ... softmax is invariant to addition of a constant: Daw article and
102 # http://www.faqs.org/faqs/ai-faq/neural-nets/part2/section-12.html#b
103 # noinspection PyUnresolvedReferences
105 if products.max() > MAX_E_EXPONENT:
106 log.warning(
107 f"OVERFLOW in softmax(): x = {x}, b = {b}, constant = {constant}, "
108 f"x*b - constant = {products}")
109 # map the maximum to 1, other things to zero
110 index_of_max = np.argmax(products)
111 answer = np.zeros(n)
112 answer[index_of_max] = 1.0
113 return answer
115 # noinspection PyUnresolvedReferences
116 exponented = np.exp(products)
117 sum_exponented = np.sum(exponented)
118 if sum_exponented == 0.0:
119 # ... avoid division by zero
120 return np.full(n, 1 / n)
122 return exponented / np.sum(exponented)
125def pick_from_probabilities(probabilities: Union[List[float],
126 np.ndarray]) -> int:
127 """
128 Given a list of probabilities like ``[0.1, 0.3, 0.6]``, returns the index
129 of the probabilistically chosen item. In this example, we would return
130 ``0`` with probability 0.1; ``1`` with probability 0.3; and ``2`` with
131 probability 0.6.
133 Args:
134 probabilities: list of probabilities, which should sum to 1
136 Returns:
137 the index of the chosen option
139 Raises:
140 :exc:`ValueError` if a random number in the range [0, 1) is greater
141 than or equal to the cumulative sum of the supplied probabilities (i.e.
142 if you've specified probabilities adding up to less than 1)
144 Does not object if you supply e.g. ``[1, 1, 1]``; it'll always pick the
145 first in this example.
147 """
148 r = random.random() # range [0.0, 1.0), i.e. 0 <= r < 1
149 cs = np.cumsum(probabilities) # e.g. [0.1, 0.4, 1] in this example
150 for idx in range(len(cs)):
151 if r < cs[idx]:
152 return idx
153 raise ValueError(
154 f"Probabilities sum to <1: probabilities = {probabilities!r}, "
155 f"cumulative sum = {cs!r}"
156 )
159# =============================================================================
160# Logistic
161# =============================================================================
163def logistic(x: Union[float, np.ndarray],
164 k: float,
165 theta: float) -> Optional[float]:
166 r"""
167 Standard logistic function.
169 .. math::
171 y = \frac {1} {1 + e^{-k (x - \theta)}}
173 Args:
174 x: :math:`x`
175 k: :math:`k`
176 theta: :math:`\theta`
178 Returns:
179 :math:`y`
181 """
182 # https://www.sharelatex.com/learn/List_of_Greek_letters_and_math_symbols
183 if x is None or k is None or theta is None:
184 return None
185 # noinspection PyUnresolvedReferences
186 return 1 / (1 + np.exp(-k * (x - theta)))
189def inv_logistic(y: Union[float, np.ndarray],
190 k: float,
191 theta: float) -> Optional[float]:
192 r"""
193 Inverse standard logistic function:
195 .. math::
197 x = ( log( \frac {1} {y} - 1) / -k ) + \theta
199 Args:
200 y: :math:`y`
201 k: :math:`k`
202 theta: :math:`\theta`
204 Returns:
205 :math:`x`
207 """
208 if y is None or k is None or theta is None:
209 return None
210 # noinspection PyUnresolvedReferences
211 return (np.log((1 / y) - 1) / -k) + theta
214# =============================================================================
215# Testing
216# =============================================================================
218def _test_softmax() -> None:
219 """
220 Tests the :func`softmax` function.
221 """
222 arrays = [
223 [1, 1],
224 [1, 1, 2],
225 [1, 1, 1, 1, 1.01],
226 [1, 2, 3, 4, 5],
227 [1, 2, 3, 4, 1000],
228 [1, 2, 3, 4, 5.0 ** 400],
229 ]
230 betas = [
231 0, 0.5, 1, 2, 10, 100,
232 ]
233 for x in arrays:
234 for b in betas:
235 y = softmax(np.array(x), b=b)
236 print(f"softmax({x!r}, b={b}) -> {y}")
239def _test_pick_from_probabilities() -> None:
240 """
241 Tests the :func:`pick_from_probabilities` function.
242 """
243 probabilities = [
244 [0.0, 1.0],
245 [0.1, 0.3, 0.6],
246 [0.25, 0.25, 5],
247 # [0.25], # crashes (appropriately)
248 [1, 1, 1], # silly, but doesn't crash
249 ]
250 n_values = [10, 1000, 100000]
251 for p in probabilities:
252 for n in n_values:
253 c = Counter()
254 c.update(pick_from_probabilities(p) for _ in range(n))
255 sc = sorted(c.items(), key=lambda kv: kv[0])
256 print(f"_test_pick_from_probabilities: p = {p}, n = {n} -> {sc}")
259if __name__ == '__main__':
260 _test_softmax()
261 _test_pick_from_probabilities()