Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/tensorflow/python/ops/distributions/special_math.py: 19%

110 statements  

« prev     ^ index     » next       coverage.py v7.4.0, created at 2024-01-03 07:57 +0000

1# Copyright 2016 The TensorFlow Authors. All Rights Reserved. 

2# 

3# Licensed under the Apache License, Version 2.0 (the "License"); 

4# you may not use this file except in compliance with the License. 

5# You may obtain a copy of the License at 

6# 

7# http://www.apache.org/licenses/LICENSE-2.0 

8# 

9# Unless required by applicable law or agreed to in writing, software 

10# distributed under the License is distributed on an "AS IS" BASIS, 

11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 

12# See the License for the specific language governing permissions and 

13# limitations under the License. 

14# ============================================================================== 

15 

16# Functions "ndtr" and "ndtri" are derived from calculations made in: 

17# https://root.cern.ch/doc/v608/SpecFuncCephesInv_8cxx_source.html 

18# In the following email exchange, the author gives his consent to redistribute 

19# derived works under an Apache 2.0 license. 

20# 

21# From: Stephen Moshier <steve@moshier.net> 

22# Date: Sat, Jun 9, 2018 at 2:36 PM 

23# Subject: Re: Licensing cephes under Apache (BSD-like) license. 

24# To: rif <rif@google.com> 

25# 

26# 

27# 

28# Hello Rif, 

29# 

30# Yes, Google may distribute Cephes files under the Apache 2 license. 

31# 

32# If clarification is needed, I do not favor BSD over other free licenses. 

33# I would agree that Apache 2 seems to cover the concern you mentioned 

34# about sublicensees. 

35# 

36# Best wishes for good luck with your projects! 

37# Steve Moshier 

38# 

39# 

40# 

41# On Thu, 31 May 2018, rif wrote: 

42# 

43# > Hello Steve. 

44# > My name is Rif. I work on machine learning software at Google. 

45# > 

46# > Your cephes software continues to be incredibly useful and widely used. I 

47# > was wondering whether it would be permissible for us to use the Cephes code 

48# > under the Apache 2.0 license, which is extremely similar in permissions to 

49# > the BSD license (Wikipedia comparisons). This would be quite helpful to us 

50# > in terms of avoiding multiple licenses on software. 

51# > 

52# > I'm sorry to bother you with this (I can imagine you're sick of hearing 

53# > about this by now), but I want to be absolutely clear we're on the level and 

54# > not misusing your important software. In former conversation with Eugene 

55# > Brevdo (ebrevdo@google.com), you wrote "If your licensing is similar to BSD, 

56# > the formal way that has been handled is simply to add a statement to the 

57# > effect that you are incorporating the Cephes software by permission of the 

58# > author." I wanted to confirm that (a) we could use the Apache license, (b) 

59# > that we don't need to (and probably you don't want to) keep getting 

60# > contacted about individual uses, because your intent is generally to allow 

61# > this software to be reused under "BSD-like" license, and (c) you're OK 

62# > letting incorporators decide whether a license is sufficiently BSD-like? 

63# > 

64# > Best, 

65# > 

66# > rif 

67# > 

68# > 

69# > 

70 

71"""Special Math Ops.""" 

72 

73import numpy as np 

74 

75from tensorflow.python.framework import constant_op 

76from tensorflow.python.framework import ops 

77from tensorflow.python.ops import array_ops 

78from tensorflow.python.ops import math_ops 

79 

80__all__ = [ 

81 "erfinv", 

82 "ndtr", 

83 "ndtri", 

84 "log_ndtr", 

85 "log_cdf_laplace", 

86] 

87 

88 

89# log_ndtr uses different functions over the ranges 

90# (-infty, lower](lower, upper](upper, infty) 

91# Lower bound values were chosen by examining where the support of ndtr 

92# appears to be zero, relative to scipy's (which is always 64bit). They were 

93# then made more conservative just to be safe. (Conservative means use the 

94# expansion more than we probably need to.) See `NdtrTest` in 

95# special_math_test.py. 

96LOGNDTR_FLOAT64_LOWER = np.array(-20, np.float64) 

97LOGNDTR_FLOAT32_LOWER = np.array(-10, np.float32) 

98 

99# Upper bound values were chosen by examining for which values of 'x' 

100# Log[cdf(x)] is 0, after which point we need to use the approximation 

