Coverage for src/lsqfitgp/_kernels/_spectral.py: 100%

1# lsqfitgp/_kernels/_spectral.py 

2# 

3# Copyright (c) 2023, Giacomo Petrillo 

4# 

5# This file is part of lsqfitgp. 

6# 

7# lsqfitgp is free software: you can redistribute it and/or modify 

8# it under the terms of the GNU General Public License as published by 

9# the Free Software Foundation, either version 3 of the License, or 

10# (at your option) any later version. 

11# 

12# lsqfitgp is distributed in the hope that it will be useful, 

13# but WITHOUT ANY WARRANTY; without even the implied warranty of 

14# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

15# GNU General Public License for more details. 

16# 

17# You should have received a copy of the GNU General Public License 

18# along with lsqfitgp. If not, see <http://www.gnu.org/licenses/>. 

19 

20from jax import numpy as jnp 6 ctx1feabcd

21 

22from .. import _special 6 ctx1feabcd

23from .._Kernel import stationarykernel 6 ctx1feabcd

24 

25@stationarykernel(derivable=True, maxdim=1) 6 ctx1feabcd

26def Cos(delta): 6 ctx1feabcd

27 """ 

28 Cosine kernel. 

29  

30 .. math:: 

31 k(\\Delta) = \\cos(\\Delta) 

32 = \\cos x \\cos y + \\sin x \\sin y 

33  

34 Samples from this kernel are harmonic functions. It can be multiplied with 

35 another kernel to introduce anticorrelations. 

36  

37 """ 

38 # TODO reference? 

39 return jnp.cos(delta) 5 ctx1eabcd

40 

41 

42@stationarykernel(maxdim=1, derivable=1, input='abs') 6 ctx1feabcd

43def Pink(delta, dw=1): 6 ctx1feabcd

44 """ 

45 Pink noise kernel. 

46  

47 .. math:: 

48 k(\\Delta) &= \\frac 1 {\\log(1 + \\delta\\omega)} 

49 \\int_1^{1+\\delta\\omega} \\mathrm d\\omega 

50 \\frac{\\cos(\\omega\\Delta)}\\omega = \\\\ 

51 &= \\frac { \\operatorname{Ci}(\\Delta (1 + \\delta\\omega)) 

52 - \\operatorname{Ci}(\\Delta) } 

53 {\\log(1 + \\delta\\omega)} 

54  

55 A process with power spectrum :math:`1/\\omega` truncated between 1 and 

56 :math:`1 + \\delta\\omega`. :math:`\\omega` is the angular frequency 

57 :math:`\\omega = 2\\pi f`. In the limit :math:`\\delta\\omega\\to\\infty` 

58 it becomes white noise. Derivable one time. 

59  

60 """ 

61 # TODO reference? 

62 

63 l = _special.ci(delta) 4 ctx1abcd

64 r = _special.ci(delta * (1 + dw)) 4 ctx1abcd

65 mean = delta * (1 + dw / 2) 4 ctx1abcd

66 norm = jnp.log1p(dw) 4 ctx1abcd

67 tol = jnp.sqrt(jnp.finfo(jnp.empty(0).dtype).eps) 4 ctx1abcd

68 

69 # TODO choose better this tolerance by estimating error terms 

70 

71 return jnp.where(delta * dw < tol, jnp.cos(mean), (r - l) / norm) 4 ctx1abcd

72 

73def _color_derivable(n=2): 6 ctx1feabcd

74 return n // 2 - 1 4 ctx1abcd

75 

76@stationarykernel(maxdim=1, derivable=_color_derivable, input='abs') 6 ctx1feabcd

77def Color(delta, n=2): 6 ctx1feabcd

78 """ 

79 Colored noise kernel. 

80  

81 .. math:: 

82 k(\\Delta) &= (n-1) \\Re E_n(-i\\Delta) = \\\\ 

83 &= (n-1) \\int_1^\\infty \\mathrm d\\omega 

84 \\frac{\\cos(\\omega\\Delta)}{\\omega^n}, 

85 \\quad n \\in \\mathbb N, n \\ge 2. 

86  

87 A process with power spectrum :math:`1/\\omega^n` truncated below 

88 :math:`\\omega = 1`. :math:`\\omega` is the angular frequency 

89 :math:`\\omega = 2\\pi f`. Derivable :math:`\\lfloor n/2 \\rfloor - 1` 

90 times. 

91  

92 """ 

93 # TODO reference? 

94 

95 # TODO real n > 1 => use usual series for small z and cont frac for large? 

96 

97 # TODO integration limit dw like Pink, use None to mean inf since the 

