import math, numpy as np, random
import numpy.polynomial as nppols
import functools
import sys

#pour l'affichage
import pygame


pygame.font.init()
font = pygame.font.SysFont('Lucida Console',20)
cur_anim=None

factorials = [1]

for i in range(1,15):
    factorials.append(i*factorials[-1])

#composition de fonctions (utile pour manipuler les flèches données par les objets universels)
def compose(*functions):
    def compose2(f, g):
        return lambda x: f(g(x))
    return functools.reduce(compose2, functions, lambda x: x)

def simplex(n):
    if n==-1:
        return Cpx()
    res=Cpx()
    res.max_faces.append({i for i in range(n+1)})
    res.init_data()
    return res

def sphere(n=2):
    return simplex(n+1).skeleton(n)

def torus(n=2):
    return circle()**n

def circle():
    return sphere(1)

def dot():
    return simplex(0)

def empty_space():
    return Cpx()

def segment(length=1):
    if length==0:
        return dot()
    if length==1: #non nécessaire mais plus rapide
        return simplex(1)
    prev=segment(length-1)
    return glue([prev,simplex(1)],[[prev.dots-1],[0]])[0]

def mobius_strip():
    strip = segment(3)*segment(1)
    return strip.quotient([[6,5],[7,4]])[0]     #numéros trouvés grâce à la numérotation dans l'animation

def real_projective_plane():
    C=segment()*segment()
    C=C.bary_div().bary_div()  # 2 subdivisions nécessaires pour avoir la torsion
    return C.quotient([[2,17],[1,44],[3,41],[8,42],[9,33],[20,34],[21,29],[24,30]])[0]     #numéros trouvés grâce à la numérotation dans l'animation

def klein_bottle():
    return real_projective_plane().connected_times(2)

 #todo algo plus performant pour trier, ici O(n^2)
def sort_face(face):
    res = []
    k=0

    for x in face:
        ok=False
        for j in range(k):
            if res[j]>=x:
                res.insert(j,x)
                ok=True
                break
        if not ok:
            res.append(x)
        k+=1
    return res

#image d'une face par une application
def map_face(f, face):
    return {f(i) for i in face}

#polynome
def poly(coefs):
    if len(coefs)==0:
        return poly([0])
    return nppols.Polynomial(coefs)

#fonction qui agit sur les points vers liste qui représente cette fonction, et inversement
def func_to_list(f,dots):
    return [f(x) for x in range(dots)]
def list_to_func(L):
    cL=L[:]
    return lambda x:cL[x]

#longueur simpliciale normale (permet un calcul de la géométrie des simplexes réguliers)
def delta_rec(n,dn):
    return (n+1)/(n+2)*math.sqrt(1-dn**2)

#non utilisée finalement
def spherical(phi,theta,r):
    ct=math.cos(theta)
    cp=math.cos(phi)
    st=math.sin(theta)
    sp=math.sin(phi)
    return np.array([   r*st*cp,
                        r*st*sp,
                        r*ct])

#convertit un set en entier binaire
def set_to_bits(set_):
    res=0
    for x in set_:
        res+=2**x
    return res

#réciproque
def bits_to_set(bits):
    res=set()
    b=bits
    i=0
    while b>0:
        q=b//2
        r=b-2*q
        if r==1:
            res.add(i)
        b=q
        i+=1
    return res

delta_ns=[0]

#on ne fera surement pas des complexes de dimension 50
for i in range(50):
    delta_ns.append(delta_rec(i,delta_ns[-1]))

#formule utilisee (qui se prouve par recurrence via Base*Hauteur/dimension) : 1/d! det gram(A1A2,...A1Ad+1)
def face_volume(points):
    d=len(points)-1
    if d<0:
        return 0
    origin = points[0]
    return det_gram([points[i]-origin for i in range(1,d+1)]) / factorials[d]

#coproduit (retourne le coproduit et les inclusions canoniques)
def coproduct(complexes):
    res = Cpx()
    inclusions=[]
    offset=0
    for c in complexes:
        #astuce pour pas que le offset change dans la fonction
        inc=(lambda x: lambda i:i+x)(offset)
        inclusions.append(inc)
        for f in c.max_faces:
            res.max_faces.append(map_face(inc,f))
        offset+=c.dots
    res.init_data()
    return res,inclusions