101# Log[cdf(x)] = Log[1 - cdf(-x)] approx -cdf(-x). We chose a value slightly 

102# conservative, meaning we use the approximation earlier than needed. 

103LOGNDTR_FLOAT64_UPPER = np.array(8, np.float64) 

104LOGNDTR_FLOAT32_UPPER = np.array(5, np.float32) 

105 

106 

107def ndtr(x, name="ndtr"): 

108 """Normal distribution function. 

109 

110 Returns the area under the Gaussian probability density function, integrated 

111 from minus infinity to x: 

112 

113 ``` 

114 1 / x 

115 ndtr(x) = ---------- | exp(-0.5 t**2) dt 

116 sqrt(2 pi) /-inf 

117 

118 = 0.5 (1 + erf(x / sqrt(2))) 

119 = 0.5 erfc(x / sqrt(2)) 

120 ``` 

121 

122 Args: 

123 x: `Tensor` of type `float32`, `float64`. 

124 name: Python string. A name for the operation (default="ndtr"). 

125 

126 Returns: 

127 ndtr: `Tensor` with `dtype=x.dtype`. 

128 

129 Raises: 

130 TypeError: if `x` is not floating-type. 

131 """ 

132 

133 with ops.name_scope(name, values=[x]): 

134 x = ops.convert_to_tensor(x, name="x") 

135 if x.dtype.as_numpy_dtype not in [np.float32, np.float64]: 

136 raise TypeError( 

137 "x.dtype=%s is not handled, see docstring for supported types." 

138 % x.dtype) 

139 return _ndtr(x) 

140 

141 

142def _ndtr(x): 

143 """Implements ndtr core logic.""" 

144 half_sqrt_2 = constant_op.constant( 

145 0.5 * np.sqrt(2.), dtype=x.dtype, name="half_sqrt_2") 

146 w = x * half_sqrt_2 

147 z = math_ops.abs(w) 

148 y = array_ops.where_v2( 

149 math_ops.less(z, half_sqrt_2), 1. + math_ops.erf(w), 

150 array_ops.where_v2( 

151 math_ops.greater(w, 0.), 2. - math_ops.erfc(z), math_ops.erfc(z))) 

152 return 0.5 * y 

153 

154 

155def ndtri(p, name="ndtri"): 

156 """The inverse of the CDF of the Normal distribution function. 

157 

158 Returns x such that the area under the pdf from minus infinity to x is equal 

159 to p. 

160 

161 A piece-wise rational approximation is done for the function. 

162 This is a port of the implementation in netlib. 

163 

164 Args: 

165 p: `Tensor` of type `float32`, `float64`. 

166 name: Python string. A name for the operation (default="ndtri"). 

167 

168 Returns: 

169 x: `Tensor` with `dtype=p.dtype`. 

170 

171 Raises: 

172 TypeError: if `p` is not floating-type. 

173 """ 

174 

175 with ops.name_scope(name, values=[p]): 

176 p = ops.convert_to_tensor(p, name="p") 

177 if p.dtype.as_numpy_dtype not in [np.float32, np.float64]: 

178 raise TypeError( 

179 "p.dtype=%s is not handled, see docstring for supported types." 

180 % p.dtype) 

181 return _ndtri(p) 

182 

183 

184def _ndtri(p): 

185 """Implements ndtri core logic.""" 

186 

187 # Constants used in piece-wise rational approximations. Taken from the cephes 

188 # library: 

189 # https://root.cern.ch/doc/v608/SpecFuncCephesInv_8cxx_source.html 

190 

191 p0 = [ 

192 -1.23916583867381258016E0, 1.39312609387279679503E1, 

193 -5.66762857469070293439E1, 9.80010754185999661536E1, 

194 -5.99633501014107895267E1 

195 ] 

196 q0 = [ 

197 -1.18331621121330003142E0, 1.59056225126211695515E1, 

198 -8.20372256168333339912E1, 2.00260212380060660359E2, 

199 -2.25462687854119370527E2, 8.63602421390890590575E1, 

200 4.67627912898881538453E0, 1.95448858338141759834E0, 1.0 

201 ] 

202 p1 = [ 

203 -8.57456785154685413611E-4, -3.50424626827848203418E-2, 

204 -1.40256079171354495875E-1, 2.18663306850790267539E0, 

205 1.46849561928858024014E1, 4.40805073893200834700E1, 

206 5.71628192246421288162E1, 3.15251094599893866154E1, 

207 4.05544892305962419923E0 

208 ] 

