import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import ConvexHull
import warnings
warnings.filterwarnings('ignore')

# Leer el archivo manualmente
data = []
with open('resultados1.txt', 'r', encoding='utf-8') as f:
    header = f.readline().strip().split('\t')
    for line in f:
        values = line.strip().split()  # Split por espacios/tabulaciones
        if len(values) >= 15:  # Tenemos 15 columnas
            data.append(values)

# Encontrar índices de columnas
tf_idx = 2  # tf es la columna 3
tw_idx = 3  # tw es la columna 4
flecha_idx = 8  # Flecha_Media es la columna 9

# Convertir datos a números
def convert_value(x):
    """Convierte string a float, manejando formato con coma decimal"""
    try:
        return float(x.replace(',', '.').strip())
    except:
        return np.nan

tf_list = []
tw_list = []
flecha_list = []

for row in data:
    tf_val = convert_value(row[tf_idx])
    tw_val = convert_value(row[tw_idx])
    flecha_val = convert_value(row[flecha_idx])
    
    if not (np.isnan(tf_val) or np.isnan(tw_val) or np.isnan(flecha_val)):
        tf_list.append(tf_val)
        tw_list.append(tw_val)
        flecha_list.append(flecha_val)

tf = np.array(tf_list)
tw = np.array(tw_list)
flecha_media = np.array(flecha_list)

# RESTRICCIÓN DE DISEÑO
LIMITE_FLECHA = -25.0  # Valores de flecha por debajo de este límite NO son válidos

# Separar puntos válidos e inválidos
valid_idx = flecha_media >= LIMITE_FLECHA
invalid_idx = flecha_media < LIMITE_FLECHA

tf_valid = tf[valid_idx]
tw_valid = tw[valid_idx]
flecha_valid = flecha_media[valid_idx]

tf_invalid = tf[invalid_idx]
tw_invalid = tw[invalid_idx]
flecha_invalid = flecha_media[invalid_idx]

print("="*70)
print("ANÁLISIS DE FRONTERA DE PARETO - PROBLEMA DE MINIMIZACIÓN 3D")
print("="*70)
print(f"\n✓ Datos cargados exitosamente")
print(f"  Total de puntos: {len(tf)}")
print(f"\n  Rango de tf (espesor ala):       [{tf.min():.6f}, {tf.max():.6f}]")
print(f"  Rango de tw (espesor alma):     [{tw.min():.6f}, {tw.max():.6f}]")
print(f"  Rango de Flecha_Media:          [{flecha_media.min():.2f}, {flecha_media.max():.2f}]")

print(f"\n⚠ RESTRICCIÓN DE DISEÑO: Flecha_Media >= {LIMITE_FLECHA} mm")
print(f"  Puntos VÁLIDOS:                 {len(tf_valid)} ({100*len(tf_valid)/len(tf):.1f}%)")
print(f"  Puntos INVÁLIDOS:               {len(tf_invalid)} ({100*len(tf_invalid)/len(tf):.1f}%)")
print(f"  Rango flecha válida:            [{flecha_valid.min():.2f}, {flecha_valid.max():.2f}]")
print(f"  Rango flecha inválida:          [{flecha_invalid.min():.2f}, {flecha_invalid.max():.2f}]")


# Calcular la frontera de Pareto con algoritmo mejorado (solo para puntos VÁLIDOS)
def compute_pareto_frontier(tf, tw, flecha):
    """
    Calcula los puntos que pertenecen a la frontera de Pareto.
    Objetivo: minimizar tf, tw y |flecha|
    """
    n = len(tf)
    is_pareto = np.ones(n, dtype=bool)
    
    # Usar valor absoluto para la comparación
    flecha_abs = np.abs(flecha)
    
    for i in range(n):
        if not is_pareto[i]:
            continue
        for j in range(n):
            if i != j and is_pareto[j]:
                # El punto j domina al i si: j <= i en todos los criterios
                # y j < i en al menos uno
                dom = (tf[j] <= tf[i]) and (tw[j] <= tw[i]) and (flecha_abs[j] <= flecha_abs[i])
                strict = (tf[j] < tf[i]) or (tw[j] < tw[i]) or (flecha_abs[j] < flecha_abs[i])
                
                if dom and strict:
                    is_pareto[i] = False
                    break
    
    return is_pareto

# Calcular Pareto solo para puntos VÁLIDOS
pareto_mask_valid = compute_pareto_frontier(tf_valid, tw_valid, flecha_valid)
pareto_tf_valid = tf_valid[pareto_mask_valid]
pareto_tw_valid = tw_valid[pareto_mask_valid]
pareto_flecha_valid = flecha_valid[pareto_mask_valid]

print(f"\n✓ Frontera de Pareto calculada (región válida)")
print(f"  Puntos en la frontera: {len(pareto_tf_valid)} ({100*len(pareto_tf_valid)/len(tf_valid):.1f}% de válidos)")


# Crear figura 3D de alta calidad
fig = plt.figure(figsize=(16, 12))
ax = fig.add_subplot(111, projection='3d')

# Graficar el PLANO DE RESTRICCIÓN en z = LIMITE_FLECHA
xx = np.linspace(tf.min(), tf.max(), 20)
yy = np.linspace(tw.min(), tw.max(), 20)
XX, YY = np.meshgrid(xx, yy)
ZZ = np.full_like(XX, LIMITE_FLECHA)  # Plano en z = -25

