Parametros_de_secciones/parametric_sweep.py

"""
Módulo de Barrido Paramétrico para Secciones de Acero

Permite realizar barridos paramétricos dinámicos sobre secciones IPE y Puente nuevo,
fijando algunos parámetros y variando otros sistemáticamente.
"""

import numpy as np
import pandas as pd
import itertools
from typing import Dict, List, Tuple
import math


def generate_ipe_points(H, b, tf, tw):
    """Genera los puntos de una sección IPE a partir de parámetros normalizados"""
    h_alma = H - 2*tf
    x_alma_ini = (b - tw) / 2
    x_alma_fin = (b + tw) / 2

    points = np.array([
        [0, 0],
        [b, 0],
        [b, tf],
        [x_alma_fin, tf],
        [x_alma_fin, H - tf],
        [b, H - tf],
        [b, H],
        [0, H],
        [0, H - tf],
        [x_alma_ini, H - tf],
        [x_alma_ini, tf],
        [0, tf],
        [0, 0],
    ])

    return points


def generate_puente_nuevo_points(h, b, tf, tw, ha, ta, tr, theta):
    """Genera los puntos de la sección tipo puente nuevo"""
    hr = (b-tw) / 2 * math.tan(math.radians(theta))
    hl = tr / math.cos(math.radians(theta))
    hr_ = (b/2 - ta - tw/2)* math.tan(math.radians(theta))

    points = np.array([
        [0, 0],
        [b, 0],
        [b, tf + ha],
        [(b + tw)/2, tf + ha + hr],
        [(b + tw)/2, tf + ha + hr - hl],
        [b - ta, tf + ha + hr - hl - hr_],
        [b - ta, tf],
        [(b + tw)/2, tf],
        [(b + tw)/2, h - tf],
        [b - ta, h - tf],
        [b - ta, h - tf - ha - hr + hl + hr_],
        [(b + tw)/2, h - tf - ha - hr + hl],
        [(b + tw)/2, h - tf - ha - hr],
        [b, h - tf - ha],
        [b, h],
        [0, h],
        [0, h - tf - ha],
        [(b - tw)/2, h - tf - ha - hr],
        [(b - tw)/2, h - tf - ha - hr + hl],
        [ta, h - tf - ha - hr + hl + hr_],
        [ta, h - tf],
        [(b - tw)/2, h - tf],
        [(b - tw)/2, tf],
        [ta, tf],
        [ta, tf + ha + hr - hl - hr_],
        [(b - tw)/2, tf + ha + hr - hl],
        [(b - tw)/2, tf + ha + hr],
        [0, tf + ha],
        [0, 0],
    ])

    return points