#algorithme détaillé dans le pdf, calcule le produit géométrique des complexes et retourne aussi les projections canoniques
def product(complexes):
    def bin_f(C,D):
        def recursively_list_max_faces_projecting_onto(w,h):
            if w==0 and h==0:
                return [[(0,0),(0,0)]]
            res=[]

            if w>=1:
                for F in recursively_list_max_faces_projecting_onto(w-1,h):
                    res.append(F+[(w,h)])
            if h>=1:
                for F in recursively_list_max_faces_projecting_onto(w,h-1):
                    res.append(F+[(w,h)])
            return res
        def list_faces_projecting_onto(F,G):
            F_num=list(F) #numérotations de F et G
            G_num=list(G) #
            numbering = lambda t: F_num[t[0]]*D.dots + G_num[t[1]]
            L=recursively_list_max_faces_projecting_onto(len(F)-1,len(G)-1) #calcule deux fois len(F) --> peut s'optimiser
            return [{numbering(v) for v in F} for F in L]

        #projections canoniques
        pi_C=lambda x:x//D.dots
        pi_D=lambda x:x%D.dots
        res=Cpx()
        for F in C.max_faces:
            for G in D.max_faces:
                L=list_faces_projecting_onto(F,G)
                res.max_faces+=L
        res.init_data()
        return res, (pi_C,pi_D)
    lg=len(complexes)
    if lg==0: #produit vide
        return dot(), []
    if lg==1:
        return complexes[0].clone(),[lambda x:x]
    p1,projs1=product(complexes[1:])
    p2,projs2=bin_f(complexes[0],p1)
    return p2, [projs2[0]]+[compose(projs1[i],projs2[1]) for i in range(lg-1)]


#pushout indexé par un ensemble fini
#glue_data[i][k] sera collé à glue_data[j][k], et les glue_data[i] ont tous la même taille
def glue(complexes,glue_data):
    if len(glue_data)==0:
        return Cpx()
    cpxs=len(complexes)
    dots_to_glue=len(glue_data[0])
    cop,inclusions=coproduct(complexes)
    classes=[[inclusions[i](glue_data[i][k]) for i in range(cpxs)] for k in range(dots_to_glue)]
    res,pi=cop.quotient(classes)
    return res,[compose(pi,inclusions[i]) for i in range(cpxs)]

#/!\les complexes[i] doivent être homogènes, sans bords, et de même dimension (exemple : deux sphères)
#/!\ les flèches renvoyées ne sont pas forcément des morphismes simpliciaux
def connected_sum(complexes):
    def bin_f(C,D):
        if len(C.max_faces)==0:
            return D.clone(), [lambda x:x, lambda x:x]
        F1=C.max_faces[0]
        F2=D.max_faces[0]
        G,inclusions=glue([C,D],[list(F1),list(F2)])
        #supprimer la face du collage
        to_remove = map_face(inclusions[0],F1)
        G.max_faces.remove(to_remove)
        G.nfaces[len(F1)-1].remove(to_remove)
        id_to_remove=-1
        for F in G.faces_tree.sub_faces:
            if F.face==to_remove:
                id_to_remove=F.identity
                break
        G.faces_tree.sub_faces=[F for F in G.faces_tree.sub_faces if F.face!=to_remove]
        if id_to_remove != -1:
            nodes=G.faces_tree.all_nodes()
            for F in nodes:
                if F.identity>=id_to_remove:
                    F.identity-=1

        return G,inclusions
    lg=len(complexes)
    if lg==0:
        return Cpx(), []
    if lg==1:
        return complexes[0].clone(),[lambda x:x]
    c1,incs=connected_sum(complexes[1:])
    G,inclusions=bin_f(complexes[0],c1)

    return G, [inclusions[0]]+[compose(inclusions[1],incs[i]) for i in range(lg-1)]

#algorithmes graphiques
def rotation_x(theta):
    cos=math.cos(theta)
    sin=math.sin(theta)
    return np.array([   [1,0,0],
                        [0,cos,-cos],
                        [0,sin,-sin]])
def rotation_y(theta):
    cos=math.cos(theta)
    sin=math.sin(theta)
    return np.array([   [cos,0,sin],
                        [0,1,0],
                        [-sin,0,cos]])
def rotation_z(theta):
    cos=math.cos(theta)
    sin=math.sin(theta)
    return np.array([   [cos,-sin,0],
                        [sin,cos,0],
                        [0,0,1]])

def draw_rect(surface, color, rect):
    try:
        shape_surf = pygame.Surface(pygame.Rect(rect).size, pygame.SRCALPHA)
        pygame.draw.rect(shape_surf, color, shape_surf.get_rect())
        surface.blit(shape_surf, rect)
    except:
        pass

def draw_circle(surface, color, center, radius):
    try:
        target_rect = pygame.Rect(center, (0, 0)).inflate((radius * 2, radius * 2))
        shape_surf = pygame.Surface(target_rect.size, pygame.SRCALPHA)
        pygame.draw.circle(shape_surf, color, (radius, radius), radius)
        surface.blit(shape_surf, target_rect)
    except:
        pass

def draw_polygon(surface, color, points):
    try:
        lx, ly = zip(*points)
        min_x, min_y, max_x, max_y = min(lx), min(ly), max(lx), max(ly)
        target_rect = pygame.Rect(min_x, min_y, max_x - min_x, max_y - min_y)
        shape_surf = pygame.Surface(target_rect.size, pygame.SRCALPHA)
        pygame.draw.polygon(shape_surf, color, [(x - min_x, y - min_y) for x, y in points])
        surface.blit(shape_surf, target_rect)
    except:
        pass

