# ReliaPy-Workbench: Python/Gradio Reliability Analysis Tool for Colab
# Copy this entire cell after installing gradio.
import io
import math
import warnings
import tempfile
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import optimize, stats
import gradio as gr
warnings.filterwarnings("ignore")
EPS = 1e-12
def _path(file_obj):
if file_obj is None:
return None
if isinstance(file_obj, str):
return file_obj
return getattr(file_obj, "name", None) or getattr(file_obj, "path", None)
def read_table(file_obj):
path = _path(file_obj)
if path is None:
raise ValueError("Please upload a CSV file.")
if path.lower().endswith((".xlsx", ".xls")):
return pd.read_excel(path)
return pd.read_csv(path)
def _clean_name(s):
return str(s).strip()
def infer_col(df, requested, candidates, fallback_index=0):
if requested and requested != "Auto":
if requested not in df.columns:
raise ValueError(f"Column '{requested}' not found. Available columns: {list(df.columns)}")
return requested
lower = {str(c).lower().strip(): c for c in df.columns}
for cand in candidates:
if cand.lower() in lower:
return lower[cand.lower()]
# fuzzy contains
for c in df.columns:
lc = str(c).lower()
if any(k in lc for k in candidates):
return c
if len(df.columns) > fallback_index:
return df.columns[fallback_index]
raise ValueError("Could not infer a required column.")
def as_event(series):
s = series.copy()
if s.dtype == bool:
return s.astype(int).values
if pd.api.types.is_numeric_dtype(s):
return (pd.to_numeric(s, errors="coerce").fillna(0) > 0).astype(int).values
ss = s.astype(str).str.lower().str.strip()
return ss.isin(["1", "true", "t", "yes", "y", "event", "failed", "failure", "fail"]).astype(int).values
def benard_ranks(n):
i = np.arange(1, n + 1)
return (i - 0.3) / (n + 0.4)
def life_rank_regression(times, dist="Weibull"):
t = np.asarray(times, dtype=float)
t = np.sort(t[t > 0])
n = len(t)
if n < 2:
raise ValueError("Rank regression needs at least two uncensored failure times.")
F = benard_ranks(n)
x = np.log(t)
if dist == "Weibull":
y = np.log(-np.log(1 - F))
slope, intercept, r, p, se = stats.linregress(x, y)
beta = max(slope, EPS)
eta = float(np.exp(-intercept / beta))
ll = weibull_loglik(t, np.ones_like(t), beta, eta)
return {"dist": dist, "method": "Rank regression", "beta_shape": beta, "eta_scale": eta,
"loglik_uncensored": ll, "rr_r2": r*r}
else:
y = stats.norm.ppf(F)
slope, intercept, r, p, se = stats.linregress(x, y)
sigma = max(1.0 / slope, EPS)
mu = -intercept / slope
ll = lognormal_loglik(t, np.ones_like(t), mu, sigma)
return {"dist": dist, "method": "Rank regression", "mu_log": mu, "sigma_log": sigma,
"median_life": float(np.exp(mu)), "loglik_uncensored": ll, "rr_r2": r*r}
def weibull_loglik(t, event, beta, eta):
t = np.asarray(t, dtype=float)
e = np.asarray(event, dtype=int)
z = (t / eta) ** beta
logpdf = np.log(beta) - beta * np.log(eta) + (beta - 1) * np.log(t) - z
logsf = -z
return float(np.sum(e * logpdf + (1 - e) * logsf))
def lognormal_loglik(t, event, mu, sigma):
t = np.asarray(t, dtype=float)
e = np.asarray(event, dtype=int)
z = (np.log(t) - mu) / sigma
logpdf = -np.log(t) - np.log(sigma) - 0.5*np.log(2*np.pi) - 0.5*z*z
logsf = stats.norm.logsf(z)
return float(np.sum(e * logpdf + (1 - e) * logsf))
def life_mle(times, event, dist="Weibull"):
t = np.asarray(times, dtype=float)
e = np.asarray(event, dtype=int)
if len(t) < 2 or e.sum() < 1:
raise ValueError("MLE needs at least two observations and at least one failure event.")
fail = t[e == 1]
rr = life_rank_regression(fail, dist) if len(fail) >= 2 else None
if dist == "Weibull":
if rr:
init = [np.log(rr["beta_shape"]), np.log(rr["eta_scale"])]
else:
init = [0.0, np.log(np.median(t))]
def nll(x):
beta, eta = np.exp(x[0]), np.exp(x[1])
return -weibull_loglik(t, e, beta, eta)
res = optimize.minimize(nll, init, method="Nelder-Mead", options={"maxiter": 20000})
beta, eta = np.exp(res.x[0]), np.exp(res.x[1])
ll = -res.fun
return {"dist": dist, "method": "MLE", "beta_shape": beta, "eta_scale": eta,
"mean_life": float(eta * math.gamma(1 + 1/beta)),
"median_life": float(eta * (np.log(2)) ** (1/beta)),
"loglik": ll, "AIC": 2*2 - 2*ll, "BIC": 2*np.log(len(t)) - 2*ll,
"events": int(e.sum()), "observations": len(t)}
else:
if rr:
init = [rr["mu_log"], np.log(rr["sigma_log"])]
else:
init = [np.log(np.median(fail)), 0.0]
def nll(x):
mu, sigma = x[0], np.exp(x[1])
return -lognormal_loglik(t, e, mu, sigma)
res = optimize.minimize(nll, init, method="Nelder-Mead", options={"maxiter": 20000})
mu, sigma = res.x[0], np.exp(res.x[1])
ll = -res.fun
return {"dist": dist, "method": "MLE", "mu_log": mu, "sigma_log": sigma,
"median_life": float(np.exp(mu)), "mean_life": float(np.exp(mu + 0.5*sigma*sigma)),
"loglik": ll, "AIC": 2*2 - 2*ll, "BIC": 2*np.log(len(t)) - 2*ll,
"events": int(e.sum()), "observations": len(t)}
def reliability_at_time(params, mission_time):
if mission_time is None or mission_time <= 0:
return np.nan
if params["dist"] == "Weibull":
beta, eta = params["beta_shape"], params["eta_scale"]
return float(np.exp(- (mission_time / eta) ** beta))
else:
mu, sigma = params["mu_log"], params["sigma_log"]
return float(stats.norm.sf((np.log(mission_time) - mu) / sigma))
def life_probability_plot(t, e, params):
fig, ax = plt.subplots(figsize=(7, 5))
fail = np.sort(np.asarray(t)[np.asarray(e) == 1])
if len(fail) < 2:
ax.text(0.05, 0.5, "At least two failures are needed for a probability plot.", transform=ax.transAxes)
return fig
F = benard_ranks(len(fail))
x = np.log(fail)
if params["dist"] == "Weibull":
y = np.log(-np.log(1 - F))
ax.scatter(x, y, label="Median-rank failures")
xx = np.linspace(x.min()*0.95, x.max()*1.05, 100)
beta, eta = params["beta_shape"], params["eta_scale"]
yy = beta * (xx - np.log(eta))
ax.plot(xx, yy, label="Fitted Weibull line")
ax.set_ylabel("ln[-ln(1-F)]")
ax.set_title("Weibull Probability Plot")
else:
y = stats.norm.ppf(F)
ax.scatter(x, y, label="Median-rank failures")
xx = np.linspace(x.min()*0.95, x.max()*1.05, 100)
mu, sigma = params["mu_log"], params["sigma_log"]
yy = (xx - mu) / sigma
ax.plot(xx, yy, label="Fitted Lognormal line")
ax.set_ylabel("Normal quantile")
ax.set_title("Lognormal Probability Plot")
ax.set_xlabel("ln(time)")
ax.grid(True, alpha=0.3)
ax.legend()
fig.tight_layout()
return fig
def life_cdf_plot(t, e, params):
fig, ax = plt.subplots(figsize=(7, 5))
fail = np.sort(np.asarray(t)[np.asarray(e) == 1])
F = benard_ranks(len(fail)) if len(fail) else np.array([])
if len(fail):
ax.scatter(fail, F, label="Median-rank empirical CDF")
x = np.linspace(max(min(np.asarray(t)) * 0.1, EPS), max(np.asarray(t)) * 1.15, 200)
if params["dist"] == "Weibull":
beta, eta = params["beta_shape"], params["eta_scale"]
cdf = 1 - np.exp(-(x/eta)**beta)
else:
mu, sigma = params["mu_log"], params["sigma_log"]
cdf = stats.norm.cdf((np.log(x) - mu)/sigma)
ax.plot(x, cdf, label="Fitted CDF")
ax.set_xlabel("Time")
ax.set_ylabel("F(t)")
ax.set_title("Fitted Life Distribution")
ax.grid(True, alpha=0.3)
ax.legend()
fig.tight_layout()
return fig
def life_contour_plot(t, e, params):
fig, ax = plt.subplots(figsize=(7, 5))
try:
if params["dist"] == "Weibull":
beta0, eta0 = params["beta_shape"], params["eta_scale"]
beta_grid = np.linspace(max(beta0*0.35, EPS), beta0*2.5, 70)
eta_grid = np.linspace(max(eta0*0.35, EPS), eta0*2.5, 70)
B, E = np.meshgrid(beta_grid, eta_grid)
LL = np.zeros_like(B)
for i in range(B.shape[0]):
for j in range(B.shape[1]):
LL[i, j] = weibull_loglik(t, e, B[i, j], E[i, j])
D = -2*(LL - np.max(LL))
cs = ax.contour(B, E, D, levels=[2.30, 6.18, 11.83])
ax.clabel(cs, inline=True, fontsize=8)
ax.scatter([beta0], [eta0], marker="x", s=80, label="Estimate")
ax.set_xlabel("Shape beta")
ax.set_ylabel("Scale eta")
ax.set_title("Weibull Likelihood Contours")
else:
mu0, sigma0 = params["mu_log"], params["sigma_log"]
mu_grid = np.linspace(mu0 - 2.0*sigma0, mu0 + 2.0*sigma0, 70)
sig_grid = np.linspace(max(sigma0*0.35, EPS), sigma0*2.5, 70)
M, S = np.meshgrid(mu_grid, sig_grid)
LL = np.zeros_like(M)
for i in range(M.shape[0]):
for j in range(M.shape[1]):
LL[i, j] = lognormal_loglik(t, e, M[i, j], S[i, j])
D = -2*(LL - np.max(LL))
cs = ax.contour(M, S, D, levels=[2.30, 6.18, 11.83])
ax.clabel(cs, inline=True, fontsize=8)
ax.scatter([mu0], [sigma0], marker="x", s=80, label="Estimate")
ax.set_xlabel("mu")
ax.set_ylabel("sigma")
ax.set_title("Lognormal Likelihood Contours")
ax.grid(True, alpha=0.3)
ax.legend()
except Exception as exc:
ax.text(0.05, 0.5, f"Contour plot could not be computed:\n{exc}", transform=ax.transAxes)
fig.tight_layout()
return fig
def fit_life_data(file_obj, dist, method, time_col, event_col, mission_time):
try:
df = read_table(file_obj)
tc = infer_col(df, time_col, ["time", "ttf", "failure_time", "life", "hours", "cycles"], 0)
raw_t = pd.to_numeric(df[tc], errors="coerce")
if event_col and event_col != "Auto" and event_col != "None":
ec = infer_col(df, event_col, ["event", "status", "failed", "failure", "censor", "censored"], 1)
event = as_event(df[ec])
else:
# Auto uses event/status column if present; otherwise all failures
possible = [c for c in df.columns if str(c).lower().strip() in ["event", "status", "failed", "failure"]]
event = as_event(df[possible[0]]) if possible else np.ones(len(df), dtype=int)
mask = np.isfinite(raw_t) & (raw_t > 0) & np.isfinite(event)
t = raw_t[mask].values.astype(float)
e = np.asarray(event)[mask].astype(int)
if method == "Rank Regression":
params = life_rank_regression(t[e == 1], dist)
params["observations"] = len(t)
params["events_used"] = int(e.sum())
params["note"] = "Rank regression uses uncensored failures; right-censored rows are not used in the line fit."
else:
params = life_mle(t, e, dist)
params[f"Reliability_at_t={mission_time}"] = reliability_at_time(params, mission_time)
summary = pd.DataFrame({"Metric": list(params.keys()), "Value": [params[k] for k in params.keys()]})
fig1 = life_probability_plot(t, e, params)
fig2 = life_contour_plot(t, e, params)
fig3 = life_cdf_plot(t, e, params)
return summary, fig1, fig2, fig3, f"Used time column: {tc}. Observations after cleaning: {len(t)}."
except Exception as exc:
return pd.DataFrame({"Error": [str(exc)]}), None, None, None, "Analysis failed."
def plp_fit(times, T=None):
t = np.asarray(times, dtype=float)
t = t[np.isfinite(t) & (t > 0)]
if T is None:
T = np.max(t)
T = float(T)
if len(t) < 2 or T <= 0:
raise ValueError("PLP/Crow-AMSAA needs at least two positive event times.")
denom = np.sum(np.log(T / np.clip(t, EPS, None)))
beta = len(t) / max(denom, EPS)
lam = len(t) / (T ** beta)
ll = float(np.sum(np.log(lam*beta) + (beta - 1)*np.log(np.clip(t, EPS, None))) - lam*(T**beta))
return {"lambda": float(lam), "beta": float(beta), "T": T, "n_events": len(t), "loglik": ll, "AIC": 2*2 - 2*ll}
def loglinear_nhpp_fit(times, T=None, n_systems=1):
t = np.asarray(times, dtype=float)
t = t[np.isfinite(t) & (t >= 0)]
if T is None:
T = np.max(t)
T = float(T)
n_systems = max(int(n_systems), 1)
if len(t) < 2:
raise ValueError("Log-linear NHPP needs at least two event times.")
init_rate = len(t) / max(n_systems*T, EPS)
init = [np.log(max(init_rate, EPS)), 0.0]
def cum_int(a, b):
if abs(b) < 1e-8:
return n_systems * np.exp(a) * T
return n_systems * np.exp(a) * (np.exp(b*T) - 1.0) / b
def nll(x):
a, b = x
val = np.sum(a + b*t) - cum_int(a, b)
if not np.isfinite(val):
return 1e100
return -val
res = optimize.minimize(nll, init, method="Nelder-Mead", options={"maxiter": 20000})
a, b = res.x
ll = -res.fun
return {"a_log_rate": float(a), "b_time_slope": float(b), "T": T, "n_events": len(t),
"loglik": float(ll), "AIC": 2*2 - 2*ll}
def parse_breakpoints(text):
if not text or str(text).strip() == "":
return []
vals = []
for x in str(text).replace(";", ",").split(","):
x = x.strip()
if x:
vals.append(float(x))
return sorted(vals)
def segment_plp(times, breakpoints=None, auto=False):
t = np.sort(np.asarray(times, dtype=float))
t = t[np.isfinite(t) & (t > 0)]
T = float(np.max(t))
if auto:
candidates = np.unique(t)
candidates = candidates[(candidates > np.quantile(t, 0.2)) & (candidates < np.quantile(t, 0.8))]
best = None
for bp in candidates:
left = t[t <= bp]
right = t[t > bp] - bp
try:
if len(left) < 3 or len(right) < 3:
continue
f1 = plp_fit(left, T=bp)
f2 = plp_fit(right, T=T-bp)
aic = f1["AIC"] + f2["AIC"]
if best is None or aic < best[0]:
best = (aic, bp, f1, f2)
except Exception:
pass
if best is not None:
breakpoints = [best[1]]
else:
breakpoints = []
if breakpoints is None:
breakpoints = []
cuts = [0.0] + [bp for bp in breakpoints if 0 < bp < T] + [T]
rows = []
for i in range(len(cuts)-1):
start, end = cuts[i], cuts[i+1]
seg_abs = t[(t > start) & (t <= end)]
local = seg_abs - start
if len(local) >= 2:
fit = plp_fit(np.clip(local, EPS, None), T=end-start)
rows.append({"segment": i+1, "start": start, "end": end, **fit})
else:
rows.append({"segment": i+1, "start": start, "end": end, "n_events": len(local),
"lambda": np.nan, "beta": np.nan, "loglik": np.nan, "AIC": np.nan})
return pd.DataFrame(rows), cuts
def growth_plots(times, summary_segments, cuts):
t = np.sort(np.asarray(times, dtype=float))
n = np.arange(1, len(t)+1)
T = max(t)
fig1, ax = plt.subplots(figsize=(7, 5))
ax.step(t, n, where="post", label="Observed cumulative failures")
xx_full = []
yy_full = []
cum_prev = 0.0
for _, row in summary_segments.iterrows():
start, end = row["start"], row["end"]
if np.isfinite(row.get("lambda", np.nan)) and np.isfinite(row.get("beta", np.nan)):
xs = np.linspace(start, end, 80)
local = np.clip(xs - start, 0, None)
ys = cum_prev + row["lambda"] * (local ** row["beta"])
xx_full.extend(xs.tolist())
yy_full.extend(ys.tolist())
cum_prev += row.get("n_events", 0)
if xx_full:
ax.plot(xx_full, yy_full, label="Fitted NHPP mean")
for bp in cuts[1:-1]:
ax.axvline(bp, linestyle="--", alpha=0.6)
ax.set_xlabel("Cumulative test time")
ax.set_ylabel("Cumulative failures")
ax.set_title("Reliability Growth / NHPP Plot")
ax.grid(True, alpha=0.3)
ax.legend()
fig1.tight_layout()
fig2, ax = plt.subplots(figsize=(7, 5))
cum_mtbf = t / n
ax.plot(t, cum_mtbf, marker="o", linewidth=1)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("Cumulative test time")
ax.set_ylabel("Cumulative MTBF = time / failures")
ax.set_title("Duane Plot")
ax.grid(True, alpha=0.3, which="both")
fig2.tight_layout()
return fig1, fig2
def fit_growth(file_obj, model_type, time_col, breakpoints_text):
try:
df = read_table(file_obj)
tc = infer_col(df, time_col, ["time", "event_time", "failure_time", "test_time", "hours", "cycles"], 0)
t = pd.to_numeric(df[tc], errors="coerce").dropna().values.astype(float)
t = np.sort(t[t > 0])
if len(t) < 2:
raise ValueError("At least two positive cumulative event times are required.")
if model_type == "Crow-AMSAA":
seg_df, cuts = segment_plp(t, breakpoints=[])
elif model_type == "Piecewise NHPP":
seg_df, cuts = segment_plp(t, breakpoints=parse_breakpoints(breakpoints_text), auto=False)
else:
seg_df, cuts = segment_plp(t, auto=True)
fig1, fig2 = growth_plots(t, seg_df, cuts)
interpretation = []
for _, r in seg_df.iterrows():
beta = r.get("beta", np.nan)
if np.isfinite(beta):
trend = "improving/decreasing event intensity" if beta < 1 else "deteriorating/increasing event intensity" if beta > 1 else "approximately constant intensity"
interpretation.append(f"Segment {int(r['segment'])}: beta={beta:.4g}, indicating {trend}.")
return seg_df, fig1, fig2, "\n".join(interpretation) + f"\nUsed time column: {tc}."
except Exception as exc:
return pd.DataFrame({"Error": [str(exc)]}), None, None, "Analysis failed."
def repair_plots(t, model_params, model_name, n_systems=1, sys_ids=None):
t = np.sort(np.asarray(t, dtype=float))
T = max(t)
fig1, ax = plt.subplots(figsize=(7, 5))
ax.step(t, np.arange(1, len(t)+1), where="post", label="Observed cumulative events")
xx = np.linspace(max(T*0.001, EPS), T, 200)
if model_name == "Power Law":
lam, beta = model_params["lambda"], model_params["beta"]
yy = n_systems * lam * (xx ** beta)
elif model_name == "Log-Linear":
a, b = model_params["a_log_rate"], model_params["b_time_slope"]
if abs(b) < 1e-8:
yy = n_systems * np.exp(a) * xx
else:
yy = n_systems * np.exp(a) * (np.exp(b*xx) - 1) / b
else:
yy = None
if yy is not None:
ax.plot(xx, yy, label="Fitted cumulative mean")
ax.set_xlabel("Time")
ax.set_ylabel("Cumulative events")
ax.set_title("Repairable-System Cumulative Events")
ax.grid(True, alpha=0.3)
ax.legend()
fig1.tight_layout()
fig2, ax = plt.subplots(figsize=(7, 5))
if model_name == "Power Law":
lam, beta = model_params["lambda"], model_params["beta"]
rate = lam * beta * (xx ** (beta - 1))
elif model_name == "Log-Linear":
a, b = model_params["a_log_rate"], model_params["b_time_slope"]
rate = np.exp(a + b*xx)
else:
# crude smoothed empirical rate for piecewise
bins = np.linspace(0, T, 12)
counts, edges = np.histogram(t, bins=bins)
centers = 0.5*(edges[1:]+edges[:-1])
widths = np.diff(edges)
ax.plot(centers, counts / np.maximum(widths*n_systems, EPS), marker="o")
rate = None
if rate is not None:
ax.plot(xx, rate)
ax.set_xlabel("Time")
ax.set_ylabel("Event rate per system")
ax.set_title("Estimated Event Rate")
ax.grid(True, alpha=0.3)
fig2.tight_layout()
fig3, ax = plt.subplots(figsize=(7, 5))
if sys_ids is None:
# one system MCF equals cumulative events
ax.step(t, np.arange(1, len(t)+1), where="post")
else:
tmp = pd.DataFrame({"time": t, "system": sys_ids}).sort_values("time")
counts = tmp.groupby("time").size().sort_index()
mcf = counts.cumsum() / max(n_systems, 1)
ax.step(mcf.index.values, mcf.values, where="post")
ax.set_xlabel("Time")
ax.set_ylabel("Mean cumulative function")
ax.set_title("Mean Cumulative Function (MCF)")
ax.grid(True, alpha=0.3)
fig3.tight_layout()
return fig1, fig2, fig3
def fit_repairable(file_obj, model_type, time_col, system_col, breakpoints_text):
try:
df = read_table(file_obj)
tc = infer_col(df, time_col, ["time", "event_time", "repair_time", "failure_time", "hours", "cycles"], 0)
time = pd.to_numeric(df[tc], errors="coerce")
mask = np.isfinite(time) & (time > 0)
t = time[mask].values.astype(float)
sys_ids = None
if system_col and system_col != "Auto" and system_col != "None":
sc = infer_col(df, system_col, ["system", "unit", "asset", "id"], 1)
sys_ids = df.loc[mask, sc].astype(str).values
n_systems = len(pd.unique(sys_ids))
else:
possible = [c for c in df.columns if any(k in str(c).lower() for k in ["system", "unit", "asset"])]
if possible:
sys_ids = df.loc[mask, possible[0]].astype(str).values
n_systems = len(pd.unique(sys_ids))
else:
n_systems = 1
if len(t) < 2:
raise ValueError("At least two repair/failure event times are required.")
if model_type == "Power Law":
fit = plp_fit(t, T=max(t))
summary = pd.DataFrame({"Metric": list(fit.keys()) + ["assumed_number_of_systems"],
"Value": list(fit.values()) + [n_systems]})
fig1, fig2, fig3 = repair_plots(t, fit, "Power Law", n_systems, sys_ids)
elif model_type == "Log-Linear":
fit = loglinear_nhpp_fit(t, T=max(t), n_systems=n_systems)
summary = pd.DataFrame({"Metric": list(fit.keys()) + ["assumed_number_of_systems"],
"Value": list(fit.values()) + [n_systems]})
fig1, fig2, fig3 = repair_plots(t, fit, "Log-Linear", n_systems, sys_ids)
else:
seg_df, cuts = segment_plp(t, parse_breakpoints(breakpoints_text), auto=False)
summary = seg_df
# plot observed plus piecewise fitted using growth plot
fig1, _ = growth_plots(t, seg_df, cuts)
# empirical rate and MCF
fig2, _, fig3 = repair_plots(t, {}, "Piecewise", n_systems, sys_ids)
note = f"Used time column: {tc}. Number of systems inferred: {n_systems}. If systems have different observation windows, add exposure handling before publication-grade inference."
return summary, fig1, fig2, fig3, note
except Exception as exc:
return pd.DataFrame({"Error": [str(exc)]}), None, None, None, "Analysis failed."
def alt_transform(stress, relationship):
s = np.asarray(stress, dtype=float)
if relationship == "Arrhenius; temperature in Celsius":
return 1.0 / (s + 273.15), "1 / absolute temperature (K^-1)"
if relationship == "Arrhenius; temperature in Kelvin":
return 1.0 / s, "1 / absolute temperature (K^-1)"
return np.log(s), "ln(stress)"
def alt_weibull_loglik(t, event, x, a, b, beta):
eta = np.exp(a + b*x)
z = (t / eta) ** beta
logpdf = np.log(beta) - beta*np.log(eta) + (beta-1)*np.log(t) - z
logsf = -z
return float(np.sum(event*logpdf + (1-event)*logsf))
def alt_lognormal_loglik(t, event, x, a, b, sigma):
mu = a + b*x
z = (np.log(t) - mu) / sigma
logpdf = -np.log(t) - np.log(sigma) - 0.5*np.log(2*np.pi) - 0.5*z*z
logsf = stats.norm.logsf(z)
return float(np.sum(event*logpdf + (1-event)*logsf))
def fit_alt_model(t, event, stress, dist, relationship):
x, xlab = alt_transform(stress, relationship)
# initialize by regression on log failure time
fail = event == 1
if fail.sum() >= 2 and len(np.unique(x[fail])) >= 2:
slope, intercept, *_ = stats.linregress(x[fail], np.log(t[fail]))
a0, b0 = intercept, slope
else:
a0, b0 = np.log(np.median(t)), 0.0
if dist == "Weibull":
init = [a0, b0, 0.0] # log beta
def nll(par):
a, b, logbeta = par
return -alt_weibull_loglik(t, event, x, a, b, np.exp(logbeta))
res = optimize.minimize(nll, init, method="Nelder-Mead", options={"maxiter": 30000})
a, b, logbeta = res.x
beta = np.exp(logbeta)
ll = -res.fun
out = {"distribution": "Weibull", "relationship": relationship, "a_intercept_log_eta": a,
"b_stress_slope": b, "beta_shape_common": beta,
"loglik": ll, "AIC": 2*3 - 2*ll, "BIC": 3*np.log(len(t)) - 2*ll}
else:
init = [a0, b0, 0.0] # log sigma
def nll(par):
a, b, logsig = par
return -alt_lognormal_loglik(t, event, x, a, b, np.exp(logsig))
res = optimize.minimize(nll, init, method="Nelder-Mead", options={"maxiter": 30000})
a, b, logsig = res.x
sigma = np.exp(logsig)
ll = -res.fun
out = {"distribution": "Lognormal", "relationship": relationship, "a_intercept_mu": a,
"b_stress_slope": b, "sigma_common": sigma,
"loglik": ll, "AIC": 2*3 - 2*ll, "BIC": 3*np.log(len(t)) - 2*ll}
return out, x, xlab
def alt_plots(t, event, stress, params, x, xlab):
fig1, ax = plt.subplots(figsize=(7, 5))
for s in sorted(pd.unique(stress)):
mask = (stress == s) & (event == 1)
tf = np.sort(t[mask])
if len(tf) >= 2:
F = benard_ranks(len(tf))
if params["distribution"] == "Weibull":
y = np.log(-np.log(1-F))
ax.scatter(np.log(tf), y, label=f"stress={s}")
else:
y = stats.norm.ppf(F)
ax.scatter(np.log(tf), y, label=f"stress={s}")
ax.set_xlabel("ln(time)")
ax.set_ylabel("Weibull/lognormal probability scale")
ax.set_title("ALT Probability Plot by Stress Level")
ax.grid(True, alpha=0.3)
ax.legend(fontsize=8)
fig1.tight_layout()
fig2, ax = plt.subplots(figsize=(7, 5))
raw_s_grid = np.linspace(np.min(stress), np.max(stress), 200)
xg, _ = alt_transform(raw_s_grid, params["relationship"])
if params["distribution"] == "Weibull":
life = np.exp(params["a_intercept_log_eta"] + params["b_stress_slope"]*xg)
ylab = "Characteristic life eta"
else:
life = np.exp(params["a_intercept_mu"] + params["b_stress_slope"]*xg)
ylab = "Median life"
# Observed stress-wise median failure times
med = pd.DataFrame({"time": t[event == 1], "stress": stress[event == 1]}).groupby("stress")["time"].median()
if len(med):
ax.scatter(med.index.values, med.values, label="Observed median failure time")
ax.plot(raw_s_grid, life, label="Fitted life-stress curve")
ax.set_xlabel("Stress")
ax.set_ylabel(ylab)
ax.set_yscale("log")
ax.set_title("Life-Stress Relationship")
ax.grid(True, alpha=0.3)
ax.legend()
fig2.tight_layout()
return fig1, fig2
def fit_alt(file_obj, dist, relationship, time_col, stress_col, event_col, use_stress):
try:
df = read_table(file_obj)
tc = infer_col(df, time_col, ["time", "ttf", "failure_time", "life", "hours", "cycles"], 0)
sc = infer_col(df, stress_col, ["stress", "temperature", "temp", "voltage", "load"], 1)
t = pd.to_numeric(df[tc], errors="coerce").values
stress = pd.to_numeric(df[sc], errors="coerce").values
if event_col and event_col != "Auto" and event_col != "None":
ec = infer_col(df, event_col, ["event", "status", "failed", "failure"], 2)
event = as_event(df[ec])
else:
possible = [c for c in df.columns if str(c).lower().strip() in ["event", "status", "failed", "failure"]]
event = as_event(df[possible[0]]) if possible else np.ones(len(df), dtype=int)
mask = np.isfinite(t) & np.isfinite(stress) & (t > 0) & (stress > 0) & np.isfinite(event)
t, stress, event = t[mask].astype(float), stress[mask].astype(float), np.asarray(event)[mask].astype(int)
if len(t) < 3 or event.sum() < 2:
raise ValueError("ALT fitting needs at least three rows and at least two failures.")
params, x, xlab = fit_alt_model(t, event, stress, dist, relationship)
if use_stress is not None and np.isfinite(use_stress) and use_stress > 0:
xu, _ = alt_transform(np.array([use_stress]), relationship)
if dist == "Weibull":
eta_use = float(np.exp(params["a_intercept_log_eta"] + params["b_stress_slope"]*xu[0]))
params[f"eta_at_use_stress_{use_stress}"] = eta_use
else:
med_use = float(np.exp(params["a_intercept_mu"] + params["b_stress_slope"]*xu[0]))
params[f"median_life_at_use_stress_{use_stress}"] = med_use
params["stress_transform"] = xlab
params["observations"] = len(t)
params["events"] = int(event.sum())
summary = pd.DataFrame({"Metric": list(params.keys()), "Value": [params[k] for k in params.keys()]})
fig1, fig2 = alt_plots(t, event, stress, params, x, xlab)
note = f"Used time column: {tc}; stress column: {sc}. Relationship transform: {xlab}."
return summary, fig1, fig2, note
except Exception as exc:
return pd.DataFrame({"Error": [str(exc)]}), None, None, "Analysis failed."
def sample_life_csv():
rng = np.random.default_rng(7)
beta, eta = 1.8, 500
t = eta * rng.weibull(beta, 80)
censor = rng.uniform(350, 900, 80)
obs = np.minimum(t, censor)
event = (t <= censor).astype(int)
return pd.DataFrame({"time": np.round(obs, 2), "event": event})
def sample_growth_csv():
rng = np.random.default_rng(9)
# NHPP with mean lambda*t^beta
beta, lam = 0.72, 0.45
n = 50
u = np.sort(rng.uniform(0, 1, n))
T = (n/lam)**(1/beta)
times = T * (u ** (1/beta))
return pd.DataFrame({"event_time": np.round(times, 2)})
def sample_repair_csv():
rng = np.random.default_rng(3)
rows = []
for sys in range(1, 6):
gaps = rng.exponential(80, size=10)
times = np.cumsum(gaps)
for tm in times[times < 600]:
rows.append({"system": f"S{sys}", "event_time": round(float(tm), 2)})
return pd.DataFrame(rows)
def sample_alt_csv():
rng = np.random.default_rng(5)
rows = []
for temp in [85, 105, 125]:
x = 1/(temp+273.15)
eta = np.exp(2.0 + 3800*x) # longer life at lower temperature
beta = 1.7
for i in range(30):
true_t = eta * rng.weibull(beta)
censor = rng.uniform(1000, 5000)
rows.append({"time": round(float(min(true_t, censor)), 2),
"event": int(true_t <= censor),
"temperature": temp})
return pd.DataFrame(rows)
# -----------------------------------------------------------------------------
# Presentation and export layer
# The analytical functions above are intentionally unchanged. The wrapper
# functions below only call the existing computations and save returned tables
# and figures as downloadable artifacts for the Gradio interface.
# -----------------------------------------------------------------------------
import os
import time
import uuid
# Gradio can only serve returned files safely from the current working directory,
# system temp directory, upload directory, or explicit allowed_paths.
EXPORT_DIR = os.path.join(tempfile.gettempdir(), "reliapy_exports")
os.makedirs(EXPORT_DIR, exist_ok=True)
def _export_token(prefix):
return f"{prefix}_{time.strftime('%Y%m%d_%H%M%S')}_{uuid.uuid4().hex[:6]}"
def _save_summary_csv(df, token):
path = os.path.join(EXPORT_DIR, f"{token}_summary.csv")
try:
if isinstance(df, pd.DataFrame):
df.to_csv(path, index=False)
return path
except Exception:
pass
return None
def _save_figure_png(fig, token, suffix):
if fig is None:
return None
path = os.path.join(EXPORT_DIR, f"{token}_{suffix}.png")
try:
fig.savefig(path, dpi=300, bbox_inches="tight")
return path
except Exception:
return None
def run_life_with_downloads(file_obj, dist, method, time_col, event_col, mission_time):
summary, fig1, fig2, fig3, note = fit_life_data(file_obj, dist, method, time_col, event_col, mission_time)
token = _export_token("life_data")
csv_path = _save_summary_csv(summary, token)
p1 = _save_figure_png(fig1, token, "probability_plot")
p2 = _save_figure_png(fig2, token, "likelihood_contour")
p3 = _save_figure_png(fig3, token, "fitted_cdf")
return summary, fig1, fig2, fig3, note, csv_path, p1, p2, p3
def run_growth_with_downloads(file_obj, model_type, time_col, breakpoints_text):
summary, fig1, fig2, note = fit_growth(file_obj, model_type, time_col, breakpoints_text)
token = _export_token("reliability_growth")
csv_path = _save_summary_csv(summary, token)
p1 = _save_figure_png(fig1, token, "growth_plot")
p2 = _save_figure_png(fig2, token, "duane_plot")
return summary, fig1, fig2, note, csv_path, p1, p2
def run_repairable_with_downloads(file_obj, model_type, time_col, system_col, breakpoints_text):
summary, fig1, fig2, fig3, note = fit_repairable(file_obj, model_type, time_col, system_col, breakpoints_text)
token = _export_token("repairable_system")
csv_path = _save_summary_csv(summary, token)
p1 = _save_figure_png(fig1, token, "cumulative_events")
p2 = _save_figure_png(fig2, token, "event_rate")
p3 = _save_figure_png(fig3, token, "mcf")
return summary, fig1, fig2, fig3, note, csv_path, p1, p2, p3
def run_alt_with_downloads(file_obj, dist, relationship, time_col, stress_col, event_col, use_stress):
summary, fig1, fig2, note = fit_alt(file_obj, dist, relationship, time_col, stress_col, event_col, use_stress)
token = _export_token("accelerated_life_testing")
csv_path = _save_summary_csv(summary, token)
p1 = _save_figure_png(fig1, token, "alt_probability_plot")
p2 = _save_figure_png(fig2, token, "life_stress_plot")
return summary, fig1, fig2, note, csv_path, p1, p2
APP_CSS = """
:root {
--rp-red: #d71920;
--rp-red-dark: #991b1b;
--rp-red-soft: rgba(215, 25, 32, .10);
--rp-ink: #0f172a;
--rp-slate: #334155;
--rp-muted: #64748b;
--rp-line: rgba(15, 23, 42, .10);
--rp-card: rgba(255, 255, 255, .92);
--rp-bg: #f5f7fb;
}
.gradio-container {
max-width: 1480px !important;
margin: 0 auto !important;
color: var(--rp-ink) !important;
background:
radial-gradient(circle at 8% 4%, rgba(215, 25, 32, .14), transparent 28%),
radial-gradient(circle at 92% 0%, rgba(2, 132, 199, .09), transparent 26%),
linear-gradient(180deg, #ffffff 0%, var(--rp-bg) 100%) !important;
}
#title-banner {
position: relative;
overflow: hidden;
padding: 34px 38px;
border-radius: 28px;
background:
linear-gradient(135deg, rgba(153, 27, 27, .98), rgba(215, 25, 32, .94) 45%, rgba(15, 23, 42, .96));
color: white;
margin: 14px 0 20px 0;
border: 1px solid rgba(255, 255, 255, .24);
box-shadow: 0 24px 70px rgba(15, 23, 42, .22);
}
#title-banner:before {
content: "";
position: absolute;
right: -95px;
top: -95px;
width: 340px;
height: 340px;
border-radius: 999px;
background: rgba(255, 255, 255, .11);
}
#title-banner:after {
content: "";
position: absolute;
right: 98px;
bottom: -120px;
width: 270px;
height: 270px;
border-radius: 999px;
border: 42px solid rgba(255, 255, 255, .07);
}
#title-banner h1 {
position: relative;
font-size: 2.8rem;
letter-spacing: -.05em;
margin: 0 0 9px 0;
line-height: 1.02;
}
#title-banner p {
position: relative;
max-width: 920px;
margin: 0;
font-size: 1.07rem;
color: rgba(255, 255, 255, .89);
line-height: 1.55;
}
#hero-kicker {
position: relative;
display: inline-flex;
align-items: center;
gap: 9px;
padding: 7px 13px;
border-radius: 999px;
background: rgba(255,255,255,.16);
color: rgba(255,255,255,.96);
font-size: .78rem;
font-weight: 800;
text-transform: uppercase;
letter-spacing: .09em;
margin-bottom: 14px;
}
.metric-card {
padding: 20px 20px;
min-height: 128px;
border: 1px solid var(--rp-line);
border-radius: 22px;
background: var(--rp-card);
box-shadow: 0 14px 38px rgba(15, 23, 42, .075);
transition: transform .18s ease, box-shadow .18s ease;
}
.metric-card:hover { transform: translateY(-2px); box-shadow: 0 20px 48px rgba(15, 23, 42, .10); }
.metric-card .icon {
display: inline-flex;
width: 38px;
height: 38px;
align-items: center;
justify-content: center;
border-radius: 13px;
color: #fff;
background: linear-gradient(135deg, var(--rp-red), #ef4444);
margin-bottom: 10px;
font-size: 1.1rem;
}
.metric-card h3 { margin: 0 0 7px 0; font-size: 1.02rem; color: var(--rp-ink); }
.metric-card p { margin: 0; color: var(--rp-muted); font-size: .91rem; line-height: 1.45; }
.module-intro {
padding: 18px 20px;
border: 1px solid var(--rp-line);
border-left: 7px solid var(--rp-red);
border-radius: 20px;
background: rgba(255, 255, 255, .88);
box-shadow: 0 10px 30px rgba(15,23,42,.055);
margin: 10px 0 13px 0;
}
.module-intro h2 { margin: 0 0 6px 0; font-size: 1.35rem; letter-spacing: -.02em; }
.module-intro p { margin: 0; color: var(--rp-muted); line-height: 1.5; }
.control-panel, .output-panel, .download-panel {
border: 1px solid var(--rp-line) !important;
border-radius: 22px !important;
background: rgba(255, 255, 255, .90) !important;
box-shadow: 0 12px 34px rgba(15, 23, 42, .06) !important;
padding: 14px !important;
}
.plot-card {
border: 1px solid var(--rp-line) !important;
border-radius: 20px !important;
background: rgba(255, 255, 255, .92) !important;
box-shadow: 0 10px 28px rgba(15, 23, 42, .055) !important;
padding: 9px !important;
}
.download-card {
border: 1px dashed rgba(215, 25, 32, .32) !important;
border-radius: 18px !important;
background: linear-gradient(180deg, rgba(255,255,255,.92), rgba(254,242,242,.55)) !important;
padding: 10px !important;
}
.footer-note {
margin-top: 16px;
padding: 17px 19px;
border-radius: 19px;
background: #fff7ed;
border: 1px solid rgba(249,115,22,.22);
color: #7c2d12;
line-height: 1.48;
}
.gr-button-primary {
border-radius: 13px !important;
font-weight: 800 !important;
box-shadow: 0 10px 24px rgba(215, 25, 32, .22) !important;
}
.gr-button-secondary, button { border-radius: 13px !important; }
.tabitem { padding-top: 14px !important; }
.block label, .gradio-container label { font-weight: 700 !important; color: #334155 !important; }
textarea, input, select { border-radius: 13px !important; }
.table-wrap, .dataframe { border-radius: 15px !important; }
.file-preview, .upload-container { border-radius: 15px !important; }
code { background: rgba(15,23,42,.06); padding: 2px 5px; border-radius: 6px; }
"""
INTRO_HTML = """
Reliability analysis · Python · Gradio
ReliaPy Workbench
Publication-oriented reliability analysis dashboard for life data, reliability growth, repairable systems, and accelerated life testing. The calculation layer is preserved; this version enriches the interface and adds downloadable CSV/PNG exports.
Developer: Partha Pratim Ray, Sikkim University · July 2026 · Email: parthapratimray1986@gmail.com
"""
CSV_NOTE = """
Input convention and export support
Upload CSV or Excel files with positive numeric times. Optional event/status columns accept 1/true/failure as failure and 0/false as right-censored. After every run, the result table can be downloaded as CSV and each plot can be downloaded as a 300-dpi PNG image.
"""
def download_area(csv_component, plot_components):
gr.Markdown("### ⬇️ Downloads")
with gr.Row():
csv_component.render()
with gr.Row():
for comp in plot_components:
comp.render()
with gr.Blocks(
css=APP_CSS,
theme=gr.themes.Soft(primary_hue="red", secondary_hue="slate", neutral_hue="slate"),
title="ReliaPy Workbench"
) as demo:
gr.Markdown(INTRO_HTML)
with gr.Row(equal_height=True):
gr.Markdown("""⏳
Life Data Analysis
Weibull and Lognormal models using MLE or Rank Regression with probability, CDF, and likelihood-contour plots.
↗
Reliability Growth
Crow-AMSAA, piecewise NHPP, and automatic change-point analysis with growth and Duane visualizations.
🔧
Repairable Systems
Power Law, Log-Linear, and Piecewise NHPP models with cumulative events, event rate, and MCF outputs.
⚡
Accelerated Life Testing
Weibull or Lognormal ALT with Arrhenius or Power Law stress relationships and use-condition prediction.
Life Data Analysis
Fit Weibull or Lognormal life distributions using MLE or median-rank regression. Outputs include parameter estimates, reliability at mission time, probability plot, likelihood contour, fitted CDF, and downloadable artifacts.
Reliability Growth Analysis
Model cumulative test failures using Crow-AMSAA/Power Law NHPP, manual piecewise NHPP, or automatic one-change-point detection.
Repairable Systems
Analyze recurrent repair/failure events with Power Law, Log-Linear, or Piecewise NHPP models. Includes cumulative events, event rate, Mean Cumulative Function, and downloadable exports.
Accelerated Life Testing
Fit accelerated life models under Arrhenius or Power Law stress relationships with common Weibull shape or common Lognormal sigma.