`reconst`¶

label='fast'

import dipy.reconst as dire
dire.bench()

import dipy.reconst as dire
dire.bench()

import dipy.reconst as dire
dire.bench()

>>> success = np.lib.bench() #doctest: +SKIP
Running benchmarks for numpy.lib
...
using 562341 items:
unique:
0.11
unique1d:
0.11
ratio: 1.0
nUnique: 56230 == 56230
...
OK

>>> success #doctest: +SKIP
True

>>> np.lib.test() #doctest: +SKIP

>>> result = np.lib.test() #doctest: +SKIP
Running unit tests for numpy.lib
...
Ran 976 tests in 3.933s

>>> result.errors #doctest: +SKIP
[]
>>> result.knownfail #doctest: +SKIP
[]

>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")  # doctest: +SKIP
Time for a single execution :  0.005 s

>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")  # doctest: +SKIP
Time for a single execution :  0.005 s

   y
   /               /               a/    
/                  /                   /__________          x      c     z

>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")  # doctest: +SKIP
Time for a single execution :  0.005 s

>>> from dipy.core.ndindex import ndindex
>>> shape = (3, 2, 1)
>>> for index in ndindex(shape):
...     print(index)
(0, 0, 0)
(0, 1, 0)
(1, 0, 0)
(1, 1, 0)
(2, 0, 0)
(2, 1, 0)

>>> arr = np.empty(3, dtype=object)
>>> arr.fill(2)
>>> old_squash(arr)
array([2, 2, 2])
>>> arr[0] = None
>>> old_squash(arr)
array([0, 2, 2])
>>> arr.fill(np.ones(2))
>>> r = old_squash(arr)
>>> r.shape == (3, 2)
True
>>> r.dtype
dtype('float64')

>>> arr = np.empty(3, dtype=object)
>>> arr.fill(2)
>>> quick_squash(arr)
array([2, 2, 2])
>>> arr[0] = None
>>> quick_squash(arr)
array([0, 2, 2])
>>> arr.fill(np.ones(2))
>>> r = quick_squash(arr)
>>> r.shape
(3, 2)
>>> r.dtype
dtype('float64')

>>> times = 10
>>> etime = np.testing.measure('for i in range(1000): np.sqrt(i**2)', times=times)
>>> print("Time for a single execution : ", etime / times, "s")  # doctest: +SKIP
Time for a single execution :  0.005 s

>>> np.random.randn()
2.1923875335537315 #random

>>> 2.5 * np.random.randn(2, 4) + 3
array([[-4.49401501,  4.00950034, -1.81814867,  7.29718677],  #random
       [ 0.39924804,  4.68456316,  4.99394529,  4.84057254]]) #random

>>> vecs = np.arange(9).reshape((1, 3, 3))
>>> vals = np.arange(3).reshape((1, 3))
>>> vec_val_vect(vecs, vals)
array([[[   9.,   24.,   39.],
        [  24.,   66.,  108.],
        [  39.,  108.,  177.]]])

>>> vecs = np.arange(9).reshape((3, 3))
>>> vals = np.arange(3)
>>> np.dot(vecs, np.dot(np.diag(vals), vecs.T))
array([[  9,  24,  39],
       [ 24,  66, 108],
       [ 39, 108, 177]])

class MyModel(Model, Cache):
    pass

class MyModelFit(Fit):
    pass

def odf(sphere):
    M = self.model.cache_get('odf_basis_matrix', key=sphere)

    if M is None:
        M = self._compute_basis_matrix(sphere)
        self.model.cache_set('odf_basis_matrix', key=sphere, value=M)

>>> def compute_expensive_matrix(parameters):
...     # Imagine the following computation is very expensive
...     return (p**2 for p in parameters)

>>> c = Cache()

>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)

>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)

>>> X1 is X2
True

>>> class MagicProp(object):
...     @auto_attr
...     def a(self):
...         return 99
...
>>> x = MagicProp()
>>> 'a' in x.__dict__
False
>>> x.a
99
>>> 'a' in x.__dict__
True

>>> import numpy as np
>>> from dipy.data import get_sphere
>>> sphere = get_sphere('symmetric362')
>>> verts, faces = sphere.vertices, sphere.faces
>>> verts.shape == (362, 3)
True
>>> faces.shape == (720, 3)
True
>>> verts, faces = get_sphere('not a sphere name') #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
    ...
DataError: No sphere called "not a sphere name"

>>> from dipy.core.ndindex import ndindex
>>> shape = (3, 2, 1)
>>> for index in ndindex(shape):
...     print(index)
(0, 0, 0)
(0, 1, 0)
(1, 0, 0)
(1, 1, 0)
(2, 0, 0)
(2, 1, 0)

double func(double x)
double func(double x, void *user_data)
double func(int n, double *xx)
double func(int n, double *xx, void *user_data)

>>> from scipy import integrate
>>> x2 = lambda x: x**2
>>> integrate.quad(x2, 0, 4)
(21.333333333333332, 2.3684757858670003e-13)
>>> print(4**3 / 3.)  # analytical result
21.3333333333

>>> invexp = lambda x: np.exp(-x)
>>> integrate.quad(invexp, 0, np.inf)
(1.0, 5.842605999138044e-11)

>>> f = lambda x,a : a*x
>>> y, err = integrate.quad(f, 0, 1, args=(1,))
>>> y
0.5
>>> y, err = integrate.quad(f, 0, 1, args=(3,))
>>> y
1.5

testlib.c =>
    double func(int n, double args[n]){
        return args[0]*args[0] + args[1]*args[1];}
compile to library testlib.*

from scipy import integrate
import ctypes
lib = ctypes.CDLL('/home/.../testlib.*') #use absolute path
lib.func.restype = ctypes.c_double
lib.func.argtypes = (ctypes.c_int,ctypes.c_double)
integrate.quad(lib.func,0,1,(1))
#(1.3333333333333333, 1.4802973661668752e-14)
print((1.0**3/3.0 + 1.0) - (0.0**3/3.0 + 0.0)) #Analytic result
# 1.3333333333333333

>>> y = lambda x: 1 if x<=0 else 0
>>> integrate.quad(y, -1, 1)
(1.0, 1.1102230246251565e-14)
>>> integrate.quad(y, -1, 100)
(1.0000000002199108, 1.0189464580163188e-08)
>>> integrate.quad(y, -1, 10000)
(0.0, 0.0)

>>> import numpy as np
>>> from dipy.core.geometry import vec2vec_rotmat
>>> u=np.array([1,0,0])
>>> v=np.array([0,1,0])
>>> R=vec2vec_rotmat(u,v)
>>> np.dot(R,u)
array([ 0.,  1.,  0.])
>>> np.dot(R.T,v)
array([ 1.,  0.,  0.])

ADC =       heta Q  heta^T

>>> import numpy as np
>>> from dipy.data import get_sphere
>>> sphere = get_sphere('symmetric362')
>>> verts, faces = sphere.vertices, sphere.faces
>>> verts.shape == (362, 3)
True
>>> faces.shape == (720, 3)
True
>>> verts, faces = get_sphere('not a sphere name') #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
    ...
DataError: No sphere called "not a sphere name"

>>> from dipy.core.ndindex import ndindex
>>> shape = (3, 2, 1)
>>> for index in ndindex(shape):
...     print(index)
(0, 0, 0)
(0, 1, 0)
(1, 0, 0)
(1, 1, 0)
(2, 0, 0)
(2, 1, 0)

>>> vecs = np.arange(9).reshape((1, 3, 3))
>>> vals = np.arange(3).reshape((1, 3))
>>> vec_val_vect(vecs, vals)
array([[[   9.,   24.,   39.],
        [  24.,   66.,  108.],
        [  39.,  108.,  177.]]])

>>> vecs = np.arange(9).reshape((3, 3))
>>> vals = np.arange(3)
>>> np.dot(vecs, np.dot(np.diag(vals), vecs.T))
array([[  9,  24,  39],
       [ 24,  66, 108],
       [ 39, 108, 177]])

>>> import numpy as np
>>> from dipy.data import get_sphere
>>> sphere = get_sphere('symmetric362')
>>> verts, faces = sphere.vertices, sphere.faces
>>> verts.shape == (362, 3)
True
>>> faces.shape == (720, 3)
True
>>> verts, faces = get_sphere('not a sphere name') #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
    ...
DataError: No sphere called "not a sphere name"

>>> from dipy.core.ndindex import ndindex
>>> shape = (3, 2, 1)
>>> for index in ndindex(shape):
...     print(index)
(0, 0, 0)
(0, 1, 0)
(1, 0, 0)
(1, 1, 0)
(2, 0, 0)
(2, 1, 0)

>>> vecs = np.arange(9).reshape((1, 3, 3))
>>> vals = np.arange(3).reshape((1, 3))
>>> vec_val_vect(vecs, vals)
array([[[   9.,   24.,   39.],
        [  24.,   66.,  108.],
        [  39.,  108.,  177.]]])

>>> vecs = np.arange(9).reshape((3, 3))
>>> vals = np.arange(3)
>>> np.dot(vecs, np.dot(np.diag(vals), vecs.T))
array([[  9,  24,  39],
       [ 24,  66, 108],
       [ 39, 108, 177]])

class MyModel(Model, Cache):
    pass

class MyModelFit(Fit):
    pass

def odf(sphere):
    M = self.model.cache_get('odf_basis_matrix', key=sphere)

    if M is None:
        M = self._compute_basis_matrix(sphere)
        self.model.cache_set('odf_basis_matrix', key=sphere, value=M)

>>> def compute_expensive_matrix(parameters):
...     # Imagine the following computation is very expensive
...     return (p**2 for p in parameters)

>>> c = Cache()

>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)

>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)

>>> X1 is X2
True

>>> import warnings
>>> from dipy.data import dsi_voxels, default_sphere
>>> data, gtab = dsi_voxels()
>>> from dipy.reconst.dsi import DiffusionSpectrumModel
>>> ds = DiffusionSpectrumModel(gtab)
>>> dsfit = ds.fit(data)
>>> from dipy.reconst.odf import gfa
>>> np.round(gfa(dsfit.odf(default_sphere))[0, 0, 0], 2)
0.11

y[j_1,..,j_d] = sum[k_1=0..n_1-1, ..., k_d=0..n_d-1]
   x[k_1,..,k_d] * prod[i=1..d] exp(-sqrt(-1)*2*pi/n_i * j_i * k_i)

>>> from scipy.fftpack import fftn, ifftn
>>> y = (-np.arange(16), 8 - np.arange(16), np.arange(16))
>>> np.allclose(y, fftn(ifftn(y)))
True

>>> freqs = np.fft.fftfreq(10, 0.1)
>>> freqs
array([ 0.,  1.,  2.,  3.,  4., -5., -4., -3., -2., -1.])
>>> np.fft.fftshift(freqs)
array([-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.])

>>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
>>> freqs
array([[ 0.,  1.,  2.],
       [ 3.,  4., -4.],
       [-3., -2., -1.]])
>>> np.fft.fftshift(freqs, axes=(1,))
array([[ 2.,  0.,  1.],
       [-4.,  3.,  4.],
       [-1., -3., -2.]])

>>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3)
>>> freqs
array([[ 0.,  1.,  2.],
       [ 3.,  4., -4.],
       [-3., -2., -1.]])
>>> np.fft.ifftshift(np.fft.fftshift(freqs))
array([[ 0.,  1.,  2.],
       [ 3.,  4., -4.],
       [-3., -2., -1.]])

