import numpy as np

# Definición del número áureo
phi = (1 + np.sqrt(5)) / 2

# --- Utilidades de reflexión ---
def reflection(v, n):
    """
    Refleja el vector v en el hiperplano ortogonal a n.
    """
    n = n / np.linalg.norm(n)
    return v - 2 * np.dot(v, n) * n

# --- Definir algunos generadores (planos de reflexión) para H4 ---
# (Estos vectores normales son escogidos de forma que, aplicando las reflexiones, se recupere la simetría de H4.
#  NOTA: Los valores aquí presentados son una aproximación; para una generación exacta se debe derivar
#  cuidadosamente a partir del diagrama de Coxeter de H4.)
normals = np.array([
    [0, 0, 0, 1],
    [0, 0, -phi, 1],
    [0, -phi, 1, 0],
    [-phi, 1, 0, 0]
], dtype=float)

# --- Generar vértices mediante el método de "reflejar y acumular" ---
# Usamos un vértice semilla conocido (por ejemplo, uno de la forma (0,0,2,2); se sabe que
# los vértices del 120-cell pueden estar escalados para que todos tengan la misma norma).
v0 = np.array([0, 0, 2, 2], dtype=float)
# Inicializamos el conjunto (usamos tuplas redondeadas para evitar duplicados numéricos)
vertices = {tuple(np.round(v0, 8))}

# Iteramos aplicando las reflexiones generadoras hasta que el conjunto se estabilice.
changed = True
while changed:
    changed = False
    new_vs = set()
    for v in list(vertices):
        v_arr = np.array(v)
        for n in normals:
            r = reflection(v_arr, n)
            # Redondeamos para evitar problemas de precisión
            r_tuple = tuple(np.round(r, 8))
            if r_tuple not in vertices:
                new_vs.add(r_tuple)
    if new_vs:
        vertices.update(new_vs)
        changed = True

# Convertir a arreglo numpy
vertices = np.array(list(vertices))
print("Número de vértices generados:", len(vertices))

# --- Ajustar la escala (opcional) ---
# Muchos autores escalan el 120-cell para que todos sus vértices estén sobre una esfera de radio √8.
# Si lo deseas, puedes reescalar:
target_radius = np.sqrt(8)
current_radii = np.linalg.norm(vertices, axis=1)
scale = target_radius / np.mean(current_radii)
vertices *= scale
print("Radio aproximado (debe ser ~√8):", np.mean(np.linalg.norm(vertices, axis=1)))

# --- Calcular la lista de aristas ---
# Se asume que dos vértices están conectados por una arista si su distancia es igual a la longitud de arista.
# Para determinarla, calculamos el mínimo valor no nulo de las distancias entre vértices.
tol = 1e-6
n = len(vertices)
d_min = None
# Primero, recolectamos las distancias (evitando la distancia 0)
for i in range(n):
    for j in range(i+1, n):
        d = np.linalg.norm(vertices[i] - vertices[j])
        if d > tol:
            if d_min is None or d < d_min:
                d_min = d

print("Longitud de arista aproximada:", d_min)

# Ahora, generamos la lista de aristas considerando como arista toda pareja cuya distancia difiera del mínimo
# por menos de una tolerancia
edges = []
for i in range(n):
    for j in range(i+1, n):
        d = np.linalg.norm(vertices[i] - vertices[j])
        if abs(d - d_min) < tol:
            edges.append((i, j))

print("Número de aristas generadas:", len(edges))

# Para comprobar, en un 120-cell se esperan 600 vértices y 1200 aristas.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def generate_120_cell_projection():
    # Generamos vértices de un 120-cell en 4D (simplificado para este ejemplo)
    vertices_4d = []
    phi = (1 + np.sqrt(5)) / 2  # Proporción áurea
    
    # Coordenadas base para un 120-cell (Esta es una versión muy simplificada)
    for i in [-1, 1]:
        for j in [-1, 1]:
            for k in [-1, 1]:
                for l in [-1, 1]:
                    vertices_4d.append([i, j, k, l])
    
    # Proyección estereográfica a 3D
    vertices_3d = []
    for v in vertices_4d:
        # Proyección simple de 4D a 3D
        factor = 1 / (2 - v[3])
        x = v[0] * factor
        y = v[1] * factor
        z = v[2] * factor
        vertices_3d.append([x, y, z])
    
    return np.array(vertices_3d)

# Visualización
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
vertices = generate_120_cell_projection()
ax.scatter(vertices[:, 0], vertices[:, 1], vertices[:, 2], c='purple', alpha=0.6, s=20)
plt.title("Proyección 3D de un 120-cell (simplificada)")
plt.savefig('120_cell_projection.png', dpi=300)
plt.show()