\begin{lstlisting}[style=pythonstyle,caption={Python code},label={lst:example}]


#!/usr/bin/env python3
"""
Supplementary script for:
"The 'b's' Knees: Logistic Stage-Transition Rules as Percentile Conventions on a Standardised Time Scale"

This script reproduces the synthetic logistic example and the free segmented-regression
breakpoint search described in the paper.

It:
1. Generates the exact synthetic logistic with:
      b = 0.35, t0 = 2012, K = 95
   on the symmetric window [t0 - 6/b, t0 + 6/b] at intervals of 0.25 time units.

2. Computes analytic stage-transition dates for:
      - third-derivative rule
      - tangent-at-inflection rule
      - 10% to 90% window
      - 15% to 85% window

3. Fits a continuous three-segment linear model
      y_hat(t) = alpha + beta1*t + beta2*max(0, t - tau1) + beta3*max(0, t - tau2),
   with tau1 < tau2,
   by brute-force search over candidate breakpoint pairs on the observation grid.

4. Reports the least-squares breakpoint pair and plots the synthetic example.

Dependencies:
    numpy
    matplotlib

Usage:
    python logistic_knees_si.py
"""

import math
import numpy as np
import matplotlib.pyplot as plt


# ---------------------------------------------------------------------
# Core functions
# ---------------------------------------------------------------------

def logistic(t, K, b, t0):
    """Standard logistic curve."""
    return K / (1.0 + np.exp(-b * (t - t0)))


def hinge_design_matrix(t, tau1, tau2):
    """
    Design matrix for the continuous three-segment linear model:
        y_hat(t) = alpha + beta1*t + beta2*max(0, t - tau1) + beta3*max(0, t - tau2)
    """
    return np.column_stack([
        np.ones_like(t),
        t,
        np.maximum(0.0, t - tau1),
        np.maximum(0.0, t - tau2),
    ])


def fit_segmented_least_squares(t, y, tau1, tau2):
    """
    Fit the continuous three-segment linear model by OLS for fixed tau1, tau2.
    Returns coefficients, fitted values, and residual sum of squares.
    """
    X = hinge_design_matrix(t, tau1, tau2)
    beta, residuals, rank, s = np.linalg.lstsq(X, y, rcond=None)
    y_hat = X @ beta
    rss = float(np.sum((y - y_hat) ** 2))
    return beta, y_hat, rss


def brute_force_breakpoints(t, y, min_frac=0.10):
    """
    Brute-force search over breakpoint pairs (tau1, tau2) on the observation grid.

    Constraint:
        At least min_frac of observations must lie in each of the three segments:
            [start, tau1], (tau1, tau2], (tau2, end]

    Returns:
        best_tau1, best_tau2, best_beta, best_yhat, best_rss
    """
    n = len(t)
    min_count = math.ceil(min_frac * n)

    best = None

    # Breakpoints restricted to observed time grid
    # Use indices i < j so tau1 = t[i], tau2 = t[j]
    for i in range(min_count - 1, n - 2 * min_count + 1):
        for j in range(i + min_count, n - min_count):
            left_count = i + 1
            middle_count = j - i
            right_count = n - (j + 1)

            if left_count < min_count or middle_count < min_count or right_count < min_count:
                continue

            tau1 = t[i]
            tau2 = t[j]

            beta, y_hat, rss = fit_segmented_least_squares(t, y, tau1, tau2)

            if (best is None) or (rss < best["rss"]):
                best = {
                    "tau1": float(tau1),
                    "tau2": float(tau2),
                    "beta": beta,
                    "y_hat": y_hat,
                    "rss": rss,
                }

    if best is None:
        raise RuntimeError("No feasible breakpoint pair found under the segment-size constraint.")

    return best["tau1"], best["tau2"], best["beta"], best["y_hat"], best["rss"]


def logistic_time_for_percentile(p, b, t0):
    """
    Return t such that logistic(t)/K = p.
    Since p = 1 / (1 + exp(-b(t - t0))),
        t = t0 + log(p / (1 - p)) / b
    """
    if not (0.0 < p < 1.0):
        raise ValueError("p must lie strictly between 0 and 1.")
    return t0 + math.log(p / (1.0 - p)) / b


# ---------------------------------------------------------------------
# Parameters from the manuscript
# ---------------------------------------------------------------------

K = 95.0
b = 0.35
t0 = 2012.0

# Observation window: [t0 - 6/b, t0 + 6/b] at intervals of 0.25
t_min = t0 - 6.0 / b
t_max = t0 + 6.0 / b
dt = 0.25

# Include endpoint with small tolerance
t = np.arange(t_min, t_max + 0.5 * dt, dt)
y = logistic(t, K=K, b=b, t0=t0)


# ---------------------------------------------------------------------
# Analytic rules
# ---------------------------------------------------------------------

# Third-derivative rule
c_third = math.log(2.0 + math.sqrt(3.0))
t_third_L = t0 - c_third / b
t_third_R = t0 + c_third / b
p_third_L = (3.0 - math.sqrt(3.0)) / 6.0
p_third_R = (3.0 + math.sqrt(3.0)) / 6.0

