circust.fitting

FMM

circust/fitting/fmm.py

Frecuency Modulated Mobius (FMM) — ajuste ritmico de un solo componente.

Modelo matematico

y(t) = M + A * cos(beta + 2*atan(omega*tan((t - alpha)/2)))

El termino interno Phi(t) = beta + 2*atan(omega*tan((t - alpha)/2)) es una transformacion de Mobius del circulo. Cuando omega = 1 se reduce a un simple desfase. Cuando omega -> 0 la forma de onda se aproxima a un coseno (Cosinor). La ventaja clave sobre Cosinor es que omega controla la asimetria de la forma de onda: el modelo puede representar picos asimetricos presentes en datos circadianos reales.

Parametros

M mesor — nivel medio A amplitud — siempre > 0 alpha — localizacion de la region del pico en [0, 2*pi) beta — parametro de forma en [0, 2*pi) omega — asimetria en (0, omegaMax]

Algoritmo de estimacion

Paso 1 — Busqueda en rejilla

Para cada combinacion (alpha, omega) en una rejilla, se fijan esos dos parametros y se estiman M, A, beta por OLS sobre el eje temporal transformado por Mobius. Se selecciona el candidato con menor RSS que tambien satisfaga las restricciones de estabilidad de amplitud.

Paso 2 — Refinamiento Nelder-Mead

Se pule la estimacion del Paso 1 minimizando el RSS sobre los cinco parametros simultaneamente, sujeto a las mismas restricciones de estabilidad (puntos infactibles devuelven infinito).

Los Pasos 1+2 se repiten num_reps veces, cada vez estrechando la rejilla alrededor de la mejor estimacion anterior.

Devuelve

FitResult

Claves de params: M, A, alpha, beta, omega

class circust.fitting.fmm.FMMModel(length_alpha_grid=48, length_omega_grid=24, omega_max=1.0, num_reps=3)[source]

Bases: RhythmModel

Ajustador ritmico FMM (Frecuency Modulated Mobius).

Ajusta y(t) = M + A*cos(beta + 2*atan(omega*tan((t - alpha)/2))) usando un algoritmo de dos pasos: busqueda en rejilla + Nelder-Mead.

Parameters:

  • length_alpha_grid (int) – Numero de valores de alpha en la rejilla. Valor R por defecto: 48.

  • length_omega_grid (int) – Numero de valores de omega en la rejilla. Valor R por defecto: 24.

  • omega_max (float) – Cota superior de omega. Valor R por defecto: 1.0.

  • num_reps (int) – Numero de iteraciones de refinamiento de rejilla. Valor R por defecto: 3.

Example

>>> import numpy as np
>>> from circust.fitting.fmm import FMMModel
>>>
>>> t    = np.linspace(0, 2*np.pi, 50, endpoint=False)
>>> data = 0.3 + 0.6 * np.cos(t - 1.2) + np.random.normal(0, 0.05, 50)
>>> result = FMMModel().fit(data, t)
>>> print(result.summary())
static peak_time(alpha, beta, omega)[source]

Tiempo de pico del modelo FMM — formula compUU de R.

t_U = alpha + 2*atan2((1/omega)*sin(-beta/2), cos(-beta/2)) mod 2*pi

Return type:

float

Parameters:

  • alpha (float)

  • beta (float)

  • omega (float)

fit(data, time_points)[source]

Ajusta el modelo FMM a data.

Parameters:

  • data (ndarray) – Valores de expresion observados (normalizados a [-1, 1]).

  • time_points (ndarray) – Eje temporal circular en [0, 2*pi), tipicamente CPCAResult.circular_scale.

Returns:

Claves de params: M, A, alpha, beta, omega

Return type:

FitResult

Cosinor

circust/fitting/cosinor.py

Modelo ritmico Cosinor — ajuste sinusoidal por minimos cuadrados.

Modelo matematico

y(t) = M + A * cos(t - phi)

donde

M mesor — nivel medio de la senal A amplitud — mitad del rango pico-a-valle phi acrofase — tiempo de pico en [0, 2*pi)

Esta es la formulacion estandar en la literatura circadiana: phi es directamente el tiempo de pico (t_peak = phi), lo que hace la interpretacion de los parametros inmediata.

El modelo se ajusta por OLS reescribiendolo como regresion lineal:

y = M + b*cos(t) + g*sin(t)

= M + A*cos(t - phi)

con b = A*cos(phi), g = A*sin(phi), de modo que

A = sqrt(b**2 + g**2) phi = atan2(g, b) mod 2*pi

Nota: la formulacion anterior de CIRCUST usaba cos(t + phi), donde el pico era -phi mod 2*pi. Esta formulacion es equivalente matematicamente (mismos valores ajustados, mismo R2) pero phi tiene signo opuesto.

class circust.fitting.cosinor.CosinorModel[source]

Bases: RhythmModel

Modelo Cosinor de un solo componente.

Ajusta y(t) = M + A*cos(t - phi) por minimos cuadrados ordinarios.