>>> from scipy import ndimage
>>> a = np.arange(12.).reshape((4, 3))
>>> a
array([[  0.,   1.,   2.],
       [  3.,   4.,   5.],
       [  6.,   7.,   8.],
       [  9.,  10.,  11.]])
>>> ndimage.map_coordinates(a, [[0.5, 2], [0.5, 1]], order=1)
array([ 2.,  7.])

>>> inds = np.array([[0.5, 2], [0.5, 4]])
>>> ndimage.map_coordinates(a, inds, order=1, cval=-33.3)
array([  2. , -33.3])
>>> ndimage.map_coordinates(a, inds, order=1, mode='nearest')
array([ 2.,  8.])
>>> ndimage.map_coordinates(a, inds, order=1, cval=0, output=bool)
array([ True, False], dtype=bool)

ADC =       heta Q  heta^T

>>> class MagicProp(object):
...     @auto_attr
...     def a(self):
...         return 99
...
>>> x = MagicProp()
>>> 'a' in x.__dict__
False
>>> x.a
99
>>> 'a' in x.__dict__
True

>>> import numpy as np
>>> from dipy.data import get_sphere
>>> sphere = get_sphere('symmetric362')
>>> verts, faces = sphere.vertices, sphere.faces
>>> verts.shape == (362, 3)
True
>>> faces.shape == (720, 3)
True
>>> verts, faces = get_sphere('not a sphere name') #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
    ...
DataError: No sphere called "not a sphere name"

>>> from dipy.core.gradients import gradient_table
>>> bvals = 1500 * np.ones(7)
>>> bvals[0] = 0
>>> sq2 = np.sqrt(2) / 2
>>> bvecs = np.array([[0, 0, 0],
...                   [1, 0, 0],
...                   [0, 1, 0],
...                   [0, 0, 1],
...                   [sq2, sq2, 0],
...                   [sq2, 0, sq2],
...                   [0, sq2, sq2]])
>>> gt = gradient_table(bvals, bvecs)
>>> gt.bvecs.shape == bvecs.shape
True
>>> gt = gradient_table(bvals, bvecs.T)
>>> gt.bvecs.shape == bvecs.T.shape
False

>>> vecs = np.arange(9).reshape((1, 3, 3))
>>> vals = np.arange(3).reshape((1, 3))
>>> vec_val_vect(vecs, vals)
array([[[   9.,   24.,   39.],
        [  24.,   66.,  108.],
        [  39.,  108.,  177.]]])

>>> vecs = np.arange(9).reshape((3, 3))
>>> vals = np.arange(3)
>>> np.dot(vecs, np.dot(np.diag(vals), vecs.T))
array([[  9,  24,  39],
       [ 24,  66, 108],
       [ 39, 108, 177]])

>>> import numpy as np
>>> vec = [[8, 15, 0], [0, 36, 77]]
>>> vector_norm(vec)
array([ 17.,  85.])
>>> vector_norm(vec, keepdims=True)
array([[ 17.],
       [ 85.]])
>>> vector_norm(vec, axis=0)
array([  8.,  39.,  77.])

class MyModel(Model, Cache):
    pass

class MyModelFit(Fit):
    pass

def odf(sphere):
    M = self.model.cache_get('odf_basis_matrix', key=sphere)

    if M is None:
        M = self._compute_basis_matrix(sphere)
        self.model.cache_set('odf_basis_matrix', key=sphere, value=M)

>>> def compute_expensive_matrix(parameters):
...     # Imagine the following computation is very expensive
...     return (p**2 for p in parameters)

>>> c = Cache()

>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)

>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)

>>> X1 is X2
True

>>> from dipy.data import default_sphere, get_3shell_gtab
>>> gtab = get_3shell_gtab()
>>> from dipy.sims.voxel import multi_tensor
>>> mevals = np.array(([0.0017, 0.0003, 0.0003], 
...                    [0.0017, 0.0003, 0.0003]))
>>> angl = [(0, 0), (60, 0)]
>>> data, sticks = multi_tensor(gtab,
...                             mevals,
...                             S0=100.0,
...                             angles=angl,
...                             fractions=[50, 50],
...                             snr=None)
>>> from dipy.reconst.forecast import ForecastModel
>>> fm = ForecastModel(gtab, sh_order=6)
>>> f_fit = fm.fit(data)
>>> d_par = f_fit.dpar
>>> d_perp = f_fit.dperp
>>> fodf = f_fit.odf(default_sphere)

x = arg min(sum(func(y)**2,axis=0))
         y

func(params) = ydata - f(xdata, params)

  min   sum((ydata - f(xdata, params))**2, axis=0)
params

>>> from dipy.utils.optpkg import optional_package
>>> pkg, have_pkg, setup_module = optional_package('not_a_package')

>>> have_pkg
False

>>> pkg.some_function() #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
    ...
TripWireError: We need package not_a_package for these functions, but
``import not_a_package`` raised an ImportError

>>> pkg, _, _ = optional_package('os')
>>> hasattr(pkg, 'path')
True

>>> subpkg, _, _ = optional_package('os.path')
>>> hasattr(subpkg, 'dirname')
True

ADC =       heta Q  heta^T

>>> from dipy.core.ndindex import ndindex
>>> shape = (3, 2, 1)
>>> for index in ndindex(shape):
...     print(index)
(0, 0, 0)
(0, 1, 0)
(1, 0, 0)
(1, 1, 0)
(2, 0, 0)
(2, 1, 0)

>>> vecs = np.arange(9).reshape((1, 3, 3))
>>> vals = np.arange(3).reshape((1, 3))
>>> vec_val_vect(vecs, vals)
array([[[   9.,   24.,   39.],
        [  24.,   66.,  108.],
        [  39.,  108.,  177.]]])

>>> vecs = np.arange(9).reshape((3, 3))
>>> vals = np.arange(3)
>>> np.dot(vecs, np.dot(np.diag(vals), vecs.T))
array([[  9,  24,  39],
       [ 24,  66, 108],
       [ 39, 108, 177]])

class MyModel(Model, Cache):
    pass

class MyModelFit(Fit):
    pass

def odf(sphere):
    M = self.model.cache_get('odf_basis_matrix', key=sphere)

    if M is None:
        M = self._compute_basis_matrix(sphere)
        self.model.cache_set('odf_basis_matrix', key=sphere, value=M)

>>> def compute_expensive_matrix(parameters):
...     # Imagine the following computation is very expensive
...     return (p**2 for p in parameters)

>>> c = Cache()

>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)

>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)

>>> X1 is X2
True

>>> from dipy.data import dsi_voxels
>>> data, gtab = dsi_voxels()
>>> from dipy.core.subdivide_octahedron import create_unit_sphere
>>> sphere = create_unit_sphere(5)
>>> from dipy.reconst.gqi import GeneralizedQSamplingModel
>>> gq = GeneralizedQSamplingModel(gtab, 'gqi2', 1.1)
>>> voxel_signal = data[0, 0, 0]
>>> odf = gq.fit(voxel_signal).odf(sphere)

\sqrt{\frac{n \sum_i{(\Psi_i - <\Psi>)^2}}{(n-1) \sum{\Psi_i ^ 2}}}

(signal -  S)^T(signal -  S)

>>> from scipy.optimize import rosen, differential_evolution
>>> bounds = [(0,2), (0, 2), (0, 2), (0, 2), (0, 2)]
>>> result = differential_evolution(rosen, bounds)
>>> result.x, result.fun
(array([1., 1., 1., 1., 1.]), 1.9216496320061384e-19)

>>> from scipy.optimize import differential_evolution
>>> import numpy as np
>>> def ackley(x):
...     arg1 = -0.2 * np.sqrt(0.5 * (x[0] ** 2 + x[1] ** 2))
...     arg2 = 0.5 * (np.cos(2. * np.pi * x[0]) + np.cos(2. * np.pi * x[1]))
...     return -20. * np.exp(arg1) - np.exp(arg2) + 20. + np.e
>>> bounds = [(-5, 5), (-5, 5)]
>>> result = differential_evolution(ackley, bounds)
>>> result.x, result.fun
(array([ 0.,  0.]), 4.4408920985006262e-16)

minimize F(x) = 0.5 * sum(rho(f_i(x)**2), i = 0, ..., m - 1)
subject to lb <= x <= ub

>>> def fun_rosenbrock(x):
...     return np.array([10 * (x[1] - x[0]**2), (1 - x[0])])

>>> from scipy.optimize import least_squares
>>> x0_rosenbrock = np.array([2, 2])
>>> res_1 = least_squares(fun_rosenbrock, x0_rosenbrock)
>>> res_1.x
array([ 1.,  1.])
>>> res_1.cost
9.8669242910846867e-30
>>> res_1.optimality
8.8928864934219529e-14

>>> def jac_rosenbrock(x):
...     return np.array([
...         [-20 * x[0], 10],
...         [-1, 0]])

>>> res_2 = least_squares(fun_rosenbrock, x0_rosenbrock, jac_rosenbrock,
...                       bounds=([-np.inf, 1.5], np.inf))
>>> res_2.x
array([ 1.22437075,  1.5       ])
>>> res_2.cost
0.025213093946805685
>>> res_2.optimality
1.5885401433157753e-07

>>> def fun_broyden(x):
...     f = (3 - x) * x + 1
...     f[1:] -= x[:-1]
...     f[:-1] -= 2 * x[1:]
...     return f

>>> from scipy.sparse import lil_matrix
>>> def sparsity_broyden(n):
...     sparsity = lil_matrix((n, n), dtype=int)
...     i = np.arange(n)
...     sparsity[i, i] = 1
...     i = np.arange(1, n)
...     sparsity[i, i - 1] = 1
...     i = np.arange(n - 1)
...     sparsity[i, i + 1] = 1
...     return sparsity
...
>>> n = 100000
>>> x0_broyden = -np.ones(n)
...
>>> res_3 = least_squares(fun_broyden, x0_broyden,
...                       jac_sparsity=sparsity_broyden(n))
>>> res_3.cost
4.5687069299604613e-23
>>> res_3.optimality
1.1650454296851518e-11

>>> def gen_data(t, a, b, c, noise=0, n_outliers=0, random_state=0):
...     y = a + b * np.exp(t * c)
...
...     rnd = np.random.RandomState(random_state)
...     error = noise * rnd.randn(t.size)
...     outliers = rnd.randint(0, t.size, n_outliers)
...     error[outliers] *= 10
...
...     return y + error
...
>>> a = 0.5
>>> b = 2.0
>>> c = -1
>>> t_min = 0
>>> t_max = 10
>>> n_points = 15
...
>>> t_train = np.linspace(t_min, t_max, n_points)
>>> y_train = gen_data(t_train, a, b, c, noise=0.1, n_outliers=3)

>>> def fun(x, t, y):
...     return x[0] + x[1] * np.exp(x[2] * t) - y
...
>>> x0 = np.array([1.0, 1.0, 0.0])

>>> res_lsq = least_squares(fun, x0, args=(t_train, y_train))

>>> res_soft_l1 = least_squares(fun, x0, loss='soft_l1', f_scale=0.1,
...                             args=(t_train, y_train))
>>> res_log = least_squares(fun, x0, loss='cauchy', f_scale=0.1,
...                         args=(t_train, y_train))