def draw_text(surface,color,pos,text):
    try:
        surface.blit(font.render(text,False,color),pos)
    except:
        pass


#complexe simplicial ordonné, implémentation par les points
#les indices des points doivent former [[0,nombre de points -1]]
class Cpx:
    def __init__(self):
        #les faces sont ici des ensembles de points (type ensemble d'entiers) PAS DE FACE DE DIMENSION -1 DANS L'IMPLEMENTATION
        self.max_faces=[]

        #membres auxilliaires, pour rendre les calculs plus rapides
        self.dim=-1
        self.nfaces=[]
        self.faces_tree=FacesTree()
        self.dots=0
        self.faces_saved=0
        self.faces_dic={}



    #remplit les membres auxilliares
    def init_data(self):
        dico={}  #face --> facestree
        for face in self.max_faces:
            self.add_face_to_list(face,self.faces_tree,dico)
        self.dim=len(self.nfaces)-1
        self.dots=len(self.get_faces(0))

    #liste des faces de dimension d
    def get_faces(self,d):
        if d>self.dim or d<0:
            return []
        return self.nfaces[d]

    def animation(self,ambient=3):
        return EmbeddedManifold(self,ambient).animation()

    #lance une animation dans un espace ambient de la dimension donnée
    def animate(self,ambient=3):
        self.animation(ambient).start()

    #classes : les classes d'équivalence (potentiellement sans les classes triviales), sous la forme d'une liste de listes
    #/!\ modifie la liste classes en y ajoutant les classes triviales
    def quotient(self,classes):
        #surjection canonique vers les classes d'équivalence
        pi1=[None]*self.dots
        for c in classes:
            for x in c:
                pi1[x]=c
        for x in range(self.dots): #classes triviales
            if pi1[x]==None:
                sing=[x]
                classes.append(sing)
                pi1[x]=sing
        #surjection canonique vers le nouveau complexe
        pi=[-1]*self.dots
        curIndic=0
        for c in classes:
            for i in c:
                pi[i]=curIndic
            curIndic+=1

        res = Cpx()

        #comme 2^pi croissante pour l'inclusion (2^* est le foncteur ensemble des parties),
        #il suffit de le faire pour les faces maximales

        pi_function=lambda x:pi[x]

        for F in self.max_faces:
            new_face = map_face(pi_function,F)
            if not new_face in res.max_faces:
                res.max_faces.append(new_face)

        res.init_data()

        return res,pi_function

    #crée les bases pour le calcul de l'homologie
    def init_bases(self):
        self.faces_dic={}
        for i in range(self.dim+1):
            Fi=self.get_faces(i)
            for j in range(len(Fi)):
                b=set_to_bits(Fi[j])
                self.faces_dic[b]=j

    #retourne un tableau d'objets à dim+1 éléments, l'élément i est le vecteur des facteurs invariants de Ker Di / Im Di+1
    def Z_homology(self):
        if self.dim==-1:
            return np.array([[]])
        res=np.array([None]*(self.dim+1),dtype=object)
        self.init_bases()
        border_matrices = [self.border_matrix(i)  for i in range(self.dim+1)]

        for i in range(0,self.dim):
            res[i]=homology_factors(border_matrices[i],border_matrices[i+1])

        res[self.dim] = np.array([0] * base_of_ker(border_matrices[self.dim]).shape[1])  #Ker dd/ Im dd+1 = Ker dd --> libre, on trouve une base et on compte

        return res

    #caractéristique d'euler
    def euler(self):
        return self.poincare()(-1)

    #affiche toutes les informations d'homologie
    def disp_homology(self):
        hom=self.Z_homology()
        for i in range(hom.shape[0]):
            print("H"+str(i)+" : " + disp_group(hom[i]))

        print()

        p = self.poincare()

        print_pol(p)

        print()
        print("Carac :")
        print(p(-1))

    #polynome de poincaré
    def poincare(self):
        if self.dim==-1:
            return poly([0])
        ranks=[0]*(self.dim+1)   #rank[i] est le rang de d_i
        bettis=[0]*(self.dim+1)
        self.init_bases()
        for i in range(1,self.dim+1):
            M=self.border_matrix(i)
            ranks[i]=gauss(M)
        for i in range(self.dim):
            bettis[i]=len(self.get_faces(i))-ranks[i]-ranks[i+1]
        bettis[-1]=len(self.get_faces(self.dim))-ranks[self.dim]
        return poly(bettis)

    #matrice de dn
    def border_matrix(self,n):
        if n==0:
            return np.zeros((0,len(self.get_faces(0))))

        Fn=self.get_faces(n)
        Fn_1 = self.get_faces(n-1)
        sizeCn=len(Fn)
        sizeCn_1=len(Fn_1)
        #SI PROBLEMES : mettre un tableau d'entiers (mais la fonction full gauss ne fonctionne pas avec les tableaux d'entiers)
        res = np.zeros((sizeCn_1,sizeCn))

        for i in range(sizeCn):
            f=Fn[i]                    #orientée par l'ordre des entiers qui représentent les points
            orientation = sort_face(f)
            b=set_to_bits(f)
            for j in range(n+1):
                sub_face_bits=b-2**orientation[j]
                res[self.faces_dic[sub_face_bits],i]=1 if (j%2)==0 else -1
        return res


    #écrase des points
    def simple_quotient(self,cl):
        return self.quotient([cl])

    #/!\ self doit être une variété
    #retourne les points appartenant au bord de la variété
    def border_dots(self):
        res=[]
        border_faces=[]
        for F in self.faces_tree.sub_faces:
            for G in F.sub_faces:
                new_face=G.face
                if new_face in border_faces:
                    border_faces.remove(new_face)  #la face doit appartenir à une seule cellule, et self est une variété donc il y en a max 2
                else:
                    border_faces.append(new_face)
        for F in border_faces:
            for d in F:
                if not d in res:
                    res.append(d)
        return res

    #copie le complexe simplicial
    def clone(self):
        res=Cpx()
        res.max_faces=[f.copy() for f in self.max_faces]
        res.init_data()
        return res

    #cone (contractile) au dessus du complexe
    def cone(self):
        res=Cpx()
        vertex=self.dots #vertex index
        for F in self.max_faces:
            res.max_faces.append(F.union({vertex}))
        res.init_data()
        return res,lambda x:x

    #à appeler avant le init data
    def relabel_dots(self):
        self.nfaces=[]
        self.faces_tree=FacesTree()
        self.dots=0
        self.faces_saved=0
        self.dim=-1

        dots_list=[]
        phi_inv={}
        for F in self.max_faces:
            for d in F:
                if not d in dots_list:
                    phi_inv[d]=len(dots_list)
                    dots_list.append(d)

        phi=list_to_func(dots_list)

        for k in range(len(self.max_faces)):
            self.max_faces[k]={phi_inv[d] for d in self.max_faces[k]}

        return phi

    #/!\ self doit être une variété
    #calcule la variété bord
    def manifold_border(self):
        res=Cpx()
        for F in self.faces_tree.sub_faces:
            for G in F.sub_faces:
                new_face=G.face.copy()
                if new_face in res.max_faces:
                    res.max_faces.remove(new_face)  #la face doit appartenir à une seule cellule, et self est une variété donc il y en a max 2
                else:
                    res.max_faces.append(new_face)
        phi=res.relabel_dots()
        res.init_data()
        return res, phi

    #retourne la liste des faces de dimension d-1 qui sont des bords et leur informe de qui est la cellule qui les contient
    #pas du tout optimal
    def border_info(self):
        res=[]

        for F in self.faces_tree.sub_faces:
            for G in F.sub_faces:
                if G in res:
                    res.remove(G)  #la face doit appartenir à une seule cellule (et elle appartient à au plus 2 cellules pour une variété
                    G.parent = None
                else:
                    res.append(G)
                    G.parent = F

        return res

    #squelette de dimension k (sous complexe dont les faces sont celles de dimension <= k)
    def skeleton(self,k):
        res = Cpx()
        for f in self.get_faces(k):
            res.max_faces.append(f.copy())
        res.init_data()
        return res

    #subdivision barycentrique
    def bary_div(self):
        res=Cpx()

        roads=self.faces_tree.roads()       #liste des chaînes totalement ordonnées maximales (du tronc de l'arbre jusqu'à une feuille)
                                            #cela constitue les faces maximales de la subdivision (voir théorie)
        for r in roads:
            res.max_faces.append({f.identity for f in r if f.identity>=0})

        res.init_data()
        return res

    #efface les infos
    def clear_additional_data(self):
        self.dim=-1
        self.nfaces=[]
        self.faces_tree=FacesTree()
        self.dots=0
        self.faces_saved=0

    #produit géométrique et puissances
    def __mul__(self,other):
        if isinstance(other,int):
            return coproduct([self for i in range(other)])[0]
        return product([self,other])[0]
    def __rmul__(self,other):
        return self.__mul__(other)
    def __pow__(self,n):
        if n==0:
            return simplex(0)
        else:
            return self*(self**(n-1))

    #coproduit
    def __add__(self,other):
        return coproduct([self,other])[0]

    #/!\ self doit être homogène
    #somme connexe avec soi même k fois
    def connected_times(self,k):
        res,incs = connected_sum([self]*k)
        return res

    #ajout d'une face proprement
    def add_face_to_list(self,face,tree_node,dico):
        d=len(face)-1

        nothing=False
        while len(self.nfaces)<=d:
            nothing=True
            self.nfaces.append([])

        if nothing or (not face in self.nfaces[d]):
            self.nfaces[d].append(face)
            new_node=FacesTree()
            new_node.face=face
            dico[set_to_bits(face)]=new_node
            new_node.identity=self.faces_saved
            self.faces_saved+=1
            tree_node.sub_faces.append(new_node)
            if d>=1:
                for v in face:
                    self.add_face_to_list(face-{v},new_node,dico)
        else:
            cur_node=dico[set_to_bits(face)]
            if cur_node!=None:
                tree_node.sub_faces.append(cur_node)