98 # derivability changes and I can not use values for conditions due to the 

99 # jit, delete Pink 

100 

101 # TODO The most general possible version I can imagine is int_a^b dw e^iwx 

102 # w^n with n in R. With n < -1 b=None is allowed, and a=None with n < 1, 

103 # meaning inf and 0. Since the general implementation would be slower, I can 

104 # detect n is an integer (at compile time) to use expn_imag. Afterward 

105 # delete Sinc and call this Power. 

106 # 

107 # This is the generalized incomplete gamma function with real argument and 

108 # imaginary parameters. DLMF says there is no known fixed-precision 

109 # implementation! 

110 # 

111 # int_l^h dw w^n e^-ixw = t = ixw 

112 # = (ix)^-(n+1) int_ixl^ixh dt t^n e^-t = 

113 # = (ix)^-(n+1) gammainc(n + 1, ixl, ixh) = z = ix, a = n+1 

114 # = z^-a gammainc(a, lz, hz) = 

115 # = z^-a (gamma(a, hz) - gamma(a, lz)) = 

116 # = z^-a (Gamma(a, lz) - Gamma(a, hz)) 

117 # 

118 # int_0^1 dw w^n e^-ixw = 

119 # = Gamma(n + 1) gamma*(n + 1, ix) 

120 # https://www.boost.org/doc/libs/1_79_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html boost impl of gamma(a,x) and Gamma(a,x) for real x 

121 # https://dlmf.nist.gov/8.2.E7 gamma*(a, z) 

122 # https://dlmf.nist.gov/8.8.E14 d/da gamma*(a, z) 

123 # https://dlmf.nist.gov/8.8.E15 d^n/dz^n z^-a gamma(a, z) 

124 # https://dlmf.nist.gov/8.9 continued fractions for gamma* and Gamma 

125 # https://dlmf.nist.gov/8.11.E2 asymp series for Gamma(a, z) 

126 # https://dlmf.nist.gov/8.12 series for a->∞ in terms of erfc (do I have complex erfc? => yes scipy, no jax) 

127 # https://dlmf.nist.gov/8.19.E17 continued fraction for E_n(z) 

128 # https://dlmf.nist.gov/8.20.E6 series of E_n(z) for large n 

129 # https://dlmf.nist.gov/8.27.i reference to Padé of gammainc for complex z 

130 # 683.f impl of E_n(z) for complex z (license?) (amos1990) 

131 # https://specialfunctions.juliamath.org/stable/functions_list/#SpecialFunctions.expint implementation in julia for arbitrary complex arguments 

132 # https://github.com/JuliaMath/SpecialFunctions.jl/blob/36c547b4a270b6089b1baf7bec05707e9fb8c7f9/src/expint.jl#L446-L503 

133 

134 # E_p(z) 

135 # following scipy's implementation, generalized to noninteger arguments: 

136 # if p > 50, use https://dlmf.nist.gov/8.20.ii (the dlmf reports it only 

137 # for real z, need to follow the references) 

138 # elif |z| < 1 and p integer, use https://dlmf.nist.gov/8.19.E8 

139 # p noninteger, use https://dlmf.nist.gov/8.19.E10 (need to 

140 # sum the external term accurately) 

141 # else (|z| >= 1), use https://dlmf.nist.gov/8.19.E17 

142 

143 # https://github.com/Radonirinaunimi/cmpx-spfunc incomplete gamma for 

144 # complex arguments (not well tested yet) 

145 

146 # Parametrize with n = 1 + 2 nu like Matérn. 

147 

148 # Bartosch, L. (2001). "Generation of colored noise". International Journal of Modern Physics C. 12 (6): 851–855. Bibcode:2001IJMPC..12..851B. doi:10.1142/S0129183101002012. S2CID 54500670. 

149 

150 assert int(n) == n and n >= 2, n 4 ctx1abcd

151 return (n - 1) * _special.expn_imag(n, delta).real 4 ctx1abcd

152 

153@stationarykernel(derivable=True, input='posabs', maxdim=1) 6 ctx1feabcd

154def Sinc(delta): 6 ctx1feabcd

155 """ 

156 Sinc kernel. 

157  

158 .. math:: k(\\Delta) = \\operatorname{sinc}(\\Delta) = 

159 \\frac{\\sin(\\pi\\Delta)}{\\pi\\Delta}. 

160  

161 Reference: Tobar (2019). 

162 """ 

163 return _special.sinc(delta) 4 ctx1abcd

164 

165 # TODO is this isotropic? My current guess is that it works up to some 

166 # dimension due to a coincidence but is not in general isotropic.