>>> t_test = np.linspace(t_min, t_max, n_points * 10)
>>> y_true = gen_data(t_test, a, b, c)
>>> y_lsq = gen_data(t_test, *res_lsq.x)
>>> y_soft_l1 = gen_data(t_test, *res_soft_l1.x)
>>> y_log = gen_data(t_test, *res_log.x)
...
>>> import matplotlib.pyplot as plt
>>> plt.plot(t_train, y_train, 'o')
>>> plt.plot(t_test, y_true, 'k', linewidth=2, label='true')
>>> plt.plot(t_test, y_lsq, label='linear loss')
>>> plt.plot(t_test, y_soft_l1, label='soft_l1 loss')
>>> plt.plot(t_test, y_log, label='cauchy loss')
>>> plt.xlabel("t")
>>> plt.ylabel("y")
>>> plt.legend()
>>> plt.show()

>>> def f(z):
...     return z - (0.5 + 0.5j)

>>> def f_wrap(x):
...     fx = f(x[0] + 1j*x[1])
...     return np.array([fx.real, fx.imag])

>>> from scipy.optimize import least_squares
>>> res_wrapped = least_squares(f_wrap, (0.1, 0.1), bounds=([0, 0], [1, 1]))
>>> z = res_wrapped.x[0] + res_wrapped.x[1]*1j
>>> z
(0.49999999999925893+0.49999999999925893j)

>>> from dipy.utils.optpkg import optional_package
>>> pkg, have_pkg, setup_module = optional_package('not_a_package')

>>> have_pkg
False

>>> pkg.some_function() #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
    ...
TripWireError: We need package not_a_package for these functions, but
``import not_a_package`` raised an ImportError

>>> pkg, _, _ = optional_package('os')
>>> hasattr(pkg, 'path')
True

>>> subpkg, _, _ = optional_package('os.path')
>>> hasattr(subpkg, 'dirname')
True

class MyModel(Model, Cache):
    pass

class MyModelFit(Fit):
    pass

def odf(sphere):
    M = self.model.cache_get('odf_basis_matrix', key=sphere)

    if M is None:
        M = self._compute_basis_matrix(sphere)
        self.model.cache_set('odf_basis_matrix', key=sphere, value=M)

>>> def compute_expensive_matrix(parameters):
...     # Imagine the following computation is very expensive
...     return (p**2 for p in parameters)

>>> c = Cache()

>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)

>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)

>>> X1 is X2
True

>>> from dipy.data import dsi_voxels, default_sphere
>>> from dipy.core.gradients import gradient_table
>>> _, gtab_ = dsi_voxels()
>>> gtab = gradient_table(gtab_.bvals, gtab_.bvecs,
...                       b0_threshold=gtab_.bvals.min())
>>> from dipy.sims.voxel import sticks_and_ball
>>> data, golden_directions = sticks_and_ball(gtab, d=0.0015, S0=1,
...                                           angles=[(0, 0),
...                                                   (90, 0)],
...                                           fractions=[50, 50],
...                                           snr=None)
>>> from dipy.reconst.mapmri import MapmriModel
>>> radial_order = 4
>>> map_model = MapmriModel(gtab, radial_order=radial_order)
>>> mapfit = map_model.fit(data)
>>> odf = mapfit.odf(default_sphere)

n!! = special.gamma(n/2+1)*2**((m+1)/2)/sqrt(pi)  n odd
    = 2**(n/2) * (n/2)!                           n even

>>> from scipy.special import factorial2
>>> factorial2(7, exact=False)
array(105.00000000000001)
>>> factorial2(7, exact=True)
105L

>>> from dipy.core.gradients import gradient_table
>>> bvals = 1500 * np.ones(7)
>>> bvals[0] = 0
>>> sq2 = np.sqrt(2) / 2
>>> bvecs = np.array([[0, 0, 0],
...                   [1, 0, 0],
...                   [0, 1, 0],
...                   [0, 0, 1],
...                   [sq2, sq2, 0],
...                   [sq2, 0, sq2],
...                   [0, sq2, sq2]])
>>> gt = gradient_table(bvals, bvecs)
>>> gt.bvecs.shape == bvecs.shape
True
>>> gt = gradient_table(bvals, bvecs.T)
>>> gt.bvecs.shape == bvecs.T.shape
False

>>> from dipy.utils.optpkg import optional_package
>>> pkg, have_pkg, setup_module = optional_package('not_a_package')

>>> have_pkg
False

>>> pkg.some_function() #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
    ...
TripWireError: We need package not_a_package for these functions, but
``import not_a_package`` raised an ImportError

>>> pkg, _, _ = optional_package('os')
>>> hasattr(pkg, 'path')
True

>>> subpkg, _, _ = optional_package('os.path')
>>> hasattr(subpkg, 'dirname')
True

n! = n * (n - 1) * (n - 2) * ... * 1

>>> from scipy.special import factorial
>>> arr = np.array([3, 4, 5])
>>> factorial(arr, exact=False)
array([   6.,   24.,  120.])
>>> factorial(arr, exact=True)
array([  6,  24, 120])
>>> factorial(5, exact=True)
120L

>>> from dipy.utils.optpkg import optional_package
>>> pkg, have_pkg, setup_module = optional_package('not_a_package')

>>> have_pkg
False

>>> pkg.some_function() #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
    ...
TripWireError: We need package not_a_package for these functions, but
``import not_a_package`` raised an ImportError

>>> pkg, _, _ = optional_package('os')
>>> hasattr(pkg, 'path')
True

>>> subpkg, _, _ = optional_package('os.path')
>>> hasattr(subpkg, 'dirname')
True

>>> class MagicProp(object):
...     @auto_attr
...     def a(self):
...         return 99
...
>>> x = MagicProp()
>>> 'a' in x.__dict__
False
>>> x.a
99
>>> 'a' in x.__dict__
True

>>> from dipy.core.ndindex import ndindex
>>> shape = (3, 2, 1)
>>> for index in ndindex(shape):
...     print(index)
(0, 0, 0)
(0, 1, 0)
(1, 0, 0)
(1, 1, 0)
(2, 0, 0)
(2, 1, 0)

>>> from dipy.core.ndindex import ndindex
>>> shape = (3, 2, 1)
>>> for index in ndindex(shape):
...     print(index)
(0, 0, 0)
(0, 1, 0)
(1, 0, 0)
(1, 1, 0)
(2, 0, 0)
(2, 1, 0)

\sqrt{\frac{n \sum_i{(\Psi_i - <\Psi>)^2}}{(n-1) \sum{\Psi_i ^ 2}}}

class MyModel(Model, Cache):
    pass

class MyModelFit(Fit):
    pass

def odf(sphere):
    M = self.model.cache_get('odf_basis_matrix', key=sphere)

    if M is None:
        M = self._compute_basis_matrix(sphere)
        self.model.cache_set('odf_basis_matrix', key=sphere, value=M)

>>> def compute_expensive_matrix(parameters):
...     # Imagine the following computation is very expensive
...     return (p**2 for p in parameters)

>>> c = Cache()

>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)

>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)

>>> X1 is X2
True

n! = n * (n - 1) * (n - 2) * ... * 1

>>> from scipy.special import factorial
>>> arr = np.array([3, 4, 5])
>>> factorial(arr, exact=False)
array([   6.,   24.,  120.])
>>> factorial(arr, exact=True)
array([  6,  24, 120])
>>> factorial(5, exact=True)
120L

n!! = special.gamma(n/2+1)*2**((m+1)/2)/sqrt(pi)  n odd
    = 2**(n/2) * (n/2)!                           n even

>>> from scipy.special import factorial2
>>> factorial2(7, exact=False)
array(105.00000000000001)
>>> factorial2(7, exact=True)
105L

>>> from dipy.core.gradients import (
...    gradient_table_from_gradient_strength_bvecs)
>>> gradient_strength = .03e-3 * np.ones(7)  # clinical strength at 30 mT/m
>>> big_delta = .03  # pulse separation of 30ms
>>> small_delta = 0.01  # pulse duration of 10ms
>>> gradient_strength[0] = 0
>>> sq2 = np.sqrt(2) / 2
>>> bvecs = np.array([[0, 0, 0],
...                   [1, 0, 0],
...                   [0, 1, 0],
...                   [0, 0, 1],
...                   [sq2, sq2, 0],
...                   [sq2, 0, sq2],
...                   [0, sq2, sq2]])
>>> gt = gradient_table_from_gradient_strength_bvecs(
...     gradient_strength, bvecs, big_delta, small_delta)

>>> from dipy.utils.optpkg import optional_package
>>> pkg, have_pkg, setup_module = optional_package('not_a_package')

>>> have_pkg
False

>>> pkg.some_function() #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
    ...
TripWireError: We need package not_a_package for these functions, but
``import not_a_package`` raised an ImportError

>>> pkg, _, _ = optional_package('os')
>>> hasattr(pkg, 'path')
True

>>> subpkg, _, _ = optional_package('os.path')
>>> hasattr(subpkg, 'dirname')
True

class MyModel(Model, Cache):
    pass

class MyModelFit(Fit):
    pass

def odf(sphere):
    M = self.model.cache_get('odf_basis_matrix', key=sphere)

    if M is None:
        M = self._compute_basis_matrix(sphere)
        self.model.cache_set('odf_basis_matrix', key=sphere, value=M)

>>> def compute_expensive_matrix(parameters):
...     # Imagine the following computation is very expensive
...     return (p**2 for p in parameters)

>>> c = Cache()

>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)

>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)

>>> X1 is X2
True

>>> class MagicProp(object):
...     @auto_attr
...     def a(self):
...         return 99
...
>>> x = MagicProp()
>>> 'a' in x.__dict__
False
>>> x.a
99
>>> 'a' in x.__dict__
True

>>> a = np.array([[1, np.nan], [3, 4]])
>>> np.nanmean(a)
2.6666666666666665
>>> np.nanmean(a, axis=0)
array([ 2.,  4.])
>>> np.nanmean(a, axis=1)
array([ 1.,  3.5])

>>> from dipy.utils.optpkg import optional_package
>>> pkg, have_pkg, setup_module = optional_package('not_a_package')

>>> have_pkg
False

>>> pkg.some_function() #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
    ...
TripWireError: We need package not_a_package for these functions, but
``import not_a_package`` raised an ImportError

>>> pkg, _, _ = optional_package('os')
>>> hasattr(pkg, 'path')
True

>>> subpkg, _, _ = optional_package('os.path')
>>> hasattr(subpkg, 'dirname')
True

>>> import dipy.data as dpd
>>> data, gtab = dpd.dsi_voxels()
>>> sphere = dpd.get_sphere()
>>> from dipy.reconst.sfm import sfm_design_matrix

>>> tensor_matrix = sfm_design_matrix(gtab, sphere,
...                                   [0.0015, 0.0005, 0.0005])

>>> stick_matrix = sfm_design_matrix(gtab, sphere, [0.001, 0, 0])

class MyModel(Model, Cache):
    pass

class MyModelFit(Fit):
    pass

def odf(sphere):
    M = self.model.cache_get('odf_basis_matrix', key=sphere)

    if M is None:
        M = self._compute_basis_matrix(sphere)
        self.model.cache_set('odf_basis_matrix', key=sphere, value=M)

>>> def compute_expensive_matrix(parameters):
...     # Imagine the following computation is very expensive
...     return (p**2 for p in parameters)

>>> c = Cache()

>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)

>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)

>>> X1 is X2
True

>>> class MagicProp(object):
...     @auto_attr
...     def a(self):
...         return 99
...
>>> x = MagicProp()
>>> 'a' in x.__dict__
False
>>> x.a
99
>>> 'a' in x.__dict__
True

>>> a = np.array([[1, 2], [3, 4]])
>>> b = np.array([[5, 6]])
>>> np.concatenate((a, b), axis=0)
array([[1, 2],
       [3, 4],
       [5, 6]])
