#!/usr/bin/env python
import numpy as np
from scipy import interpolate, special
import mpmath as mp
[docs]def log_rician(observed, signal, noise_var, derivative="None", log_bessel='approximate'):
"""
Models the noise as a Rician distribution
:param observed: observed signal
:param signal: model signal
:param noise_var: variance of the noise
:param derivative: set to one of:
- "None": no derivative
- "signal": add derivatives wrt signal
- "noise_var": add derivative wrt noise variance
- "both": add derivatives wrt signal and noise variance
:param log_bessel: method to evaluate the log-bessel function; one of:
- "approximate": Uses scipy if possible and approximations if too small or large (error of up to 1%, but typically much smaller)
- "scipy": uses the default scipy interpretation, which will return -np.inf or np.inf if too small/large
- "precise": Uses mpmath to precisely evaluate the value. This is very slow
- "interpolate": use cubic interpolation between pre-computed "precise" values. Will become slow if values vary too much.
:return: log(P) of the Rician distribution
"""
return log_non_central_chi(observed, signal, noise_var, 1, derivative=derivative, log_bessel=log_bessel)
[docs]def log_non_central_chi(observed, signal, noise_var, ncoils, derivative="None", log_bessel='approximate'):
"""
Models the noise as a non-central chi distribution
:param observed: observed signal
:param signal: model signal
:param noise_var: variance of the noise
:param ncoils: number of coils (equal to degrees of freedom divided by 2)
:param derivative: set to one of:
- "None": no derivative
- "signal": add derivatives wrt signal
- "noise_var": add derivative wrt noise variance
- "both": add derivatives wrt signal and noise variance
:param log_bessel: method to evaluate the log-bessel function; one of:
- "approximate": Uses scipy if possible and approximations if too small or large (error of up to 1%, but typically much smaller)
- "scipy": uses the default scipy interpretation, which will return -np.inf or np.inf if too small/large
- "precise": Uses mpmath to precisely evaluate the value. This is very slow
- "interpolate": use cubic interpolation between pre-computed "precise" values. Will become slow if values vary too much.
:return: log(P) of the non-central chi distribution and its derivatives if set
"""
observed = np.asarray(observed).astype('float')
bessel_arg = signal * observed / noise_var
log_bessel_val = LogBessel.eval(ncoils -1, bessel_arg, method=log_bessel)
logP = (
ncoils * np.log(observed) - np.log(noise_var) - (ncoils - 1) * np.log(signal) -
(signal ** 2 + observed ** 2) / (2 * noise_var) +
log_bessel_val
)
if derivative.lower() == "none":
return logP
if derivative.lower() not in ('signal', 'noise_var', 'both'):
raise ValueError("derivative should be one of ('none', 'signal', 'noise_var', 'both'")
log_bessel_alt = LogBessel.eval(ncoils, bessel_arg, method=log_bessel)
res = [logP]
if derivative.lower() in ("signal", "both"):
res.append(
(np.exp(log_bessel_alt - log_bessel_val) * observed - signal) / noise_var
)
if derivative.lower() in ("noise_var", "both"):
res.append(
(signal ** 2 + observed ** 2) / (2 * noise_var ** 2) -
np.exp(log_bessel_alt - log_bessel_val) * signal * observed / noise_var ** 2 -
ncoils / noise_var
)
return res
[docs]class LogBessel(object):
"""
Pre-compute the log-bessel function
"""
pre_computed = {}
precise_log_bessel = np.vectorize(lambda nu, y: float(mp.log(mp.besseli(float(nu), float(y)))))
[docs] @classmethod
def interpolate(cls, nu, value):
if nu in cls.pre_computed:
minval, maxval, interp = cls.pre_computed[nu]
if (value < maxval).all() and (value > minval).all():
return interp(value)
max_order = int(np.maximum(6, np.ceil(np.log10(np.amax(value)))))
min_order = int(np.minimum(-6, np.floor(np.log10(np.amin(value)))))
xval = np.logspace(min_order, max_order, max_order * 1000 + 1)
interp = interpolate.interp1d(xval, cls.precise_log_bessel(nu, xval), bounds_error=True, kind='cubic')
cls.pre_computed[nu] = (xval.min(), xval.max(), interp)
return interp(value)
[docs] @classmethod
def approximate(cls, nu, value):
lb = np.log(special.iv(nu, value))
lower = lb == -np.inf
if lower.any():
lb[lower] = (-special.loggamma(nu + 1) + nu * np.log(value / 2))[lower]
upper = lb == np.inf
if upper.any():
lb[upper] = np.broadcast_to(-np.log(2 * np.pi * value) / 2 + value, upper.shape)[upper]
return lb
[docs] @classmethod
def eval(cls, nu, value, method='approximate'):
"""
Evaluates the log of the modified Bessel function of the first kind using chosen method
:param nu: order of the Bessel
:param value: scalar/array to evaluate
:param method: method to evaluate the log-bessel function; one of:
- "approximate": Uses scipy if possible and approximations if too small or large (error of up to 1%, but typically much smaller)
- "scipy": uses the default scipy interpretation, which will return -np.inf or np.inf if too small/large
- "precise": Uses mpmath to precisely evaluate the value. This is very slow
- "interpolate": use cubic interpolation between pre-computed "precise" values. Will become slow if values vary too much.
:return: log(Bessel_<nu>(<value>))
"""
if method == 'approximate':
return cls.approximate(nu, value)
if method == 'scipy':
return np.log(special.iv(nu, value))
if method == 'precise':
return cls.precise_log_bessel(nu, value)
if method == 'interpolate':
return cls.interpolate(nu, value)
raise ValueError(f"Method {method} not in ('approximate', 'scipy', 'precise', 'interpolate')")