#copyright Igor Rivin, 2005, all rights reserved.

from numarray import *
import numarray.linear_algebra as la
import math
from igorutils import *

#returns the gram matrix of a list of vectors.

def gram(veclist):
    return outerprod(innerproduct, veclist, veclist)


# the problem is this: given a (not necessarily orthonormal) basis
# v1, ..., vn, and a collection of numbers t1,...tn,
# find a vector x, such that <x, vi> = ti. Supposing that
# x = sum ai vi, we see that tj = sum ai <vi, vj>, so to find the
# coefficients ai, we must solve G A = T, where G is the gram matrix
# of v1,...,vn. In the orthonormal case, G = I, so we recover the
# well-known observation that that ai=ti, for all i.

def findcoeffs(veclist, innerprodlist):
    Gmat = a64(gram(veclist))
    rhs = a64(innerprodlist)
    try:
        soln = list(la.solve_linear_equations(Gmat, rhs))
        weighted = map(lambda x, y: x * y, soln, veclist)
        theans = reduce(operator.add, weighted)
        return theans
    except:
        return None

#given k points in general position in  Rn returns the circumradius
#and the circumcenter of
#the k-1 dimensional sphere circumscribing the (k-1 dimensional)
#simplex spanning those points. If the simplex happens to be
#degenerate, will bomb out and return the tuple (None, None).

def circumcenter(pointlist):
    veclist = [i - pointlist[0] for i in pointlist[1:]]
    rhs = [innerproduct(i, i)/2.0 for i in veclist]
    thediff = findcoeffs(veclist, rhs)
    if thediff != None:
        thecenter = thediff + pointlist[0]
        circumradius = math.sqrt(innerproduct(thediff, thediff))
        return (thecenter, circumradius)
    else:
        return (None, None)

#given k points in general position in Rn returns the volume of the
#(k-1)-dimensional simplex spanned by the points.

def simpvol(pointlist):
    veclist = [i - pointlist[0] for i in pointlist[1:]]
    gmat = gram(veclist)
    truedim = len(gmat)
    Gmat = a64(gmat)
    thedet = la.determinant(Gmat)
    return math.sqrt(thedet)/factorial(truedim)

    
def surfarea(pointlist):
    thefaces = nless1subs(pointlist)
    areas = [simpvol(i) for i in thefaces]
    total = reduce(operator.add, areas)
    return (areas, total)

def inradius(pointlist):
    thevol = simpvol(pointlist)
    thearea = surfarea(pointlist)[1]
    thedim = len(pointlist) - 1
    return thedim * thevol/thearea