>>> np.concatenate((a, b.T), axis=1)
array([[1, 2, 5],
       [3, 4, 6]])
>>> np.concatenate((a, b), axis=None)
array([1, 2, 3, 4, 5, 6])

>>> a = np.ma.arange(3)
>>> a[1] = np.ma.masked
>>> b = np.arange(2, 5)
>>> a
masked_array(data=[0, --, 2],
             mask=[False,  True, False],
       fill_value=999999)
>>> b
array([2, 3, 4])
>>> np.concatenate([a, b])
masked_array(data=[0, 1, 2, 2, 3, 4],
             mask=False,
       fill_value=999999)
>>> np.ma.concatenate([a, b])
masked_array(data=[0, --, 2, 2, 3, 4],
             mask=[False,  True, False, False, False, False],
       fill_value=999999)

>>> x = np.arange(9).reshape((3,3))
>>> x
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

>>> np.diag(x)
array([0, 4, 8])
>>> np.diag(x, k=1)
array([1, 5])
>>> np.diag(x, k=-1)
array([3, 7])

>>> np.diag(np.diag(x))
array([[0, 0, 0],
       [0, 4, 0],
       [0, 0, 8]])

>>> u8_arr = np.array([1, 0], dtype=np.uint8)
>>> np.diff(u8_arr)
array([255], dtype=uint8)
>>> u8_arr[1,...] - u8_arr[0,...]
array(255, np.uint8)

>>> i16_arr = u8_arr.astype(np.int16)
>>> np.diff(i16_arr)
array([-1], dtype=int16)

>>> x = np.array([1, 2, 4, 7, 0])
>>> np.diff(x)
array([ 1,  2,  3, -7])
>>> np.diff(x, n=2)
array([  1,   1, -10])

>>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]])
>>> np.diff(x)
array([[2, 3, 4],
       [5, 1, 2]])
>>> np.diff(x, axis=0)
array([[-1,  2,  0, -2]])

>>> x = np.arange('1066-10-13', '1066-10-16', dtype=np.datetime64)
>>> np.diff(x)
array([1, 1], dtype='timedelta64[D]')

dot(a, b)[i,j,k,m] = sum(a[i,j,:] * b[k,:,m])

>>> np.dot(3, 4)
12

>>> np.dot([2j, 3j], [2j, 3j])
(-13+0j)

>>> a = [[1, 0], [0, 1]]
>>> b = [[4, 1], [2, 2]]
>>> np.dot(a, b)
array([[4, 1],
       [2, 2]])

>>> a = np.arange(3*4*5*6).reshape((3,4,5,6))
>>> b = np.arange(3*4*5*6)[::-1].reshape((5,4,6,3))
>>> np.dot(a, b)[2,3,2,1,2,2]
499128
>>> sum(a[2,3,2,:] * b[1,2,:,2])
499128

>>> np.empty([2, 2])
array([[ -9.74499359e+001,   6.69583040e-309],
       [  2.13182611e-314,   3.06959433e-309]])         #random

>>> np.empty([2, 2], dtype=int)
array([[-1073741821, -1067949133],
       [  496041986,    19249760]])                     #random

>>> np.eye(2, dtype=int)
array([[1, 0],
       [0, 1]])
>>> np.eye(3, k=1)
array([[ 0.,  1.,  0.],
       [ 0.,  0.,  1.],
       [ 0.,  0.,  0.]])

>>> a = np.random.randn(9, 6)
>>> B = np.linalg.pinv(a)
>>> np.allclose(a, np.dot(a, np.dot(B, a)))
True
>>> np.allclose(B, np.dot(B, np.dot(a, B)))
True

>>> np.random.randint(2, size=10)
array([1, 0, 0, 0, 1, 1, 0, 0, 1, 0])
>>> np.random.randint(1, size=10)
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

>>> np.random.randint(5, size=(2, 4))
array([[4, 0, 2, 1],
       [3, 2, 2, 0]])

Real(:math:`Y^m_n`)       if m > 0
:math:`Y^0_n`             if m = 0
Imag(:math:`Y^|m|_n`)     if m < 0

>>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
>>> b = np.random.randn(2, 7, 8, 3) + 1j*np.random.randn(2, 7, 8, 3)

>>> u, s, vh = np.linalg.svd(a, full_matrices=True)
>>> u.shape, s.shape, vh.shape
((9, 9), (6,), (6, 6))
>>> np.allclose(a, np.dot(u[:, :6] * s, vh))
True
>>> smat = np.zeros((9, 6), dtype=complex)
>>> smat[:6, :6] = np.diag(s)
>>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
True

>>> u, s, vh = np.linalg.svd(a, full_matrices=False)
>>> u.shape, s.shape, vh.shape
((9, 6), (6,), (6, 6))
>>> np.allclose(a, np.dot(u * s, vh))
True
>>> smat = np.diag(s)
>>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
True

>>> u, s, vh = np.linalg.svd(b, full_matrices=True)
>>> u.shape, s.shape, vh.shape
((2, 7, 8, 8), (2, 7, 3), (2, 7, 3, 3))
>>> np.allclose(b, np.matmul(u[..., :3] * s[..., None, :], vh))
True
>>> np.allclose(b, np.matmul(u[..., :3], s[..., None] * vh))
True

>>> u, s, vh = np.linalg.svd(b, full_matrices=False)
>>> u.shape, s.shape, vh.shape
((2, 7, 8, 3), (2, 7, 3), (2, 7, 3, 3))
>>> np.allclose(b, np.matmul(u * s[..., None, :], vh))
True
>>> np.allclose(b, np.matmul(u, s[..., None] * vh))
True

>>> np.unique([1, 1, 2, 2, 3, 3])
array([1, 2, 3])
>>> a = np.array([[1, 1], [2, 3]])
>>> np.unique(a)
array([1, 2, 3])

>>> a = np.array([[1, 0, 0], [1, 0, 0], [2, 3, 4]])
>>> np.unique(a, axis=0)
array([[1, 0, 0], [2, 3, 4]])

>>> a = np.array(['a', 'b', 'b', 'c', 'a'])
>>> u, indices = np.unique(a, return_index=True)
>>> u
array(['a', 'b', 'c'],
       dtype='|S1')
>>> indices
array([0, 1, 3])
>>> a[indices]
array(['a', 'b', 'c'],
       dtype='|S1')

>>> a = np.array([1, 2, 6, 4, 2, 3, 2])
>>> u, indices = np.unique(a, return_inverse=True)
>>> u
array([1, 2, 3, 4, 6])
>>> indices
array([0, 1, 4, 3, 1, 2, 1])
>>> u[indices]
array([1, 2, 6, 4, 2, 3, 2])

class MyModel(Model, Cache):
    pass

class MyModelFit(Fit):
    pass

def odf(sphere):
    M = self.model.cache_get('odf_basis_matrix', key=sphere)

    if M is None:
        M = self._compute_basis_matrix(sphere)
        self.model.cache_set('odf_basis_matrix', key=sphere, value=M)

>>> def compute_expensive_matrix(parameters):
...     # Imagine the following computation is very expensive
...     return (p**2 for p in parameters)

>>> c = Cache()

>>> parameters = (1, 2, 3)
>>> X1 = compute_expensive_matrix(parameters)

>>> c.cache_set('expensive_matrix', parameters, X1)
>>> X2 = c.cache_get('expensive_matrix', parameters)

>>> X1 is X2
True

>>> from dipy.utils.optpkg import optional_package
>>> pkg, have_pkg, setup_module = optional_package('not_a_package')

>>> have_pkg
False

>>> pkg.some_function() #doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
    ...
TripWireError: We need package not_a_package for these functions, but
``import not_a_package`` raised an ImportError

>>> pkg, _, _ = optional_package('os')
>>> hasattr(pkg, 'path')
True

>>> subpkg, _, _ = optional_package('os.path')
>>> hasattr(subpkg, 'dirname')
True

:nowrap:
    \begin{equation}
        \textbf{E}(q\textbf{u})=\sum_{l=0, even}^{N_{max}}
                                \sum_{n=l}^{(N_{max}+l)/2}
                                \sum_{m=-l}^l c_{nlm}
                                \phi_{nlm}(q\textbf{u})
    \end{equation}

:nowrap:
    \begin{equation}
        \phi_{nlm}^{SHORE}(q\textbf{u})=\Biggl[\dfrac{2(n-l)!}
            {\zeta^{3/2} \Gamma(n+3/2)} \Biggr]^{1/2}
            \Biggl(\dfrac{q^2}{\zeta}\Biggr)^{l/2}
            exp\Biggl(\dfrac{-q^2}{2\zeta}\Biggr)
            L^{l+1/2}_{n-l} \Biggl(\dfrac{q^2}{\zeta}\Biggr)
            Y_l^m(\textbf{u}).
    \end{equation}

:nowrap:
    \begin{equation}
        \phi_{nlm}^{SHORE}(q\textbf{u})=\Biggl[\dfrac{2(n-l)!}
            {\zeta^{3/2} \Gamma(n+3/2)} \Biggr]^{1/2}
            \Biggl(\dfrac{q^2}{\zeta}\Biggr)^{l/2}
            exp\Biggl(\dfrac{-q^2}{2\zeta}\Biggr)
            L^{l+1/2}_{n-l} \Biggl(\dfrac{q^2}{\zeta}\Biggr)
            Y_l^m(\textbf{u}).
    \end{equation}

`bench`([label, verbose, extra_argv])	Run benchmarks for module using nose.
`test`([label, verbose, extra_argv, doctests, …])	Run tests for module using nose.

`ReconstFit`(model, data)	Abstract class which holds the fit result of ReconstModel
`ReconstModel`(gtab)	Abstract class for signal reconstruction models

`bench_bounding_box`()
`bounding_box`(vol)	Compute the bounding box of nonzero intensity voxels in the volume.
`measure`(code_str[, times, label])	Return elapsed time for executing code in the namespace of the caller.

`ConstrainedSphericalDeconvModel`(gtab, response)	Methods
`GradientTable`(gradients[, big_delta, …])	Diffusion gradient information
`bench_csdeconv`([center, width])
`num_grad`(gtab)
`read_stanford_labels`()	Read stanford hardi data and label map

`bench_local_maxima`()
`local_maxima`	Local maxima of a function evaluated on a discrete set of points.
`measure`(code_str[, times, label])	Return elapsed time for executing code in the namespace of the caller.
`unique_edges`(faces[, return_mapping])	Extract all unique edges from given triangular faces.

`bench_quick_squash`()
`measure`(code_str[, times, label])	Return elapsed time for executing code in the namespace of the caller.
`ndindex`(shape)	An N-dimensional iterator object to index arrays.
`old_squash`(arr[, mask, fill])	Try and make a standard array from an object array
`quick_squash`	Try and make a standard array from an object array
`reduce`(function, sequence[, initial])	Apply a function of two arguments cumulatively to the items of a sequence, from left to right, so as to reduce the sequence to a single value.

`bench_vec_val_vect`()
`measure`(code_str[, times, label])	Return elapsed time for executing code in the namespace of the caller.
`randn`(d0, d1, …, dn)	Return a sample (or samples) from the “standard normal” distribution.
`vec_val_vect`	Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs
`with_einsum`(f)

`Cache`	Cache values based on a key object (such as a sphere or gradient table).
`auto_attr`(func)	Decorator to create OneTimeProperty attributes.

`coeff_of_determination`(data, model[, axis])	Calculate the coefficient of determination for a model prediction, relative
`kfold_xval`(model, data, folds, *model_args, …)	Perform k-fold cross-validation to generate out-of-sample predictions for each measurement.