def calculate_section_properties(points):
    """Calcula propiedades de la sección a partir de puntos"""
    # Constantes de material
    DENS = 7850  # kg/m3
    FY = 355  # MPa
    GAMMAS = 1.05
    EYOUNG = 210000  # MPa
    NU = 0.3
    FYD = FY * 10**6 / GAMMAS

    npuntos = len(points)
    px = points[:, 0]
    py = points[:, 1]

    bmax = np.amax(px) - np.amin(px)
    hmax = np.amax(py) - np.amin(py)

    # Perímetro
    long_i = np.zeros(npuntos - 1)
    for i in range(npuntos - 1):
        long_i[i] = ((points[i+1, 0] - points[i, 0])**2 + (points[i+1, 1] - points[i, 1])**2)**(1/2)
    perimetro = abs(sum(long_i))

    # Área
    area_i = np.zeros(npuntos - 1)
    for i in range(npuntos - 1):
        area_i[i] = (points[i+1, 0] - points[i, 0]) * (points[i+1, 1] + points[i, 1]) / 2
    area = abs(sum(area_i))

    # Peso
    peso = area * DENS

    # Centro de gravedad
    cdg_i = np.zeros([npuntos, 2])
    for i in range(npuntos - 1):
        h1 = points[i, 1]
        h2 = points[i+1, 1]
        b = points[i+1, 0] - points[i, 0]
        d = points[i, 0]

        if h1 + h2 == 0:
            cdg_i[i, 1] = 0
        else:
            cdg_i[i, 1] = 1/3 * (h1*h1 + h1*h2 + h2*h2) / (h1 + h2)

        if h1 + h2 == 0:
            cdg_i[i, 0] = d + b/2
        else:
            cdg_i[i, 0] = d + b/3 * (h1 + 2*h2) / (h1 + h2)

    statico_i = np.zeros([npuntos, 2])
    for i in range(npuntos - 1):
        statico_i[i, 1] = area_i[i] * cdg_i[i, 1]
        statico_i[i, 0] = area_i[i] * cdg_i[i, 0]

    cdg = sum(statico_i) / sum(area_i)
    xg = cdg[0]
    yg = cdg[1]

    # Fibras más alejadas
    v1y = np.amax(py) - yg
    v2y = np.amin(py) - yg
    v1x = np.amax(px) - xg
    v2x = np.amin(px) - xg

    # Momentos de inercia
    inercia_i = np.zeros([npuntos, 3])
    for i in range(npuntos - 1):
        h1 = points[i, 1]
        h2 = points[i+1, 1]
        b = points[i+1, 0] - points[i, 0]
        d = points[i, 0]
        xgi = cdg_i[i, 0]
        ygi = cdg_i[i, 1]
        ai = area_i[i]

        if h2 >= h1:
            ixcuad_G_local = 1/12 * b * (h1**3) + b * h1 * (h1/2 - ygi)**2
            ixtriang_G_loc = 1/36 * b * (h2-h1)**3 + 1/2 * b * (h2-h1) * ((2*h1+h2)/3 - ygi)**2
        else:
            ixcuad_G_local = 1/12 * b * (h2**3) + b * h2 * (h2/2 - ygi)**2
            ixtriang_G_loc = 1/36 * b * (h1-h2)**3 + 1/2 * b * (h1-h2) * ((2*h2+h1)/3 - ygi)**2
        inercia_i[i, 0] = ixcuad_G_local + ixtriang_G_loc + ai * (yg - ygi)**2

        if h2 >= h1:
            iycuad = 1/12 * h1 * b**3 + h1 * b * (b/2 + d - xgi)**2
            iytrian = 1/36 * (h2-h1) * b**3 + 1/2 * b * (h2-h1) * (2/3*b + d - xgi)**2
        else:
            iycuad = 1/12 * h2 * b**3 + h2 * b * (b/2 + d - xgi)**2
            iytrian = 1/36 * (h1-h2) * b**3 + 1/2 * b * (h1-h2) * (1/3*b + d - xgi)**2
        inercia_i[i, 1] = iycuad + iytrian + ai * (xg - xgi)**2

        if h2 >= h1:
            pxygcuadrado = b * h1 * (-h1/2 + ygi) * (-d - b/2 + xgi)
            pxytriangulo = b*b * (h2-h1)**2 / 72 + b * (h2-h1) / 2 * (-(h2-h1)/3 - h1 + ygi) * (-d - 2/3*b + xgi)
        else:
            pxygcuadrado = b * h2 * (-h2/2 + ygi) * (-d - b/2 + xgi)
            pxytriangulo = -b*b * (h1-h2)**2 / 72 + b * (h1-h2) / 2 * (-(h1-h2)/3 - h2 + ygi) * (-d - 1/3*b + xgi)
        inercia_i[i, 2] = pxygcuadrado + pxytriangulo + ai * (-xg + xgi) * (-yg + ygi)

    ig = sum(inercia_i)
    ixg = abs(ig[0])
    iyg = abs(ig[1])
    pxyg = ig[2]
    if sum(area_i) >= 0:
        pxyg = pxyg
    else:
        pxyg = -pxyg

    # Radios de giro
    rx = (ixg / area)**0.5
    ry = (iyg / area)**0.5

    # Ejes principales
    ic = (ixg + iyg) / 2
    ir = (((ixg - iyg)/2)**2 + pxyg**2)**0.5
    imax = ic + ir
    imin = ic - ir
    rmax = (imax / area)**0.5
    rmin = (imin / area)**0.5

    # Conversión a unidades prácticas (cm, cm2, cm3, cm4, kN)
    pot = 2
    area_cm2 = area * 10**(pot*2)
    ixg_cm4 = ixg * 10**(pot*4)
    iyg_cm4 = iyg * 10**(pot*4)
    imax_cm4 = imax * 10**(pot*4)
    imin_cm4 = imin * 10**(pot*4)

    return {
        'area': area_cm2,
        'ixg': ixg_cm4,
        'iyg': iyg_cm4,
        'imax': imax_cm4,
        'imin': imin_cm4,
        'peso': peso,
    }


