# SPDX-License-Identifier: AGPL-3.0-or-later
# Copyright (C) 2025-2026 Djura | Risk - Data - Engineering S.r.l.
from numpy import pi, append
from functools import cached_property
from pydantic import BaseModel
from typing import List, Optional, Dict, Union
[docs]
class SPOModel(BaseModel):
displacement: Union[List[float], List[List[float]]]
force: List[float]
[docs]
class DispForceModel(BaseModel):
displacement: Optional[List[float]] = None
force: Optional[List[float]] = None
mstar: Optional[float] = None
gamma: Optional[float] = None
sdof: Optional[bool] = True
spo: Optional[Dict] = SPOModel
damping: Optional[float] = 0.05
residual: Optional[float] = 0.0
[docs]
class Backbone:
def __init__(
self, data, method: str = "backbone", model: str = "bilin"
) -> None:
self.model = model
self.method = method
self.data = data
@cached_property
def backbone(self):
if self.model.lower() == "isol":
return self.data
if self.method == "backbone":
return self.data
elif self.method == "idealized":
return self._estimate_backbone(self.data)
elif self.method == "spo":
return self._read_from_spo(self.data)
else:
raise ValueError(
f"Unsupported method '{self.method}'. Expected 'backbone',"
"'spo', or 'idealized")
def _estimate_backbone_bilin(self, data: DispForceModel):
# for SDOF: force = cy*g
# this is not exactly Force, which would be cy*g*mstar
sdof = data.sdof
force = data.force
disp = data.displacement
if force is None or disp is None:
return None
if not sdof:
# Perform MDOF to ESDOF conversion first
mstar = data.mstar
if mstar is None:
return None
# this is without inclusion of Gamma
# for truly converting MDOF to ESDOF, the forces and displacements
# need to be divided by Gamma too, whith is the
# first mode participation factor
force[0] /= mstar
force[1] /= mstar
dy = disp[0]
fy = force[0]
ductility = disp[1] / dy
# T*, period of ESDOF
period = 2 * pi * (dy / fy) ** 0.5
# Hardening ratio of Ductility vs Strength Ratio curve
ah = (force[1] / force[0] - 1) / (ductility - 1)
return {
"damping": data.damping,
"ductility": ductility,
"hardening_ratio": ah,
"period": period,
"ductility_f": disp[2] / disp[0],
}
def _estimate_backbone_infill(self, data: DispForceModel):
sdof = data.sdof
force = data.force
disp = data.displacement
if force is None or disp is None:
return None
if not sdof:
# Perform MDOF to ESDOF conversion first
mstar = data.mstar
gamma = data.gamma
if mstar is None or gamma is None:
return None
# this is without inclusion of Gamma
# for truly converting MDOF to ESDOF, the forces and displacements
# need to be divided by Gamma too, whith is the
# first mode participation factor
force[0] /= mstar
force[-2] /= mstar
dy = disp[0]
fy = force[0]
# T*, period of ESDOF
period = 2 * pi * (dy / fy) ** 0.5
if not sdof and gamma is not None:
fy /= gamma
force[-2] /= gamma
return {
"period": period,
"c_y": fy,
"c_rp": force[-2],
"mu_h": disp[2] / dy,
"mu_s": disp[-3] / dy,
"mu_rp": disp[-2] / dy,
"mu_ult": disp[-1] / dy,
}
def _estimate_backbone(self, data):
"""
Estimate backbone parameters from idealized curve
"""
data = DispForceModel.model_validate(data)
if self.model.lower() == "bilin":
return self._estimate_backbone_bilin(data)
elif self.model.lower() == "infill":
return self._estimate_backbone_infill(data)
return None
def _read_from_spo(self, data: DispForceModel):
"""
Fit an idealized curve to the static pushover curve
then estimate the backbone parameters
"""
from numpy import asarray
data = DispForceModel.model_validate(data)
mstar = data.mstar
gamma = data.gamma
spo = data.spo
sdof = data.sdof
spo = SPOModel.model_validate(spo)
residual = data.residual
if (mstar is None or spo is None) and not sdof:
return None
if sdof:
mstar = 1
gamma = 1
v = asarray(spo.force)
d = asarray(spo.displacement)
if d.ndim > 1:
# Base shear and all floor displacements
d = d[:, -1]
# else, Base shear and top displacement
v[v < 0] = 0.0
if d[0] != 0:
d = append(0, d)
v = append(0, v)
if self.model.lower() == "bilin":
out = self._fit_spo_bilin(mstar, d, v, residual)
elif self.model.lower() == "infill":
if gamma is None:
return None
out = self._fit_spo_infill(mstar, gamma, d, v, residual)
else:
return None
out["damping"] = data.damping
return out
def _fit_spo_bilin(self, mstar, d, v, residual):
disp, force = self._piecewise_fit(v, d, 4)
if force is None or disp is None:
return None
# Count indices backwards as last points are more stable
f_i = -1
s_i = -2
y_i = -3
v_res = residual * force[y_i]
d_res = disp[s_i] + (disp[f_i] - disp[s_i]) * \
force[s_i] / (force[s_i] - force[f_i])
disp[f_i] = d_res
force[f_i] = v_res
force[y_i] /= mstar
force[s_i] /= mstar
dy = disp[y_i]
fy = force[y_i]
ductility = disp[s_i] / dy
# T*, period of ESDOF
period = 2 * pi * (dy / fy) ** 0.5
# Hardening ratio of Ductility vs Strength Ratio curve
ah = (force[s_i] / force[y_i] - 1) / (ductility - 1)
return {
"ductility": ductility,
"hardening_ratio": ah,
"period": period,
"ductility_f": disp[f_i] / disp[y_i],
"dy_sdof_wo_gamma": dy,
"fy_sdof_wo_gamma": fy,
}
def _fit_spo_infill(self, mstar, gamma, d, v, residual):
from numpy import argmax
kel = v[1] / d[1]
# First point
idx = argmax(v)
force = [v[idx], v[idx]]
disp = [force[0] / kel]
d_fit, _ = self._piecewise_fit(v[idx + 1:], d[idx + 1:], 4)
disp.append(d_fit[1])
idx1 = argmax(d > d_fit[1])
d_fit, f_fit = self._piecewise_fit(v[idx1 + 1:], d[idx1 + 1:], 3)
disp = disp + d_fit
force = force + f_fit
# This is a better fitting function, but it does not fit
# Nafeh and O'Reilly model (2022)
# disp, force = self._piecewise_fit(v, d, 5)
if force is None or disp is None:
return None
v_res = residual * force[0]
force[-1] = v_res
dy = disp[0]
fy = force[0] / mstar
# T*, period of ESDOF
period = 2 * pi * (dy / fy) ** 0.5
# Idealised displacement vs shear curve
_d = [0, dy, disp[2], disp[-3], disp[-2], disp[-1]]
_v = [0, fy, fy, force[-2] / gamma /
mstar, force[-2] / gamma / mstar, 0]
return {
"period": period,
"c_y": force[0] / gamma / mstar,
"c_rp": force[-2] / gamma / mstar,
"mu_h": disp[2] / dy,
"mu_s": disp[-3] / dy,
"mu_rp": disp[-2] / dy,
"mu_ult": disp[-1] / dy,
"d": _d,
"v": _v,
"dy_sdof_wo_gamma": dy,
"fy_sdof_wo_gamma": fy,
}
def _piecewise_fit(self, v, d, npoints):
from pwlf import PiecewiseLinFit
pwlf_model = PiecewiseLinFit(d, v)
breaks = pwlf_model.fit(npoints, x_c=[0], y_c=[0])
breaks[-1] = d[-1]
d_ideal = []
v_ideal = []
for x_break in breaks:
y_break = pwlf_model.predict([x_break])[0]
if y_break < 0:
y_break = 0
d_ideal.append(x_break)
v_ideal.append(y_break)
return d_ideal, v_ideal
[docs]
def get_infill_spo(data):
c_y = data["c_y"]
c_rp = data["c_rp"]
mu_h = data["mu_h"]
mu_s = data["mu_s"]
mu_rp = data["mu_rp"]
mu_ult = data["mu_ult"]
r = [0, 1, 1, c_rp / c_y, c_rp / c_y, 0]
m = [0, 1, mu_h, mu_s, mu_rp, mu_ult]
return r, m