# ----------------------------------------------------------------------------
# Copyright (c) 2013--, scikit-bio development team.
#
# Distributed under the terms of the Modified BSD License.
#
# The full license is in the file COPYING.txt, distributed with this software.
# ----------------------------------------------------------------------------

from __future__ import absolute_import, division, print_function
import unittest

import numpy as np
import numpy.testing as npt

from skbio.stats.distance import mantel
from skbio.stats.evolve import hommola_cospeciation
from skbio.stats.evolve._hommola import _get_dist, _gen_lists


class HommolaCospeciationTests(unittest.TestCase):
    def setUp(self):
        # Test matrices, as presented in original paper by Hommola et al.
        self.hdist = np.array([[0, 3, 8, 8, 9], [3, 0, 7, 7, 8], [
            8, 7, 0, 6, 7], [8, 7, 6, 0, 3], [9, 8, 7, 3, 0]])
        self.pdist = np.array([[0, 5, 8, 8, 8], [5, 0, 7, 7, 7], [
            8, 7, 0, 4, 4], [8, 7, 4, 0, 2], [8, 7, 4, 2, 0]])
        self.interact = np.array([[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [
            0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 1, 1]])

        # Reduced-size host matrix for testing asymmetric interaction matrix
        self.hdist_4x4 = np.array([[0, 3, 8, 8], [3, 0, 7, 7], [8, 7, 0, 6],
                                  [8, 7, 6, 0]])
        self.interact_5x4 = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0],
                                      [0, 0, 0, 1], [0, 0, 0, 1]])

        # One to one interaction matrix for comparing against Mantel output
        self.interact_1to1 = np.array([[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [
            0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1]])

        # interaction matrix yielding non-significant results.
        # this matrix was picked because it will generate an r value that's
        # less than a standard deviation away from the mean of the normal
        # distribution of r vals
        self.interact_ns = np.array(
            [[0, 0, 0, 1, 0], [0, 0, 0, 0, 1], [1, 0, 0, 0, 0],
             [1, 0, 0, 0, 0], [0, 0, 0, 0, 1]])

        # minimal size matrices for sanity checks of inputs
        self.h_dist_3x3 = np.array([[0, 1, 2], [1, 0, 1], [2, 1, 0]])
        self.h_dist_2x2 = np.array([[0, 3], [3, 0]])
        self.p_dist_3x3 = np.array([[0, 3, 2], [3, 0, 1], [2, 1, 0]])
        self.interact_3x3 = np.array([[0, 1, 1], [1, 0, 1], [0, 0, 1]])
        self.interact_3x2 = np.array([[0, 1], [1, 0], [1, 1]])
        self.interact_2x3 = np.array([[0, 1, 1], [1, 0, 1]])
        self.interact_zero = np.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]])

    def test_hommola_cospeciation_sig(self):
        np.random.seed(1)

        obs_r, obs_p, obs_perm_stats = hommola_cospeciation(
            self.hdist, self.pdist, self.interact, 9)
        exp_p = .1
        exp_r = 0.83170965463247915
        exp_perm_stats = np.array([-0.14928122, 0.26299538, -0.21125858,
                                   0.24143838, 0.61557855, -0.24258293,
                                   0.09885203, 0.02858, 0.42742399])
        self.assertAlmostEqual(obs_p, exp_p)
        self.assertAlmostEqual(obs_r, exp_r)

        npt.assert_allclose(obs_perm_stats, exp_perm_stats)

    def test_hommola_cospeciation_asymmetric(self):
        np.random.seed(1)

        obs_r, obs_p, obs_perm_stats = hommola_cospeciation(
            self.hdist_4x4, self.pdist, self.interact_5x4, 9)
        exp_p = 0.2
        exp_r = 0.85732140997411233
        exp_perm_stats = np.array([-0.315244162496, -0.039405520312,
                                   0.093429386594, -0.387835875941,
                                   0.183711730709,  0.056057631956,
                                   0.945732487487,  0.056057631956,
                                   -0.020412414523])
        self.assertAlmostEqual(obs_p, exp_p)
        self.assertAlmostEqual(obs_r, exp_r)

        npt.assert_allclose(obs_perm_stats, exp_perm_stats)

    def test_hommola_cospeciation_no_sig(self):
        np.random.seed(1)

        obs_r, obs_p, obs_perm_stats = hommola_cospeciation(
            self.hdist, self.pdist, self.interact_ns, 9)
        exp_p = .6
        exp_r = -0.013679391379114569
        exp_perm_stats = np.array([-0.22216543, -0.14836061, -0.04434843,
                                   0.1478281, -0.29105645, 0.56395839,
                                   0.47304992, 0.79125657, 0.06804138])
        self.assertAlmostEqual(obs_p, exp_p)
        self.assertAlmostEqual(obs_r, exp_r)
        npt.assert_allclose(obs_perm_stats, exp_perm_stats)

    def test_hommola_vs_mantel(self):
        # we don't compare p-values because the two methods use different
        # permutation strategies
        r_mantel, p_mantel, _ = mantel(
            self.hdist, self.pdist, method='pearson', permutations=0,
            alternative='greater'
        )
        r_hommola, p_hommola, _ = hommola_cospeciation(
            self.hdist, self.pdist, self.interact_1to1, permutations=0
        )

        self.assertAlmostEqual(r_hommola, r_mantel)
        npt.assert_equal(p_hommola, p_mantel)

    def test_zero_permutations(self):
        obs_r, obs_p, obs_perm_stats = hommola_cospeciation(
            self.hdist, self.pdist, self.interact, 0)

        exp_p = np.nan
        exp_r = 0.83170965463247915
        exp_perm_stats = np.array([])

        npt.assert_equal(obs_p, exp_p)
        self.assertAlmostEqual(obs_r, exp_r)
        npt.assert_equal(obs_perm_stats, exp_perm_stats)

    def test_get_dist(self):
        labels = np.array([0, 1, 1, 2, 3])
        k_labels, t_labels = _gen_lists(labels)
        dists = np.array([[0, 2, 6, 3], [2, 0, 5, 4], [6, 5, 0, 7],
                          [3, 4, 7, 0]])
        index = np.array([2, 3, 1, 0])

        expected_vec = np.array([7, 7, 5, 6, 0, 4, 3, 4, 3, 2])
        actual_vec = _get_dist(k_labels, t_labels, dists, index)

        npt.assert_allclose(actual_vec, expected_vec)

    def test_gen_lists(self):
        exp_pars_k_labels = np.array([0, 0, 0, 0, 0, 1, 1, 1,
                                      1, 2, 2, 2, 3, 3, 4])
        exp_pars_t_labels = np.array([1, 2, 3, 4, 4, 2, 3, 4,
                                      4, 3, 4, 4, 4, 4, 4])
        exp_host_k_labels = np.array([0, 0, 0, 0, 0, 1, 1, 1,
                                      1, 2, 2, 2, 3, 3, 3])
        exp_host_t_labels = np.array([1, 2, 3, 3, 4, 2, 3, 3,
                                      4, 3, 3, 4, 3, 4, 4])

        pars, hosts = np.nonzero(self.interact)

        obs_pars_k_labels, obs_pars_t_labels = _gen_lists(pars)
        obs_hosts_k_labels, obs_hosts_t_labels = _gen_lists(hosts)

        npt.assert_allclose(exp_pars_k_labels, obs_pars_k_labels)
        npt.assert_allclose(exp_pars_t_labels, obs_pars_t_labels)
        npt.assert_allclose(exp_host_k_labels, obs_hosts_k_labels)
        npt.assert_allclose(exp_host_t_labels, obs_hosts_t_labels)

    def test_dm_too_small(self):
        with self.assertRaises(ValueError):
            hommola_cospeciation(self.h_dist_2x2, self.p_dist_3x3,
                                 self.interact_3x3)

    def test_host_interaction_not_equal(self):
        with self.assertRaises(ValueError):
            hommola_cospeciation(self.h_dist_3x3, self.p_dist_3x3,
                                 self.interact_2x3)

    def test_par_interaction_not_equal(self):
        with self.assertRaises(ValueError):
            hommola_cospeciation(self.h_dist_3x3, self.p_dist_3x3,
                                 self.interact_3x2)

    def test_interaction_too_few(self):
        with self.assertRaises(ValueError):
            hommola_cospeciation(self.h_dist_3x3, self.p_dist_3x3,
                                 self.interact_zero)

    def test_permutations_too_few(self):
        with self.assertRaises(ValueError):
            hommola_cospeciation(self.h_dist_3x3, self.p_dist_3x3,
                                 self.interact_3x3, -1)


if __name__ == '__main__':
    unittest.main()