class ParametricSweep:
    """Ejecuta barridos paramétricos de secciones"""

    def __init__(self, section_type: str, fixed_params: Dict, sweep_configs: Dict):
        """
        Args:
            section_type: "ipe" o "puente_nuevo"
            fixed_params: {param_name: value, ...}
            sweep_configs: {param_name: {"min": v, "max": v, "steps": n}, ...}
        """
        self.section_type = section_type
        self.fixed_params = fixed_params.copy()
        self.sweep_configs = sweep_configs.copy()
        self.results_df = None
        self.pareto_front = None

    def _is_valid_geometry(self, params: Dict) -> bool:
        """Valida que la geometría sea válida para IPE"""
        if self.section_type == "ipe":
            H = params.get('H', 0)
            b = params.get('b', 0)
            tf = params.get('tf', 0)
            tw = params.get('tw', 0)

            # Validaciones básicas
            if H <= 0 or b <= 0 or tf <= 0 or tw <= 0:
                return False
            # Espesor de flange no puede ser más de la mitad de altura
            if tf > H / 2:
                return False
            # Espesor de alma no puede ser mayor que ancho
            if tw > b:
                return False
            # Altura del alma positiva
            if H - 2*tf <= 0:
                return False

        return True

    def generate_sweep_matrix(self) -> List[Dict]:
        """Genera todas las combinaciones de parámetros a barrer"""
        # Obtener parámetros a barrer
        sweep_params = list(self.sweep_configs.keys())

        # Generar valores para cada parámetro
        param_values = {}
        for param in sweep_params:
            config = self.sweep_configs[param]
            values = np.linspace(config['min'], config['max'], config['steps'])
            param_values[param] = values

        # Generar todas las combinaciones
        combinations = []
        for combo in itertools.product(*[param_values[p] for p in sweep_params]):
            params = self.fixed_params.copy()
            for i, param in enumerate(sweep_params):
                params[param] = combo[i]

            if self._is_valid_geometry(params):
                combinations.append(params)

        return combinations

    def execute_sweep(self) -> pd.DataFrame:
        """Ejecuta el barrido y retorna DataFrame con resultados"""
        combinations = self.generate_sweep_matrix()

        if not combinations:
            raise ValueError("No hay combinaciones válidas en el barrido")

        results = []

        for params in combinations:
            # Generar puntos según tipo
            if self.section_type == "ipe":
                points = generate_ipe_points(
                    params['H'], params['b'], params['tf'], params['tw']
                )
            elif self.section_type == "puente_nuevo":
                points = generate_puente_nuevo_points(
                    params['h'], params['b'], params['tf'], params['tw'],
                    params['ha'], params['ta'], params['tr'], params['theta']
                )
            else:
                continue

            # Calcular propiedades
            props = calculate_section_properties(points)

            # Construir fila de resultados
            row = params.copy()
            row.update(props)

            results.append(row)

        self.results_df = pd.DataFrame(results)

        # Normalizar valores para cálculo de Pareto
        self.results_df['peso_norm'] = (
            self.results_df['peso'] / self.results_df['peso'].max()
        )
        self.results_df['ixg_norm'] = (
            self.results_df['ixg'] / self.results_df['ixg'].max()
        )
        self.results_df['iyg_norm'] = (
            self.results_df['iyg'] / self.results_df['iyg'].max()
        )

        return self.results_df

    def find_pareto_front(self) -> pd.DataFrame:
        """Ordena soluciones por eficiencia: (Ix+Iy) / peso
        Mayor eficiencia = mejor inercia con menos peso
        """
        if self.results_df is None:
            raise ValueError("Ejecutar execute_sweep() primero")

        # Calcular eficiencia para cada solución
        self.results_df['eficiencia'] = (
            (self.results_df['ixg'] + self.results_df['iyg']) /
            (self.results_df['peso'] * 1000)
        )

        # Ordenar por eficiencia descendente (mejor primero)
        self.pareto_front = self.results_df.sort_values('eficiencia', ascending=False)

        return self.pareto_front

    def get_results_dataframe(self) -> pd.DataFrame:
        """Retorna DataFrame con todos los resultados"""
        return self.results_df

    def get_pareto_front(self) -> pd.DataFrame:
        """Retorna DataFrame con Pareto front"""
        return self.pareto_front