"""Test spherical harmonic models and the tools associated with those models"""
import numpy as np
import numpy.linalg as npl

from nose.tools import assert_equal, assert_raises, assert_true
from numpy.testing import assert_array_equal, assert_array_almost_equal
import numpy.testing as npt
from scipy.special import sph_harm as sph_harm_sp

from dipy.core.sphere import hemi_icosahedron
from dipy.core.gradients import gradient_table
from dipy.sims.voxel import single_tensor
from dipy.direction.peaks import peak_directions
from dipy.reconst.shm import sf_to_sh, sh_to_sf
from dipy.reconst.interpolate import NearestNeighborInterpolator
from dipy.sims.voxel import multi_tensor_odf
from dipy.data import mrtrix_spherical_functions
from dipy.reconst import odf


from dipy.reconst.shm import (real_sph_harm, real_sym_sh_basis,
                              real_sym_sh_mrtrix, sph_harm_ind_list,
                              order_from_ncoef,
                              OpdtModel, normalize_data, hat, lcr_matrix,
                              smooth_pinv, bootstrap_data_array,
                              bootstrap_data_voxel, ResidualBootstrapWrapper,
                              CsaOdfModel, QballModel, SphHarmFit,
                              spherical_harmonics, anisotropic_power,
                              calculate_max_order)


def test_order_from_ncoeff():
    """

    """
    # Just try some out:
    for sh_order in [2, 4, 6, 8, 12, 24]:
        m, n = sph_harm_ind_list(sh_order)
        n_coef = m.shape[0]
        npt.assert_equal(order_from_ncoef(n_coef), sh_order)


def test_sph_harm_ind_list():
    m_list, n_list = sph_harm_ind_list(8)
    assert_equal(m_list.shape, n_list.shape)
    assert_equal(m_list.shape, (45,))
    assert_true(np.all(np.abs(m_list) <= n_list))
    assert_array_equal(n_list % 2, 0)
    assert_raises(ValueError, sph_harm_ind_list, 1)


def test_real_sph_harm():
    # Tests derived from tables in
    # http://en.wikipedia.org/wiki/Table_of_spherical_harmonics
    # where real spherical harmonic $Y^m_n$ is defined to be:
    #    Real($Y^m_n$) * sqrt(2) if m > 0
    #    $Y^m_n$                 if m == 0
    #    Imag($Y^m_n$) * sqrt(2) if m < 0

    rsh = real_sph_harm
    pi = np.pi
    exp = np.exp
    sqrt = np.sqrt
    sin = np.sin
    cos = np.cos
    assert_array_almost_equal(rsh(0, 0, 0, 0),
                              0.5 / sqrt(pi))
    assert_array_almost_equal(rsh(-2, 2, pi / 5, pi / 3),
                              0.25 * sqrt(15. / (2. * pi)) *
                              (sin(pi / 5.)) ** 2. * cos(0 + 2. * pi / 3) *
                              sqrt(2))
    assert_array_almost_equal(rsh(2, 2, pi / 5, pi / 3),
                              -1 * 0.25 * sqrt(15. / (2. * pi)) *
                              (sin(pi / 5.)) ** 2. * sin(0 - 2. * pi / 3) *
                              sqrt(2))
    assert_array_almost_equal(rsh(-2, 2, pi / 2, pi),
                              0.25 * sqrt(15 / (2. * pi)) *
                              cos(2. * pi) * sin(pi / 2.) ** 2. * sqrt(2))
    assert_array_almost_equal(rsh(2, 4, pi / 3., pi / 4.),
                              -1 * (3. / 8.) * sqrt(5. / (2. * pi)) *
                              sin(0 - 2. * pi / 4.) *
                              sin(pi / 3.) ** 2. *
                              (7. * cos(pi / 3.) ** 2. - 1) * sqrt(2))
    assert_array_almost_equal(rsh(-4, 4, pi / 6., pi / 8.),
                              (3. / 16.) * sqrt(35. / (2. * pi)) *
                              cos(0 + 4. * pi / 8.) * sin(pi / 6.) ** 4. *
                              sqrt(2))
    assert_array_almost_equal(rsh(4, 4, pi / 6., pi / 8.),
                              -1 * (3. / 16.) * sqrt(35. / (2. * pi)) *
                              sin(0 - 4. * pi / 8.) * sin(pi / 6.) ** 4. *
                              sqrt(2))

    aa = np.ones((3, 1, 1, 1))
    bb = np.ones((1, 4, 1, 1))
    cc = np.ones((1, 1, 5, 1))
    dd = np.ones((1, 1, 1, 6))
    assert_equal(rsh(aa, bb, cc, dd).shape, (3, 4, 5, 6))


def test_real_sym_sh_mrtrix():
    coef, expected, sphere = mrtrix_spherical_functions()
    basis, m, n = real_sym_sh_mrtrix(8, sphere.theta, sphere.phi)
    func = np.dot(coef, basis.T)
    assert_array_almost_equal(func, expected, 4)