# Tangent-at-inflection rule
c_tangent = 2.0
t_tangent_L = t0 - c_tangent / b
t_tangent_R = t0 + c_tangent / b
p_tangent_L = 1.0 / (1.0 + math.exp(c_tangent))
p_tangent_R = 1.0 / (1.0 + math.exp(-c_tangent))

# 10% to 90%
c_10_90 = math.log(0.90 / 0.10)
t_10_90_L = logistic_time_for_percentile(0.10, b=b, t0=t0)
t_10_90_R = logistic_time_for_percentile(0.90, b=b, t0=t0)

# 15% to 85%
c_15_85 = math.log(0.85 / 0.15)
t_15_85_L = logistic_time_for_percentile(0.15, b=b, t0=t0)
t_15_85_R = logistic_time_for_percentile(0.85, b=b, t0=t0)


# ---------------------------------------------------------------------
# Segmented regression by brute-force search
# ---------------------------------------------------------------------

tau1_hat, tau2_hat, beta_hat, y_hat, rss_hat = brute_force_breakpoints(t, y, min_frac=0.10)


# ---------------------------------------------------------------------
# Console output
# ---------------------------------------------------------------------

print("\nSynthetic logistic example")
print("--------------------------")
print(f"K  = {K}")
print(f"b  = {b}")
print(f"t0 = {t0}")
print(f"Observation window = [{t_min:.6f}, {t_max:.6f}]")
print(f"Number of observations = {len(t)}")
print(f"Time step = {dt}")

print("\nAnalytic stage-transition dates")
print("-------------------------------")
print(f"Third-derivative rule:   c = {c_third:.6f}, dates = ({t_third_L:.2f}, {t_third_R:.2f}), "
      f"window = ({100*p_third_L:.2f}%, {100*p_third_R:.2f}%)")
print(f"Tangent-at-inflection:   c = {c_tangent:.6f}, dates = ({t_tangent_L:.2f}, {t_tangent_R:.2f}), "
      f"window = ({100*p_tangent_L:.2f}%, {100*p_tangent_R:.2f}%)")
print(f"10% to 90% convention:   c = {c_10_90:.6f}, dates = ({t_10_90_L:.2f}, {t_10_90_R:.2f})")
print(f"15% to 85% convention:   c = {c_15_85:.6f}, dates = ({t_15_85_L:.2f}, {t_15_85_R:.2f})")

print("\nSegmented regression results")
print("----------------------------")
print(f"Estimated breakpoints: tau1 = {tau1_hat:.2f}, tau2 = {tau2_hat:.2f}")
print(f"Residual sum of squares: {rss_hat:.10f}")
print("Estimated coefficients (alpha, beta1, beta2, beta3):")
print(beta_hat)

# Expected values close to those reported in the manuscript:
# tau1 approx 2006.36, tau2 approx 2017.61


# ---------------------------------------------------------------------
# Plot
# ---------------------------------------------------------------------

fig, ax = plt.subplots(figsize=(9, 5.5))

# Main curves
ax.plot(t, y, label="Logistic curve", linewidth=2.5)
ax.plot(t, y_hat, label="Free segmented regression", linewidth=2.5)

# Vertical markers
ax.axvline(t_third_L, linestyle="--", linewidth=1.5, label="Third-derivative")
ax.axvline(t_third_R, linestyle="--", linewidth=1.5)

ax.axvline(t_tangent_L, linestyle="--", linewidth=1.5, label="Tangent-at-inflection")
ax.axvline(t_tangent_R, linestyle="--", linewidth=1.5)

ax.axvline(tau1_hat, linestyle=":", linewidth=2.0, label="Segmented-regression breakpoints")
ax.axvline(tau2_hat, linestyle=":", linewidth=2.0)

# Guide lines
ax.axvline(t0, color="gray", linewidth=1.0, alpha=0.8)
ax.axhline(K / 2.0, color="gray", linewidth=1.0, alpha=0.8)
ax.axhline(0.10 * K, color="gray", linewidth=0.8, alpha=0.6)
ax.axhline(0.15 * K, color="gray", linewidth=0.8, alpha=0.6)
ax.axhline(0.85 * K, color="gray", linewidth=0.8, alpha=0.6)
ax.axhline(0.90 * K, color="gray", linewidth=0.8, alpha=0.6)

# Annotations
ax.text(t0 + 0.15, K / 2.0 + 1.5, r"$K/2$", fontsize=10)
ax.text(t0 + 0.15, 0.10 * K + 1.5, r"$0.10K$", fontsize=10)
ax.text(t0 + 0.15, 0.15 * K + 1.5, r"$0.15K$", fontsize=10)
ax.text(t0 + 0.15, 0.85 * K + 1.5, r"$0.85K$", fontsize=10)
ax.text(t0 + 0.15, 0.90 * K + 1.5, r"$0.90K$", fontsize=10)
ax.text(t0 + 0.15, K - 6, r"$t_0$", fontsize=10)

ax.set_xlabel("Time")
ax.set_ylabel(r"$y(t)$")
ax.set_title("Synthetic logistic example and stage-transition estimators")
ax.set_xlim(t_min, t_max)
ax.set_ylim(min(-2, np.min(y_hat) - 1), max(K + 5, np.max(y_hat) + 1))
ax.legend(frameon=False)
ax.grid(False)

plt.tight_layout()
plt.show()

\end{lstlisting}