`AxSymShResponse`(S0, dwi_response[, bvalue])	A simple wrapper for response functions represented using only axially symmetric, even spherical harmonic functions (ie, m == 0 and n even).
`ConstrainedSDTModel`(gtab, ratio[, …])	Methods
`ConstrainedSphericalDeconvModel`(gtab, response)	Methods
`SphHarmFit`(model, shm_coef, mask)	Diffusion data fit to a spherical harmonic model
`SphHarmModel`(gtab)	To be subclassed by all models that return a SphHarmFit when fit.
`TensorModel`(gtab[, fit_method, return_S0_hat])	Diffusion Tensor
`auto_response`(gtab, data[, roi_center, …])	Automatic estimation of response function using FA.
`cart2sphere`(x, y, z)	Return angles for Cartesian 3D coordinates x, y, and z
`csdeconv`(dwsignal, X, B_reg[, tau, …])	Constrained-regularized spherical deconvolution (CSD) [1]
`estimate_response`(gtab, evals, S0)	Estimate single fiber response function
`fa_inferior`(FA, fa_thr)	Check that the FA is lower than the FA threshold
`fa_superior`(FA, fa_thr)	Check that the FA is greater than the FA threshold
`fa_trace_to_lambdas`([fa, trace])
`forward_sdeconv_mat`(r_rh, n)	Build forward spherical deconvolution matrix
`forward_sdt_deconv_mat`(ratio, n[, r2_term])	Build forward sharpening deconvolution transform (SDT) matrix
`fractional_anisotropy`(evals[, axis])	Fractional anisotropy (FA) of a diffusion tensor.
`get_sphere`([name])	provide triangulated spheres
`lazy_index`(index)	Produces a lazy index
`lpn`(n, z)	Legendre function of the first kind.
`multi_voxel_fit`(single_voxel_fit)	Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition
`ndindex`(shape)	An N-dimensional iterator object to index arrays.
`odf_deconv`(odf_sh, R, B_reg[, lambda_, tau, …])	ODF constrained-regularized spherical deconvolution using the Sharpening Deconvolution Transform (SDT) [1], [2].
`odf_sh_to_sharp`(odfs_sh, sphere[, basis, …])	Sharpen odfs using the sharpening deconvolution transform [2]
`peaks_from_model`(model, data, sphere, …[, …])	Fit the model to data and computes peaks and metrics
`quad`(func, a, b[, args, full_output, …])	Compute a definite integral.
`real_sph_harm`(m, n, theta, phi)	Compute real spherical harmonics.
`real_sym_sh_basis`(sh_order, theta, phi)	Samples a real symmetric spherical harmonic basis at point on the sphere
`recursive_response`(gtab, data[, mask, …])	Recursive calibration of response function using peak threshold
`response_from_mask`(gtab, data, mask)	Estimate the response function from a given mask.
`sh_to_rh`(r_sh, m, n)	Spherical harmonics (SH) to rotational harmonics (RH)
`single_tensor`(gtab[, S0, evals, evecs, snr])	Simulated Q-space signal with a single tensor.
`sph_harm_ind_list`(sh_order)	Returns the degree (n) and order (m) of all the symmetric spherical harmonics of degree less then or equal to sh_order.
`vec2vec_rotmat`(u, v)	rotation matrix from 2 unit vectors

`DiffusionKurtosisFit`(model, model_params)	Class for fitting the Diffusion Kurtosis Model
`DiffusionKurtosisModel`(gtab[, fit_method])	Class for the Diffusion Kurtosis Model
`ReconstModel`(gtab)	Abstract class for signal reconstruction models
`TensorFit`(model, model_params[, model_S0])	Attributes
`Wcons`(k_elements)	Construct the full 4D kurtosis tensors from its 15 independent elements
`Wrotate`(kt, Basis)	Rotate a kurtosis tensor from the standard Cartesian coordinate system to another coordinate system basis
`Wrotate_element`(kt, indi, indj, indk, indl, B)	Computes the the specified index element of a kurtosis tensor rotated to the coordinate system basis B.
`apparent_kurtosis_coef`(dki_params, sphere[, …])	Calculates the apparent kurtosis coefficient (AKC) in each direction of a sphere [1].
`axial_kurtosis`(dki_params[, min_kurtosis, …])	Computes axial Kurtosis (AK) from the kurtosis tensor.
`carlson_rd`(x, y, z[, errtol])	Computes the Carlson’s incomplete elliptic integral of the second kind defined as:
`carlson_rf`(x, y, z[, errtol])	Computes the Carlson’s incomplete elliptic integral of the first kind defined as:
`cart2sphere`(x, y, z)	Return angles for Cartesian 3D coordinates x, y, and z
`check_multi_b`(gtab, n_bvals[, non_zero, bmag])	Check if you have enough different b-values in your gradient table
`decompose_tensor`(tensor[, min_diffusivity])	Returns eigenvalues and eigenvectors given a diffusion tensor
`design_matrix`(gtab)	Construct B design matrix for DKI.
`directional_diffusion`(dt, V[, min_diffusivity])	Calculates the apparent diffusion coefficient (adc) in each direction of a sphere for a single voxel [1].
`directional_diffusion_variance`(kt, V[, …])	Calculates the apparent diffusion variance (adv) in each direction of a sphere for a single voxel [1].
`directional_kurtosis`(dt, md, kt, V[, …])	Calculates the apparent kurtosis coefficient (akc) in each direction of a sphere for a single voxel [1].
`dki_prediction`(dki_params, gtab[, S0])	Predict a signal given diffusion kurtosis imaging parameters.
`from_lower_triangular`(D)	Returns a tensor given the six unique tensor elements
`get_sphere`([name])	provide triangulated spheres
`kurtosis_maximum`(dki_params[, sphere, gtol, …])	Computes kurtosis maximum value
`local_maxima`	Local maxima of a function evaluated on a discrete set of points.
`lower_triangular`(tensor[, b0])	Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None
`mean_diffusivity`(evals[, axis])	Mean Diffusivity (MD) of a diffusion tensor.
`mean_kurtosis`(dki_params[, min_kurtosis, …])	Computes mean Kurtosis (MK) from the kurtosis tensor [1].
`ndindex`(shape)	An N-dimensional iterator object to index arrays.
`ols_fit_dki`(design_matrix, data)	Computes ordinary least squares (OLS) fit to calculate the diffusion tensor and kurtosis tensor using a linear regression diffusion kurtosis model 1.
`radial_kurtosis`(dki_params[, min_kurtosis, …])	Radial Kurtosis (RK) of a diffusion kurtosis tensor [1].
`sphere2cart`(r, theta, phi)	Spherical to Cartesian coordinates
`split_dki_param`(dki_params)	Extract the diffusion tensor eigenvalues, the diffusion tensor eigenvector matrix, and the 15 independent elements of the kurtosis tensor from the model parameters estimated from the DKI model
`vec_val_vect`	Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs
`wls_fit_dki`(design_matrix, data)	Computes weighted linear least squares (WLS) fit to calculate the diffusion tensor and kurtosis tensor using a weighted linear regression diffusion kurtosis model 1.

`Cache`	Cache values based on a key object (such as a sphere or gradient table).
`DiffusionSpectrumDeconvFit`(model, data)	Methods
`DiffusionSpectrumDeconvModel`(gtab[, …])	Methods
`DiffusionSpectrumFit`(model, data)	Methods
`DiffusionSpectrumModel`(gtab[, qgrid_size, …])	Methods
`OdfFit`(model, data)	Methods
`OdfModel`(gtab)	An abstract class to be sub-classed by specific odf models
`LR_deconv`(prop, psf[, numit, acc_factor])	Perform Lucy-Richardson deconvolution algorithm on a 3D array.
`create_qspace`(gtab, origin)	create the 3D grid which holds the signal values (q-space)
`create_qtable`(gtab, origin)	create a normalized version of gradients
`fftn`(x[, shape, axes, overwrite_x])	Return multidimensional discrete Fourier transform.
`fftshift`(x[, axes])	Shift the zero-frequency component to the center of the spectrum.
`gen_PSF`(qgrid_sampling, siz_x, siz_y, siz_z)	Generate a PSF for DSI Deconvolution by taking the ifft of the binary q-space sampling mask and truncating it to keep only the center.
`half_to_full_qspace`(data, gtab)	Half to full Cartesian grid mapping
`hanning_filter`(gtab, filter_width, origin)	create a hanning window
`ifftshift`(x[, axes])	The inverse of fftshift.
`map_coordinates`(input, coordinates[, …])	Map the input array to new coordinates by interpolation.
`multi_voxel_fit`(single_voxel_fit)	Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition
`pdf_interp_coords`(sphere, rradius, origin)	Precompute coordinates for ODF calculation from the PDF
`pdf_odf`(Pr, rradius, interp_coords)	Calculates the real ODF from the diffusion propagator(PDF) Pr
`project_hemisph_bvecs`(gtab)	Project any near identical bvecs to the other hemisphere
`threshold_propagator`(P[, estimated_snr])	Applies hard threshold on the propagator to remove background noise for the deconvolution.

`Cache`	Cache values based on a key object (such as a sphere or gradient table).
`ForecastFit`(model, data, sh_coef, d_par, d_perp)	Attributes
`ForecastModel`(gtab[, sh_order, lambda_lb, …])	Fiber ORientation Estimated using Continuous Axially Symmetric Tensors (FORECAST) [1,2,3]_.
`OdfFit`(model, data)	Methods
`OdfModel`(gtab)	An abstract class to be sub-classed by specific odf models
`cart2sphere`(x, y, z)	Return angles for Cartesian 3D coordinates x, y, and z
`csdeconv`(dwsignal, X, B_reg[, tau, …])	Constrained-regularized spherical deconvolution (CSD) [1]
`find_signal_means`(b_unique, data_norm, …)	Calculate the mean signal for each shell.
`forecast_error_func`(x, b_unique, E)	Calculates the difference between the mean signal calculated using the parameter vector x and the average signal E using FORECAST and SMT
`forecast_matrix`(sh_order, d_par, d_perp, bvals)	Compute the FORECAST radial matrix
`lb_forecast`(sh_order)	Returns the Laplace-Beltrami regularization matrix for FORECAST
`leastsq`(func, x0[, args, Dfun, full_output, …])	Minimize the sum of squares of a set of equations.
`multi_voxel_fit`(single_voxel_fit)	Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition
`optional_package`(name[, trip_msg])	Return package-like thing and module setup for package name
`psi_l`(l, b)
`real_sph_harm`(m, n, theta, phi)	Compute real spherical harmonics.
`rho_matrix`(sh_order, vecs)	Compute the SH matrix \(\rho\)
`warn`	Issue a warning, or maybe ignore it or raise an exception.

