import numpy as np
import numpy.linalg as L
import numpy.random as R

from nipy.testing import assert_equal, assert_almost_equal, dec

from nipy.algorithms.statistics import intvol, utils


def symnormal(p=10):
    M = R.standard_normal((p,p))
    return (M + M.T) / np.sqrt(2)


def randorth(p=10):
    """
    A random orthogonal matrix.
    """
    A = symnormal(p)
    return L.eig(A)[1]
                

def box(shape, edges):
    data = np.zeros(shape)
    sl = []
    for i in range(len(shape)):
        sl.append(slice(edges[i][0], edges[i][1],1))
    data[sl] = 1
    return data.astype(np.int)


def randombox(shape):
    """
    Generate a random box, returning the box and the edge lengths
    """
    edges = [R.random_integers(0, shape[j], size=(2,)) 
             for j in range(len(shape))]

    for j in range(len(shape)):
        edges[j].sort()
        if edges[j][0] == edges[j][1]:
            edges[j][0] = 0; edges[j][1] = shape[j]/2+1
    return edges, box(shape, edges)
        

def elsym(edgelen, order=1):
    """
    Elementary symmetric polynomial of a given order
    """
    l = len(edgelen)
    if order == 0:
        return 1
    r = 0
    for v in utils.combinations(range(l), order):
        r += np.product([edgelen[vv] for vv in v])
    return r


def nonintersecting_boxes(shape):
    """
    The Lips's are supposed to be additive, so disjoint things
    should be additive. But, if they ALMOST intersect, different
    things get added to the triangulation.

    >>> b1 = np.zeros(40, np.int)
    >>> b1[:11] = 1
    >>> b2 = np.zeros(40, np.int)
    >>> b2[11:] = 1
    >>> (b1*b2).sum()
    0
    >>> c = np.indices((40,)).astype(np.float)
    >>> intvol.Lips1d(c, b1)
    10.0
    >>> intvol.Lips1d(c, b2)
    28.0
    >>> intvol.Lips1d(c, b1+b2)
    39.0

    The function creates two boxes such that
    the 'dilated' box1 does not intersect with box2.
    Additivity works in this case.
    """
    while True:
        edge1, box1 = randombox(shape)
        edge2, box2 = randombox(shape)

        diledge1 = [[max(ed[0]-1, 0), min(ed[1]+1, sh)] 
                    for ed, sh in zip(edge1, box1.shape)]

        dilbox1 = box(box1.shape, diledge1)

        if set(np.unique(dilbox1 + box2)).issubset([0,1]):
            break
    return box1, box2, edge1, edge2


def test_mu3tet():
    assert_equal(intvol.mu3_tet(0,0,0,0,1,0,0,1,0,1), 1./6)


def test_mu2tri():
    assert_equal(intvol.mu2_tri(0,0,0,1,0,1), 1./2)


def test_mu1tri():
    assert_equal(intvol.mu1_tri(0,0,0,1,0,1), 1+np.sqrt(2)/2)


def test_mu2tet():
    assert_equal(intvol.mu2_tet(0,0,0,0,1,0,0,1,0,1), (3./2 + np.sqrt(3./4))/2)