#/!\ pas un vrai arbre : certaines faces apparaissent plusieurs fois
class FacesTree:
    def __init__(self):
        #none pour le tronc
        self.face=None
        self.identity=-1 #un entier qui représente cette face en tant que point pour la subdivision barycentrique (-1 pour le tronc)
        self.sub_faces=[]

        #info en plus pour les faces qui sont des bords
        self.parent=None

    #chemins allant du tronc à une feuille
    def roads(self):
        if len(self.sub_faces)==0:
            return [[self]]
        return [[self]+r for f in self.sub_faces for r in f.roads()]

    #exécute la fonction f une et une seule fois pour chaque élément de l'arbre
    def propagate(self,f):
        prop=[]
        self.rec_prop(f,prop)

    def rec_prop(self,f,propagated):
        if not (self.identity in propagated):
            propagated.append(self.identity)
            f(self)
            for F in self.sub_faces:
                F.rec_prop(f,propagated)

    #liste des noeuds
    def all_nodes(self):
        res=[]
        def action(f):
            res.append(f)
        self.propagate(action)
        return res

#pour l'interface graphique
class Button:
    def __init__(self,x,y,w,h,text,ev):
        self.x=x
        self.y=y
        self.w=w
        self.h=h
        self.rect=(x,y,w,h)
        self.pos=(x,y)
        self.ev=ev
        self.text=text
    def draw(self,wind):
        draw_rect(wind,(30,30,30),self.rect)
        draw_text(wind,(255,255,255),(self.x+5,self.y+5),self.text)