`FreeWaterTensorFit`(model, model_params)	Class for fitting the Free Water Tensor Model
`FreeWaterTensorModel`(gtab[, fit_method])	Class for the Free Water Elimination Diffusion Tensor Model
`ReconstModel`(gtab)	Abstract class for signal reconstruction models
`TensorFit`(model, model_params[, model_S0])	Attributes
`check_multi_b`(gtab, n_bvals[, non_zero, bmag])	Check if you have enough different b-values in your gradient table
`cholesky_to_lower_triangular`(R)	Convert Cholesky decompostion elements to the diffusion tensor elements
`decompose_tensor`(tensor[, min_diffusivity])	Returns eigenvalues and eigenvectors given a diffusion tensor
`design_matrix`(gtab[, dtype])	Constructs design matrix for DTI weighted least squares or least squares fitting.
`from_lower_triangular`(D)	Returns a tensor given the six unique tensor elements
`fwdti_prediction`(params, gtab[, S0, Diso])	Signal prediction given the free water DTI model parameters.
`lower_triangular`(tensor[, b0])	Returns the six lower triangular values of the tensor and a dummy variable if b0 is not None
`lower_triangular_to_cholesky`(tensor_elements)	Performs Cholesky decomposition of the diffusion tensor
`multi_voxel_fit`(single_voxel_fit)	Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition
`ndindex`(shape)	An N-dimensional iterator object to index arrays.
`nls_fit_tensor`(gtab, data[, mask, Diso, …])	Fit the water elimination tensor model using the non-linear least-squares.
`nls_iter`(design_matrix, sig, S0[, Diso, …])	Applies non linear least squares fit of the water free elimination model to single voxel signals.
`vec_val_vect`	Vectorize vecs.diag(vals).`vecs`.T for last 2 dimensions of vecs
`wls_fit_tensor`(gtab, data[, Diso, mask, …])	Computes weighted least squares (WLS) fit to calculate self-diffusion tensor using a linear regression model [1].
`wls_iter`(design_matrix, sig, S0[, Diso, …])	Applies weighted linear least squares fit of the water free elimination model to single voxel signals.

`Cache`	Cache values based on a key object (such as a sphere or gradient table).
`GeneralizedQSamplingFit`(model, data)	Methods
`GeneralizedQSamplingModel`(gtab[, method, …])	Methods
`OdfFit`(model, data)	Methods
`OdfModel`(gtab)	An abstract class to be sub-classed by specific odf models
`equatorial_maximum`(vertices, odf, pole, width)
`equatorial_zone_vertices`(vertices, pole[, width])	finds the ‘vertices’ in the equatorial zone conjugate to ‘pole’ with width half ‘width’ degrees
`gfa`(samples)	The general fractional anisotropy of a function evaluated on the unit sphere
`local_maxima`	Local maxima of a function evaluated on a discrete set of points.
`multi_voxel_fit`(single_voxel_fit)	Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition
`normalize_qa`(qa[, max_qa])	Normalize quantitative anisotropy.
`npa`(self, odf[, width])	non-parametric anisotropy
`odf_sum`(odf)
`patch_maximum`(vertices, odf, pole, width)
`patch_sum`(vertices, odf, pole, width)
`patch_vertices`(vertices, pole, width)	find ‘vertices’ within the cone of ‘width’ degrees around ‘pole’
`polar_zone_vertices`(vertices, pole[, width])	finds the ‘vertices’ in the equatorial band around the ‘pole’ of radius ‘width’ degrees
`remove_similar_vertices`	Remove vertices that are less than theta degrees from any other
`squared_radial_component`(x[, tol])	Part of the GQI2 integral
`triple_odf_maxima`(vertices, odf, width)
`upper_hemi_map`(v)	maps a 3-vector into the z-upper hemisphere

`cache_clear`()	Clear the cache.
`cache_get`(tag, key[, default])	Retrieve a value from the cache.
`cache_set`(tag, key, value)	Store a value in the cache.
`fit`(data[, mask])	Fit method for every voxel in data
`predict`(sh_coeff[, gtab, S0])	Compute a signal prediction given spherical harmonic coefficients for the provided GradientTable class instance.
`sampling_matrix`(sphere)	The matrix needed to sample ODFs from coefficients of the model.

b0s_mask
bvals
bvecs
gradient_strength
qvals
tau

`cache_clear`()	Clear the cache.
`cache_get`(tag, key[, default])	Retrieve a value from the cache.
`cache_set`(tag, key, value)	Store a value in the cache.

`cache_clear`()	Clear the cache.
`cache_get`(tag, key[, default])	Retrieve a value from the cache.
`cache_set`(tag, key, value)	Store a value in the cache.
`fit`(data[, mask])	Fit method for every voxel in data
`sampling_matrix`(sphere)	The matrix needed to sample ODFs from coefficients of the model.

`IvimFit`(model, model_params)	Attributes
`IvimModelTRR`(gtab[, split_b_D, split_b_S0, …])	Ivim model
`IvimModelVP`(gtab[, bounds, maxiter, xtol])	Methods
`ReconstModel`(gtab)	Abstract class for signal reconstruction models
`IvimModel`(gtab[, fit_method])	Selector function to switch between the 2-stage Trust-Region Reflective based NLLS fitting method (also containing the linear fit): trr and the Variable Projections based fitting method: varpro.
`differential_evolution`(func, bounds[, args, …])	Finds the global minimum of a multivariate function.
`f_D_star_error`(params, gtab, signal, S0, D)	Error function used to fit f and D_star keeping S0 and D fixed
`f_D_star_prediction`(params, gtab, S0, D)	Function used to predict IVIM signal when S0 and D are known by considering f and D_star as the unknown parameters.
`ivim_model_selector`(gtab[, fit_method])	Selector function to switch between the 2-stage Trust-Region Reflective based NLLS fitting method (also containing the linear fit): trr and the Variable Projections based fitting method: varpro.
`ivim_prediction`(params, gtab)	The Intravoxel incoherent motion (IVIM) model function.
`least_squares`(fun, x0[, jac, bounds, …])	Solve a nonlinear least-squares problem with bounds on the variables.
`multi_voxel_fit`(single_voxel_fit)	Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition
`optional_package`(name[, trip_msg])	Return package-like thing and module setup for package name

`Cache`	Cache values based on a key object (such as a sphere or gradient table).
`MapmriFit`(model, mapmri_coef, mu, R, lopt[, …])	Attributes
`MapmriModel`(gtab[, radial_order, …])	Mean Apparent Propagator MRI (MAPMRI) [1] of the diffusion signal.
`Optimizer`(fun, x0[, args, method, jac, …])	Attributes
`ReconstFit`(model, data)	Abstract class which holds the fit result of ReconstModel
`ReconstModel`(gtab)	Abstract class for signal reconstruction models
`b_mat`(index_matrix)	Calculates the B coefficients from [1] Eq.
`b_mat_isotropic`(index_matrix)	Calculates the isotropic B coefficients from [1] Fig 8.
`binomialfloat`(n, k)	Custom Binomial function
`cart2sphere`(x, y, z)	Return angles for Cartesian 3D coordinates x, y, and z
`create_rspace`(gridsize, radius_max)	Create the real space table, that contains the points in which to compute the pdf.
`delta`(n, m)
`factorial2`(n[, exact])	Double factorial.
`gcv_cost_function`(weight, args)	The GCV cost function that is iterated [4].
`generalized_crossvalidation`(data, M, LR[, …])	Generalized Cross Validation Function [Rb690cd738504-1] eq.
`generalized_crossvalidation_array`(data, M, LR)	Generalized Cross Validation Function 1 eq.
`genlaguerre`(n, alpha[, monic])	Generalized (associated) Laguerre polynomial.
`gradient_table`(bvals[, bvecs, big_delta, …])	A general function for creating diffusion MR gradients.
`hermite`(n[, monic])	Physicist’s Hermite polynomial.
`isotropic_scale_factor`(mu_squared)	Estimated isotropic scaling factor _[1] Eq.
`map_laplace_s`(n, m)	R(m,n) static matrix for Laplacian regularization [R932dd40ca52e-1] eq.
`map_laplace_t`(n, m)	L(m, n) static matrix for Laplacian regularization [Reb78d789d6c4-1] eq.
`map_laplace_u`(n, m)	S(n, m) static matrix for Laplacian regularization [Rb93dd9dab8c9-1] eq.
`mapmri_STU_reg_matrices`(radial_order)	Generate the static portions of the Laplacian regularization matrix according to [R1d585103467a-1] eq.
`mapmri_index_matrix`(radial_order)	Calculates the indices for the MAPMRI [1] basis in x, y and z.
`mapmri_isotropic_K_mu_dependent`(…)	Computes mu dependent part of M.
`mapmri_isotropic_K_mu_independent`(…)	Computes mu independent part of K.
`mapmri_isotropic_M_mu_dependent`(…)	Computed the mu dependent part of the signal design matrix.
`mapmri_isotropic_M_mu_independent`(…)	Computed the mu independent part of the signal design matrix.
`mapmri_isotropic_index_matrix`(radial_order)	Calculates the indices for the isotropic MAPMRI basis [1] Fig 8.
`mapmri_isotropic_laplacian_reg_matrix`(…)	Computes the Laplacian regularization matrix for MAP-MRI’s isotropic implementation [R156f27ca005f-1] eq.
`mapmri_isotropic_laplacian_reg_matrix_from_index_matrix`(…)	Computes the Laplacian regularization matrix for MAP-MRI’s isotropic implementation [Rdcc29394f577-1] eq.
`mapmri_isotropic_odf_matrix`(radial_order, …)	Compute the isotropic MAPMRI ODF matrix [1] Eq.
`mapmri_isotropic_odf_sh_matrix`(radial_order, …)	Compute the isotropic MAPMRI ODF matrix [1] Eq.
`mapmri_isotropic_phi_matrix`(radial_order, mu, q)	Three dimensional isotropic MAPMRI signal basis function from [1] Eq.
`mapmri_isotropic_psi_matrix`(radial_order, …)	Three dimensional isotropic MAPMRI propagator basis function from [1] Eq.
`mapmri_isotropic_radial_pdf_basis`(j, l, mu, r)	Radial part of the isotropic 1D-SHORE propagator basis [1] eq.
`mapmri_isotropic_radial_signal_basis`(j, l, …)	Radial part of the isotropic 1D-SHORE signal basis [1] eq.
`mapmri_laplacian_reg_matrix`(ind_mat, mu, …)	Put the Laplacian regularization matrix together [Rc66aaccd07c1-1] eq.
`mapmri_odf_matrix`(radial_order, mu, s, vertices)	Compute the MAPMRI ODF matrix [1] Eq.
`mapmri_phi_1d`(n, q, mu)	One dimensional MAPMRI basis function from [1] Eq.
`mapmri_phi_matrix`(radial_order, mu, q_gradients)	Compute the MAPMRI phi matrix for the signal [1] eq.
`mapmri_psi_1d`(n, x, mu)	One dimensional MAPMRI propagator basis function from [1] Eq.
`mapmri_psi_matrix`(radial_order, mu, rgrad)	Compute the MAPMRI psi matrix for the propagator [1] eq.
`mfactorial`	factorial(x) -> Integral
`multi_voxel_fit`(single_voxel_fit)	Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition
`optional_package`(name[, trip_msg])	Return package-like thing and module setup for package name
`real_sph_harm`(m, n, theta, phi)	Compute real spherical harmonics.
`sfactorial`(n[, exact])	The factorial of a number or array of numbers.
`sph_harm_ind_list`(sh_order)	Returns the degree (n) and order (m) of all the symmetric spherical harmonics of degree less then or equal to sh_order.
`warn`	Issue a warning, or maybe ignore it or raise an exception.

`MSDeconvFit`(model, coeff, mask)	Attributes
`MultiShellDeconvModel`(gtab, response[, …])	Methods
`MultiShellResponse`(response, sh_order, shells)	Attributes
`QpFitter`(X, reg)	Methods
`multi_tissue_basis`(gtab, sh_order, iso_comp)	Builds a basis for multi-shell multi-tissue CSD model.
`multi_voxel_fit`(single_voxel_fit)	Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition
`optional_package`(name[, trip_msg])	Return package-like thing and module setup for package name
`solve_qp`(P, Q, G, H)	Helper function to set up and solve the Quadratic Program (QP) in CVXPY.