209 q1 = [ 

210 -9.33259480895457427372E-4, -3.80806407691578277194E-2, 

211 -1.42182922854787788574E-1, 2.50464946208309415979E0, 

212 1.50425385692907503408E1, 4.13172038254672030440E1, 

213 4.53907635128879210584E1, 1.57799883256466749731E1, 1.0 

214 ] 

215 p2 = [ 

216 6.23974539184983293730E-9, 2.65806974686737550832E-6, 

217 3.01581553508235416007E-4, 1.23716634817820021358E-2, 

218 2.01485389549179081538E-1, 1.33303460815807542389E0, 

219 3.93881025292474443415E0, 6.91522889068984211695E0, 

220 3.23774891776946035970E0 

221 ] 

222 q2 = [ 

223 6.79019408009981274425E-9, 2.89247864745380683936E-6, 

224 3.28014464682127739104E-4, 1.34204006088543189037E-2, 

225 2.16236993594496635890E-1, 1.37702099489081330271E0, 

226 3.67983563856160859403E0, 6.02427039364742014255E0, 1.0 

227 ] 

228 

229 def _create_polynomial(var, coeffs): 

230 """Compute n_th order polynomial via Horner's method.""" 

231 coeffs = np.array(coeffs, var.dtype.as_numpy_dtype) 

232 if not coeffs.size: 

233 return array_ops.zeros_like(var) 

234 return coeffs[0] + _create_polynomial(var, coeffs[1:]) * var 

235 

236 maybe_complement_p = array_ops.where_v2(p > -np.expm1(-2.), 1. - p, p) 

237 # Write in an arbitrary value in place of 0 for p since 0 will cause NaNs 

238 # later on. The result from the computation when p == 0 is not used so any 

239 # number that doesn't result in NaNs is fine. 

240 sanitized_mcp = array_ops.where_v2( 

241 maybe_complement_p <= 0., 

242 array_ops.fill(array_ops.shape(p), np.array(0.5, p.dtype.as_numpy_dtype)), 

243 maybe_complement_p) 

244 

245 # Compute x for p > exp(-2): x/sqrt(2pi) = w + w**3 P0(w**2)/Q0(w**2). 

246 w = sanitized_mcp - 0.5 

247 ww = w ** 2 

248 x_for_big_p = w + w * ww * (_create_polynomial(ww, p0) 

249 / _create_polynomial(ww, q0)) 

250 x_for_big_p *= -np.sqrt(2. * np.pi) 

251 

252 # Compute x for p <= exp(-2): x = z - log(z)/z - (1/z) P(1/z) / Q(1/z), 

253 # where z = sqrt(-2. * log(p)), and P/Q are chosen between two different 

254 # arrays based on whether p < exp(-32). 

255 z = math_ops.sqrt(-2. * math_ops.log(sanitized_mcp)) 

256 first_term = z - math_ops.log(z) / z 

257 second_term_small_p = ( 

258 _create_polynomial(1. / z, p2) / 

259 _create_polynomial(1. / z, q2) / z) 

260 second_term_otherwise = ( 

261 _create_polynomial(1. / z, p1) / 

262 _create_polynomial(1. / z, q1) / z) 

263 x_for_small_p = first_term - second_term_small_p 

264 x_otherwise = first_term - second_term_otherwise 

265 

266 x = array_ops.where_v2( 

267 sanitized_mcp > np.exp(-2.), x_for_big_p, 

268 array_ops.where_v2(z >= 8.0, x_for_small_p, x_otherwise)) 

269 

270 x = array_ops.where_v2(p > 1. - np.exp(-2.), x, -x) 

271 infinity_scalar = constant_op.constant(np.inf, dtype=p.dtype) 

272 infinity = array_ops.fill(array_ops.shape(p), infinity_scalar) 

273 x_nan_replaced = array_ops.where_v2(p <= 0.0, -infinity, 

274 array_ops.where_v2(p >= 1.0, infinity, x)) 

275 return x_nan_replaced 

276 

277 

278def log_ndtr(x, series_order=3, name="log_ndtr"): 