Es el mas rapido de los modelos ritmicos — O(n) — y se usa tanto como test de ritmicidad independiente como estimador de parametros iniciales dentro de FMM.

Example

>>> import numpy as np
>>> from circust.fitting.cosinor import CosinorModel
>>>
>>> t    = np.linspace(0, 2*np.pi, 50, endpoint=False)
>>> data = 0.3 + 0.6 * np.cos(t - 1.2) + np.random.normal(0, 0.05, 50)
>>> result = CosinorModel().fit(data, t)
>>> print(result.summary())
static peak_time(phi)[source]

Tiempo de pico del modelo Cosinor.

Return type:

float

Parameters:

phi (float)

El maximo de M + A*cos(t - phi) ocurre cuando t - phi = 0:

t_peak = phi

fit(data, time_points)[source]

Ajusta el modelo Cosinor a data.

Parameters:

  • data (ndarray) – Valores de expresion observados (normalizados a [-1, 1]).

  • time_points (ndarray) – Eje temporal circular en [0, 2*pi).

Returns:

Claves de params: M, A, phi phi es directamente el tiempo de pico en [0, 2*pi).

Return type:

FitResult

ORI (Isotonic Regression Circular)

circust/fitting/ori.py

Modelo ritmico ORI — regresion isotonica circular unimodal (no parametrica).

Referencia

Larriba, Rueda & Fernandez (2016). «Identification of rhythmic genes in complex biological datasets and application to the study of human circadian system gene expression.» Statistical Applications in Genetics and Molecular Biology, 15(3), 205-223.

Modelo matematico

Dado un perfil de expresion genica ordenado por tiempo circular, el modelo ORI busca el mejor ajuste de forma subida-bajada (un pico, un valle) mediante el Algoritmo Pool-Adjacent-Violators (PAVA).

Para cada par candidato (L=valle, U=pico) de extremos locales se ajusta:

  • Un PAVA creciente de L a U (circularmente).

  • Un PAVA decreciente de U a L (circularmente).

Se selecciona el par (L_opt, U_opt) que minimiza el MSE del ajuste unimodal resultante. La puntuacion de ritmicidad es:

R2_ORI = 1 - MSE_unimodal / MSE_plano

donde MSE_plano = Var(datos) es el error del modelo nulo (media constante).

Algoritmo PAVA (Pool-Adjacent-Violators)

Complejidad: O(n) tiempo, O(n) espacio.

Dado un vector de observaciones y con pesos w, produce el vector y_hat que minimiza sum w_i (y_i - y_hat_i)^2 sujeto a restricciones de monotonia (creciente o decreciente).

Aceleracion

Cuando Numba esta disponible, las funciones criticas (PAVA, busqueda de pares L-U) se compilan JIT a codigo maquina, obteniendo ~20-50x de aceleracion dado que el bucle (L, U) x PAVA es O(n^3) por gen.

circust.fitting.ori.pava_increasing(y, w=None)[source]

Regresion isotonica creciente mediante PAVA.

Complejidad: O(n) tiempo, O(n) espacio.

Parameters:

  • y (ndarray) – Valores observados.

  • w (ndarray | None) – Pesos positivos. Si es None, pesos uniformes.

Return type:

ndarray

circust.fitting.ori.pava_decreasing(y, w=None)[source]

Regresion isotonica decreciente.

Equivalente a negar y, ejecutar pava_increasing, y negar de vuelta.

Return type:

ndarray

Parameters:

  • y (ndarray)

  • w (ndarray | None)

circust.fitting.ori.find_local_extrema(v)[source]

Encuentra minimos locales (valles) y maximos locales (picos) usando vecinos circulares.

Un punto i es un maximo local si v[i-1] <= v[i] Y v[i] >= v[i+1], con indices envolventes circulares. Minimos definidos simetricamente.

Parameters:

v (ndarray)

Return type:

tuple[ndarray, ndarray]

Returns:

  • candL (array int — indices de minimos locales (valles).)

  • candU (array int — indices de maximos locales (picos).)

circust.fitting.ori.circular_unimodal_fit(v)[source]

Ajusta una regresion isotonica circular unimodal (subida-bajada).

Busca sobre todos los pares candidatos (valle, pico) — minimos y maximos locales de v — para encontrar el par (L, U) que produce el menor MSE al ajustar un PAVA creciente de L->U y un PAVA decreciente de U->L (circularmente).

Parameters:

v (ndarray) – Valores de expresion genica ordenados por tiempo circular.

Return type:

tuple[ndarray, float, int, int] | None

Returns:

  • fitted (array (n,) — valores ajustados.)

  • mse (float — error cuadratico medio.)

  • L_opt (int — indice del valle optimo (base 0).)

  • U_opt (int — indice del pico optimo (base 0).)

  • None si no se encuentra un ajuste valido.

class circust.fitting.ori.ORIModel[source]

Bases: RhythmModel

Modelo de ritmicidad no parametrica basado en regresion isotonica circular unimodal (ORI / ORIOS).