`MeanDiffusionKurtosisFit`(model, model_params)	Attributes
`MeanDiffusionKurtosisModel`(gtab[, bmag, …])	Mean signal Diffusion Kurtosis Model
`ReconstModel`(gtab)	Abstract class for signal reconstruction models
`auto_attr`(func)	Decorator to create OneTimeProperty attributes.
`check_multi_b`(gtab, n_bvals[, non_zero, bmag])	Check if you have enough different b-values in your gradient table
`design_matrix`(ubvals)	Constructs design matrix for the mean signal diffusion kurtosis model
`mean_signal_bvalue`(data, gtab[, bmag])	Computes the average signal across different diffusion directions for each unique b-value
`msdki_prediction`(msdki_params, gtab[, S0])	Predict the mean signal given the parameters of the mean signal DKI, an GradientTable object and S0 signal.
`ndindex`(shape)	An N-dimensional iterator object to index arrays.
`round_bvals`(bvals[, bmag])	“This function rounds the b-values
`unique_bvals`(bvals[, bmag, rbvals])	This function gives the unique rounded b-values of the data
`wls_fit_msdki`(design_matrix, msignal, ng[, …])	Fits the mean signal diffusion kurtosis imaging based on a weighted least square solution [1].

`CallableArray`	An array which can be called like a function
`MultiVoxelFit`(model, fit_array, mask)	Holds an array of fits and allows access to their attributes and methods
`ReconstFit`(model, data)	Abstract class which holds the fit result of ReconstModel
`as_strided`(x[, shape, strides, subok, writeable])	Create a view into the array with the given shape and strides.
`multi_voxel_fit`(single_voxel_fit)	Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition
`ndindex`(shape)	An N-dimensional iterator object to index arrays.

`Cache`	Cache values based on a key object (such as a sphere or gradient table).
`QtdmriFit`(model, qtdmri_coef, us, ut, …)	Methods
`QtdmriModel`(gtab[, radial_order, …])	The q:math:tau-dMRI model [1] to analytically and continuously represent the q:math:tau diffusion signal attenuation over diffusion sensitization q and diffusion time \(\tau\).
`GCV_cost_function`(weight, arguments)	Generalized Cross Validation Function that is iterated [1].
`H`(value)	Step function of H(x)=1 if x>=0 and zero otherwise.
`angular_basis_EAP_opt`(j, l, m, r, theta, phi)
`angular_basis_opt`(l, m, q, theta, phi)	Angular basis independent of spatial scaling factor us.
`cart2sphere`(x, y, z)	Return angles for Cartesian 3D coordinates x, y, and z
`create_rt_space_grid`(grid_size_r, …)	Generates EAP grid (for potential positivity constraint).
`design_matrix_spatial`(bvecs, qvals[, dtype])	Constructs design matrix for DTI weighted least squares or least squares fitting.
`elastic_crossvalidation`(b0s_mask, E, M, L, lopt)	cross-validation function to find the optimal weight of alpha for sparsity regularization when also Laplacian regularization is used.
`factorial`(n[, exact])	The factorial of a number or array of numbers.
`factorial2`(n[, exact])	Double factorial.
`fmin_l_bfgs_b`(func, x0[, fprime, args, …])	Minimize a function func using the L-BFGS-B algorithm.
`generalized_crossvalidation`(data, M, LR[, …])	Generalized Cross Validation Function [1].
`genlaguerre`(n, alpha[, monic])	Generalized (associated) Laguerre polynomial.
`gradient_table_from_gradient_strength_bvecs`(…)	A general function for creating diffusion MR gradients.
`l1_crossvalidation`(b0s_mask, E, M[, …])	cross-validation function to find the optimal weight of alpha for sparsity regularization
`multi_voxel_fit`(single_voxel_fit)	Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition
`optional_package`(name[, trip_msg])	Return package-like thing and module setup for package name
`part1_reg_matrix_tau`(ind_mat, ut)	Partial temporal Laplacian regularization matrix following Appendix B in [1].
`part23_iso_reg_matrix_q`(ind_mat, us)	Partial spherical spatial Laplacian regularization matrix following the equation below Eq.
`part23_reg_matrix_q`(ind_mat, U_mat, T_mat, us)	Partial cartesian spatial Laplacian regularization matrix following second line of Eq.
`part23_reg_matrix_tau`(ind_mat, ut)	Partial temporal Laplacian regularization matrix following Appendix B in [1].
`part4_iso_reg_matrix_q`(ind_mat, us)	Partial spherical spatial Laplacian regularization matrix following the equation below Eq.
`part4_reg_matrix_q`(ind_mat, U_mat, us)	Partial cartesian spatial Laplacian regularization matrix following equation Eq.
`part4_reg_matrix_tau`(ind_mat, ut)	Partial temporal Laplacian regularization matrix following Appendix B in [1].
`qtdmri_anisotropic_scaling`(data, q, bvecs, tau)	Constructs design matrix for fitting an exponential to the diffusion time points.
`qtdmri_eap_matrix`(radial_order, time_order, …)	Constructs the design matrix as a product of 3 separated radial, angular and temporal design matrices.
`qtdmri_eap_matrix_`(radial_order, time_order, …)
`qtdmri_index_matrix`(radial_order, time_order)	Computes the SHORE basis order indices according to [1].
`qtdmri_isotropic_eap_matrix`(radial_order, …)	Constructs the design matrix as a product of 3 separated radial, angular and temporal design matrices.
`qtdmri_isotropic_eap_matrix_`(radial_order, …)
`qtdmri_isotropic_index_matrix`(radial_order, …)	Computes the SHORE basis order indices according to [1].
`qtdmri_isotropic_laplacian_reg_matrix`(…[, …])	Computes the spherical qt-dMRI Laplacian regularization matrix.
`qtdmri_isotropic_scaling`(data, q, tau)	Constructs design matrix for fitting an exponential to the diffusion time points.
`qtdmri_isotropic_signal_matrix`(radial_order, …)
`qtdmri_isotropic_signal_matrix_`(…[, …])
`qtdmri_isotropic_to_mapmri_matrix`(…)	Generates the matrix that maps the spherical qtdmri coefficients to MAP-MRI coefficients.
`qtdmri_laplacian_reg_matrix`(ind_mat, us, ut)	Computes the cartesian qt-dMRI Laplacian regularization matrix.
`qtdmri_mapmri_isotropic_normalization`(j, l, u0)	Normalization factor for Spherical MAP-MRI basis.
`qtdmri_mapmri_normalization`(mu)	Normalization factor for Cartesian MAP-MRI basis.
`qtdmri_number_of_coefficients`(radial_order, …)	Computes the total number of coefficients of the qtdmri basis given a radial and temporal order.
`qtdmri_signal_matrix`(radial_order, …)	Constructs the design matrix as a product of 3 separated radial, angular and temporal design matrices.
`qtdmri_signal_matrix_`(radial_order, …[, …])	Function to generate the qtdmri signal basis.
`qtdmri_temporal_normalization`(ut)	Normalization factor for the temporal basis
`qtdmri_to_mapmri_matrix`(radial_order, …)	Generates the matrix that maps the qtdmri coefficients to MAP-MRI coefficients.
`radial_basis_EAP_opt`(j, l, us, r)
`radial_basis_opt`(j, l, us, q)	Spatial basis dependent on spatial scaling factor us
`real_sph_harm`(m, n, theta, phi)	Compute real spherical harmonics.
`temporal_basis`(o, ut, tau)	Temporal basis dependent on temporal scaling factor ut
`visualise_gradient_table_G_Delta_rainbow`(gtab)	This function visualizes a q-tau acquisition scheme as a function of gradient strength and pulse separation (big_delta).
`warn`	Issue a warning, or maybe ignore it or raise an exception.

`Cache`	Cache values based on a key object (such as a sphere or gradient table).
`ExponentialIsotropicFit`(model, params)	A fit to the ExponentialIsotropicModel object, based on data.
`ExponentialIsotropicModel`(gtab)	Representing the isotropic signal as a fit to an exponential decay function with b-values
`IsotropicFit`(model, params)	A fit object for representing the isotropic signal as the mean of the diffusion-weighted signal.
`IsotropicModel`(gtab)	A base-class for the representation of isotropic signals.
`ReconstFit`(model, data)	Abstract class which holds the fit result of ReconstModel
`ReconstModel`(gtab)	Abstract class for signal reconstruction models
`SparseFascicleFit`(model, beta, S0, iso)	Methods
`SparseFascicleModel`(gtab[, sphere, …])	Methods
`auto_attr`(func)	Decorator to create OneTimeProperty attributes.
`nanmean`(a[, axis, dtype, out, keepdims])	Compute the arithmetic mean along the specified axis, ignoring NaNs.
`optional_package`(name[, trip_msg])	Return package-like thing and module setup for package name
`sfm_design_matrix`(gtab, sphere, response[, mode])	Construct the SFM design matrix

`Cache`	Cache values based on a key object (such as a sphere or gradient table).
`CsaOdfModel`(gtab, sh_order[, smooth, …])	Implementation of Constant Solid Angle reconstruction method.
`LooseVersion`([vstring])	Version numbering for anarchists and software realists.
`OdfFit`(model, data)	Methods
`OdfModel`(gtab)	An abstract class to be sub-classed by specific odf models
`OpdtModel`(gtab, sh_order[, smooth, …])	Implementation of Orientation Probability Density Transform reconstruction method.
`QballBaseModel`(gtab, sh_order[, smooth, …])	To be subclassed by Qball type models.
`QballModel`(gtab, sh_order[, smooth, …])	Implementation of regularized Qball reconstruction method.
`ResidualBootstrapWrapper`(signal_object, B, …)	Returns a residual bootstrap sample of the signal_object when indexed
`SphHarmFit`(model, shm_coef, mask)	Diffusion data fit to a spherical harmonic model
`SphHarmModel`(gtab)	To be subclassed by all models that return a SphHarmFit when fit.
`anisotropic_power`(sh_coeffs[, norm_factor, …])	Calculates anisotropic power map with a given SH coefficient matrix
`auto_attr`(func)	Decorator to create OneTimeProperty attributes.
`bootstrap_data_array`(data, H, R[, permute])	Applies the Residual Bootstraps to the data given H and R
`bootstrap_data_voxel`(data, H, R[, permute])	Like bootstrap_data_array but faster when for a single voxel
`calculate_max_order`(n_coeffs)	Calculate the maximal harmonic order, given that you know the
`cart2sphere`(x, y, z)	Return angles for Cartesian 3D coordinates x, y, and z
`concatenate`([axis, out])	Join a sequence of arrays along an existing axis.
`diag`(v[, k])	Extract a diagonal or construct a diagonal array.
`diff`(a[, n, axis, prepend, append])	Calculate the n-th discrete difference along the given axis.
`dot`(a, b[, out])	Dot product of two arrays.
`empty`(shape[, dtype, order])	Return a new array of given shape and type, without initializing entries.
`eye`(N[, M, k, dtype, order])	Return a 2-D array with ones on the diagonal and zeros elsewhere.
`forward_sdeconv_mat`(r_rh, n)	Build forward spherical deconvolution matrix
`gen_dirac`(m, n, theta, phi)	Generate Dirac delta function orientated in (theta, phi) on the sphere
`hat`(B)	Returns the hat matrix for the design matrix B
`lazy_index`(index)	Produces a lazy index
`lcr_matrix`(H)	Returns a matrix for computing leveraged, centered residuals from data
`lpn`(n, z)	Legendre function of the first kind.
`normalize_data`(data, where_b0[, min_signal, out])	Normalizes the data with respect to the mean b0
`order_from_ncoef`(ncoef)	Given a number n of coefficients, calculate back the sh_order
`pinv`(a[, rcond])	Compute the (Moore-Penrose) pseudo-inverse of a matrix.
`randint`(low[, high, size, dtype])	Return random integers from low (inclusive) to high (exclusive).
`real_sph_harm`(m, n, theta, phi)	Compute real spherical harmonics.
`real_sym_sh_basis`(sh_order, theta, phi)	Samples a real symmetric spherical harmonic basis at point on the sphere
`real_sym_sh_mrtrix`(sh_order, theta, phi)	Compute real spherical harmonics as in Tournier 2007 [2], where the real harmonic \(Y^m_n\) is defined to be.
`sf_to_sh`(sf, sphere[, sh_order, basis_type, …])	Spherical function to spherical harmonics (SH).
`sh_to_rh`(r_sh, m, n)	Spherical harmonics (SH) to rotational harmonics (RH)
`sh_to_sf`(sh, sphere, sh_order[, basis_type])	Spherical harmonics (SH) to spherical function (SF).
`sh_to_sf_matrix`(sphere, sh_order[, …])	Matrix that transforms Spherical harmonics (SH) to spherical function (SF).
`smooth_pinv`(B, L)	Regularized pseudo-inverse
`sph_harm_ind_list`(sh_order)	Returns the degree (n) and order (m) of all the symmetric spherical harmonics of degree less then or equal to sh_order.
`spherical_harmonics`(m, n, theta, phi)	Compute spherical harmonics
`svd`(a[, full_matrices, compute_uv])	Singular Value Decomposition.
`unique`(ar[, return_index, return_inverse, …])	Find the unique elements of an array.

