#!/usr/bin/env python

from __future__ import division
from other.core import *
import numpy as np
from numpy import linalg as LA
import sys
#import other.cgal as cg

octa_base_vec = {
'cube':[1,1,1],  
'cuboctahedron':[1,0,1],
'octahedron':[1,0,0],
'small-rhombicuboctahedron': [1,1,sqrt(2)+1],
'truncated-cube':[1,1,sqrt(2)-1],
'truncated-octahedron':[1,0,2],
  }

phi = (1+sqrt(5))/2
icos_base_vec = {
'dodecahedron':[1+phi,1,0],
'icosahedron':[0,1,phi],
'icosidodecahedron':[0,0,1],
'small-rhombicosidodecahedron':[-1,phi,2],
'truncated-dodecahedron':[2,phi,1+phi],
'truncated-icosahedron':[1+2*phi,phi,2]
}



def close_find(l,x,tol=1e-5):
    for i,y in enumerate(l):
        if relative_error(x,y)<tol:
            return i
    return -1

def close_finds(l,x,tol=1e-5):
    inds = []
    for i,y in enumerate(l):
        if relative_error(x,y)<tol:
            inds.append(i)
    return inds

def close_mins(l,tol=1e-5):
    mins = [sys.float_info.max]
    inds = [-1]
    for i,y in enumerate(l):
        if relative_error(y,mins[0])<tol:
            mins.append(y)
            inds.append(i)
        elif y < mins[0]:
            mins = [y]
            inds = [i]
    return mins,inds


#given generators, close the group under multiplication
def close_group(G,tol=1e-5):
    G = list(G)
    grew = 1
    while grew: 
        grew = 0
        for a in G:
            for b in G:
                ab = a*b
                if close_find(G,ab)>=0:
                    break
                else:
                    G.append(ab)
                    grew = 1
    return G

def real_eigvec(m):
    val,vec = LA.eig(m)
    i = close_find(val,1.)
    return vec[i]


def icosahedral_group():
    P = Matrix([[0,1,0],[0,0,1],[1,0,0]],dtype=float)
    p = (+1+sqrt(5))/2
    q = (-1+sqrt(5))/2
    A = Matrix([[q,-p,1],[p,1,q],[-1,q,p]])/2
    assert allclose(dot(A.T,A),eye(3))
    G = close_group([P,A])
    assert len(G)==60
    return G

def octahedral_group():
    R1 = Matrix([[0,0,1],[1,0,0],[0,1,0]],dtype=float)
    R2 = Matrix([[0,1,0],[-1,0,0],[0,0,1]])
    G = close_group([R1,R2])
    assert len(G)==24
    return G

def orbit(G,v):
    v = np.array(v)
    X = v.reshape(1,3)
    close_points = np.array([])
    for g in G:
        y = g*v
        if close_find(X,y)<0:
            X=vstack((X,y))            
    return X

def find_char_length(G,v):
    o = orbit(G,v)
    l = [magnitude(x-v) for x in o]
    mins,ind = close_mins(l[1:],1e-5)
    return mins[0]

def eig(G):
    return np.linalg.eig(G)

def gv(name):
    if name in icos_base_vec.keys():
        g = icosahedral_group()
        v = icos_base_vec[name]
    elif name in octa_base_vec.keys():
        g = octahedral_group()
        v = octa_base_vec[name]
    else:
        print "polyhedron name not found"
        return -1
    return g,v

def find_edges(pts,l,tol=1e-5):
    edges = []
    for i,x in enumerate(pts):
        mags = [magnitude(x-y) for y in pts]
        inds = close_finds(mags,l)
        for j in inds:
            if j>i: edges.append([i,j])
    return edges


def m_rep(p):
    return "{%(x)04f,%(y)04f,%(z)04f}" % {"x":p[0],"y":p[1],"z":p[2]}

def e_rep(p1,p2):
    return "%(x1)05f,%(y1)05f,%(z1)05f,%(x2)05f,%(y2)05f,%(z2)05f" % {"x1":p1[0],"y1":p1[1],"z1":p1[2],"x2":p2[0],"y2":p2[1],"z2":p2[2]}
def v_rep(p1):
    return "%(x1)05f,%(y1)05f,%(z1)05f" % {"x1":p1[0],"y1":p1[1],"z1":p1[2]}

if __name__ == "__main__":
    name = 'small-rhombicosidodecahedron'
    g,v = gv(name)
    v = asarray(v)/magnitude(v)
    l = find_char_length(g,v)
    pts = orbit(g,v)
    edges = find_edges(pts,l)
 #   print edges
 #   print pts
    print "n_verts=",len(pts)," n_edges=",len(edges)," char_length=",l

    f = open(name+'.verts','w')
    for pt in pts:
        f.write(v_rep(pt)+'\n')

    f = open(name+'.edges','w')
    for e in edges:
        f.write(e_rep(pts[e[0]],pts[e[1]])+'\n')

 #   for e in edges:
 #       print "Line[{"+m_rep(pts[e[0]])+','+m_rep(pts[e[1]])+'}],'



#    for poly,vector in icos_base_vec.items():
#        print poly
#    	print find_tri(icosahedral_group(), vector)
#    for poly,vector in octa_base_vec.items():
#        print poly
#        print find_tri(octahedral_group(), vector)