279 """Log Normal distribution function. 

280 

281 For details of the Normal distribution function see `ndtr`. 

282 

283 This function calculates `(log o ndtr)(x)` by either calling `log(ndtr(x))` or 

284 using an asymptotic series. Specifically: 

285 - For `x > upper_segment`, use the approximation `-ndtr(-x)` based on 

286 `log(1-x) ~= -x, x << 1`. 

287 - For `lower_segment < x <= upper_segment`, use the existing `ndtr` technique 

288 and take a log. 

289 - For `x <= lower_segment`, we use the series approximation of erf to compute 

290 the log CDF directly. 

291 

292 The `lower_segment` is set based on the precision of the input: 

293 

294 ``` 

295 lower_segment = { -20, x.dtype=float64 

296 { -10, x.dtype=float32 

297 upper_segment = { 8, x.dtype=float64 

298 { 5, x.dtype=float32 

299 ``` 

300 

301 When `x < lower_segment`, the `ndtr` asymptotic series approximation is: 

302 

303 ``` 

304 ndtr(x) = scale * (1 + sum) + R_N 

305 scale = exp(-0.5 x**2) / (-x sqrt(2 pi)) 

306 sum = Sum{(-1)^n (2n-1)!! / (x**2)^n, n=1:N} 

307 R_N = O(exp(-0.5 x**2) (2N+1)!! / |x|^{2N+3}) 

308 ``` 

309 

310 where `(2n-1)!! = (2n-1) (2n-3) (2n-5) ... (3) (1)` is a 

311 [double-factorial](https://en.wikipedia.org/wiki/Double_factorial). 

312 

313 

314 Args: 

315 x: `Tensor` of type `float32`, `float64`. 

316 series_order: Positive Python `integer`. Maximum depth to 

317 evaluate the asymptotic expansion. This is the `N` above. 

318 name: Python string. A name for the operation (default="log_ndtr"). 

319 

320 Returns: 

321 log_ndtr: `Tensor` with `dtype=x.dtype`. 

322 

323 Raises: 

324 TypeError: if `x.dtype` is not handled. 

325 TypeError: if `series_order` is a not Python `integer.` 

326 ValueError: if `series_order` is not in `[0, 30]`. 

327 """ 

328 if not isinstance(series_order, int): 

329 raise TypeError("series_order must be a Python integer.") 

330 if series_order < 0: 

331 raise ValueError("series_order must be non-negative.") 

332 if series_order > 30: 

333 raise ValueError("series_order must be <= 30.") 

334 

335 with ops.name_scope(name, values=[x]): 

336 x = ops.convert_to_tensor(x, name="x") 

337 

338 if x.dtype.as_numpy_dtype == np.float64: 

339 lower_segment = LOGNDTR_FLOAT64_LOWER 

340 upper_segment = LOGNDTR_FLOAT64_UPPER 

341 elif x.dtype.as_numpy_dtype == np.float32: 

342 lower_segment = LOGNDTR_FLOAT32_LOWER 

343 upper_segment = LOGNDTR_FLOAT32_UPPER 

344 else: 

345 raise TypeError("x.dtype=%s is not supported." % x.dtype) 

346 

347 # The basic idea here was ported from: 

348 # https://root.cern.ch/doc/v608/SpecFuncCephesInv_8cxx_source.html 

349 # We copy the main idea, with a few changes 

350 # * For x >> 1, and X ~ Normal(0, 1), 

351 # Log[P[X < x]] = Log[1 - P[X < -x]] approx -P[X < -x], 

352 # which extends the range of validity of this function. 

353 # * We use one fixed series_order for all of 'x', rather than adaptive. 

354 # * Our docstring properly reflects that this is an asymptotic series, not a 

355 # Taylor series. We also provided a correct bound on the remainder. 

356 # * We need to use the max/min in the _log_ndtr_lower arg to avoid nan when 

357 # x=0. This happens even though the branch is unchosen because when x=0 

358 # the gradient of a select involves the calculation 1*dy+0*(-inf)=nan 

359 # regardless of whether dy is finite. Note that the minimum is a NOP if 

360 # the branch is chosen. 

361 return array_ops.where_v2( 

362 math_ops.greater(x, upper_segment), 

363 -_ndtr(-x), # log(1-x) ~= -x, x << 1 # pylint: disable=invalid-unary-operand-type 

364 array_ops.where_v2( 

365 math_ops.greater(x, lower_segment), 

366 math_ops.log(_ndtr(math_ops.maximum(x, lower_segment))), 

367 _log_ndtr_lower(math_ops.minimum(x, lower_segment), series_order))) 

368 

369 