`Cache`	Cache values based on a key object (such as a sphere or gradient table).
`ShoreFit`(model, shore_coef)	Attributes
`ShoreModel`(gtab[, radial_order, zeta, …])	Simple Harmonic Oscillator based Reconstruction and Estimation (SHORE) [1] of the diffusion signal.
`cart2sphere`(x, y, z)	Return angles for Cartesian 3D coordinates x, y, and z
`create_rspace`(gridsize, radius_max)	Create the real space table, that contains the points in which
`factorial`(x)	Find x!.
`genlaguerre`(n, alpha[, monic])	Generalized (associated) Laguerre polynomial.
`l_shore`(radial_order)	Returns the angular regularisation matrix for SHORE basis
`multi_voxel_fit`(single_voxel_fit)	Method decorator to turn a single voxel model fit definition into a multi voxel model fit definition
`n_shore`(radial_order)	Returns the angular regularisation matrix for SHORE basis
`optional_package`(name[, trip_msg])	Return package-like thing and module setup for package name
`real_sph_harm`(m, n, theta, phi)	Compute real spherical harmonics.
`shore_indices`(radial_order, index)	Given the basis order and the index, return the shore indices n, l, m for modified Merlet’s 3D-SHORE ..math:: :nowrap: begin{equation} textbf{E}(qtextbf{u})=sum_{l=0, even}^{N_{max}} sum_{n=l}^{(N_{max}+l)/2} sum_{m=-l}^l c_{nlm} phi_{nlm}(qtextbf{u}) end{equation}
`shore_matrix`(radial_order, zeta, gtab[, tau])	Compute the SHORE matrix for modified Merlet’s 3D-SHORE [1]
`shore_matrix_odf`(radial_order, zeta, …)	Compute the SHORE ODF matrix [1]”
`shore_matrix_pdf`(radial_order, zeta, rtab)	Compute the SHORE propagator matrix [1]”
`shore_order`(n, l, m)	Given the indices (n,l,m) of the basis, return the minimum order for those indices and their index for modified Merlet’s 3D-SHORE.
`warn`	Issue a warning, or maybe ignore it or raise an exception.

`basis`(sphere)	A basis that maps the response coefficients onto a sphere.
`on_sphere`(sphere)	Evaluates the response function on sphere.

reconst¶

Module: reconst.base¶

Module: reconst.benchmarks¶

Module: reconst.benchmarks.bench_bounding_box¶

Module: reconst.benchmarks.bench_csd¶

Module: reconst.benchmarks.bench_peaks¶

Module: reconst.benchmarks.bench_squash¶

Module: reconst.benchmarks.bench_vec_val_sum¶

Module: reconst.cache¶

Module: reconst.cross_validation¶

Module: reconst.csdeconv¶

Module: reconst.dki¶

Module: reconst.dki_micro¶

Module: reconst.dsi¶

Module: reconst.dti¶

Module: reconst.forecast¶

Module: reconst.fwdti¶

Module: reconst.gqi¶

Module: reconst.ivim¶

Module: reconst.mapmri¶

Module: reconst.mcsd¶

Module: reconst.msdki¶

Module: reconst.multi_voxel¶

Module: reconst.odf¶

Module: reconst.qtdmri¶

Module: reconst.sfm¶

Module: reconst.shm¶

References¶

Module: reconst.shore¶

Module: reconst.utils¶

bench¶

test¶

bench_bounding_box¶

bounding_box¶

measure¶

bench_csdeconv¶

num_grad¶

read_stanford_labels¶

bench_local_maxima¶

local_maxima¶

measure¶

unique_edges¶

bench_quick_squash¶

measure¶

ndindex¶

old_squash¶

quick_squash¶

reduce¶

bench_vec_val_vect¶

measure¶

randn¶

vec_val_vect¶

with_einsum¶

auto_attr¶

coeff_of_determination¶

kfold_xval¶

auto_response¶

cart2sphere¶

csdeconv¶

estimate_response¶

fa_inferior¶

fa_superior¶

fa_trace_to_lambdas¶

forward_sdeconv_mat¶

forward_sdt_deconv_mat¶

fractional_anisotropy¶

get_sphere¶

lazy_index¶

lpn¶

multi_voxel_fit¶

ndindex¶

odf_deconv¶

odf_sh_to_sharp¶

peaks_from_model¶

quad¶

real_sph_harm¶

real_sym_sh_basis¶

recursive_response¶

response_from_mask¶

sh_to_rh¶

`reconst`¶

Module: `reconst.base`¶

Module: `reconst.benchmarks`¶

Module: `reconst.benchmarks.bench_bounding_box`¶

Module: `reconst.benchmarks.bench_csd`¶

Module: `reconst.benchmarks.bench_peaks`¶

Module: `reconst.benchmarks.bench_squash`¶

Module: `reconst.benchmarks.bench_vec_val_sum`¶

Module: `reconst.cache`¶

Module: `reconst.cross_validation`¶

Module: `reconst.csdeconv`¶

Module: `reconst.dki`¶

Module: `reconst.dki_micro`¶

Module: `reconst.dsi`¶

Module: `reconst.dti`¶

Module: `reconst.forecast`¶

Module: `reconst.fwdti`¶

Module: `reconst.gqi`¶

Module: `reconst.ivim`¶

Module: `reconst.mapmri`¶

Module: `reconst.mcsd`¶

Module: `reconst.msdki`¶

Module: `reconst.multi_voxel`¶

Module: `reconst.odf`¶

Module: `reconst.qtdmri`¶

Module: `reconst.sfm`¶

Module: `reconst.shm`¶

Module: `reconst.shore`¶

Module: `reconst.utils`¶

`odf`(sphere)	Samples the odf function on the points of a sphere
`predict`([gtab, S0])	Predict the diffusion signal from the model coefficients.

`fit`(data[, mask])	Fit method of the DTI model class
`predict`(dti_params[, S0])	Predict a signal for this TensorModel class instance given parameters.

`weight`	Weight function used	`wvar`
‘cos’	cos(w*x)	wvar = w
‘sin’	sin(w*x)	wvar = w
‘alg’	g(x) = ((x-a)*alpha)((b-x)**beta)	wvar = (alpha, beta)
‘alg-loga’	g(x)*log(x-a)	wvar = (alpha, beta)
‘alg-logb’	g(x)*log(b-x)	wvar = (alpha, beta)
‘alg-log’	g(x)log(x-a)log(b-x)	wvar = (alpha, beta)
‘cauchy’	1/(x-c)	wvar = c

`ad`()	Axial diffusivity (AD) calculated from cached eigenvalues.
`adc`(sphere)	Calculate the apparent diffusion coefficient (ADC) in each direction on
`ak`([min_kurtosis, max_kurtosis])	Axial Kurtosis (AK) of a diffusion kurtosis tensor [1].
`akc`(sphere)	Calculates the apparent kurtosis coefficient (AKC) in each direction on the sphere for each voxel in the data
`color_fa`()	Color fractional anisotropy of diffusion tensor
`fa`()	Fractional anisotropy (FA) calculated from cached eigenvalues.
`ga`()	Geodesic anisotropy (GA) calculated from cached eigenvalues.
`kmax`([sphere, gtol, mask])	Computes the maximum value of a single voxel kurtosis tensor
`linearity`()	Returns
`md`()	Mean diffusivity (MD) calculated from cached eigenvalues.
`mk`([min_kurtosis, max_kurtosis])	Computes mean Kurtosis (MK) from the kurtosis tensor.
`mode`()	Tensor mode calculated from cached eigenvalues.
`odf`(sphere)	The diffusion orientation distribution function (dODF).
`planarity`()	Returns
`predict`(gtab[, S0])	Given a DKI model fit, predict the signal on the vertices of a gradient table
`rd`()	Radial diffusivity (RD) calculated from cached eigenvalues.
`rk`([min_kurtosis, max_kurtosis])	Radial Kurtosis (RK) of a diffusion kurtosis tensor [1].
`sphericity`()	Returns
`trace`()	Trace of the tensor calculated from cached eigenvalues.

`fit`(data[, mask])	Fit method of the DKI model class
`predict`(dki_params[, S0])	Predict a signal for this DKI model class instance given parameters.

`fit`(data[, mask, sphere, gtol, awf_only])	Fit method of the Diffusion Kurtosis Microstructural Model
`predict`(params[, S0])	Predict a signal for the DKI microstructural model class instance given parameters.

`msd_discrete`([normalized])	Calculates the mean squared displacement on the discrete propagator
`odf`(sphere)	Calculates the real discrete odf for a given discrete sphere
`pdf`()	Applies the 3D FFT in the q-space grid to generate the DSI diffusion propagator, remove the background noise with a hard threshold and then deconvolve the propagator with the Lucy-Richardson deconvolution algorithm
`rtop_pdf`([normalized])	Calculates the return to origin probability from the propagator, which is the propagator evaluated at zero (see Descoteaux et Al.
`rtop_signal`([filtering])	Calculates the return to origin probability (rtop) from the signal rtop equals to the sum of all signal values

`fractional_anisotropy`()	Calculates the fractional anisotropy.
`mean_diffusivity`()	Calculates the mean diffusivity.
`odf`(sphere[, clip_negative])	Calculates the fODF for a given discrete sphere.
`predict`([gtab, S0])	Calculates the fODF for a given discrete sphere.

`fit`(data[, mask])	Fit method for every voxel in data
`predict`(fwdti_params[, S0])	Predict a signal for this TensorModel class instance given parameters.

`estimate_f_D_star`(params_f_D_star, data, S0, D)	Estimate f and D_star using the values of all the other parameters obtained from a linear fit.
`estimate_linear_fit`(data, split_b[, less_than])	Estimate a linear fit by taking log of data.
`fit`(data[, mask])	Fit method for every voxel in data
`predict`(ivim_params, gtab[, S0])	Predict a signal for this IvimModel class instance given parameters.