AMIDIQ/src/analysis.py

import  numpy   as  np

def calcular_pendiente(df, x_col, y_col):
    """
    Calcular pendiente de ajuste lineal 
    """
    x   =   df[x_col].values
    y   =   df[y_col].values

    m, b    =   np.polyfit(x, y, 1)

    return m, b

def filtrar_por_limites(df, x_col, v_min=None, v_max=None):

    out =   df.copy()

    if v_min is not None:
        out =   out[out[x_col] >= v_min]
    if v_max is not None:
        out =   out[out[x_col] <= v_max]
    
    return out.reset_index(drop=True)

def calcular_r2(df, x_col, y_col):
    """
    Calcular coeficiente de determinación R^2.
    """

    x   =   df[x_col].values
    y   =   df[y_col].values

    m, b    =   np.polyfit(x, y, 1)

    y_pred  =   m * x + b

    ss_res  =   np.sum((y - y_pred) ** 2)
    ss_tot  =   np.sum((y - np.mean(y)) ** 2)

    r2  =   1 - (ss_res / ss_tot)

    return r2

def derivada_numerica(df, x_col, y_col):
    """
    Calcula primera derivada númerica dy/dx
    """

    x   =   np.asanyarray(df[x_col].values, dtype=float)
    y   =   np.asanyarray(df[y_col].values, dtype=float)

    idx     =   np.argsort(x)
    x       =   x[idx]
    y       =   y[idx]

    return x, y, np.gradient(y,x)

def segunda_derivada(df, x_col, y_col):
    """
    Calcula segunda derivada númerica d^2y/d^2x
    """

    x   =   np.asanyarray(df[x_col].values, dtype=float)
    y   =   np.asanyarray(df[y_col].values, dtype=float)

    idx     =   np.argsort(x)
    x       =   x[idx]
    y       =   x[idx]

    dy_dx   =   np.gradient(y, x)
    d2y_dx2 =   np.gradient(dy_dx, x)

    return x, y, dy_dx, d2y_dx2

def find_local_peaks(signal):
    """
    Encuentra indices de maximos locales simples
    """
    peaks   =   []

    for i in range(1, len(signal) - 1):
        if signal[i] > signal[i - 1] and signal[i] > signal[i + 1]:
            peaks.append(i)

    return np.array(peaks, dtype=int)

def variación_resistencia(df, col_res="Resistencia"):
    """
    Calcula variación relativa de resistencia.
    """
    R   =   df[col_res].values

    R0  =   R[0]

    delta_R =   (R - R0) / R0

    return delta_R

def comparar_pendientes(data, rondas, experimentos, x_col, y_col):
    """
    Compara pendientes entre diferentes experimentos o rondas
    """

    resultados = {}

    for r in rondas:
        for e in experimentos:
            df      =   data[(r, e)]
            m, b    =   calcular_pendiente(df, x_col, y_col)
            resultados[(r, e)]  =   m
    
    return resultados

def zonas_manual(df, x_col, limites):
    """
    Divide el dataset en zonas usando límites definidos manualmente
    """ 
    zonas   =   {}

    for i in range(len(limites) - 1):

        v_min   =   limites[i]
        v_max   =   limites[i + 1]

        zona_df =   df[(df[x_col] >= v_min) & (df[x_col] < v_max)]

        zonas[f"Zona_{i+1}"]    =   zona_df
    
    return  zonas

def detectar_zonas_auto(df, x_col, y_col, n_ruido=8):

    """
    Detecta automaticamente los limites de 4 zonas:

    1) No conduccion 
    2) Conduccion
    3) Generacion H
    4) Saturacion
    """

    x, y, dy_dx, d2y_d2x = segunda_derivada(df, x_col=x_col, y_col=y_col)

    #===============================
    #Estimacion automatica del ruido 
    #===============================
    n_ruido     =   min(n_ruido, len(y))
    ruido_base  =   np.std(y[:n_ruido])

    #Umbral fisico basado en ruido experimental
    tol_corriente   =   max(3 * ruido_base, 1e-12)

    #=============================
    # Fin de no conduccion
    #=============================
    idx_no_conduccion = None
    for i in range(len(y)):
        if np.abs(y[i]) > tol_corriente:
            idx_no_conduccion = i
            break

    if idx_no_conduccion is None:
        return{
            "Zona no Conduccion": None,
            "Zona Conduccion": None,
            "Generacion": None,
            "Saturacion": None,
        }
    
    #=====================
    # Picos de Curvatura
    #=====================
    abs_d2  =   np.abs(d2y_d2x)
    peaks   =   find_local_peaks(abs_d2)

    # Nos quedamos solo con picos posteriores a no conduccion
    peaks   =   peaks[peaks > idx_no_conduccion]

    # Si no hay suficientes picos, usamos la estrategia de respaldo
    if len(peaks) < 2:
        idx_conduccion = idx_no_conduccion

        # Generación: punto donde la pendiente alcanza su maximo crecimiento relativo
        idx_generacion  =   np.argmax(np.abs(dy_dx[idx_no_conduccion:])) + idx_no_conduccion
        # Saturación: máximo de segunda derivada despues de generacion
        idx_saturacion  =   np.argmax(abs_d2[idx_generacion]) + idx_generacion
    else:
        # Primer pico importante de curvatura - transicion a conduccion
        idx_conduccion  =   peaks[0]

        # Segundo pico importante - inicio de generacion 
        idx_generacion  =   peaks[1]

        # Saturacion: buscar el mayro pico de curvatura despues de generacion 
        peaks_after_gen = peaks[peaks > idx_generacion]
        if len(peaks_after_gen) > 0:
            idx_saturacion  =   peaks_after_gen[np.argmax(abs_d2[peaks_after_gen])]
        else:
            idx_saturacion  =   np.argmax(abs_d2[idx_generacion:]) + idx_generacion
    
    limites =   {
        "Zona_no_conduccion":   float(x[idx_no_conduccion]),
        "Zona_conduccion"   :   float(x[idx_conduccion]),
        "Zona_Generacion"   :   float(x[idx_generacion]),
        "Zona_Saturacion"   :   float(x[idx_saturacion]),
    }

    return limites