from __future__ import division, print_function, absolute_import

import os
import copy
import pytest

import numpy as np
from numpy.testing import (assert_equal, assert_almost_equal,
                           assert_, assert_allclose, assert_array_equal)
import pytest
from pytest import raises as assert_raises
from scipy._lib.six import xrange

import scipy.spatial.qhull as qhull
from scipy.spatial import cKDTree as KDTree
from scipy.spatial import Voronoi

import itertools

def sorted_tuple(x):
    return tuple(sorted(x))


def sorted_unique_tuple(x):
    return tuple(np.unique(x))


def assert_unordered_tuple_list_equal(a, b, tpl=tuple):
    if isinstance(a, np.ndarray):
        a = a.tolist()
    if isinstance(b, np.ndarray):
        b = b.tolist()
    a = list(map(tpl, a))
    a.sort()
    b = list(map(tpl, b))
    b.sort()
    assert_equal(a, b)


np.random.seed(1234)

points = [(0,0), (0,1), (1,0), (1,1), (0.5, 0.5), (0.5, 1.5)]

pathological_data_1 = np.array([
    [-3.14,-3.14], [-3.14,-2.36], [-3.14,-1.57], [-3.14,-0.79],
    [-3.14,0.0], [-3.14,0.79], [-3.14,1.57], [-3.14,2.36],
    [-3.14,3.14], [-2.36,-3.14], [-2.36,-2.36], [-2.36,-1.57],
    [-2.36,-0.79], [-2.36,0.0], [-2.36,0.79], [-2.36,1.57],
    [-2.36,2.36], [-2.36,3.14], [-1.57,-0.79], [-1.57,0.79],
    [-1.57,-1.57], [-1.57,0.0], [-1.57,1.57], [-1.57,-3.14],
    [-1.57,-2.36], [-1.57,2.36], [-1.57,3.14], [-0.79,-1.57],
    [-0.79,1.57], [-0.79,-3.14], [-0.79,-2.36], [-0.79,-0.79],
    [-0.79,0.0], [-0.79,0.79], [-0.79,2.36], [-0.79,3.14],
    [0.0,-3.14], [0.0,-2.36], [0.0,-1.57], [0.0,-0.79], [0.0,0.0],
    [0.0,0.79], [0.0,1.57], [0.0,2.36], [0.0,3.14], [0.79,-3.14],
    [0.79,-2.36], [0.79,-0.79], [0.79,0.0], [0.79,0.79],
    [0.79,2.36], [0.79,3.14], [0.79,-1.57], [0.79,1.57],
    [1.57,-3.14], [1.57,-2.36], [1.57,2.36], [1.57,3.14],
    [1.57,-1.57], [1.57,0.0], [1.57,1.57], [1.57,-0.79],
    [1.57,0.79], [2.36,-3.14], [2.36,-2.36], [2.36,-1.57],
    [2.36,-0.79], [2.36,0.0], [2.36,0.79], [2.36,1.57],
    [2.36,2.36], [2.36,3.14], [3.14,-3.14], [3.14,-2.36],
    [3.14,-1.57], [3.14,-0.79], [3.14,0.0], [3.14,0.79],
    [3.14,1.57], [3.14,2.36], [3.14,3.14],
])

pathological_data_2 = np.array([
    [-1, -1], [-1, 0], [-1, 1],
    [0, -1], [0, 0], [0, 1],
    [1, -1 - np.finfo(np.float_).eps], [1, 0], [1, 1],
])

bug_2850_chunks = [np.random.rand(10, 2),
                   np.array([[0,0], [0,1], [1,0], [1,1]])  # add corners
                   ]

# same with some additional chunks
bug_2850_chunks_2 = (bug_2850_chunks +
                     [np.random.rand(10, 2),
                      0.25 + np.array([[0,0], [0,1], [1,0], [1,1]])])

DATASETS = {
    'some-points': np.asarray(points),
    'random-2d': np.random.rand(30, 2),
    'random-3d': np.random.rand(30, 3),
    'random-4d': np.random.rand(30, 4),
    'random-5d': np.random.rand(30, 5),
    'random-6d': np.random.rand(10, 6),
    'random-7d': np.random.rand(10, 7),
    'random-8d': np.random.rand(10, 8),
    'pathological-1': pathological_data_1,
    'pathological-2': pathological_data_2
}

INCREMENTAL_DATASETS = {
    'bug-2850': (bug_2850_chunks, None),
    'bug-2850-2': (bug_2850_chunks_2, None),
}


def _add_inc_data(name, chunksize):
    """
    Generate incremental datasets from basic data sets
    """
    points = DATASETS[name]
    ndim = points.shape[1]

    opts = None
    nmin = ndim + 2

    if name == 'some-points':
        # since Qz is not allowed, use QJ
        opts = 'QJ Pp'
    elif name == 'pathological-1':
        # include enough points so that we get different x-coordinates
        nmin = 12

    chunks = [points[:nmin]]
    for j in xrange(nmin, len(points), chunksize):
        chunks.append(points[j:j+chunksize])

    new_name = "%s-chunk-%d" % (name, chunksize)
    assert new_name not in INCREMENTAL_DATASETS
    INCREMENTAL_DATASETS[new_name] = (chunks, opts)


