Try a small prediction example: source code
This source belongs to the new teaching example, not a recovered employer artifact or the original research implementation.
Open raw demo.py"""Synthetic lead-pipe evaluation companion, Python 3.10+, standard library only.
Newly authored educational code, not the original utility implementation or
field performance. Never use its outputs for material or public-health decisions.
"""
from __future__ import annotations
import argparse
import hashlib
import json
import math
import random
from dataclasses import dataclass
from html import escape
from pathlib import Path
DATA_SEED = 2026
SPLIT_SEED = 4173
THRESHOLD = 0.5
EPOCHS = 900
RATE = 0.4
PENALTY = 0.002
DISCLOSURE = (
"Entirely synthetic educational companion, not the original TF-IDF/gradient-boosting "
"prototype, field validation, verified material classifications, or evidence of savings."
)
@dataclass(frozen=True)
class Record:
property_id: str
group_id: str
installation_age: float | None
keyword_evidence: int
conflicting_records: int
target: int
def sigmoid(value: float) -> float:
if value >= 0:
return 1.0 / (1.0 + math.exp(-value))
exp_value = math.exp(value)
return exp_value / (1.0 + exp_value)
def make_records(seed: int = DATA_SEED, groups: int = 160, per_group: int = 5) -> list[Record]:
"""Invent records with noisy evidence and an intentionally learnable signal."""
if groups < 4 or per_group < 1:
raise ValueError("Use at least four groups and one property per group.")
rng = random.Random(seed)
records = []
for group in range(groups):
group_id = f"synthetic-group-{group:03d}"
cohort_age = rng.uniform(10, 110)
unobserved_group_effect = rng.gauss(0, 0.6)
for property_number in range(per_group):
age = min(120, max(0, cohort_age + rng.gauss(0, 12)))
material_probability = sigmoid(-3.5 + age / 19 + unobserved_group_effect)
target = int(rng.random() < material_probability)
# Simulated evidence observes the state noisily; it is not a label rule.
keyword = int(rng.random() < (0.72 if target else 0.23))
conflict = int(rng.random() < 0.15)
if conflict:
keyword = int(rng.random() < 0.5)
observed_age = None if rng.random() < 0.12 else round(age, 1)
records.append(Record(
f"{group_id}-property-{property_number:03d}", group_id,
observed_age, keyword, conflict, target
))
return records
def grouped_split(records: list[Record], seed: int = SPLIT_SEED,
test_fraction: float = 0.25) -> tuple[list[Record], list[Record]]:
"""Whole groups held out by stable hash, without inspecting predictors or labels."""
if not 0 < test_fraction < 1:
raise ValueError("test_fraction must lie strictly between zero and one.")
if len({row.property_id for row in records}) != len(records):
raise ValueError("Expected exactly one prepared row per property.")
groups = sorted({row.group_id for row in records}, key=lambda group: (
hashlib.sha256(f"{seed}:{group}".encode("utf-8")).hexdigest(), group
))
if len(groups) < 2:
raise ValueError("At least two groups are required.")
count = min(len(groups) - 1, max(1, round(len(groups) * test_fraction)))
held_out = set(groups[:count])
train = sorted((row for row in records if row.group_id not in held_out),
key=lambda row: row.property_id)
test = sorted((row for row in records if row.group_id in held_out),
key=lambda row: row.property_id)
return train, test
def features(record: Record) -> tuple[float, ...]:
"""Fixed transforms; no test-fitted scaling, identifiers, or target features."""
age = record.installation_age if record.installation_age is not None else 60.0
return (1.0, (age - 60.0) / 60.0, float(record.keyword_evidence),
float(record.conflicting_records), float(record.installation_age is None))
def predict(weights: list[float], record: Record) -> float:
if len(weights) != len(features(record)):
raise ValueError("Expected one weight per feature.")
return sigmoid(sum(weight * value for weight, value in zip(weights, features(record))))
def fit(train: list[Record], epochs: int = EPOCHS, rate: float = RATE,
penalty: float = PENALTY) -> list[float]:
"""Full-batch L2 logistic regression; no test-driven tuning or calibration."""
if not train or epochs < 1 or rate <= 0 or penalty < 0:
raise ValueError("Training data and positive optimization settings are required.")
if not math.isfinite(rate) or not math.isfinite(penalty):
raise ValueError("Optimization settings must be finite.")
weights = [0.0] * 5
matrix = [(features(row), row.target) for row in train]
for _ in range(epochs):
gradient = [0.0] * len(weights)
for values, target in matrix:
error = sigmoid(sum(w * x for w, x in zip(weights, values))) - target
for index, value in enumerate(values):
gradient[index] += error * value
weights = [weight - rate * (gradient[index] / len(train) +
(penalty * weight if index else 0.0))
for index, weight in enumerate(weights)]
return weights
def validate_predictions(targets: list[int], probabilities: list[float]) -> None:
if not targets or len(targets) != len(probabilities):
raise ValueError("Provide equally sized nonempty target and probability lists.")
if any(target not in (0, 1) for target in targets):
raise ValueError("Targets must be binary.")
if any(not math.isfinite(p) or not 0 <= p <= 1 for p in probabilities):
raise ValueError("Probabilities must be finite and between zero and one.")
def evaluate(targets: list[int], probabilities: list[float], threshold: float = THRESHOLD) -> dict:
validate_predictions(targets, probabilities)
if not math.isfinite(threshold) or not 0 <= threshold <= 1:
raise ValueError("Threshold must lie between zero and one.")
clipped = [min(1 - 1e-12, max(1e-12, p)) for p in probabilities]
confusion = {"true_negative": 0, "false_positive": 0,
"false_negative": 0, "true_positive": 0}
names = {(0, False): "true_negative", (0, True): "false_positive",
(1, False): "false_negative", (1, True): "true_positive"}
for target, probability in zip(targets, probabilities):
confusion[names[target, probability >= threshold]] += 1
size = len(targets)
return {
"n": size, "positive_labels": sum(targets), "threshold": threshold,
"log_loss": sum(-y * math.log(p) - (1 - y) * math.log(1 - p)
for y, p in zip(targets, clipped)) / size,
"brier_score": sum((y - p) ** 2 for y, p in zip(targets, probabilities)) / size,
"accuracy": (confusion["true_positive"] + confusion["true_negative"]) / size,
"confusion": confusion,
}
def calibration_bins(targets: list[int], probabilities: list[float], bins: int = 5) -> list[dict]:
"""Descriptive equal-width reliability bins, not fitted probability calibration."""
validate_predictions(targets, probabilities)
if not isinstance(bins, int) or bins < 1:
raise ValueError("Use a positive integer number of bins.")
grouped = [[] for _ in range(bins)]
for target, probability in zip(targets, probabilities):
grouped[min(bins - 1, int(probability * bins))].append((target, probability))
return [{
"interval": f"[{i / bins:.1f}, {(i + 1) / bins:.1f}{']' if i == bins - 1 else ')'}",
"count": len(items),
"mean_prediction": sum(p for _, p in items) / len(items) if items else None,
"observed_fraction": sum(y for y, _ in items) / len(items) if items else None,
} for i, items in enumerate(grouped)]
def error_slices(records: list[Record], probabilities: list[float]) -> list[dict]:
validate_predictions([row.target for row in records], probabilities)
predicates = [
("Age present", lambda row: row.installation_age is not None),
("Age missing", lambda row: row.installation_age is None),
("No conflicting record", lambda row: row.conflicting_records == 0),
("Conflicting record", lambda row: row.conflicting_records == 1),
]
results = []
for name, predicate in predicates:
selected = [(row.target, p) for row, p in zip(records, probabilities) if predicate(row)]
metrics = evaluate([y for y, _ in selected], [p for _, p in selected]) if selected else None
results.append({"slice": name, "metrics": metrics})
return results
def fmt(value: float | None) -> str:
return "n/a" if value is None else f"{value:.4f}"
def run(seed: int = DATA_SEED) -> dict:
train, test = grouped_split(make_records(seed=seed))
shared_groups = {row.group_id for row in train} & {row.group_id for row in test}
shared_properties = {row.property_id for row in train} & {row.property_id for row in test}
assert not shared_groups and not shared_properties
weights = fit(train)
probabilities = [predict(weights, row) for row in test]
prevalence = sum(row.target for row in train) / len(train)
targets = [row.target for row in test]
model = evaluate(targets, probabilities)
baseline = evaluate(targets, [prevalence] * len(test))
reliability = calibration_bins(targets, probabilities)
slices = error_slices(test, probabilities)
errors = sorted(((row, p) for row, p in zip(test, probabilities)
if (p >= THRESHOLD) != bool(row.target)), key=lambda item: (
-abs(item[1] - item[0].target), item[0].property_id
))[:5]
review = sorted(zip(test, probabilities), key=lambda item: (
abs(item[1] - THRESHOLD), item[0].property_id
))[:5]
notes = [
"All 800 records, group IDs, evidence and target labels are invented; signal is intentionally learnable.",
"Fixed seeds: data " + str(seed) + ", split 4173. Whole-group 75/25 split; no shared groups or properties.",
"Features and training settings were fixed before evaluation. No tuning, threshold choice or calibration uses test labels.",
"Reliability bins describe this one synthetic holdout; they do not establish calibrated probabilities.",
"The 0.5 threshold is illustrative, not an approved material classification or inspection rule.",
"Age and record-conflict slices overlap; small slices have unstable estimates. No confidence intervals are claimed.",
"Real evaluation needs verified labels and spatial, temporal and cross-utility holdouts; this example has none.",
"This companion is separate from the original professional OCR effort and any historical cost savings.",
]
tables = [
{"heading": "Held-out model comparison (lower loss is better)",
"headers": ["Predictor", "Brier score", "Log loss", "Accuracy at 0.5"],
"rows": [["Logistic model", fmt(model["brier_score"]), fmt(model["log_loss"]), fmt(model["accuracy"])],
["Training-prevalence baseline", fmt(baseline["brier_score"]), fmt(baseline["log_loss"]), fmt(baseline["accuracy"])]]},
{"heading": "Confusion counts at illustrative threshold 0.5",
"headers": ["Predictor", "True negative", "False positive", "False negative", "True positive"],
"rows": [[name, *[metrics["confusion"][key] for key in (
"true_negative", "false_positive", "false_negative", "true_positive")]]
for name, metrics in (("Logistic model", model), ("Prevalence baseline", baseline))]},
{"heading": "Reliability bins (descriptive, not fitted calibration)",
"headers": ["Probability interval", "Records", "Mean prediction", "Observed positive fraction"],
"rows": [[item["interval"], item["count"], fmt(item["mean_prediction"]),
fmt(item["observed_fraction"])] for item in reliability]},
{"heading": "Error slices (overlapping, descriptive)",
"headers": ["Slice", "Records", "Brier score", "False positive", "False negative"],
"rows": [[item["slice"], item["metrics"]["n"] if item["metrics"] else 0,
fmt(item["metrics"]["brier_score"]) if item["metrics"] else "n/a",
item["metrics"]["confusion"]["false_positive"] if item["metrics"] else "n/a",
item["metrics"]["confusion"]["false_negative"] if item["metrics"] else "n/a"] for item in slices]},
{"heading": "Five largest probability errors (invented held-out properties)",
"headers": ["Synthetic property", "Simulated target", "Prediction", "Age missing", "Record conflict"],
"rows": [[row.property_id, row.target, fmt(p), "Yes" if row.installation_age is None else "No",
"Yes" if row.conflicting_records else "No"] for row, p in errors]},
]
return {
"title": "Lead pipe prediction: synthetic evaluation",
"disclosure": DISCLOSURE,
"metrics": [{"label": "Held-out properties", "value": len(test)},
{"label": "Held-out groups", "value": len({row.group_id for row in test})},
{"label": "Shared groups", "value": len(shared_groups)},
{"label": "Model Brier score", "value": fmt(model["brier_score"])}],
"tables": tables, "notes": notes, "seed": seed,
"settings": {"split_seed": SPLIT_SEED, "test_fraction": 0.25, "threshold": THRESHOLD,
"epochs": EPOCHS, "learning_rate": RATE, "l2_penalty": PENALTY,
"calibration_bins": 5, "optimizer": "Full-batch gradient descent"},
"split": {"train_properties": len(train), "test_properties": len(test),
"train_groups": len({row.group_id for row in train}),
"test_groups": len({row.group_id for row in test}),
"shared_groups": len(shared_groups), "shared_properties": len(shared_properties)},
"learned_weights": weights, "training_prevalence": prevalence,
"logistic_model": model, "training_prevalence_baseline": baseline,
"reliability_bins": reliability, "error_slices": slices,
"illustrative_uncertainty_review": [
{"property_id": row.property_id, "synthetic_probability": p} for row, p in review
],
}
def chart_svg(report: dict) -> str:
"""Produce a self-contained, accessible vector chart from the run's raw metrics."""
model = report["logistic_model"]
baseline = report["training_prevalence_baseline"]
split = report["split"]
loss_max = max(1.0, math.ceil(max(metrics[key] for metrics in (model, baseline)
for key in ("brier_score", "log_loss")) * 4) / 4)
loss_width = 300
parts = [
'<svg xmlns="http://www.w3.org/2000/svg" width="1200" height="866" viewBox="0 0 1200 866" role="img" aria-labelledby="title desc">',
'<title id="title">Synthetic lead-pipe model evaluation</title>',
'<desc id="desc">Brier and log-loss comparison against training prevalence, plus a reliability plot. '
f'Generated from {split["test_properties"]} held-out invented records in {split["test_groups"]} disjoint groups. '
'Not real-world performance. Both loss metrics use one common zero-based scale; '
'the reliability plot uses equal zero-to-one axes.</desc>',
'<rect width="1200" height="866" fill="#fff"/>',
'<g font-family="Arial, sans-serif" fill="#253746">',
'<text x="40" y="47" font-size="16" fill="#5f6b73">Synthetic example</text>',
'<text x="40" y="94" font-size="32" font-weight="700">Lead-pipe model evaluation</text>',
'<text x="40" y="125" font-size="18" fill="#5f6b73">A reproducible comparison on invented records, with a fully separate holdout.</text>',
f'<text x="40" y="174" font-size="18">Holdout: {split["test_properties"]} properties in {split["test_groups"]} groups</text>',
f'<text x="40" y="203" font-size="17" fill="#5f6b73">Shared groups: {split["shared_groups"]} | Data seed: {report["seed"]}</text>',
'<line x1="40" y1="229" x2="1160" y2="229" stroke="#d9dfe3"/>',
'<text x="60" y="280" font-size="23" font-weight="700">Held-out losses</text>',
'<text x="60" y="309" font-size="16" fill="#5f6b73">Lower is better. Both metrics share the same scale.</text>',
'<circle cx="68" cy="338" r="5" fill="#147d86"/>',
'<text x="82" y="344" font-size="15">Logistic model</text>',
'<circle cx="250" cy="338" r="5" fill="#87939b"/>',
'<text x="264" y="344" font-size="15">Training-prevalence baseline</text>',
]
for block, (key, label) in enumerate((("brier_score", "Brier score"), ("log_loss", "Log loss"))):
y = 390 + block * 149
parts.append(f'<text x="60" y="{y}" font-size="19" font-weight="700">{label}</text>')
for index, (metrics, color, name) in enumerate(((model, "#147d86", "Model"), (baseline, "#87939b", "Baseline"))):
line_y = y + 19 + index * 43
width = metrics[key] / loss_max * loss_width
parts.extend([
f'<text x="60" y="{line_y + 20}" font-size="16" fill="#5f6b73">{name}</text>',
f'<rect x="151" y="{line_y}" width="{width:.3f}" height="28" fill="{color}" '
f'data-metric="{key}" data-predictor="{name.lower()}"><title>{label}, {name}: {metrics[key]:.4f}</title></rect>',
f'<text x="465" y="{line_y + 20}" font-size="17" font-weight="700">{metrics[key]:.4f}</text>',
])
for tick in range(5):
x = 151 + tick * loss_width / 4
parts.extend([
f'<line x1="{x:g}" y1="644" x2="{x:g}" y2="651" stroke="#87939b"/>',
f'<text x="{x:g}" y="673" text-anchor="middle" font-size="14" fill="#5f6b73">{tick * loss_max / 4:g}</text>',
])
parts.extend([
'<line x1="151" y1="644" x2="451" y2="644" stroke="#87939b"/>',
f'<text x="60" y="711" font-size="15" fill="#5f6b73">Common loss scale: 0 to {loss_max:g}. Bar lengths are proportional.</text>',
'<text x="630" y="280" font-size="23" font-weight="700">Reliability by probability bin</text>',
'<text x="630" y="309" font-size="16" fill="#5f6b73">Descriptive check; no calibrator was fitted.</text>',
])
plot_left, plot_bottom, plot_size = 724, 654, 326
for tick in range(6):
x = plot_left + tick * plot_size / 5
y = plot_bottom - tick * plot_size / 5
parts.extend([
f'<line x1="724" y1="{y:g}" x2="1050" y2="{y:g}" stroke="#e5e9ec"/>',
f'<line x1="{x:g}" y1="328" x2="{x:g}" y2="654" stroke="#e5e9ec"/>',
f'<text x="710" y="{y + 5:g}" text-anchor="end" font-size="14" fill="#5f6b73">{tick / 5:.1f}</text>',
f'<text x="{x:g}" y="676" text-anchor="middle" font-size="14" fill="#5f6b73">{tick / 5:.1f}</text>',
])
parts.extend([
'<path d="M724 328V654H1050" fill="none" stroke="#87939b" stroke-width="1.5"/>',
'<line x1="724" y1="654" x2="1050" y2="328" stroke="#87939b" stroke-width="2" stroke-dasharray="7 6"/>',
])
for item in report["reliability_bins"]:
if not item["count"]:
continue
x = plot_left + item["mean_prediction"] * plot_size
y = plot_bottom - item["observed_fraction"] * plot_size
label = escape(f'{item["interval"]}: n={item["count"]}, mean prediction {item["mean_prediction"]:.4f}, observed {item["observed_fraction"]:.4f}')
parts.append(f'<circle cx="{x:.3f}" cy="{y:.3f}" r="7" fill="#147d86" stroke="#fff" stroke-width="2" data-reliability-bin="{escape(item["interval"])}"><title>{label}</title></circle>')
parts.append(f'<text x="{x + 12:.3f}" y="{y - 12:.3f}" font-size="14" font-weight="700">n={item["count"]}</text>')
parts.extend([
'<text x="887" y="700" text-anchor="middle" font-size="16">Mean predicted probability</text>',
'<text transform="translate(658 491) rotate(-90)" text-anchor="middle" font-size="16">Observed positive fraction</text>',
'<line x1="726" y1="721" x2="752" y2="721" stroke="#87939b" stroke-width="2" stroke-dasharray="7 6"/>',
'<text x="761" y="725" font-size="12" fill="#5f6b73">Perfect agreement</text>',
'<circle cx="922" cy="721" r="4" fill="#147d86"/>',
'<text x="933" y="725" font-size="12" fill="#5f6b73">Bin mean; n = records</text>',
'<line x1="40" y1="755" x2="1160" y2="755" stroke="#d9dfe3"/>',
'<text x="54" y="783" font-size="17" font-weight="700">Synthetic-only results, not field-validated utility performance.</text>',
f'<text x="54" y="806" font-size="15" fill="#5f6b73">Split seed {report["settings"]["split_seed"]}. Threshold {report["settings"]["threshold"]:g} is illustrative. Test labels never fit or tune the model.</text>',
'<text x="54" y="827" font-size="15" fill="#5f6b73">The report includes confusion counts, overlapping error slices and the largest prediction errors.</text>',
'</g></svg>',
])
return "\n".join(parts) + "\n"
def main() -> None:
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--seed", type=int, default=DATA_SEED, help="Synthetic-data seed (default: 2026)")
parser.add_argument("--report", type=Path, help="Optional output JSON path; overwrites that file")
parser.add_argument("--chart", type=Path, help="Optional output SVG path; overwrites that file")
args = parser.parse_args()
report = run(args.seed)
serialized = json.dumps(report, indent=2, allow_nan=False) + "\n"
for path, content in ((args.report, serialized), (args.chart, chart_svg(report))):
if path:
path.parent.mkdir(parents=True, exist_ok=True)
path.write_text(content, encoding="utf-8", newline="\n")
print(serialized, end="")
if __name__ == "__main__":
main()