A diferencia de Cosinor y FMM, este modelo no asume una forma de onda parametrica. Busca la mejor forma subida-bajada (unimodal) compatible con restricciones de monotonia impuestas por PAVA.

Los “parametros” del modelo son los indices del valle (L_opt) y pico (U_opt) optimos; el R2 compara contra el modelo plano.

Example

>>> import numpy as np
>>> from circust.fitting.ori import ORIModel
>>>
>>> t    = np.linspace(0, 2*np.pi, 50, endpoint=False)
>>> data = 0.3 + 0.6 * np.cos(t - 1.2) + np.random.normal(0, 0.05, 50)
>>> result = ORIModel().fit(data, t)
>>> print(result.r2)
>>> print(result.summary())
static peak_time(U_opt, time_points)[source]

Tiempo de pico del modelo ORI.

El pico se localiza en la posicion U_opt del vector ordenado circularmente. El tiempo de pico es time_points[U_opt].

Return type:

float

Parameters:

  • U_opt (int)

  • time_points (ndarray)

fit(data, time_points)[source]

Ajusta el modelo ORI a data.

Parameters:

  • data (ndarray) – Valores de expresion ordenados por el tiempo circular preliminar.

  • time_points (ndarray) – Eje temporal circular en [0, 2*pi), producido por CPCA.

Returns:

Claves de params: L_opt, U_opt.

  • L_opt : indice del valle optimo (base 0)

  • U_opt : indice del pico optimo (base 0)

peak_time = time_points[U_opt]. r2 = R2_ORI = 1 - MSE_unimodal / MSE_plano.

Return type:

FitResult

Clase base RhythmModel

Interfaz compartida para todos los modelos de ajuste ritmico del pipeline CIRCUST.

Nota de diseno

FitResult es un contenedor de datos puro: almacena todo lo que fit() produce, incluido el tiempo de pico precalculado por el modelo. Asi el resultado es autosuficiente — nadie necesita llamar al modelo despues de obtenerlo.

La logica especifica de cada modelo vive exclusivamente en RhythmModel y sus subclases:

fit() — abstracto, produce un FitResult con peak_time ya calculado.

Para el calculo del pico con floats sueltos (reparametrizacion manual) cada modelo expone un metodo estatico auxiliar, p.ej. FMMModel.peak_time().

class circust.fitting.rhythm_model.FitResult(fitted, params, peak_time, r2, residuals, residuals_std, sse, model_name)[source]

Bases: object

Resultado de cualquier modelo de ajuste ritmico.

Variables:

  • fitted (np.ndarray, dim (n_muestras,)) – Valores predichos por el modelo en cada punto temporal.

  • params (dict) –

    Parametros estimados del modelo.

    Claves CosinorM, A, phi

    • M : mesor (nivel medio)

    • A : amplitud

    • phi : acrofase en [0, 2*pi)

    Claves FMMM, A, alpha, beta, omega

    • M : mesor

    • A : amplitud

    • alpha : parametro de localizacion en [0, 2*pi)

    • beta : parametro de forma en [0, 2*pi)

    • omega : parametro de asimetria en (0, 1]

  • peak_time (float) – Tiempo de pico en [0, 2*pi), precalculado por el modelo al hacer fit(). Cosinor: phi mod 2*pi (phi es directamente el tiempo de pico). FMM: formula compUU — alpha + 2*atan2(…) mod 2*pi.

  • r2 (float) – Coeficiente de determinacion R-cuadrado = 1 - SS_res / SS_tot.

  • residuals (np.ndarray, forma (n_muestras,)) – Residuos brutos: datos - ajustado.

  • residuals_std (np.ndarray, forma (n_muestras,)) – Residuos estandarizados: (residuos - media) / desv_std.

  • sse (float) – Suma de errores al cuadrado.

  • model_name (str) – "cosinor" o "fmm".

Parameters:

  • fitted (ndarray)

  • params (dict)

  • peak_time (float)

  • r2 (float)

  • residuals (ndarray)

  • residuals_std (ndarray)

  • sse (float)

  • model_name (str)

fitted: ndarray
params: dict
peak_time: float
r2: float
residuals: ndarray
residuals_std: ndarray
sse: float
model_name: str
summary()[source]
Return type:

str

class circust.fitting.rhythm_model.RhythmModel[source]

Bases: ABC

Interfaz comun para todos los modelos de ajuste de ritmos circadianos.

Las subclases deben implementar fit().

Los parametros pasados al constructor son hiperparametros que permanecen constantes entre llamadas a fit() (p.ej. tamanos de rejilla para FMM).

abstractmethod fit(data, time_points)[source]

Ajusta el modelo a data observados en time_points.

El FitResult devuelto incluye peak_time ya calculado.

Parameters:

  • data (ndarray) – Valores de expresion normalizados (tipicamente en [-1, 1]).

  • time_points (ndarray) – Eje temporal circular en [0, 2*pi), producido por CPCA.

Return type:

FitResult