for name in DATASETS:
    for chunksize in 1, 4, 16:
        _add_inc_data(name, chunksize)


class Test_Qhull(object):
    def test_swapping(self):
        # Check that Qhull state swapping works

        x = qhull._Qhull(b'v',
                         np.array([[0,0],[0,1],[1,0],[1,1.],[0.5,0.5]]),
                         b'Qz')
        xd = copy.deepcopy(x.get_voronoi_diagram())

        y = qhull._Qhull(b'v',
                         np.array([[0,0],[0,1],[1,0],[1,2.]]),
                         b'Qz')
        yd = copy.deepcopy(y.get_voronoi_diagram())

        xd2 = copy.deepcopy(x.get_voronoi_diagram())
        x.close()
        yd2 = copy.deepcopy(y.get_voronoi_diagram())
        y.close()

        assert_raises(RuntimeError, x.get_voronoi_diagram)
        assert_raises(RuntimeError, y.get_voronoi_diagram)

        assert_allclose(xd[0], xd2[0])
        assert_unordered_tuple_list_equal(xd[1], xd2[1], tpl=sorted_tuple)
        assert_unordered_tuple_list_equal(xd[2], xd2[2], tpl=sorted_tuple)
        assert_unordered_tuple_list_equal(xd[3], xd2[3], tpl=sorted_tuple)
        assert_array_equal(xd[4], xd2[4])

        assert_allclose(yd[0], yd2[0])
        assert_unordered_tuple_list_equal(yd[1], yd2[1], tpl=sorted_tuple)
        assert_unordered_tuple_list_equal(yd[2], yd2[2], tpl=sorted_tuple)
        assert_unordered_tuple_list_equal(yd[3], yd2[3], tpl=sorted_tuple)
        assert_array_equal(yd[4], yd2[4])

        x.close()
        assert_raises(RuntimeError, x.get_voronoi_diagram)
        y.close()
        assert_raises(RuntimeError, y.get_voronoi_diagram)

    def test_issue_8051(self):
        points = np.array([[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2],[2, 0], [2, 1], [2, 2]])
        Voronoi(points)