def test_real_sym_sh_basis():
    # This test should do for now
    # The mrtrix basis should be the same as re-ordering and re-scaling the
    # fibernav basis
    new_order = [0, 5, 4, 3, 2, 1, 14, 13, 12, 11, 10, 9, 8, 7, 6]
    sphere = hemi_icosahedron.subdivide(2)
    basis, m, n = real_sym_sh_mrtrix(4, sphere.theta, sphere.phi)
    expected = basis[:, new_order]
    expected *= np.where(m == 0, 1., np.sqrt(2))

    fibernav_basis, m, n = real_sym_sh_basis(4, sphere.theta, sphere.phi)
    assert_array_almost_equal(fibernav_basis, expected)


def test_smooth_pinv():
    hemi = hemi_icosahedron.subdivide(2)
    m, n = sph_harm_ind_list(4)
    B = real_sph_harm(m, n, hemi.theta[:, None], hemi.phi[:, None])

    L = np.zeros(len(m))
    C = smooth_pinv(B, L)
    D = np.dot(npl.inv(np.dot(B.T, B)), B.T)
    assert_array_almost_equal(C, D)

    L = n * (n + 1) * .05
    C = smooth_pinv(B, L)
    L = np.diag(L)
    D = np.dot(npl.inv(np.dot(B.T, B) + L * L), B.T)

    assert_array_almost_equal(C, D)

    L = np.arange(len(n)) * .05
    C = smooth_pinv(B, L)
    L = np.diag(L)
    D = np.dot(npl.inv(np.dot(B.T, B) + L * L), B.T)
    assert_array_almost_equal(C, D)


def test_normalize_data():

    sig = np.arange(1, 66)[::-1]

    where_b0 = np.zeros(65, 'bool')
    where_b0[0] = True
    d = normalize_data(sig, where_b0, 1)
    assert_raises(ValueError, normalize_data, sig, where_b0, out=sig)

    norm_sig = normalize_data(sig, where_b0, min_signal=1)
    assert_array_almost_equal(norm_sig, sig / 65.)
    norm_sig = normalize_data(sig, where_b0, min_signal=5)
    assert_array_almost_equal(norm_sig[-5:], 5 / 65.)

    where_b0[[0, 1]] = [True, True]
    norm_sig = normalize_data(sig, where_b0, min_signal=1)
    assert_array_almost_equal(norm_sig, sig / 64.5)
    norm_sig = normalize_data(sig, where_b0, min_signal=5)
    assert_array_almost_equal(norm_sig[-5:], 5 / 64.5)

    sig = sig * np.ones((2, 3, 1))

    where_b0[[0, 1]] = [True, False]
    norm_sig = normalize_data(sig, where_b0, min_signal=1)
    assert_array_almost_equal(norm_sig, sig / 65.)
    norm_sig = normalize_data(sig, where_b0, min_signal=5)
    assert_array_almost_equal(norm_sig[..., -5:], 5 / 65.)

    where_b0[[0, 1]] = [True, True]
    norm_sig = normalize_data(sig, where_b0, min_signal=1)
    assert_array_almost_equal(norm_sig, sig / 64.5)
    norm_sig = normalize_data(sig, where_b0, min_signal=5)
    assert_array_almost_equal(norm_sig[..., -5:], 5 / 64.5)


def make_fake_signal():
    hemisphere = hemi_icosahedron.subdivide(2)
    bvecs = np.concatenate(([[0, 0, 0]], hemisphere.vertices))
    bvals = np.zeros(len(bvecs)) + 2000
    bvals[0] = 0
    gtab = gradient_table(bvals, bvecs)

    evals = np.array([[2.1, .2, .2], [.2, 2.1, .2]]) * 10 ** -3
    evecs0 = np.eye(3)
    sq3 = np.sqrt(3) / 2.
    evecs1 = np.array([[sq3, .5, 0],
                       [.5, sq3, 0],
                       [0, 0, 1.]])
    evecs1 = evecs0
    a = evecs0[0]
    b = evecs1[1]
    S1 = single_tensor(gtab, .55, evals[0], evecs0)
    S2 = single_tensor(gtab, .45, evals[1], evecs1)
    return S1 + S2, gtab, np.vstack([a, b])


class TestQballModel(object):

    model = QballModel

    def test_single_voxel_fit(self):
        signal, gtab, expected = make_fake_signal()
        sphere = hemi_icosahedron.subdivide(4)

        model = self.model(gtab, sh_order=4, min_signal=1e-5,
                           assume_normed=True)
        fit = model.fit(signal)
        odf = fit.odf(sphere)
        assert_equal(odf.shape, sphere.phi.shape)
        directions, _, _ = peak_directions(odf, sphere)
        # Check the same number of directions
        n = len(expected)
        assert_equal(len(directions), n)
        # Check directions are unit vectors
        cos_similarity = (directions * directions).sum(-1)
        assert_array_almost_equal(cos_similarity, np.ones(n))
        # Check the directions == expected or -expected
        cos_similarity = (directions * expected).sum(-1)
        assert_array_almost_equal(abs(cos_similarity), np.ones(n))

        # Test normalize data
        model = self.model(gtab, sh_order=4, min_signal=1e-5,
                           assume_normed=False)
        fit = model.fit(signal * 5)
        odf_with_norm = fit.odf(sphere)