#complexe plongé dans Rn
class EmbeddedManifold:
    def __init__(self,c,ambient=3):
        self.cpx=c
        self.ambient=ambient
        self.positions=[]
        self.speeds=[]
        self.friction=0.1
        self.dt=0.02

        #dictionnaire qui à une face de dimension d-1 associe la liste des cellules voisines (1 ou 2 pour une variété)
        self.links = {}

        #pression intérieur * n-volume (constant)
        self.pv = 1

        #pression extérieure
        self.p_ext = 1

        self.infl=False
        self.spring=0.5 #raideur du ressort

        self.shuffle()
        self.force_center()

        #temporaire
        self.vec_o = np.zeros(self.ambient)
        self.vec_v = np.zeros(self.ambient)

        #liste de couples : face et à quelle unique cellule elle appartient
        self.border_info = self.cpx.border_info()

    #isobarycentre géométrique
    def center(self):
        res=np.zeros(self.ambient)
        if len(self.positions)==0:
            return res
        for p in self.positions:
            res+=p
        return res/len(self.positions)

    #centre d'une face
    def face_center(self,face):
        res=np.zeros(self.ambient)
        if len(face)==0:
            return res
        for p in face:
            res+=self.positions[p]
        return res/len(face)

    #centre la caméra sur le complexe
    def force_center(self):
        c=self.center()
        for p in self.positions:
            p-=c

    #place les points à des lieux aléatoires
    def shuffle(self):
        self.positions=[]
        self.speeds=[]
        for i in range(self.cpx.dots):
            self.positions.append(np.random.uniform(-2,2,size=self.ambient))
            self.speeds.append(np.zeros(self.ambient))

    #actualise
    def update(self,times=1):
        for j in range(times):
            for i in range(self.cpx.dots):
                self.positions[i]+=self.speeds[i]*self.dt
                self.speeds[i]*=self.friction**self.dt
            if self.infl:
                self.update_dynamics()

    #place les points de sorte que i soit une inclusion géométrique
    def force_inclusion(self,i,parent):
        for d in range(self.cpx.dots):
            self.positions[d]=parent.positions[i(d)]

    #bord
    def manifold_border(self):
        b,i=self.cpx.manifold_border()
        res=EmbeddedManifold(b,self.ambient)
        res.force_inclusion(i,self)
        return res

    #subdivision barycentrique
    def bary_div(self):
        res=EmbeddedManifold(self.cpx.bary_div(),self.ambient)
        def action(F):
            if F.identity>=0:
                res.positions[F.identity]=self.face_center(F.face)
        self.cpx.faces_tree.propagate(action)
        return res

    #copie
    def clone(self):
        res=EmbeddedManifold(self.cpx.clone(),self.ambient)
        for d in self.dots:
            res.positions[d]=self.positions[d].copy()
            res.speeds[d]=self.speeds[d].copy()

        res.friction=self.friction
        res.dt=self.dt

        res.infl=self.infl
        res.spring=self.spring #raideur du ressort
        return res

    def animation(self):
        return Animation(self)
    def animate(self):
        self.animation().start()


    #description
    def update_dynamics(self):
        if self.cpx.dim<=0: #évite les exceptions
            return

        if len(self.links) == 0:
            #construction du dictionnaire links
            for F in self.cpx.get_faces(self.cpx.dim-1):
                self.links[F]=[]

            for F in self.cpx.get_faces(self.cpx.dim):
                for x in F:
                    G=F-{x}
                    if not F in self.links[G]:
                        self.links[G].append(F)

        V=self.volume()

        pulled=[]

        #tension superficielle à un simplexe du bord
        def pull(g):
            if g in pulled:
                return
            pulled.append(g)
            G=g.face
            d=len(G)-1
            if d>=1:
                delta=delta_ns[d]
                center=self.face_center(G)

                for p in G:
                    pos=self.positions[p]
                    self.speeds[p]-=self.spring*self.dt*(np.linalg.norm(pos-center)-delta)*(pos-center)

        #tension superficielle
        for f in self.border_info:
            f.propagate(pull)

        for l in links:
            cells = links[l]
            S=self.face_vol(l)
            C=self.face_center(l)
            basis = self.ortho_basis(l)

            if len(cells)==1: #bord
                cell=cells[0]
                #point opposé dans cell
                pt = (cell.face - f.face).pop()
                n=(C-self.positions[pt])

                for u in basis:
                    n = n - np.dot(n,u) * u

                nrm = np.linalg.norm(n)
                if abs(nrm)>=0.0000001:
                    n=n/nrm

                Vint = self.face_vol(cell)  #peut être optimisé en calculant S*h/n

                if Vint>0.00000001:
                    p_int = self.pv / Vint
                    #calcul du vecteur force
                    force = (p_int-self.p_ext)*S*n

                    #répartition de la force aux sommets de l
                    for x in l:
                        self.speeds[x]+=force/self.cpx.dim*self.dt #rajouter un coefficient de masse ?

            elif len(cells)==2:  #connexion
                cell1,cell2=cells[0],cells[1]
                #point opposé dans cell1
                pt = (cell1.face - f.face).pop()

                #vecteur normal de cell1 vers cell2
                n=(C-self.positions[pt])

                for u in basis:
                    n = n - np.dot(n,u) * u

                nrm = np.linalg.norm(n)
                if abs(nrm)>=0.0000001:
                    n=n/nrm

                V1 = self.face_vol(cell1)  #peut être optimisé en calculant S*h/n
                V2 = self.face_vol(cell2)

                if V1>0.00000001 and V2>0.000000001:
                    P1 = self.pv / V1
                    P2 = self.pv / V2
                    #calcul du vecteur force
                    force = (P1-P2)*S*n

                    #répartition de la force aux sommets de l
                    for x in l:
                        self.speeds[x]+=force/self.cpx.dim*self.dt #rajouter un coefficient de masse ?

    #volume d'une face
    def face_vol(self, F):
        return face_volume([self.positions[p] for p in F])

    #calcule le n-volume
    #peu de sens pour des complexes non connexes ou non homogenes
    def volume(self):
        res=0
        for F in self.cpx.nfaces[self.cpx.dim]:
            res+=self.face_vol(F)
        return res

    #fournit une base orthonormée de la (direction de la) face
    def ortho_basis(self,F):
        if len(F)==0:
            return []
        clone_F = F.copy()
        origin = self.positions[clone_F.pop()]
        vecs = []
        for x in clone_F:
            vecs.append(self.positions[x]-origin)
        return gram_process(vecs)



