Coverage for physioblocks / utils / math_utils.py: 94%

1# SPDX-FileCopyrightText: Copyright INRIA 

2# 

3# SPDX-License-Identifier: LGPL-3.0-only 

4# 

5# Copyright INRIA 

6# 

7# This file is part of PhysioBlocks, a library mostly developed by the 

8# [Ananke project-team](https://team.inria.fr/ananke) at INRIA. 

9# 

10# Authors: 

11# - Colin Drieu 

12# - Dominique Chapelle 

13# - François Kimmig 

14# - Philippe Moireau 

15# 

16# PhysioBlocks is free software: you can redistribute it and/or modify it under the 

17# terms of the GNU Lesser General Public License as published by the Free Software 

18# Foundation, version 3 of the License. 

19# 

20# PhysioBlocks is distributed in the hope that it will be useful, but WITHOUT ANY 

21# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 

22# PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. 

23# 

24# You should have received a copy of the GNU Lesser General Public License along with 

25# PhysioBlocks. If not, see <https://www.gnu.org/licenses/>. 

26 

27"""Defines some math functions to help compute quantities""" 

28 

29from typing import Any 

30 

31import numpy as np 

32 

33 

34def exp_diff(a: float, b: float, diff: float) -> Any: 

35 """ 

36 Computes exponential differences with more numerical accuracy when a and b are close 

37 using the relationship: 

38 

39 .. math:: 

40 

41 e^a - e^b = 2\\ e^{\\frac{a + b}{2}}\\ \\text{sinh}(\\frac{a + b}{2}) 

42 

43 :param a: value of a 

44 :type a: np.float64 

45 

46 :param b: value of b 

47 :type b: np.float64 

48 

49 :param diff: value of a - b 

50 :type diff: np.float64 

51 

52 :return: exponential difference value 

53 :rtype: np.float64 

54 """ 

55 return np.exp(0.5 * (a + b)) * 2.0 * np.sinh(0.5 * diff) 

56 

57 

58def power_diff(a: float, b: float, diff: float, n: int) -> Any: 

59 """ 

60 Computes the a and b power n difference with more numerical accuracy when a and 

61 b are close by using the following relationship: 

62 

63 .. math:: 

64 

65 a^n - b^n = (a-b) * \\sum_{i = 0}^{n -1}{a^i b^{n-1-i}} 

66 

67 When n is negative, the following transformation is used: 

68 

69 .. math:: 

70 

71 \\frac{1}{a^n} - \\frac{1}{b^n} = -\\frac{a^n-b^n}{{ab}^n} 

72 

73 :param a: value of a 

74 :type a: np.float64 

75 

76 :param b: value of b 

77 :type b: np.float64 

78 

79 :param diff: value of a - b 

80 :type diff: np.float64 

81 

82 :return: the a and b power n difference value 

83 :rtype: np.float64 

84 """ 

85 

86 if n < 0: 

87 return -power_diff(a, b, diff, -n) / np.pow(a * b, -n) 

88 elif n > 0: 

89 return diff * _power_dec(a, b, n - 1, 0) 

90 else: 

91 return 0.0 

92 

93 

94def _power_dec(a: float, b: float, n1: int, n2: int) -> Any: 

95 """ 

96 Helper recursive function used by :func:`power_diff` to compute the term: 

97 

98 .. math:: 

99 

100 \\sum_{i = 0}^{n -1}{a^i b^{n-1-i}} 

101 """ 

102 

103 if n1 != 0 or n2 != 0: 

104 if n1 != 0: 

105 return np.pow(a, n1) * np.pow(b, n2) + _power_dec(a, b, n1 - 1, n2 + 1) 

106 else: 

107 return np.pow(b, n2) 

108 

109 return 1.0