class TestUtilities(object):
    """
    Check that utility functions work.

    """

    def test_find_simplex(self):
        # Simple check that simplex finding works
        points = np.array([(0,0), (0,1), (1,1), (1,0)], dtype=np.double)
        tri = qhull.Delaunay(points)

        # +---+
        # |\ 0|
        # | \ |
        # |1 \|
        # +---+

        assert_equal(tri.vertices, [[1, 3, 2], [3, 1, 0]])

        for p in [(0.25, 0.25, 1),
                  (0.75, 0.75, 0),
                  (0.3, 0.2, 1)]:
            i = tri.find_simplex(p[:2])
            assert_equal(i, p[2], err_msg='%r' % (p,))
            j = qhull.tsearch(tri, p[:2])
            assert_equal(i, j)

    def test_plane_distance(self):
        # Compare plane distance from hyperplane equations obtained from Qhull
        # to manually computed plane equations
        x = np.array([(0,0), (1, 1), (1, 0), (0.99189033, 0.37674127),
                      (0.99440079, 0.45182168)], dtype=np.double)
        p = np.array([0.99966555, 0.15685619], dtype=np.double)

        tri = qhull.Delaunay(x)

        z = tri.lift_points(x)
        pz = tri.lift_points(p)

        dist = tri.plane_distance(p)

        for j, v in enumerate(tri.vertices):
            x1 = z[v[0]]
            x2 = z[v[1]]
            x3 = z[v[2]]

            n = np.cross(x1 - x3, x2 - x3)
            n /= np.sqrt(np.dot(n, n))
            n *= -np.sign(n[2])

            d = np.dot(n, pz - x3)

            assert_almost_equal(dist[j], d)

    def test_convex_hull(self):
        # Simple check that the convex hull seems to works
        points = np.array([(0,0), (0,1), (1,1), (1,0)], dtype=np.double)
        tri = qhull.Delaunay(points)

        # +---+
        # |\ 0|
        # | \ |
        # |1 \|
        # +---+

        assert_equal(tri.convex_hull, [[3, 2], [1, 2], [1, 0], [3, 0]])

    def test_volume_area(self):
        #Basic check that we get back the correct volume and area for a cube
        points = np.array([(0, 0, 0), (0, 1, 0), (1, 0, 0), (1, 1, 0),
                           (0, 0, 1), (0, 1, 1), (1, 0, 1), (1, 1, 1)])
        hull = qhull.ConvexHull(points)

        assert_allclose(hull.volume, 1., rtol=1e-14,
                        err_msg="Volume of cube is incorrect")
        assert_allclose(hull.area, 6., rtol=1e-14,
                        err_msg="Area of cube is incorrect")

    def test_random_volume_area(self):
        #Test that the results for a random 10-point convex are
        #coherent with the output of qconvex Qt s FA
        points = np.array([(0.362568364506, 0.472712355305, 0.347003084477),
                           (0.733731893414, 0.634480295684, 0.950513180209),
                           (0.511239955611, 0.876839441267, 0.418047827863),
                           (0.0765906233393, 0.527373281342, 0.6509863541),
                           (0.146694972056, 0.596725793348, 0.894860986685),
                           (0.513808585741, 0.069576205858, 0.530890338876),
                           (0.512343805118, 0.663537132612, 0.037689295973),
                           (0.47282965018, 0.462176697655, 0.14061843691),
                           (0.240584597123, 0.778660020591, 0.722913476339),
                           (0.951271745935, 0.967000673944, 0.890661319684)])

        hull = qhull.ConvexHull(points)
        assert_allclose(hull.volume, 0.14562013, rtol=1e-07,
                        err_msg="Volume of random polyhedron is incorrect")
        assert_allclose(hull.area, 1.6670425, rtol=1e-07,
                        err_msg="Area of random polyhedron is incorrect")

    def test_incremental_volume_area_random_input(self):
        """Test that incremental mode gives the same volume/area as
        non-incremental mode and incremental mode with restart"""
        nr_points = 20
        dim = 3
        points = np.random.random((nr_points, dim))
        inc_hull = qhull.ConvexHull(points[:dim+1, :], incremental=True)
        inc_restart_hull = qhull.ConvexHull(points[:dim+1, :], incremental=True)
        for i in range(dim+1, nr_points):
            hull = qhull.ConvexHull(points[:i+1, :])
            inc_hull.add_points(points[i:i+1, :])
            inc_restart_hull.add_points(points[i:i+1, :], restart=True)
            assert_allclose(hull.volume, inc_hull.volume, rtol=1e-7)
            assert_allclose(hull.volume, inc_restart_hull.volume, rtol=1e-7)
            assert_allclose(hull.area, inc_hull.area, rtol=1e-7)
            assert_allclose(hull.area, inc_restart_hull.area, rtol=1e-7)

    def _check_barycentric_transforms(self, tri, err_msg="",
                                      unit_cube=False,
                                      unit_cube_tol=0):
        """Check that a triangulation has reasonable barycentric transforms"""
        vertices = tri.points[tri.vertices]
        sc = 1/(tri.ndim + 1.0)
        centroids = vertices.sum(axis=1) * sc

        # Either: (i) the simplex has a `nan` barycentric transform,
        # or, (ii) the centroid is in the simplex

        def barycentric_transform(tr, x):
            ndim = tr.shape[1]
            r = tr[:,-1,:]
            Tinv = tr[:,:-1,:]
            return np.einsum('ijk,ik->ij', Tinv, x - r)

        eps = np.finfo(float).eps

        c = barycentric_transform(tri.transform, centroids)
        olderr = np.seterr(invalid="ignore")
        try:
            ok = np.isnan(c).all(axis=1) | (abs(c - sc)/sc < 0.1).all(axis=1)
        finally:
            np.seterr(**olderr)

        assert_(ok.all(), "%s %s" % (err_msg, np.nonzero(~ok)))

        # Invalid simplices must be (nearly) zero volume
        q = vertices[:,:-1,:] - vertices[:,-1,None,:]
        volume = np.array([np.linalg.det(q[k,:,:])
                           for k in range(tri.nsimplex)])
        ok = np.isfinite(tri.transform[:,0,0]) | (volume < np.sqrt(eps))
        assert_(ok.all(), "%s %s" % (err_msg, np.nonzero(~ok)))

        # Also, find_simplex for the centroid should end up in some
        # simplex for the non-degenerate cases
        j = tri.find_simplex(centroids)
        ok = (j != -1) | np.isnan(tri.transform[:,0,0])
        assert_(ok.all(), "%s %s" % (err_msg, np.nonzero(~ok)))

        if unit_cube:
            # If in unit cube, no interior point should be marked out of hull
            at_boundary = (centroids <= unit_cube_tol).any(axis=1)
            at_boundary |= (centroids >= 1 - unit_cube_tol).any(axis=1)

            ok = (j != -1) | at_boundary
            assert_(ok.all(), "%s %s" % (err_msg, np.nonzero(~ok)))

    def test_degenerate_barycentric_transforms(self):
        # The triangulation should not produce invalid barycentric
        # transforms that stump the simplex finding
        data = np.load(os.path.join(os.path.dirname(__file__), 'data',
                                    'degenerate_pointset.npz'))
        points = data['c']