class Animation:
    def __init__(self,manifold):
        global cur_anim
        cur_anim=self
        self.manifold=manifold
        self.cpx=manifold.cpx
        self.background=(20,20,20)
        self.windw=1000
        self.windh=700
        self.active=False
        self.pressed_keys=[]
        self.old_pressed=[]
        self.window = pygame.display.set_mode((self.windw,self.windh))
        pygame.display.set_caption("Visualisation")

        self.zoom=100.0

        #pour la 3d
        self.origin=manifold.center()
        self.dir=np.array([1.0,0.0,0.0])
        self.eye_norm=10.0
        self.eye=spherical(0,math.pi/2,self.eye_norm)
        #vecteur unitaire tangent à la sphère qui contient l'oeil au point de l'oeil qui indique la droite de l'écran
        self.local_right=np.array([0.0,1.0,0.0])            #base orthonormale du plan focal
        self.local_top=np.cross(self.local_right,self.dir)  #
        self.focal=0.1
        #vecteur unitaire qui représente la direction de la vue
        #recalculés à chaque update
        self.focal_vector = np.zeros(3)

        self.disp_dots=False
        self.disp_lines=True
        self.disp_faces=True

        #point en cours de déplacement, -1 pour aucun point
        self.selected_dot=-1

        self.disp_num=False

        self.faces_alpha=50 #entre 0 et 255

        self.view_sub_complex=None  #complexe x morphisme simplicial

        def tog_dots(): self.disp_dots=not self.disp_dots
        def tog_lines(): self.disp_lines=not self.disp_lines
        def tog_faces(): self.disp_faces=not self.disp_faces
        def shuf(): self.manifold.shuffle()
        def tog_infl(): self.manifold.infl=not self.manifold.infl
        def cent(): self.manifold.force_center()
        def tog_num(): self.disp_num=not self.disp_num
        def div(): self.bary_div()
        def bord(): self.manifold_border()
        def yz_plane():
            self.dir=np.array([1.0,0.0,0.0])
            self.eye=self.origin + np.array([-self.eye_norm,0.0,0.0])
            self.local_right=np.array([0.0,1.0,0.0])
            self.local_top=np.cross(self.local_right,self.dir)
        def dsp_hom():
            self.cpx.disp_homology()

        self.buttons=[  Button(0,0,100,30,"Points",tog_dots),
                        Button(100,0,100,30,"Lignes",tog_lines),
                        Button(200,0,100,30,"Faces",tog_faces),
                        Button(300,0,120,30,"Mélanger",shuf),
                        Button(420,0,120,30,"Gonflage",tog_infl),
                        Button(520,0,120,30,"Plan yz",yz_plane),
                        Button(620,0,200,30,"Numérotation",tog_num),
                        Button(820,0,200,30,"Division",div),
                        Button(0,30,100,30,"Bord",bord),
                        Button(100,30,130,30,"Homologie", dsp_hom)
                        ]

    #ne retourne rien
    def bary_div(self):
        C=self.manifold.bary_div()
        C.infl=False
        self.manifold=C
        self.cpx=C.cpx
        self.origin=C.center()

    #ne retourne rien /!\ self.manifold doit être une variété plongée
    def manifold_border(self):
        B=self.manifold.manifold_border()
        self.set_manifold(B)

    def set_manifold(self,M,pos_morphism=None):
        manif=M if isinstance(M, EmbeddedManifold) else EmbeddedManifold(M,self.manifold.ambient)

        if pos_morphism != None:
            for i in range(manif.cpx.dots):
                manif.positions[i]=self.manifold.positions[pos_morphism(i)]
                manif.speeds[i]=np.zeros(manif.ambient)
        self.manifold=manif
        manif.infl=False
        self.cpx=manif.cpx
        self.origin=manif.center()

    #dessin du complexe
    def draw(self):
        #todo anti alias
        pygame.draw.rect(self.window, self.background,(0,0,self.windw,self.windh))
        self.compute_focal()
        screen_positions=[self.screen_pos(p) for p in self.manifold.positions]

        for b in self.buttons:
            b.draw(self.window)

        def dot_size(i):
            if self.selected_dot==i:
                return 10
            return 6

        if self.disp_dots:
            for i in range(self.cpx.dots):
                pos=screen_positions[i]
                if pos.size==2: #image réelle
                    draw_circle(self.window,(80,150,255),tuple(pos),dot_size(i))

        if self.disp_lines:
            for l in self.cpx.get_faces(1):
                if all([screen_positions[p].size==2 for p in l]):  #image réelle
                    points=[tuple(screen_positions[p]) for p in l]
                    pygame.draw.line(self.window,(80,150,255),points[0],points[1],2)
        if self.disp_faces:
            for l in self.cpx.get_faces(2):
                if all([screen_positions[p].size==2 for p in l]):  #image réelle
                    points=[tuple(screen_positions[p]) for p in l]
                    draw_polygon(self.window,(80,150,255,self.faces_alpha),[points[0],points[1],points[2]])

        if self.view_sub_complex!=None:
            K,f=self.view_sub_complex
            if self.disp_dots:
                for i in range(K.dots):
                    pos=screen_positions[f(i)]
                    if pos.size==2: #image réelle
                        draw_circle(self.window,(47,202,13),tuple(pos),dot_size(i))
            if self.disp_lines:
                for l in K.get_faces(1):
                    if all([screen_positions[f(p)].size==2 for p in l]):  #image réelle
                        points=[tuple(screen_positions[f(p)]) for p in l]
                        pygame.draw.line(self.window,(47,202,13),points[0],points[1],2)
            if self.disp_faces:
                for l in K.get_faces(2):
                    if all([screen_positions[f(p)].size==2 for p in l]):  #image réelle
                        points=[tuple(screen_positions[f(p)]) for p in l]
                        draw_polygon(self.window,(47,202,13,self.faces_alpha//2),[points[0],points[1],points[2]])

        #temporaire
        vec_A = self.screen_pos(self.manifold.vec_o)
        vec_B = self.screen_pos(self.manifold.vec_o + self.manifold.vec_v)
        draw_circle(self.window,(255,0,0),tuple(vec_A),5)
        draw_circle(self.window,(255,0,0),tuple(vec_B),5)
        try:
            pygame.draw.line(self.window,(255,0,0),tuple(vec_A),tuple(vec_B),2)
        except:
            pass

        if self.disp_num:
            for i in range(self.cpx.dots):
                pos=screen_positions[i]
                if pos.size==2: #image réelle
                    draw_text(self.window,(255,0,0),tuple(pos),str(i))

    def to_3d(self,u):
        if self.manifold.ambient>=3:
            return np.array([u[i] for i in range(0,self.manifold.ambient)])
        else:
            return np.array([0.0 for i in range(3-self.manifold.ambient)]+u.tolist())

    #actualise la cinétique
    def update(self):
        self.manifold.update()

        if self.selected_dot==-1:
            self.origin=self.to_3d(self.manifold.center())

        self.compute_focal()
        eye_speed=0.02
        dtheta=0.01

        if self.selected_dot>=0:
            if self.cpx.dots>=self.selected_dot+1:
                self.manifold.positions[self.selected_dot]=self.screen_to_space(pygame.mouse.get_pos())

        #position de l'oeil
        if 100 in self.pressed_keys:
            self.eye+=eye_speed*self.eye_norm*self.local_right #déplacement quasi linéaire
        if 113 in self.pressed_keys:
            self.eye-=eye_speed*self.eye_norm*self.local_right
        if 122 in self.pressed_keys:
            self.eye+=eye_speed*self.eye_norm*self.local_top
        if 115 in self.pressed_keys:
            self.eye-=eye_speed*self.eye_norm*self.local_top
        if 101 in self.pressed_keys:   #commande
            try:
                exec(input())
            except Exception as e:
                print(e)

        self.eye=self.origin+(self.eye-self.origin)/np.linalg.norm(self.eye-self.origin)*self.eye_norm  #force l'oeil à rester à la bonne distance de l'origine

        #direction du regard vers l'origine
        self.dir=(self.origin-self.eye)/np.linalg.norm(self.origin-self.eye)

        #projette la base locale sur le plan tangent à la sphère qui contient l'oeil au point de l'oeil, et le force à être unitaire
        self.local_right=self.local_right-np.dot(self.local_right,self.dir)*self.dir
        self.local_right=self.local_right/np.linalg.norm(self.local_right)
    def start(self):
        self.active=True
        while self.active:
            for event in pygame.event.get():
                if event.type==pygame.KEYDOWN:
                    if not event.key in self.pressed_keys:
                        self.pressed_keys.append(event.key)
                    #print(event.key)
                    if event.key==27:
                        self.active=False
                        pygame.display.quit()
                        pygame.quit()
                        exit()
                        break
                elif event.type==pygame.KEYUP:
                    if event.key in self.pressed_keys:
                        self.pressed_keys.remove(event.key)
                elif event.type == pygame.MOUSEBUTTONUP:
                    if event.button==1:
                        self.selected_dot=-1
                        pos = pygame.mouse.get_pos()
                        for b in self.buttons:
                            if b.x<=pos[0] and b.y<=pos[1] and b.x+b.w>=pos[0] and b.y+b.h>=pos[1]:
                                b.ev()
                elif event.type == pygame.MOUSEBUTTONDOWN:
                    if event.button == 4:
                        self.eye_norm-=0.5
                        self.eye_norm=max(self.eye_norm,0.3)
                    elif event.button == 5:
                        self.eye_norm+=0.5
                    elif event.button == 1:
                        for d in range(self.cpx.dots):
                            sp=self.screen_pos(self.manifold.positions[d])
                            if np.linalg.norm(sp-np.array(pygame.mouse.get_pos()))<=6:
                                self.selected_dot=d
                                break


            if self.active:
                self.draw()
                self.update()
                self.old_pressed=self.pressed_keys.copy()
                pygame.display.update()

    #calculs dédiés à la 3D
    def compute_focal(self):
        self.focal_vector=self.focal*self.dir
        self.local_top=np.cross(self.local_right,self.dir)

    #coordonnées centrées vers coordonnées écran
    def mid_to_screen(self,point):
        return np.array([self.zoom*point[0]+self.windw/2,-self.zoom*point[1]+self.windh/2])

    #position à l'écran
    def screen_pos(self, point):
        pt=self.to_3d(point)
        vec=pt-self.eye
        scal=np.dot(vec,self.focal_vector)
        if scal<0:
            return np.array([0]) #image virtuelle
        eye_to_intersection=self.focal**2/(scal)*vec
        xfoc = np.dot(eye_to_intersection,self.local_right)*self.zoom
        yfoc = np.dot(eye_to_intersection, self.local_top)*self.zoom
        return self.mid_to_screen([xfoc,yfoc])

    #une section de screen_pos
    def screen_to_space(self,point):
        foc=np.array([(point[0]-self.windw/2),-(point[1]-self.windh/2)])*(self.zoom**(-2))
        #car base orthonormée -> composantes covariantes = composantes contravariantes
        eye_to_target = foc[0]*self.local_right + foc[1]*self.local_right + self.eye_norm/2*self.dir
        return self.eye+eye_to_target

    def just_pressed(self,k):
        return k in self.pressed_keys and (not k in self.old_pressed)
