## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4
)

## ----mixed-types--------------------------------------------------------------
library(serieshaz)

disk <- dfr_weibull(shape = 2, scale = 500)
mem  <- dfr_exponential(0.001)
psu  <- dfr_gompertz(a = 0.0001, b = 0.02)

server <- dfr_dist_series(list(disk, mem, psu))

## ----hazard-decomp------------------------------------------------------------
h_disk <- component_hazard(server, 1)
h_mem  <- component_hazard(server, 2)
h_psu  <- component_hazard(server, 3)
h_sys  <- hazard(server)

t_grid <- seq(10, 300, by = 10)
cat("Time | Disk     | Memory   | PSU      | System\n")
cat("-----|----------|----------|----------|--------\n")
for (t0 in c(50, 100, 150, 200, 250)) {
    cat(sprintf("%4d | %.6f | %.6f | %.6f | %.6f\n",
                t0, h_disk(t0), h_mem(t0), h_psu(t0), h_sys(t0)))
}

## ----nested-------------------------------------------------------------------
# CPU subsystem: two exponential failure modes
cpu <- dfr_dist_series(list(
    dfr_exponential(0.001),  # thermal failure
    dfr_exponential(0.0005)  # electrical failure
))

# Storage subsystem: Weibull wear + exponential shock
storage <- dfr_dist_series(list(
    dfr_weibull(shape = 2, scale = 1000),
    dfr_exponential(0.0002)
))

# Full system: CPU + storage + power supply
full <- dfr_dist_series(list(
    cpu,
    storage,
    dfr_gompertz(a = 0.0001, b = 0.01)
))

cat("Full system components:", ncomponents(full), "\n")
cat("Total parameters:", length(params(full)), "\n")

## ----nested-verify------------------------------------------------------------
h_full <- hazard(full)

# Manually sum all leaf hazards
h_cpu1 <- component_hazard(cpu, 1)
h_cpu2 <- component_hazard(cpu, 2)
h_stor1 <- component_hazard(storage, 1)
h_stor2 <- component_hazard(storage, 2)
h_psu <- component_hazard(full, 3)

t0 <- 100
nested_sum <- h_cpu1(t0) + h_cpu2(t0) + h_stor1(t0) + h_stor2(t0) + h_psu(t0)
system_h <- h_full(t0)

cat(sprintf("Sum of leaf hazards: %.8f\n", nested_sum))
cat(sprintf("System hazard:       %.8f\n", system_h))
cat(sprintf("Difference:          %.2e\n", abs(nested_sum - system_h)))

## ----custom-component---------------------------------------------------------
# Bathtub: h(t) = a/t^b + c*t^d  (infant mortality + wear-out)
bathtub <- dfr_dist(
    rate = function(t, par, ...) {
        a <- par[1]; b <- par[2]; c <- par[3]; d <- par[4]
        a * t^(-b) + c * t^d
    },
    par = c(0.5, 0.5, 0.0001, 2)
)

# Series system with bathtub + exponential random shock
sys <- dfr_dist_series(list(bathtub, dfr_exponential(0.001)))

h_sys <- hazard(sys)
for (t0 in c(1, 10, 50, 100, 200)) {
    cat(sprintf("t=%3d: h_sys = %.6f\n", t0, h_sys(t0)))
}

## ----custom-cumhaz------------------------------------------------------------
cat("Analytical cum_haz available:", !is.null(sys$cum_haz_rate), "\n")

# The system still computes survival correctly via numerical integration
S <- surv(sys)
for (t0 in c(10, 50, 100)) {
    cat(sprintf("t=%3d: S(t) = %.6f\n", t0, S(t0)))
}

## ----failure-attribution------------------------------------------------------
server <- dfr_dist_series(list(
    dfr_weibull(shape = 2, scale = 200),     # mechanical wear
    dfr_exponential(0.005),            # random shock
    dfr_gompertz(a = 0.0001, b = 0.02)       # degradation
))

set.seed(42)
n <- 10000
mat <- sample_components(server, n = n)

# System lifetime = minimum across components
t_sys <- apply(mat, 1, min)

# Identify the failing component for each observation
failing <- apply(mat, 1, which.min)

## ----attribution-props--------------------------------------------------------
props <- table(failing) / n
names(props) <- c("Weibull (wear)", "Exp (shock)", "Gompertz (degrad.)")
cat("Failure attribution:\n")
print(round(props, 3))

## ----attribution-time---------------------------------------------------------
# Split into early, middle, and late failures
breaks <- quantile(t_sys, probs = c(0, 1/3, 2/3, 1))
period <- cut(t_sys, breaks = breaks, labels = c("Early", "Middle", "Late"),
              include.lowest = TRUE)

cat("\nAttribution by time period:\n")
for (p in levels(period)) {
    idx <- which(period == p)
    props_p <- table(factor(failing[idx], levels = 1:3)) / length(idx)
    cat(sprintf("  %s (n=%d): Wear=%.1f%% Shock=%.1f%% Degrad=%.1f%%\n",
                p, length(idx),
                100 * props_p[1], 100 * props_p[2], 100 * props_p[3]))
}

## ----sensitivity--------------------------------------------------------------
server <- dfr_dist_series(list(
    dfr_weibull(shape = 2, scale = 100),
    dfr_exponential(0.01)
))

S_fn <- surv(server)
base_par <- params(server)  # c(2, 100, 0.01)
t0 <- 50

cat("Sensitivity of S(50) to Weibull scale parameter:\n")
for (scale_val in c(50, 75, 100, 125, 150)) {
    par_mod <- base_par
    par_mod[2] <- scale_val  # layout tells us par[2] is Weibull scale
    cat(sprintf("  scale = %3d: S(50) = %.4f\n", scale_val, S_fn(t0, par = par_mod)))
}

cat("\nSensitivity of S(50) to exponential rate:\n")
for (rate_val in c(0.005, 0.01, 0.02, 0.05, 0.1)) {
    par_mod <- base_par
    par_mod[3] <- rate_val  # layout tells us par[3] is exp rate
    cat(sprintf("  rate = %.3f: S(50) = %.4f\n", rate_val, S_fn(t0, par = par_mod)))
}