# Plotear el plano de restricción
plane = ax.plot_surface(XX, YY, ZZ, alpha=0.3, color='orange', label='Plano de restricción')
# Agregar wireframe al plano para hacerlo más visible
ax.plot_wireframe(XX, YY, ZZ, color='darkorange', alpha=0.5, linewidth=0.5)

# Graficar puntos INVÁLIDOS (debajo del plano)
scatter_invalid = ax.scatter(tf_invalid, tw_invalid, flecha_invalid, 
                             c='red', 
                             marker='x', 
                             s=100, 
                             alpha=0.3, 
                             label=f'Puntos inválidos (n={len(tf_invalid)})',
                             linewidth=2)

# Graficar puntos VÁLIDOS (encima del plano)
scatter_valid = ax.scatter(tf_valid, tw_valid, flecha_valid, 
                           c='lightgreen', 
                           marker='o', 
                           s=80, 
                           alpha=0.6, 
                           label=f'Puntos válidos (n={len(tf_valid)})',
                           edgecolors='darkgreen',
                           linewidth=0.5)

# Intentar graficar la frontera de Pareto solo para puntos válidos
if len(tf_valid) >= 4:
    try:
        # Crear matriz de coordenadas para la frontera de Pareto (solo válidos)
        points_pareto = np.column_stack((tf_valid, tw_valid, np.abs(flecha_valid)))
        
        # Calcular la envolvente convexa (convex hull)
        hull = ConvexHull(points_pareto)
        
        # Graficar la superficie
        for simplex in hull.simplices:
            triangle = points_pareto[simplex]
            # Graficar los bordes del triángulo
            for i in range(3):
                p1 = triangle[i]
                p2 = triangle[(i+1)%3]
                ax.plot([p1[0], p2[0]], 
                       [p1[1], p2[1]], 
                       [p1[2], p2[2]], 
                       'b-', alpha=0.2, linewidth=0.5)
        
        print(f"  Envolvente convexa (región válida): {len(hull.simplices)} triángulos")
    except Exception as e:
        print(f"  ⚠ Envolvente convexa: No se pudo calcular ({str(e)[:50]})")

# Etiquetas de ejes
ax.set_xlabel('tf - Espesor del ala (mm)', fontsize=12, fontweight='bold')
ax.set_ylabel('tw - Espesor del alma (mm)', fontsize=12, fontweight='bold')
ax.set_zlabel('Flecha_Media - Deflexión (mm)', fontsize=12, fontweight='bold')
ax.set_title(f'Frontera de Pareto con Restricción: Flecha_Media ≥ {LIMITE_FLECHA} mm\ntf vs tw vs Deflexión Media', 
             fontsize=14, fontweight='bold', pad=20)

# Leyenda mejorada
ax.legend(loc='upper left', fontsize=11, framealpha=0.95)

# Mejorar vista
ax.view_init(elev=25, azim=45)
ax.grid(True, alpha=0.3)

# Agregar línea de referencia en el plano de restricción
ax.text2D(0.02, 0.98, f'Plano de restricción: Flecha = {LIMITE_FLECHA} mm\n(Naranja)', 
          transform=ax.transAxes, fontsize=10, verticalalignment='top',
          bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))

plt.tight_layout()
plt.savefig('grafico_3d_pareto_con_restriccion.png', dpi=300, bbox_inches='tight')
print(f"\n✓ Gráfico guardado: 'grafico_3d_pareto_con_restriccion.png'")

# Mostrar estadísticas de los puntos de Pareto
print("\n" + "="*70)
print("ESTADÍSTICAS DE LA FRONTERA DE PARETO (REGIÓN VÁLIDA)")
print("="*70)
print(f"{'Parámetro':<20} {'Mínimo':<15} {'Máximo':<15} {'Promedio':<15}")
print("-"*70)
print(f"{'tf':<20} {pareto_tf_valid.min():<15.6f} {pareto_tf_valid.max():<15.6f} {pareto_tf_valid.mean():<15.6f}")
print(f"{'tw':<20} {pareto_tw_valid.min():<15.6f} {pareto_tw_valid.max():<15.6f} {pareto_tw_valid.mean():<15.6f}")
print(f"{'Flecha_Media':<20} {pareto_flecha_valid.min():<15.2f} {pareto_flecha_valid.max():<15.2f} {pareto_flecha_valid.mean():<15.2f}")

# Mostrar detalles de los puntos de Pareto VÁLIDOS
print("\n" + "="*70)
print(f"PUNTOS EN LA FRONTERA DE PARETO - REGIÓN VÁLIDA (primeros 30 de {len(pareto_tf_valid)})")
print("="*70)
print(f"{'#':<4} {'tf':<12} {'tw':<12} {'Flecha_Media':<15}")
print("-"*70)
for i in range(min(30, len(pareto_tf_valid))):
    print(f"{i+1:<4} {pareto_tf_valid[i]:<12.6f} {pareto_tw_valid[i]:<12.6f} {pareto_flecha_valid[i]:<15.2f}")

if len(pareto_tf_valid) > 30:
    print(f"\n... y {len(pareto_tf_valid) - 30} puntos más")

print("\n" + "="*70)
plt.show()