370def _log_ndtr_lower(x, series_order): 

371 """Asymptotic expansion version of `Log[cdf(x)]`, appropriate for `x<<-1`.""" 

372 x_2 = math_ops.square(x) 

373 # Log of the term multiplying (1 + sum) 

374 log_scale = -0.5 * x_2 - math_ops.log(-x) - 0.5 * np.log(2. * np.pi) 

375 return log_scale + math_ops.log(_log_ndtr_asymptotic_series(x, series_order)) 

376 

377 

378def _log_ndtr_asymptotic_series(x, series_order): 

379 """Calculates the asymptotic series used in log_ndtr.""" 

380 dtype = x.dtype.as_numpy_dtype 

381 if series_order <= 0: 

382 return np.array(1, dtype) 

383 x_2 = math_ops.square(x) 

384 even_sum = array_ops.zeros_like(x) 

385 odd_sum = array_ops.zeros_like(x) 

386 x_2n = x_2 # Start with x^{2*1} = x^{2*n} with n = 1. 

387 for n in range(1, series_order + 1): 

388 y = np.array(_double_factorial(2 * n - 1), dtype) / x_2n 

389 if n % 2: 

390 odd_sum += y 

391 else: 

392 even_sum += y 

393 x_2n *= x_2 

394 return 1. + even_sum - odd_sum 

395 

396 

397def erfinv(x, name="erfinv"): 

398 """The inverse function for erf, the error function. 

399 

400 Args: 

401 x: `Tensor` of type `float32`, `float64`. 

402 name: Python string. A name for the operation (default="erfinv"). 

403 

404 Returns: 

405 x: `Tensor` with `dtype=x.dtype`. 

406 

407 Raises: 

408 TypeError: if `x` is not floating-type. 

409 """ 

410 

411 with ops.name_scope(name, values=[x]): 

412 x = ops.convert_to_tensor(x, name="x") 

413 if x.dtype.as_numpy_dtype not in [np.float32, np.float64]: 

414 raise TypeError( 

415 "x.dtype=%s is not handled, see docstring for supported types." 

416 % x.dtype) 

417 return ndtri((x + 1.0) / 2.0) / np.sqrt(2) 

418 

419 

420def _double_factorial(n): 

421 """The double factorial function for small Python integer `n`.""" 

422 return np.prod(np.arange(n, 1, -2)) 

423 

424 

425def log_cdf_laplace(x, name="log_cdf_laplace"): 

426 """Log Laplace distribution function. 

427 

428 This function calculates `Log[L(x)]`, where `L(x)` is the cumulative 

429 distribution function of the Laplace distribution, i.e. 

430 

431 ```L(x) := 0.5 * int_{-infty}^x e^{-|t|} dt``` 

432 

433 For numerical accuracy, `L(x)` is computed in different ways depending on `x`, 

434 

435 ``` 

436 x <= 0: 

437 Log[L(x)] = Log[0.5] + x, which is exact 

438 

439 0 < x: 

440 Log[L(x)] = Log[1 - 0.5 * e^{-x}], which is exact 

441 ``` 

442 

443 Args: 

444 x: `Tensor` of type `float32`, `float64`. 

445 name: Python string. A name for the operation (default="log_ndtr"). 

446 

447 Returns: 

448 `Tensor` with `dtype=x.dtype`. 

449 

450 Raises: 

451 TypeError: if `x.dtype` is not handled. 

452 """ 

453 

454 with ops.name_scope(name, values=[x]): 

455 x = ops.convert_to_tensor(x, name="x") 

456 

457 # For x < 0, L(x) = 0.5 * exp{x} exactly, so Log[L(x)] = log(0.5) + x. 

458 lower_solution = -np.log(2.) + x 

459 

460 # safe_exp_neg_x = exp{-x} for x > 0, but is 

461 # bounded above by 1, which avoids 

462 # log[1 - 1] = -inf for x = log(1/2), AND 

463 # exp{-x} --> inf, for x << -1 

464 safe_exp_neg_x = math_ops.exp(-math_ops.abs(x)) 

465 

466 # log1p(z) = log(1 + z) approx z for |z| << 1. This approximation is used 

467 # internally by log1p, rather than being done explicitly here. 

468 upper_solution = math_ops.log1p(-0.5 * safe_exp_neg_x) 

469 

470 return array_ops.where_v2(x < 0., lower_solution, upper_solution)