Demographic predictors of load profile mix
-We begin by aggregating load profile clusters at the block group level. If \(j\) is the block group, and \(q\) is the cluster (\(q\in\{1,\ldots,k\}\)), then
+We begin by aggregating load profile clusters at the block group level. If \(j\) is the block group, and \(q\) is the cluster (\(q\in\{0,1,2,3\}\)), then
\[
C_{jq} = \sum_{i\in \mathrm{bg}(j)} \mathrm{I}(c_i = q)
\]
is the count of household-day observations assigned to cluster \(q\) in block group \(j\), where \(c_i\) denotes the cluster assignment for household-day observation \(i\). Under this aggregation, a single household with multiple sampled days contributes multiple observations to the block group totals.
We further normalize these to \(\pi_{jq} = \frac{C_{jq}}{\sum_q C_{jq}}\), the proportion of cluster assignments in each block group.
-Our aim is to understand how block group level demographics predict the load cluster mix. Because \((\pi_{j1},\ldots,\pi_{jk})\) are compositional proportions that must sum to 1, we model log-ratios of cluster proportions. To prevent numerical instability from zero proportions, we apply Laplace smoothing (\(\alpha = 0.5\) pseudocount) before computing log-ratios. Specifically, we define smoothed proportions
+Our aim is to understand how block group level demographics predict the load cluster mix. Because \((\pi_{j0},\pi_{j1},\pi_{j2},\pi_{j3})\) are compositional proportions that must sum to 1, we use multinomial logistic regression to model the distribution of cluster assignments. This approach properly accounts for the compositional constraint, handles zero counts naturally through the multinomial likelihood, and provides statistically correct standard errors and hypothesis tests.
+
+Model specification
+We model the counts \((C_{j0}, C_{j1}, C_{j2}, C_{j3})\) as following a multinomial distribution:
\[
-\tilde{\pi}_{jq} = \frac{C_{jq} + \alpha}{\sum_{q=1}^{k} C_{jq} + k\alpha},
-\qquad \alpha = 0.5,
+(C_{j0}, C_{j1}, C_{j2}, C_{j3}) \sim \mathrm{Multinomial}\left(\sum_q C_{jq}; \pi_{j0}, \pi_{j1}, \pi_{j2}, \pi_{j3}\right)
\]
-which ensures \(\tilde{\pi}_{jq} > 0\) for all \(q\) and prevents \(\log(0) = -\infty\).
-We fit one model per cluster outcome (except the baseline), modeling the log-ratio of each cluster’s share relative to baseline cluster 1. Specifically, for each \(q\in\{2,\ldots,k\}\), we fit:
+We fit a single unified multinomial logit model with Cluster 2 as the baseline (reference) category. This cluster contains the largest share of household-day observations and exhibits the most stable and representative load-shape pattern among the clusters, providing a stable comparison point for the log-ratio regressions. The model estimates coefficients for Clusters 0, 1, and 3 relative to this baseline, parameterizing the log-odds as:
\[
-\log\left(\frac{\tilde{\pi}_{jq}}{\tilde{\pi}_{j1}}\right)
-= \beta_{0q}
-+ \beta_{1q} X_{1j}
-+ \beta_{2q} X_{2j}
-+ \cdots
-+ \varepsilon_{jq}.
+\log\left(\frac{\pi_{jq}}{\pi_{j2}}\right) = \beta_{0q} + \beta_{1q} X_{1j} + \beta_{2q} X_{2j} + \cdots + \beta_{pq} X_{pj}
\]
-We estimate these as separate weighted least squares (WLS) regressions, weighting each block group \(j\) by \(\mathrm{total\_obs}_j\), the total number of household-day observations available for that block group. This weighting gives greater influence to block groups with more observed data. As a robustness check, we also estimate unweighted OLS versions of the same log-ratio regressions; large differences between WLS and OLS estimates would indicate results are disproportionately driven by high-observation block groups.
-Coefficients are interpreted on the log-ratio scale: \(\beta_{pq}\) is the expected change in \(\log(\tilde{\pi}_{jq}/\tilde{\pi}_{j1})\) for a one-unit increase in predictor \(X_{pj}\), holding other predictors constant. Equivalently, \(\exp(\beta_{pq})\) is the multiplicative effect on the proportion ratio \(\tilde{\pi}_{jq}/\tilde{\pi}_{j1}\) for a one-unit increase in \(X_{pj}\).
+for \(q \in \{0, 1, 3\}\) (i.e., all non-baseline clusters).
+This parameterization ensures that:
+
+- All predicted probabilities are positive: \(\pi_{jq} > 0\)
+- Probabilities sum to 1: \(\sum_q \pi_{jq} = 1\)
+- The baseline cluster probability is determined by the constraint: \(\pi_{j2} = \frac{1}{1 + \sum_{q \in \{0,1,3\}} \exp(\beta_{0q} + \beta_{1q} X_{1j} + \cdots)}\)
+
+We estimate this model via maximum likelihood using the VGAM package in R, which implements iteratively reweighted least squares (IRLS) optimization. The multinomial likelihood naturally weights observations by their sample size—block groups with more household-day observations contribute more to the likelihood and thus have greater influence on parameter estimates.
+
+
+Why multinomial logit?
+An alternative approach would be to compute empirical log-ratios \(\log(C_{jq}/C_{j2})\) and model them using separate weighted ordinary least squares (OLS) regressions. However, this approach has several statistical problems:
+
+- Zero counts: Log-ratios are undefined when counts are zero, requiring arbitrary smoothing constants
+- Dependent errors: Log-ratios for different clusters share a common denominator and are thus correlated, violating OLS independence assumptions
+- Heteroskedasticity: The variance of log-ratios depends on unknown proportions, creating heteroskedasticity that weighted OLS only partially addresses
+- Invalid predictions: Predicted proportions may not sum to 1
+
+Multinomial logit avoids all these issues by using the correct likelihood for multinomial count data.
+
+
+Predictor variables and rank deficiency
+Our analysis includes 47 demographic predictor variables across five categories (spatial, economic, housing, household, and demographic characteristics). However, several of these variables form compositional groups where predictors sum to a constant (e.g., income brackets sum to 100% of households). This creates linear dependencies (rank deficiency) in the design matrix.
+To handle this, we:
+
+- Detect rank deficiency via QR decomposition of the predictor matrix
+- Drop redundant predictors that are exact linear combinations of others
+- In practice, this reduces our predictor set from 47 to approximately 42 active predictors (plus an intercept term)
+
+The final model includes 43 terms per cluster equation (42 predictors + 1 intercept), yielding 129 total coefficients across the 3 non-baseline cluster equations.
+
+
+
+Interpretation
+Coefficients are interpreted on the log-odds scale: \(\beta_{pq}\) is the expected change in \(\log(\pi_{jq}/\pi_{j2})\) for a one-unit increase in predictor \(X_{pj}\), holding other predictors constant.
+Equivalently, \(\exp(\beta_{pq})\) is the multiplicative effect on the odds ratio \(\pi_{jq}/\pi_{j2}\) for a one-unit increase in \(X_{pj}\). For example:
+
+- If \(\beta_{pq} = 0.02\) for median household income (in thousands of dollars), then a $10k increase in income multiplies the odds of cluster \(q\) vs. baseline by \(\exp(0.02 \times 10) = e^{0.2} \approx 1.22\) (a 22% increase in odds).
+
+To obtain predicted probabilities for a block group with covariates \(X_j\), we use the softmax transformation:
+\[
+\hat{\pi}_{jq} = \begin{cases}
+\frac{\exp(\hat{\beta}_{0q} + \hat{\beta}_{1q} X_{1j} + \cdots)}
+{1 + \sum_{m \in \{0,1,3\}} \exp(\hat{\beta}_{0m} + \hat{\beta}_{1m} X_{1j} + \cdots)}
+& \text{if } q \in \{0,1,3\} \\[1em]
+\frac{1}
+{1 + \sum_{m \in \{0,1,3\}} \exp(\hat{\beta}_{0m} + \hat{\beta}_{1m} X_{1j} + \cdots)}
+& \text{if } q = 2 \text{ (baseline)}
+\end{cases}
+\]
+These predicted probabilities automatically sum to 1 and are guaranteed to be positive.
+
+
diff --git a/docs/index.qmd b/docs/index.qmd
index 2c3b497..b47a8ea 100644
--- a/docs/index.qmd
+++ b/docs/index.qmd
@@ -60,7 +60,8 @@ We selected cluster 1 as the baseline (reference category) because it contains t
## Demographic predictors of load profile mix
-We begin by aggregating load profile clusters at the block group level. If $j$ is the block group, and $q$ is the cluster ($q\in\{1,\ldots,k\}$), then
+
+We begin by aggregating load profile clusters at the block group level. If $j$ is the block group, and $q$ is the cluster ($q\in\{0,1,2,3\}$), then
$$
C_{jq} = \sum_{i\in \mathrm{bg}(j)} \mathrm{I}(c_i = q)
@@ -70,26 +71,82 @@ is the count of household-day observations assigned to cluster $q$ in block grou
We further normalize these to $\pi_{jq} = \frac{C_{jq}}{\sum_q C_{jq}}$, the proportion of cluster assignments in each block group.
-Our aim is to understand how block group level demographics predict the load cluster mix. Because $(\pi_{j1},\ldots,\pi_{jk})$ are compositional proportions that must sum to 1, we model log-ratios of cluster proportions. To prevent numerical instability from zero proportions, we apply Laplace smoothing ($\alpha = 0.5$ pseudocount) before computing log-ratios. Specifically, we define smoothed proportions
+Our aim is to understand how block group level demographics predict the load cluster mix. Because $(\pi_{j0},\pi_{j1},\pi_{j2},\pi_{j3})$ are compositional proportions that must sum to 1, we use **multinomial logistic regression** to model the distribution of cluster assignments. This approach properly accounts for the compositional constraint, handles zero counts naturally through the multinomial likelihood, and provides statistically correct standard errors and hypothesis tests.
+
+### Model specification
+
+We model the counts $(C_{j0}, C_{j1}, C_{j2}, C_{j3})$ as following a multinomial distribution:
+
+$$
+(C_{j0}, C_{j1}, C_{j2}, C_{j3}) \sim \mathrm{Multinomial}\left(\sum_q C_{jq}; \pi_{j0}, \pi_{j1}, \pi_{j2}, \pi_{j3}\right)
+$$
+
+We fit a single unified multinomial logit model with **Cluster 2 as the baseline (reference) category**. This cluster contains the largest share of household-day observations and exhibits the most stable and representative load-shape pattern among the clusters, providing a stable comparison point for the log-ratio regressions. The model estimates coefficients for Clusters 0, 1, and 3 relative to this baseline, parameterizing the log-odds as:
$$
-\tilde{\pi}_{jq} = \frac{C_{jq} + \alpha}{\sum_{q=1}^{k} C_{jq} + k\alpha},
-\qquad \alpha = 0.5,
+\log\left(\frac{\pi_{jq}}{\pi_{j2}}\right) = \beta_{0q} + \beta_{1q} X_{1j} + \beta_{2q} X_{2j} + \cdots + \beta_{pq} X_{pj}
$$
-which ensures $\tilde{\pi}_{jq} > 0$ for all $q$ and prevents $\log(0) = -\infty$.
+for $q \in \{0, 1, 3\}$ (i.e., all non-baseline clusters).
+
+This parameterization ensures that:
+
+- All predicted probabilities are positive: $\pi_{jq} > 0$
+- Probabilities sum to 1: $\sum_q \pi_{jq} = 1$
+- The baseline cluster probability is determined by the constraint: $\pi_{j2} = \frac{1}{1 + \sum_{q \in \{0,1,3\}} \exp(\beta_{0q} + \beta_{1q} X_{1j} + \cdots)}$
+
+We estimate this model via **maximum likelihood** using the VGAM package in R, which implements iteratively reweighted least squares (IRLS) optimization. The multinomial likelihood naturally weights observations by their sample size—block groups with more household-day observations contribute more to the likelihood and thus have greater influence on parameter estimates.
+
+### Why multinomial logit?
+
+An alternative approach would be to compute empirical log-ratios $\log(C_{jq}/C_{j2})$ and model them using separate weighted ordinary least squares (OLS) regressions. However, this approach has several statistical problems:
+
+1. **Zero counts:** Log-ratios are undefined when counts are zero, requiring arbitrary smoothing constants
+2. **Dependent errors:** Log-ratios for different clusters share a common denominator and are thus correlated, violating OLS independence assumptions
+3. **Heteroskedasticity:** The variance of log-ratios depends on unknown proportions, creating heteroskedasticity that weighted OLS only partially addresses
+4. **Invalid predictions:** Predicted proportions may not sum to 1
+
+Multinomial logit avoids all these issues by using the correct likelihood for multinomial count data.
-We fit one model per cluster outcome (except the baseline), modeling the log-ratio of each cluster’s share relative to baseline cluster 1. Specifically, for each $q\in\{2,\ldots,k\}$, we fit:
+### Predictor variables and rank deficiency
+
+Our analysis includes 47 demographic predictor variables across five categories (spatial, economic, housing, household, and demographic characteristics). However, several of these variables form **compositional groups** where predictors sum to a constant (e.g., income brackets sum to 100% of households). This creates linear dependencies (rank deficiency) in the design matrix.
+
+To handle this, we:
+
+1. **Detect rank deficiency** via QR decomposition of the predictor matrix
+2. **Drop redundant predictors** that are exact linear combinations of others
+3. In practice, this reduces our predictor set from 47 to approximately **42 active predictors** (plus an intercept term)
+
+The final model includes **43 terms per cluster equation** (42 predictors + 1 intercept), yielding **129 total coefficients** across the 3 non-baseline cluster equations.
+
+::: {.callout-note}
+## Variable Selection is still in progress.
+:::
+
+### Interpretation
+
+Coefficients are interpreted on the **log-odds scale**: $\beta_{pq}$ is the expected change in $\log(\pi_{jq}/\pi_{j2})$ for a one-unit increase in predictor $X_{pj}$, holding other predictors constant.
+
+Equivalently, $\exp(\beta_{pq})$ is the **multiplicative effect** on the odds ratio $\pi_{jq}/\pi_{j2}$ for a one-unit increase in $X_{pj}$. For example:
+
+- If $\beta_{pq} = 0.02$ for median household income (in thousands of dollars), then a \$10k increase in income multiplies the odds of cluster $q$ vs. baseline by $\exp(0.02 \times 10) = e^{0.2} \approx 1.22$ (a 22% increase in odds).
+
+To obtain **predicted probabilities** for a block group with covariates $X_j$, we use the softmax transformation:
$$
-\log\left(\frac{\tilde{\pi}_{jq}}{\tilde{\pi}_{j1}}\right)
-= \beta_{0q}
-+ \beta_{1q} X_{1j}
-+ \beta_{2q} X_{2j}
-+ \cdots
-+ \varepsilon_{jq}.
+\hat{\pi}_{jq} = \begin{cases}
+\frac{\exp(\hat{\beta}_{0q} + \hat{\beta}_{1q} X_{1j} + \cdots)}
+{1 + \sum_{m \in \{0,1,3\}} \exp(\hat{\beta}_{0m} + \hat{\beta}_{1m} X_{1j} + \cdots)}
+& \text{if } q \in \{0,1,3\} \\[1em]
+\frac{1}
+{1 + \sum_{m \in \{0,1,3\}} \exp(\hat{\beta}_{0m} + \hat{\beta}_{1m} X_{1j} + \cdots)}
+& \text{if } q = 2 \text{ (baseline)}
+\end{cases}
$$
-We estimate these as separate weighted least squares (WLS) regressions, weighting each block group $j$ by $\mathrm{total\_obs}_j$, the total number of household-day observations available for that block group. This weighting gives greater influence to block groups with more observed data. As a robustness check, we also estimate unweighted OLS versions of the same log-ratio regressions; large differences between WLS and OLS estimates would indicate results are disproportionately driven by high-observation block groups.
+These predicted probabilities automatically sum to 1 and are guaranteed to be positive.
+
+### Model performance
-Coefficients are interpreted on the log-ratio scale: $\beta_{pq}$ is the expected change in $\log(\tilde{\pi}_{jq}/\tilde{\pi}_{j1})$ for a one-unit increase in predictor $X_{pj}$, holding other predictors constant. Equivalently, $\exp(\beta_{pq})$ is the multiplicative effect on the proportion ratio $\tilde{\pi}_{jq}/\tilde{\pi}_{j1}$ for a one-unit increase in $X_{pj}$.
+Our final model achieves a **pseudo-R² of 0.567** (McFadden's R²), indicating that demographic predictors explain approximately 57% of the variation in cluster distributions across block groups. The model successfully converged using VGAM's IRLS optimizer, and all coefficients have valid standard errors for hypothesis testing.
diff --git a/pyproject.toml b/pyproject.toml
index c7297bb..746b468 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -20,6 +20,7 @@ dependencies = [
"boto3>=1.40.46",
"botocore>=1.40.47",
"cenpy>=1.0.1",
+ "git-filter-repo>=2.47.0",
"ipykernel>=6.30.1",
"matplotlib>=3.9.4",
"memory-profiler>=0.61.0",
@@ -59,7 +60,7 @@ DEP002 = [
"snakeviz",
"tslearn",
]
-DEP003 = ["botocore", "analysis", "smart_meter_analysis", "pandas", "scipy"]
+DEP003 = ["botocore", "analysis", "smart_meter_analysis", "pandas", "scipy", "psutil"]
DEP004 = ["botocore"]
[dependency-groups]
diff --git a/scripts/run_comed_pipeline.py b/scripts/run_comed_pipeline.py
index 9501ef7..7d7cc25 100644
--- a/scripts/run_comed_pipeline.py
+++ b/scripts/run_comed_pipeline.py
@@ -1,16 +1,36 @@
#!/usr/bin/env python3
-"""ComEd Smart Meter Analysis Pipeline Orchestrator (Phase 3)
-
-Directory structure (required):
-data/runs/{run_name}/
-├── raw/ # Downloaded CSVs
-├── processed/ # Canonical interval parquet
-├── clustering/ # Stage 1 outputs + clustering outputs
-├── stage2/ # Stage 2 outputs
-├── logs/ # All logs
-├── run_manifest.json # Pipeline metadata
-├── download_manifest.jsonl # S3 download record
-└── processing_manifest.jsonl # CSV processing record
+"""
+ComEd Smart Meter Analysis Pipeline Orchestrator (Production-Grade)
+
+This orchestrator integrates the existing repo components to run:
+1) S3 download (optional; supports hybrid loop)
+2) CSV -> canonical parquet ingest (hybrid-safe parts promotion)
+3) Prepare clustering sample/profiles
+4) OPTIONAL: K-search (must run BEFORE single-k clustering)
+5) Single-k clustering
+6) OPTIONAL: Stage 2 regression (R multinomial blockgroup model)
+
+Non-negotiables implemented here
+--------------------------------
+- All imports are at module top (no imports inside functions).
+- Deterministic ordering (sorted inputs; fixed seeds passed through).
+- Resumable + safe against partial failures (hybrid promotion, recursive raw detection, skip flags, idempotent checks).
+- Logs under data/runs/{run_name}/logs/.
+- Never assumes "no files" from a single glob; always checks recursively.
+
+Stage 2 (Updated)
+-----------------
+Stage 2 regression is executed via:
+analysis/stage2/stage2_multinom_blockgroup_weighted.R
+
+Inputs:
+- clusters parquet: data/runs/{run_name}/clustering/cluster_assignments.parquet
+- crosswalk: data/reference/2023_comed_zip4_census_crosswalk.txt (default, overrideable)
+- census predictors parquet: data/runs/{run_name}/stage2/census_predictors_blockgroups.parquet
+ (generated by this orchestrator from census_cache + census_specs.py VARIABLE_SPECS)
+
+New capability:
+- --stage2-build-census-only: regenerate census cache + predictors parquet without running Stage 2 regression
"""
from __future__ import annotations
@@ -19,15 +39,29 @@
import json
import logging
import os
+import re
+import shutil
import subprocess
import sys
+from dataclasses import dataclass
from datetime import datetime, timezone
+from inspect import signature
from pathlib import Path
from typing import Any
import polars as pl
+import pyarrow.parquet as pq
+
+try:
+ import psutil # type: ignore[import-untyped]
+except Exception:
+ psutil = None # type: ignore[assignment]
+
+import contextlib
from smart_meter_analysis.aws_loader import download_s3_batch, list_s3_files
+from smart_meter_analysis.census import fetch_census_data
+from smart_meter_analysis.census_specs import VARIABLE_SPECS
logger = logging.getLogger(__name__)
@@ -36,25 +70,67 @@
DEFAULT_STATE_FIPS = "17"
DEFAULT_ACS_YEAR = 2023
+DEFAULT_PARTS_DIRNAME = "parts"
+DEFAULT_MIN_FREE_GB = 25.0
+DEFAULT_DOWNLOAD_CHUNK_SIZE = 50_000
+
+DEFAULT_SAMPLE_DAYS = 31
+DEFAULT_SEED = 42
+
+DEFAULT_K_MIN = 3
+DEFAULT_K_MAX = 8
+DEFAULT_K_SEARCH_DIRNAME = "k_search"
+DEFAULT_K_SEARCH_MAX_PROFILES = 500_000
+
+BATCH_BASENAME_RE = re.compile(r"^batch_\d+\.parquet$")
+
+DEFAULT_STAGE2_R_OUTDIRNAME = "stage2_R_multi"
+DEFAULT_STAGE2_CENSUS_PREDICTORS_BASENAME = "census_predictors_blockgroups.parquet"
+
+
+@dataclass(frozen=True)
+class StepPaths:
+ run_dir: Path
+ raw_dir: Path
+ processed_dir: Path
+ clustering_dir: Path
+ stage2_dir: Path
+ stage2_r_dir: Path
+ logs_dir: Path
+ run_manifest: Path
+ run_state: Path
+ download_manifest: Path
+ processing_manifest: Path
+
def _utc_now_iso() -> str:
return datetime.now(timezone.utc).replace(microsecond=0).isoformat()
-def _configure_logging(log_path: Path) -> None:
- log_path.parent.mkdir(parents=True, exist_ok=True)
+def _atomic_write_text(path: Path, text: str) -> None:
+ path.parent.mkdir(parents=True, exist_ok=True)
+ tmp = path.with_suffix(path.suffix + ".tmp")
+ tmp.write_text(text, encoding="utf-8")
+ tmp.replace(path)
+
+
+def _configure_logging(logs_dir: Path) -> Path:
+ logs_dir.mkdir(parents=True, exist_ok=True)
+ stamp = datetime.now().strftime("%Y%m%d_%H%M%S")
+ log_path = logs_dir / f"pipeline_{stamp}.log"
handlers: list[logging.Handler] = [
logging.StreamHandler(sys.stdout),
logging.FileHandler(str(log_path), encoding="utf-8"),
]
-
logging.basicConfig(
level=logging.INFO,
- format="%(asctime)s - %(levelname)s - %(message)s",
+ format="%(asctime)s %(levelname)s %(name)s %(message)s",
handlers=handlers,
force=True,
)
+ logger.info("Logging to %s", log_path)
+ return log_path
def _get_git_sha(repo_root: Path) -> str | None:
@@ -74,102 +150,349 @@ def _get_git_sha(repo_root: Path) -> str | None:
return None
-def create_directories(run_dir: Path) -> dict[str, Path]:
- paths = {
- "run_dir": run_dir,
- "raw": run_dir / "raw",
- "processed": run_dir / "processed",
- "clustering": run_dir / "clustering",
- "stage2": run_dir / "stage2",
- "logs": run_dir / "logs",
- }
- for p in paths.values():
+def _disk_free_gb(path: Path) -> float:
+ usage = shutil.disk_usage(str(path))
+ return float(usage.free) / (1024.0**3)
+
+
+def _log_system_diagnostics(*, base_path: Path) -> None:
+ free_gb = _disk_free_gb(base_path)
+ msg: dict[str, Any] = {"event": "system_diagnostics", "disk_free_gb": round(free_gb, 2)}
+ if psutil is not None:
+ try:
+ vm = psutil.virtual_memory()
+ msg["mem_total_gb"] = round(float(vm.total) / (1024.0**3), 2)
+ msg["mem_available_gb"] = round(float(vm.available) / (1024.0**3), 2)
+ msg["mem_percent"] = float(vm.percent)
+ except Exception: # noqa: S110
+ pass
+ logger.info("DIAG %s", json.dumps(msg, sort_keys=True))
+
+
+def _ensure_paths(run_dir: Path) -> StepPaths:
+ raw_dir = run_dir / "raw"
+ processed_dir = run_dir / "processed"
+ clustering_dir = run_dir / "clustering"
+ stage2_dir = run_dir / "stage2"
+ stage2_r_dir = run_dir / DEFAULT_STAGE2_R_OUTDIRNAME
+ logs_dir = run_dir / "logs"
+
+ for p in [run_dir, raw_dir, processed_dir, clustering_dir, stage2_dir, stage2_r_dir, logs_dir]:
p.mkdir(parents=True, exist_ok=True)
- return paths
+
+ return StepPaths(
+ run_dir=run_dir,
+ raw_dir=raw_dir,
+ processed_dir=processed_dir,
+ clustering_dir=clustering_dir,
+ stage2_dir=stage2_dir,
+ stage2_r_dir=stage2_r_dir,
+ logs_dir=logs_dir,
+ run_manifest=run_dir / "run_manifest.json",
+ run_state=run_dir / "run_state.json",
+ download_manifest=run_dir / "download_manifest.jsonl",
+ processing_manifest=run_dir / "processing_manifest.jsonl",
+ )
+
+
+def _validate_crosswalk(path: Path) -> Path:
+ if not path.exists():
+ raise FileNotFoundError(
+ f"ZIP+4 crosswalk not found at {path}. Provide --stage2-crosswalk or ensure the file exists."
+ )
+ return path
-def _run_subprocess(cmd: list[str], *, cwd: Path | None = None, env: dict[str, str] | None = None) -> None:
+def _run_subprocess(cmd: list[str], *, cwd: Path | None = None) -> None:
logger.info("Command: %s", " ".join(cmd))
- r = subprocess.run(cmd, cwd=str(cwd) if cwd else None, env=env, check=False)
+ r = subprocess.run(cmd, cwd=str(cwd) if cwd else None, check=False)
if r.returncode != 0:
- raise RuntimeError(f"Command failed with return code {r.returncode}: {' '.join(cmd)}")
+ raise RuntimeError(f"Command failed (rc={r.returncode}): {' '.join(cmd)}")
-def download_from_s3(*, year_month: str, num_files: int, run_dir: Path) -> dict[str, Any]:
- raw_dir = run_dir / "raw"
- raw_dir.mkdir(parents=True, exist_ok=True)
+def _raw_csvs_recursive(raw_dir: Path) -> list[Path]:
+ if not raw_dir.exists():
+ return []
+ csvs = [p for p in raw_dir.rglob("*.csv") if p.is_file()]
+ csvs.extend([p for p in raw_dir.rglob("*.CSV") if p.is_file()])
+ return sorted(set(csvs))
- manifest_path = run_dir / "download_manifest.jsonl"
- # list_s3_files MUST return keys (not URIs).
- s3_keys = list_s3_files(year_month=year_month, max_files=num_files)
+def _raw_processing_dirs(raw_dir: Path) -> list[Path]:
+ """
+ Return directories that contain top-level CSVs (non-recursive within each returned dir),
+ in deterministic (sorted) order.
+ """
+ if not raw_dir.exists():
+ return []
- logger.info("Downloading %d files for %s into %s", len(s3_keys), year_month, raw_dir)
- result = download_s3_batch(
- s3_keys=s3_keys,
- output_dir=raw_dir,
- manifest_path=manifest_path,
- max_files=num_files,
- fail_fast=True,
- max_errors=10,
- retries=3,
- backoff_factor=2.0,
- )
- logger.info(
- "Download complete: downloaded=%s failed=%s manifest=%s",
- result.get("downloaded"),
- result.get("failed"),
- result.get("manifest_path"),
- )
- return result
+ dirs: set[Path] = set()
+ top = list(raw_dir.glob("*.csv")) + list(raw_dir.glob("*.CSV"))
+ if any(p.is_file() for p in top):
+ dirs.add(raw_dir)
-def process_csvs(*, run_dir: Path, year_month: str, day_mode: str = "calendar") -> Path:
- raw_dir = run_dir / "raw"
- processed_dir = run_dir / "processed"
- processed_dir.mkdir(parents=True, exist_ok=True)
+ for d in [p for p in raw_dir.rglob("*") if p.is_dir()]:
+ files = list(d.glob("*.csv")) + list(d.glob("*.CSV"))
+ if any(p.is_file() for p in files):
+ dirs.add(d)
+
+ return sorted(dirs)
+
+
+def _delete_csvs_under(dir_path: Path) -> None:
+ if not dir_path.exists():
+ return
+ for p in list(dir_path.glob("*.csv")) + list(dir_path.glob("*.CSV")):
+ if p.is_file():
+ try:
+ p.unlink(missing_ok=True)
+ except Exception as exc:
+ logger.warning("Failed to delete raw CSV %s: %s", p, exc)
+
+
+def _prune_empty_dirs(root: Path) -> None:
+ if not root.exists():
+ return
+ for d in sorted([p for p in root.rglob("*") if p.is_dir()], reverse=True):
+ with contextlib.suppress(OSError):
+ d.rmdir()
+
+
+def _parts_root(processed_dir: Path, parts_dirname: str) -> Path:
+ return processed_dir / parts_dirname
+
+
+def _existing_cycle_ids(parts_root: Path) -> list[int]:
+ if not parts_root.exists():
+ return []
+ ids: list[int] = []
+ for d in parts_root.glob("cycle_*"):
+ if not d.is_dir():
+ continue
+ m = re.match(r"^cycle_(\d{4})$", d.name)
+ if not m:
+ continue
+ ids.append(int(m.group(1)))
+ return sorted(ids)
+
+
+def _next_cycle_id(parts_root: Path) -> int:
+ ids = _existing_cycle_ids(parts_root)
+ return (max(ids) + 1) if ids else 0
+
+
+def _collect_part_parquets(parts_root: Path) -> list[Path]:
+ if not parts_root.exists():
+ return []
+ return sorted([p for p in parts_root.rglob("*.parquet") if p.is_file()])
+
+
+def _persist_temp_batches_as_parts(
+ *,
+ processed_dir: Path,
+ parts_dirname: str,
+ cycle_id: int,
+) -> list[Path]:
+ """
+ Promote processed/temp_batches/batch_####.parquet -> processed/{parts_dirname}/cycle_####/batch_####.parquet
+ """
+ temp_dir = processed_dir / "temp_batches"
+ if not temp_dir.exists():
+ return []
+
+ candidates = sorted([p for p in temp_dir.glob("*.parquet") if p.is_file()])
+ batch_files = [p for p in candidates if BATCH_BASENAME_RE.match(p.name)]
+ sub_files = [p for p in candidates if p.name.startswith("batch_") and ("_sub_" in p.name)]
+
+ if not batch_files and not sub_files:
+ return []
+
+ parts_root = _parts_root(processed_dir, parts_dirname)
+ parts_root.mkdir(parents=True, exist_ok=True)
+ cycle_dir = parts_root / f"cycle_{cycle_id:04d}"
+ cycle_dir.mkdir(parents=True, exist_ok=True)
+
+ moved: list[Path] = []
+ for src in batch_files:
+ dst = cycle_dir / src.name
+ if dst.exists():
+ raise FileExistsError(f"Refusing to overwrite existing parts file: {dst}")
+ src.replace(dst)
+ moved.append(dst)
+
+ for s in sub_files:
+ with contextlib.suppress(Exception):
+ s.unlink(missing_ok=True)
+
+ with contextlib.suppress(Exception):
+ shutil.rmtree(temp_dir, ignore_errors=True)
+
+ logger.info("Promoted %d durable batch parquet(s) into %s", len(moved), cycle_dir)
+ return moved
+
+
+def _download_manifest_supported_kwargs() -> dict[str, bool]:
+ try:
+ params = set(signature(download_s3_batch).parameters.keys())
+ except Exception:
+ params = set()
+ return {
+ "overwrite_manifest": "overwrite_manifest" in params,
+ "skip_existing": "skip_existing" in params,
+ "log_every": "log_every" in params,
+ }
+
+
+def _s3_key_basename(s3_key: str) -> str:
+ return os.path.basename(s3_key)
+
+
+def _append_file(src: Path, dst: Path) -> None:
+ dst.parent.mkdir(parents=True, exist_ok=True)
+ if not src.exists():
+ return
+ with src.open("r", encoding="utf-8") as r, dst.open("a", encoding="utf-8") as w:
+ for line in r:
+ w.write(line)
+
+
+def _download_s3_chunk(
+ *,
+ s3_keys: list[str],
+ output_dir: Path,
+ manifest_path: Path,
+ fail_fast: bool,
+ max_errors: int,
+ retries: int,
+ backoff_factor: float,
+ log_every: int,
+ chunk_id: int,
+) -> dict[str, Any]:
+ support = _download_manifest_supported_kwargs()
+
+ output_dir.mkdir(parents=True, exist_ok=True)
+ manifest_path.parent.mkdir(parents=True, exist_ok=True)
+
+ kwargs: dict[str, Any] = {
+ "s3_keys": s3_keys,
+ "output_dir": output_dir,
+ "manifest_path": manifest_path,
+ "max_files": None,
+ "fail_fast": bool(fail_fast),
+ "max_errors": int(max_errors),
+ "retries": int(retries),
+ "backoff_factor": float(backoff_factor),
+ }
+
+ if support["log_every"]:
+ kwargs["log_every"] = int(log_every)
+
+ if support["overwrite_manifest"] and support["skip_existing"]:
+ kwargs["overwrite_manifest"] = False
+ kwargs["skip_existing"] = True
+ result = download_s3_batch(**kwargs)
+ return result if isinstance(result, dict) else {"result": result}
+
+ existing = {p.name for p in output_dir.glob("*") if p.is_file()}
+ filtered = [k for k in s3_keys if _s3_key_basename(k) not in existing]
+
+ if not filtered:
+ logger.info("Download chunk %d: all files already exist locally (fallback skip).", chunk_id)
+ return {"downloaded": 0, "skipped": len(s3_keys), "failed": 0, "manifest_path": str(manifest_path)}
+
+ tmp_manifest = manifest_path.parent / f"download_manifest_chunk_{chunk_id:06d}.jsonl"
+ kwargs["s3_keys"] = filtered
+ kwargs["manifest_path"] = tmp_manifest
+
+ result = download_s3_batch(**kwargs)
+ _append_file(tmp_manifest, manifest_path)
+ with contextlib.suppress(Exception):
+ tmp_manifest.unlink(missing_ok=True)
+
+ return result if isinstance(result, dict) else {"result": result}
+
+
+def _safe_list_s3_keys(*, year_month: str, max_files: int | None) -> list[str]:
+ try:
+ keys = list_s3_files(year_month=year_month, max_files=max_files)
+ except TypeError:
+ keys = list_s3_files(year_month=year_month, max_files=max_files if max_files is not None else 0)
+ if keys is None:
+ return []
+ return list(keys)
+
+
+def _load_run_state(path: Path) -> dict[str, Any]:
+ if not path.exists():
+ return {}
+ try:
+ return json.loads(path.read_text(encoding="utf-8"))
+ except Exception:
+ return {}
- output_file = processed_dir / f"comed_{year_month}.parquet"
- manifest_file = run_dir / "processing_manifest.jsonl"
+
+def _update_run_state(path: Path, patch: dict[str, Any]) -> None:
+ state = _load_run_state(path)
+ state.setdefault("created_utc", _utc_now_iso())
+ state["updated_utc"] = _utc_now_iso()
+ for k, v in patch.items():
+ state[k] = v
+ _atomic_write_text(path, json.dumps(state, indent=2, sort_keys=True) + "\n")
+
+
+def process_csvs(
+ *,
+ input_dir: Path,
+ processed_dir: Path,
+ run_dir: Path,
+ year_month: str,
+ day_mode: str,
+ keep_temp: bool,
+ output_path: Path,
+ processing_manifest: Path,
+) -> None:
+ input_dir.mkdir(parents=True, exist_ok=True)
+ processed_dir.mkdir(parents=True, exist_ok=True)
+ processing_manifest.parent.mkdir(parents=True, exist_ok=True)
cmd = [
sys.executable,
"scripts/process_csvs_batched_optimized.py",
"--input-dir",
- str(raw_dir),
+ str(input_dir),
"--output",
- str(output_file),
+ str(output_path),
"--processing-manifest",
- str(manifest_file),
+ str(processing_manifest),
"--day-mode",
str(day_mode),
]
+ if keep_temp:
+ cmd.append("--keep-temp")
- logger.info("Ingesting CSVs into canonical interval parquet: %s", output_file)
+ logger.info("Ingesting CSVs from %s -> %s", input_dir, output_path)
_run_subprocess(cmd)
- if not output_file.exists():
- raise FileNotFoundError(f"Expected processed parquet missing: {output_file}")
-
- return output_file
-
def prepare_clustering(
*,
- run_dir: Path,
- processed_parquet: Path,
+ clustering_dir: Path,
+ processed_inputs: list[Path],
sample_days: int,
sample_households: int | None,
seed: int,
) -> Path:
- clustering_dir = run_dir / "clustering"
clustering_dir.mkdir(parents=True, exist_ok=True)
+ if not processed_inputs:
+ raise ValueError("No processed parquet inputs provided for clustering preparation.")
+
cmd = [
sys.executable,
"analysis/clustering/prepare_clustering_data_households.py",
"--input",
- str(processed_parquet),
+ *[str(p) for p in processed_inputs],
"--output-dir",
str(clustering_dir),
"--sample-days",
@@ -178,11 +501,10 @@ def prepare_clustering(
str(int(seed)),
"--streaming",
]
-
if sample_households is not None:
cmd.extend(["--sample-households", str(int(sample_households))])
- logger.info("Preparing Stage 1 clustering data in %s", clustering_dir)
+ logger.info("Preparing clustering data (inputs=%d) -> %s", len(processed_inputs), clustering_dir)
_run_subprocess(cmd)
profiles = clustering_dir / "sampled_profiles.parquet"
@@ -191,16 +513,163 @@ def prepare_clustering(
return profiles
+def _script_help_text(script: Path) -> str:
+ r = subprocess.run([sys.executable, str(script), "--help"], capture_output=True, text=True, check=False)
+ return (r.stdout or "") + "\n" + (r.stderr or "")
+
+
+def _help_supports(help_text: str, flag: str) -> bool:
+ return flag in help_text
+
+
+def run_k_search(
+ *,
+ profiles_path: Path,
+ output_dir: Path,
+ seed: int,
+ k_values: list[int],
+ sample_households: int | None,
+ max_profiles: int | None,
+) -> int:
+ script = Path("analysis/clustering/euclidean_clustering_k_search.py")
+ if not script.exists():
+ raise FileNotFoundError(f"Missing k-search script: {script}")
+
+ output_dir.mkdir(parents=True, exist_ok=True)
+ help_text = _script_help_text(script)
+
+ cmd: list[str] = [
+ sys.executable,
+ str(script),
+ "--input",
+ str(profiles_path),
+ "--output-dir",
+ str(output_dir),
+ ]
+
+ if _help_supports(help_text, "--random-state"):
+ cmd.extend(["--random-state", str(int(seed))])
+ elif _help_supports(help_text, "--seed"):
+ cmd.extend(["--seed", str(int(seed))])
+
+ if _help_supports(help_text, "--k-values"):
+ cmd.extend(["--k-values", ",".join(str(k) for k in k_values)])
+ elif _help_supports(help_text, "--k-min") and _help_supports(help_text, "--k-max"):
+ cmd.extend(["--k-min", str(min(k_values)), "--k-max", str(max(k_values))])
+ else:
+ if _help_supports(help_text, "--k"):
+ for k in k_values:
+ cmd.extend(["--k", str(int(k))])
+ else:
+ raise RuntimeError(
+ "k-search script does not expose --k-values, (--k-min/--k-max), or --k; cannot drive it safely."
+ )
+
+ if sample_households is not None:
+ if _help_supports(help_text, "--sample-households"):
+ cmd.extend(["--sample-households", str(int(sample_households))])
+ elif _help_supports(help_text, "--k-search-sample-households"):
+ cmd.extend(["--k-search-sample-households", str(int(sample_households))])
+
+ if max_profiles is not None:
+ if _help_supports(help_text, "--max-profiles"):
+ cmd.extend(["--max-profiles", str(int(max_profiles))])
+ elif _help_supports(help_text, "--k-search-max-profiles"):
+ cmd.extend(["--k-search-max-profiles", str(int(max_profiles))])
+
+ logger.info("Running k-search -> %s", output_dir)
+ _run_subprocess(cmd)
+
+ return read_recommended_k(output_dir=output_dir, k_values=k_values)
+
+
+def read_recommended_k(*, output_dir: Path, k_values: list[int]) -> int:
+ explicit_candidates = [
+ output_dir / "recommended_k.json",
+ output_dir / "recommended_k.txt",
+ output_dir / "best_k.json",
+ output_dir / "best_k.txt",
+ output_dir / "k_recommended.json",
+ output_dir / "k_recommended.txt",
+ ]
+ for p in explicit_candidates:
+ if not p.exists():
+ continue
+ try:
+ if p.suffix.lower() == ".txt":
+ k = int(p.read_text(encoding="utf-8").strip())
+ logger.info("k-search recommended k from %s: %d", p.name, k)
+ return k
+ obj = json.loads(p.read_text(encoding="utf-8"))
+ for key in ["recommended_k", "best_k", "k", "recommended"]:
+ if key in obj:
+ k = int(obj[key])
+ logger.info("k-search recommended k from %s[%s]: %d", p.name, key, k)
+ return k
+ except Exception: # noqa: S110
+ pass
+
+ for p in sorted(output_dir.glob("*.json")):
+ try:
+ obj = json.loads(p.read_text(encoding="utf-8"))
+ for key in ["recommended_k", "best_k"]:
+ if key in obj:
+ k = int(obj[key])
+ logger.info("k-search recommended k from %s[%s]: %d", p.name, key, k)
+ return k
+ except Exception: # noqa: S112
+ continue
+
+ table_candidates = (
+ list(output_dir.glob("*.parquet")) + list(output_dir.glob("*.csv")) + list(output_dir.glob("*.tsv"))
+ )
+ for p in sorted(table_candidates):
+ try:
+ if p.suffix.lower() == ".parquet":
+ df = pl.read_parquet(p)
+ elif p.suffix.lower() == ".csv":
+ df = pl.read_csv(p)
+ elif p.suffix.lower() == ".tsv":
+ df = pl.read_csv(p, separator="\t")
+ else:
+ continue
+ except Exception: # noqa: S112
+ continue
+
+ cols = set(df.columns)
+ if "k" not in cols:
+ continue
+
+ silhouette_col = None
+ for c in ["silhouette", "silhouette_score", "mean_silhouette"]:
+ if c in cols:
+ silhouette_col = c
+ break
+
+ if silhouette_col is not None:
+ df2 = df.select(["k", silhouette_col]).drop_nulls()
+ if df2.height == 0:
+ continue
+ best = df2.sort([silhouette_col, "k"], descending=[True, False]).row(0)
+ k = int(best[0])
+ logger.info("k-search recommended k from %s by max %s: %d", p.name, silhouette_col, k)
+ return k
+
+ k = int(sorted(k_values)[0])
+ logger.warning("k-search outputs did not expose a recommended k; falling back deterministically to k=%d", k)
+ return k
+
+
def run_clustering(
*,
- run_dir: Path,
+ clustering_dir: Path,
profiles_path: Path,
k: int,
clustering_seed: int,
) -> None:
- clustering_dir = run_dir / "clustering"
clustering_dir.mkdir(parents=True, exist_ok=True)
+ # NOTE: this is where clustering standardization/normalization is applied.
cmd = [
sys.executable,
"analysis/clustering/euclidean_clustering_minibatch.py",
@@ -210,14 +679,14 @@ def run_clustering(
str(clustering_dir),
"--k",
str(int(k)),
- "--random-state",
+ "--seed",
str(int(clustering_seed)),
"--normalize",
"--normalize-method",
"minmax",
]
- logger.info("Running clustering (MiniBatchKMeans) k=%d output=%s", k, clustering_dir)
+ logger.info("Running clustering (MiniBatchKMeans) k=%d -> %s", k, clustering_dir)
_run_subprocess(cmd)
expected = clustering_dir / "cluster_assignments.parquet"
@@ -226,86 +695,343 @@ def run_clustering(
def _default_stage2_census_cache(*, stage2_dir: Path, state_fips: str, acs_year: int) -> Path:
- # Option 1: cache lives under the run directory
return stage2_dir / f"census_cache_{state_fips}_{acs_year}.parquet"
-def run_stage2_logratio(
+def _default_stage2_census_predictors(*, stage2_dir: Path) -> Path:
+ return stage2_dir / DEFAULT_STAGE2_CENSUS_PREDICTORS_BASENAME
+
+
+def _predictor_names_from_specs() -> list[str]:
+ names: list[str] = []
+ for spec in VARIABLE_SPECS:
+ nm = spec.get("name")
+ if nm and isinstance(nm, str):
+ names.append(nm)
+ seen: set[str] = set()
+ out: list[str] = []
+ for n in names:
+ if n in seen:
+ continue
+ seen.add(n)
+ out.append(n)
+ return out
+
+
+def build_stage2_census_predictors(
*,
- run_dir: Path,
- crosswalk_path: Path,
+ stage2_dir: Path,
state_fips: str,
acs_year: int,
- min_obs_per_bg: int,
- alpha: float,
- standardize: bool,
+ decennial_year: int,
fetch_census: bool,
- no_ols: bool,
- baseline_cluster: str | None,
- predictors_from: Path | None,
census_cache_path: Path | None,
-) -> None:
- stage2_dir = run_dir / "stage2"
+ census_predictors_path: Path | None,
+) -> Path:
+ """
+ Build (or reuse) the block-group census cache and derive a predictors parquet for the R stage2 script.
+
+ Source of truth for predictors:
+ smart_meter_analysis.census_specs.VARIABLE_SPECS[*]["name"]
+
+ Output:
+ stage2/census_predictors_blockgroups.parquet with:
+ - GEOID (12-digit BG GEOID)
+ - all predictor columns listed in VARIABLE_SPECS names (must exist; strict by default)
+
+ Write format:
+ Parquet with SNAPPY compression (broadly supported by R/arrow builds).
+ """
stage2_dir.mkdir(parents=True, exist_ok=True)
+ cache = census_cache_path or _default_stage2_census_cache(
+ stage2_dir=stage2_dir, state_fips=state_fips, acs_year=acs_year
+ )
+ preds_out = census_predictors_path or _default_stage2_census_predictors(stage2_dir=stage2_dir)
+
+ if fetch_census or (not cache.exists()):
+ logger.info("Building census cache via census.py -> %s", cache)
+ fetch_census_data(
+ state_fips=str(state_fips),
+ acs_year=int(acs_year),
+ decennial_year=int(decennial_year),
+ county_fips=None,
+ output_path=cache,
+ keep_raw_debug_cols=None,
+ )
+ else:
+ logger.info("Using existing census cache: %s", cache)
+
+ if not cache.exists():
+ raise FileNotFoundError(f"Census cache not found after build attempt: {cache}")
+
+ df = pl.read_parquet(str(cache))
+ cols = set(df.columns)
+
+ if "GEOID" not in cols:
+ raise ValueError(f"Census cache missing GEOID column. Columns={df.columns}")
+
+ predictor_names = _predictor_names_from_specs()
+ missing = [c for c in predictor_names if c not in cols]
+ if missing:
+ raise ValueError(
+ "Census cache is missing predictor column(s) required by census_specs.py VARIABLE_SPECS. "
+ f"Missing {len(missing)}: {missing[:50]}{' ...' if len(missing) > 50 else ''}"
+ )
+
+ out_df = df.select(["GEOID", *predictor_names]).unique(subset=["GEOID"], keep="first")
+ preds_out.parent.mkdir(parents=True, exist_ok=True)
+
+ table = out_df.to_arrow()
+ pq.write_table(table, preds_out, compression="snappy")
+ logger.info("Wrote census predictors parquet: %s (rows=%d cols=%d)", preds_out, out_df.height, out_df.width)
+ return preds_out
+
+
+def run_stage2_multinom_r(
+ *,
+ run_dir: Path,
+ crosswalk_path: Path,
+ census_predictors_path: Path,
+ out_dir: Path,
+ baseline_cluster: str | None,
+ alpha: float | None,
+ min_obs_per_bg: int,
+ allow_missing_predictors: bool,
+ standardize: int,
+ use_vgam: int,
+ verbose: bool,
+ no_emoji: bool,
+) -> None:
+ """
+ Execute the R Stage 2 multinomial BG regression script with the required inputs.
+
+ Pass-through flags:
+ --standardize (SD-only scaling; no mean centering)
+ --use-vgam (VGAM vglm vs nnet multinom)
+ """
clusters_path = run_dir / "clustering" / "cluster_assignments.parquet"
if not clusters_path.exists():
raise FileNotFoundError(f"Missing cluster assignments for Stage 2: {clusters_path}")
- if not crosswalk_path.exists():
- raise FileNotFoundError(f"Missing crosswalk for Stage 2: {crosswalk_path}")
+ crosswalk_path = _validate_crosswalk(crosswalk_path)
- script_path = Path("analysis/clustering/stage2_logratio_regression.py")
+ if not census_predictors_path.exists():
+ raise FileNotFoundError(f"Missing census predictors parquet for Stage 2: {census_predictors_path}")
+
+ script_path = Path("analysis/stage2/stage2_multinom_blockgroup_weighted.R")
if not script_path.exists():
- raise FileNotFoundError(f"Missing Stage 2 script: {script_path}")
+ raise FileNotFoundError(f"Missing Stage 2 R script: {script_path}")
- resolved_cache = census_cache_path or _default_stage2_census_cache(
- stage2_dir=stage2_dir,
- state_fips=state_fips,
- acs_year=acs_year,
- )
+ out_dir.mkdir(parents=True, exist_ok=True)
- cmd = [
- sys.executable,
+ cmd: list[str] = [
+ "Rscript",
str(script_path),
"--clusters",
str(clusters_path),
"--crosswalk",
str(crosswalk_path),
- "--output-dir",
- str(stage2_dir),
- "--census-cache",
- str(resolved_cache),
- "--state-fips",
- str(state_fips),
- "--acs-year",
- str(int(acs_year)),
+ "--census",
+ str(census_predictors_path),
+ "--out-dir",
+ str(out_dir),
"--min-obs-per-bg",
str(int(min_obs_per_bg)),
- "--alpha",
- str(float(alpha)),
+ "--allow-missing-predictors",
+ "1" if allow_missing_predictors else "0",
+ "--standardize",
+ str(int(standardize)),
+ "--use-vgam",
+ str(int(use_vgam)),
+ "--verbose",
+ "1" if verbose else "0",
+ "--no-emoji",
+ "1" if no_emoji else "0",
]
- if fetch_census:
- cmd.append("--fetch-census")
-
- if standardize:
- cmd.append("--standardize")
-
- if no_ols:
- cmd.append("--no-ols")
-
if baseline_cluster is not None:
cmd.extend(["--baseline-cluster", str(baseline_cluster)])
- if predictors_from is not None:
- cmd.extend(["--predictors-from", str(predictors_from)])
+ if alpha is not None:
+ cmd.extend(["--alpha", str(float(alpha))])
- logger.info("Running Stage 2 log-ratio regression output=%s", stage2_dir)
+ logger.info("Running Stage 2 (R multinomial) -> %s", out_dir)
_run_subprocess(cmd)
-def write_run_manifest(*, run_dir: Path, args: argparse.Namespace) -> None:
+def _hybrid_download_and_process(
+ *,
+ paths: StepPaths,
+ year_month: str,
+ day_mode: str,
+ download_chunk_size: int,
+ min_free_gb: float,
+ parts_dirname: str,
+ delete_raw_after_process: bool,
+ download_fail_fast: bool,
+ download_max_errors: int,
+ download_retries: int,
+ download_backoff_factor: float,
+ download_log_every: int,
+ skip_download: bool,
+ skip_ingest: bool,
+ hybrid_write_monolith: bool,
+) -> list[Path]:
+ raw_dir = paths.raw_dir
+ processed_dir = paths.processed_dir
+ parts_root = _parts_root(processed_dir, parts_dirname)
+ parts_root.mkdir(parents=True, exist_ok=True)
+
+ ym_raw_dir = raw_dir / year_month
+ ym_raw_dir.mkdir(parents=True, exist_ok=True)
+
+ cycle_id = _next_cycle_id(parts_root)
+ logger.info("Hybrid: next cycle_id=%d (existing cycles=%s)", cycle_id, _existing_cycle_ids(parts_root))
+
+ promoted = _persist_temp_batches_as_parts(
+ processed_dir=processed_dir, parts_dirname=parts_dirname, cycle_id=cycle_id
+ )
+ if promoted:
+ cycle_id += 1
+
+ def process_any_raw_dirs() -> None:
+ nonlocal cycle_id
+ dirs = _raw_processing_dirs(raw_dir)
+ if not dirs:
+ rec = _raw_csvs_recursive(raw_dir)
+ logger.info("Hybrid: raw recursive CSV count=%d; no processable top-level dirs found.", len(rec))
+ return
+
+ for d in dirs:
+ csv_count = len([p for p in list(d.glob("*.csv")) + list(d.glob("*.CSV")) if p.is_file()])
+ rec_total = len(_raw_csvs_recursive(d))
+ logger.info("Hybrid: process dir=%s top_level_csvs=%d recursive_csvs=%d", d, csv_count, rec_total)
+
+ if skip_ingest:
+ logger.info("Hybrid: --skip-ingest set; will not process raw dir %s", d)
+ continue
+
+ free_gb = _disk_free_gb(paths.run_dir)
+ if free_gb < float(min_free_gb):
+ raise RuntimeError(
+ f"Insufficient disk space before processing: free_gb={free_gb:.1f} < min_free_gb={min_free_gb:.1f}"
+ )
+
+ if hybrid_write_monolith:
+ out_path = processed_dir / f"comed_{year_month}_cycle_{cycle_id:04d}.parquet"
+ else:
+ out_path = processed_dir / f"_scratch_comed_{year_month}.parquet"
+
+ process_csvs(
+ input_dir=d,
+ processed_dir=processed_dir,
+ run_dir=paths.run_dir,
+ year_month=year_month,
+ day_mode=day_mode,
+ keep_temp=True,
+ output_path=out_path,
+ processing_manifest=paths.processing_manifest,
+ )
+
+ _persist_temp_batches_as_parts(
+ processed_dir=processed_dir,
+ parts_dirname=parts_dirname,
+ cycle_id=cycle_id,
+ )
+
+ if delete_raw_after_process:
+ logger.info("Hybrid: deleting raw CSVs in %s after successful processing (cycle=%d)", d, cycle_id)
+ _delete_csvs_under(d)
+
+ if not hybrid_write_monolith and out_path.exists():
+ with contextlib.suppress(Exception):
+ out_path.unlink(missing_ok=True)
+
+ cycle_id += 1
+
+ _prune_empty_dirs(raw_dir)
+
+ process_any_raw_dirs()
+
+ if skip_download:
+ logger.info("Hybrid: --skip-download set; no further S3 downloads will be performed.")
+ parts = _collect_part_parquets(parts_root)
+ if not parts:
+ raise RuntimeError(f"Hybrid mode (skip-download) found no parquet parts under: {parts_root}")
+ return parts
+
+ s3_keys_all = _safe_list_s3_keys(year_month=year_month, max_files=None)
+ logger.info("Hybrid: total S3 keys for %s: %d", year_month, len(s3_keys_all))
+
+ for i, start in enumerate(range(0, len(s3_keys_all), int(download_chunk_size))):
+ _log_system_diagnostics(base_path=paths.run_dir)
+
+ free_gb = _disk_free_gb(paths.run_dir)
+ if free_gb < float(min_free_gb):
+ logger.info(
+ "Hybrid: free space %.1f GB below threshold %.1f GB; processing existing raw before downloading more.",
+ free_gb,
+ float(min_free_gb),
+ )
+ process_any_raw_dirs()
+ free_gb = _disk_free_gb(paths.run_dir)
+ if free_gb < float(min_free_gb):
+ raise RuntimeError(
+ f"Insufficient disk after processing/cleanup: free_gb={free_gb:.1f} < min_free_gb={min_free_gb:.1f}"
+ )
+
+ end = min(len(s3_keys_all), start + int(download_chunk_size))
+ chunk = s3_keys_all[start:end]
+ logger.info("Hybrid: downloading chunk %d keys [%d, %d) size=%d -> %s", i, start, end, len(chunk), ym_raw_dir)
+
+ _download_s3_chunk(
+ s3_keys=chunk,
+ output_dir=ym_raw_dir,
+ manifest_path=paths.download_manifest,
+ fail_fast=bool(download_fail_fast),
+ max_errors=int(download_max_errors),
+ retries=int(download_retries),
+ backoff_factor=float(download_backoff_factor),
+ log_every=int(download_log_every),
+ chunk_id=i,
+ )
+
+ process_any_raw_dirs()
+
+ process_any_raw_dirs()
+
+ parts = _collect_part_parquets(parts_root)
+ if not parts:
+ raise RuntimeError(f"Hybrid mode produced no parquet parts under: {parts_root}")
+
+ logger.info("Hybrid complete: durable parts=%d under %s", len(parts), parts_root)
+ return parts
+
+
+def _parse_k_values(args: argparse.Namespace) -> list[int]:
+ if args.k_values:
+ parts = [p.strip() for p in str(args.k_values).split(",") if p.strip()]
+ ks = sorted({int(p) for p in parts})
+ if any(k < 2 for k in ks):
+ raise ValueError("--k-values must all be >= 2")
+ return ks
+ kmin = int(args.k_min)
+ kmax = int(args.k_max)
+ if kmin < 2 or kmax < 2 or kmax < kmin:
+ raise ValueError("--k-min/--k-max must be >=2 and k-max >= k-min")
+ return list(range(kmin, kmax + 1))
+
+
+def write_run_manifest(
+ *,
+ run_dir: Path,
+ args: argparse.Namespace,
+ selected_k: int | None,
+ k_search_enabled: bool,
+ k_search_output_dir: Path | None,
+ log_path: Path,
+) -> None:
repo_root = Path().resolve()
manifest: dict[str, Any] = {
@@ -315,153 +1041,527 @@ def write_run_manifest(*, run_dir: Path, args: argparse.Namespace) -> None:
"git_sha": _get_git_sha(repo_root),
"python_version": sys.version,
"polars_version": pl.__version__,
- "seeds": {
- "sampling": args.seed,
- "clustering": args.clustering_seed,
- },
- "directory_structure": "data/runs/{run_name}/",
- "pipeline_version": "v1.0-restored-baseline",
+ "selected_k": selected_k,
+ "k_search_enabled": bool(k_search_enabled),
+ "k_search_output_dir": str(k_search_output_dir) if k_search_output_dir else None,
+ "logs": {"pipeline_log": str(log_path)},
"cwd": str(Path.cwd()),
- "platform": {
- "os_name": os.name,
- "sys_platform": sys.platform,
- },
+ "platform": {"os_name": os.name, "sys_platform": sys.platform},
}
out = run_dir / "run_manifest.json"
- out.write_text(json.dumps(manifest, indent=2, sort_keys=True) + "\n", encoding="utf-8")
+ _atomic_write_text(out, json.dumps(manifest, indent=2, sort_keys=True) + "\n")
logger.info("Wrote run manifest: %s", out)
def run_pipeline(args: argparse.Namespace) -> None:
run_dir = DEFAULT_RUNS_DIR / args.run_name
- paths = create_directories(run_dir)
-
- _configure_logging(paths["logs"] / "pipeline.log")
+ paths = _ensure_paths(run_dir)
+ log_path = _configure_logging(paths.logs_dir)
+ _update_run_state(paths.run_state, {"run_name": args.run_name, "year_month": args.year_month, "status": "running"})
logger.info("Run directory: %s", run_dir)
logger.info("Year-month: %s", args.year_month)
- if args.from_s3:
- download_from_s3(year_month=args.year_month, num_files=int(args.num_files), run_dir=run_dir)
- else:
- logger.info("Skipping S3 download (--from-s3 not set). Expecting CSVs in: %s", paths["raw"])
+ raw_rec = _raw_csvs_recursive(paths.raw_dir)
+ logger.info("Raw recursive CSV count under %s: %d", paths.raw_dir, len(raw_rec))
- processed_parquet = process_csvs(run_dir=run_dir, year_month=args.year_month, day_mode=args.day_mode)
+ _log_system_diagnostics(base_path=paths.run_dir)
- profiles_path = prepare_clustering(
- run_dir=run_dir,
- processed_parquet=processed_parquet,
- sample_days=int(args.sample_days),
- sample_households=(int(args.sample_households) if args.sample_households is not None else None),
- seed=int(args.seed),
+ processed_inputs: list[Path] = []
+ selected_k: int | None = None
+ k_search_out_dir: Path | None = None
+
+ # ----------------------------
+ # Early exit: Stage 2 census build only
+ # ----------------------------
+ if bool(args.stage2_build_census_only):
+ preds_path = build_stage2_census_predictors(
+ stage2_dir=paths.stage2_dir,
+ state_fips=str(args.stage2_state_fips),
+ acs_year=int(args.stage2_acs_year),
+ decennial_year=int(args.stage2_decennial_year),
+ fetch_census=bool(args.stage2_fetch_census),
+ census_cache_path=(Path(args.stage2_census_cache) if args.stage2_census_cache is not None else None),
+ census_predictors_path=(
+ Path(args.stage2_census_predictors) if args.stage2_census_predictors is not None else None
+ ),
+ )
+ _update_run_state(
+ paths.run_state,
+ {"status": "running", "stage2_census_build": {"complete": True, "census_predictors": str(preds_path)}},
+ )
+ write_run_manifest(
+ run_dir=paths.run_dir,
+ args=args,
+ selected_k=None,
+ k_search_enabled=False,
+ k_search_output_dir=None,
+ log_path=log_path,
+ )
+ _update_run_state(paths.run_state, {"status": "complete"})
+ logger.info("Pipeline completed (stage2 census build only).")
+ return
+
+ # ----------------------------
+ # Early exit: Stage 2 only (no ingest/prepare/clustering/k-search)
+ # ----------------------------
+ do_stage2 = bool(args.run_stage2) and not bool(args.skip_stage2)
+ stage2_only = (
+ do_stage2
+ and bool(args.skip_download)
+ and bool(args.skip_ingest)
+ and bool(args.skip_prepare)
+ and bool(args.skip_k_search)
+ and bool(args.skip_cluster)
+ and (not bool(args.k_search))
)
+ if stage2_only:
+ preds_path = build_stage2_census_predictors(
+ stage2_dir=paths.stage2_dir,
+ state_fips=str(args.stage2_state_fips),
+ acs_year=int(args.stage2_acs_year),
+ decennial_year=int(args.stage2_decennial_year),
+ fetch_census=bool(args.stage2_fetch_census),
+ census_cache_path=(Path(args.stage2_census_cache) if args.stage2_census_cache is not None else None),
+ census_predictors_path=(
+ Path(args.stage2_census_predictors) if args.stage2_census_predictors is not None else None
+ ),
+ )
+ run_stage2_multinom_r(
+ run_dir=paths.run_dir,
+ crosswalk_path=Path(args.stage2_crosswalk),
+ census_predictors_path=preds_path,
+ out_dir=paths.stage2_r_dir,
+ baseline_cluster=(str(args.stage2_baseline_cluster) if args.stage2_baseline_cluster is not None else None),
+ alpha=(float(args.stage2_alpha) if args.stage2_alpha is not None else None),
+ min_obs_per_bg=int(args.stage2_min_obs_per_bg),
+ allow_missing_predictors=bool(args.stage2_allow_missing_predictors),
+ standardize=int(args.stage2_standardize),
+ use_vgam=int(args.stage2_use_vgam),
+ verbose=bool(args.stage2_verbose),
+ no_emoji=bool(args.stage2_no_emoji),
+ )
+ _update_run_state(
+ paths.run_state,
+ {"status": "running", "stage2": {"complete": True, "engine": "R", "out_dir": str(paths.stage2_r_dir)}},
+ )
+ write_run_manifest(
+ run_dir=paths.run_dir,
+ args=args,
+ selected_k=None,
+ k_search_enabled=False,
+ k_search_output_dir=None,
+ log_path=log_path,
+ )
+ _update_run_state(paths.run_state, {"status": "complete"})
+ logger.info("Pipeline completed (stage2-only).")
+ return
- run_clustering(
- run_dir=run_dir,
- profiles_path=profiles_path,
- k=int(args.k),
- clustering_seed=int(args.clustering_seed),
+ # ----------------------------
+ # Stage 0/1: Download + Ingest
+ # ----------------------------
+ if args.hybrid:
+ if not args.from_s3 and not args.skip_download:
+ raise ValueError("--hybrid requires --from-s3 unless --skip-download is used.")
+
+ if args.skip_ingest:
+ parts_root = _parts_root(paths.processed_dir, args.parts_dirname)
+ processed_inputs = _collect_part_parquets(parts_root)
+ if not processed_inputs:
+ raise RuntimeError(
+ f"--skip-ingest specified but no parts found under {parts_root}. "
+ "Either run ingest or point parts-dirname correctly."
+ )
+ else:
+ processed_inputs = _hybrid_download_and_process(
+ paths=paths,
+ year_month=args.year_month,
+ day_mode=args.day_mode,
+ download_chunk_size=int(args.download_chunk_size),
+ min_free_gb=float(args.min_free_gb),
+ parts_dirname=str(args.parts_dirname),
+ delete_raw_after_process=bool(args.delete_raw_after_process),
+ download_fail_fast=not bool(args.no_fail_fast),
+ download_max_errors=int(args.max_errors),
+ download_retries=int(args.retries),
+ download_backoff_factor=float(args.backoff_factor),
+ download_log_every=int(args.log_every),
+ skip_download=bool(args.skip_download),
+ skip_ingest=bool(args.skip_ingest),
+ hybrid_write_monolith=bool(args.hybrid_write_monolith),
+ )
+ else:
+ monolith = paths.processed_dir / f"comed_{args.year_month}.parquet"
+ if args.skip_ingest:
+ if not monolith.exists():
+ raise FileNotFoundError(f"--skip-ingest specified but monolith missing: {monolith}")
+ processed_inputs = [monolith]
+ else:
+ if args.from_s3 and not args.skip_download:
+ ym_raw_dir = paths.raw_dir / args.year_month
+ ym_raw_dir.mkdir(parents=True, exist_ok=True)
+
+ s3_keys = _safe_list_s3_keys(year_month=args.year_month, max_files=int(args.num_files))
+ logger.info("Non-hybrid: downloading %d file(s) -> %s", len(s3_keys), ym_raw_dir)
+ _download_s3_chunk(
+ s3_keys=s3_keys,
+ output_dir=ym_raw_dir,
+ manifest_path=paths.download_manifest,
+ fail_fast=not bool(args.no_fail_fast),
+ max_errors=int(args.max_errors),
+ retries=int(args.retries),
+ backoff_factor=float(args.backoff_factor),
+ log_every=int(args.log_every),
+ chunk_id=0,
+ )
+ else:
+ logger.info("Non-hybrid: skipping S3 download (expecting CSVs already under %s).", paths.raw_dir)
+
+ dirs = _raw_processing_dirs(paths.raw_dir)
+ if not dirs:
+ raise RuntimeError(
+ "Non-hybrid ingest found no processable CSV directories. "
+ "Check raw recursive CSV count and ensure CSVs exist under raw/ or raw/{YYYYMM}/."
+ )
+ if len(dirs) > 1:
+ logger.warning("Non-hybrid: multiple raw CSV dirs detected; using first deterministically: %s", dirs[0])
+
+ process_csvs(
+ input_dir=dirs[0],
+ processed_dir=paths.processed_dir,
+ run_dir=paths.run_dir,
+ year_month=args.year_month,
+ day_mode=args.day_mode,
+ keep_temp=False,
+ output_path=monolith,
+ processing_manifest=paths.processing_manifest,
+ )
+ if not monolith.exists():
+ raise FileNotFoundError(f"Non-hybrid expected monolith missing: {monolith}")
+ processed_inputs = [monolith]
+
+ _update_run_state(
+ paths.run_state,
+ {"status": "running", "ingest": {"complete": True, "inputs": [str(p) for p in processed_inputs]}},
)
- if bool(args.run_stage2):
- run_stage2_logratio(
- run_dir=run_dir,
- crosswalk_path=Path(args.stage2_crosswalk),
+ # ----------------------------
+ # Stage 1: Prepare clustering
+ # ----------------------------
+ profiles_path = paths.clustering_dir / "sampled_profiles.parquet"
+ if args.skip_prepare:
+ if not profiles_path.exists():
+ raise FileNotFoundError(f"--skip-prepare specified but missing: {profiles_path}")
+ logger.info("Skipping clustering preparation (using existing %s).", profiles_path)
+ else:
+ profiles_path = prepare_clustering(
+ clustering_dir=paths.clustering_dir,
+ processed_inputs=processed_inputs,
+ sample_days=int(args.sample_days),
+ sample_households=(int(args.sample_households) if args.sample_households is not None else None),
+ seed=int(args.seed),
+ )
+
+ _update_run_state(
+ paths.run_state, {"status": "running", "prepare": {"complete": True, "profiles": str(profiles_path)}}
+ )
+
+ # ----------------------------
+ # Stage 1.5: K-search (optional)
+ # ----------------------------
+ if args.k_search:
+ k_search_out_dir = paths.clustering_dir / str(args.k_search_output_dirname)
+ k_values = _parse_k_values(args)
+
+ if args.skip_k_search:
+ logger.info("Skipping k-search by flag.")
+ selected_k = None
+ if k_search_out_dir.exists():
+ selected_k = read_recommended_k(output_dir=k_search_out_dir, k_values=k_values)
+ logger.info("Using existing k-search recommendation: k=%d", selected_k)
+ else:
+ selected_k = run_k_search(
+ profiles_path=profiles_path,
+ output_dir=k_search_out_dir,
+ seed=int(args.seed),
+ k_values=k_values,
+ sample_households=(
+ int(args.k_search_sample_households) if args.k_search_sample_households is not None else None
+ ),
+ max_profiles=(int(args.k_search_max_profiles) if args.k_search_max_profiles is not None else None),
+ )
+
+ _update_run_state(
+ paths.run_state,
+ {
+ "status": "running",
+ "k_search": {"enabled": True, "output_dir": str(k_search_out_dir), "selected_k": selected_k},
+ },
+ )
+ else:
+ _update_run_state(paths.run_state, {"status": "running", "k_search": {"enabled": False}})
+
+ # ----------------------------
+ # Stage 1.6: Single-k clustering
+ # ----------------------------
+ assignments = paths.clustering_dir / "cluster_assignments.parquet"
+ if args.skip_cluster:
+ if not assignments.exists():
+ raise FileNotFoundError(f"--skip-cluster specified but missing: {assignments}")
+ logger.info("Skipping clustering (using existing %s).", assignments)
+ else:
+ k_final = selected_k if args.k is None or str(args.k).lower() == "auto" else int(args.k)
+
+ if k_final is None:
+ raise ValueError(
+ "No k available for clustering. Provide --k, or enable --k-search and use --k auto (or omit --k)."
+ )
+ if k_final < 2:
+ raise ValueError("--k must be >= 2")
+
+ run_clustering(
+ clustering_dir=paths.clustering_dir,
+ profiles_path=profiles_path,
+ k=int(k_final),
+ clustering_seed=int(args.clustering_seed),
+ )
+ selected_k = int(k_final)
+
+ _update_run_state(paths.run_state, {"status": "running", "cluster": {"complete": True, "selected_k": selected_k}})
+
+ # ----------------------------
+ # Stage 2: Regression (optional)
+ # ----------------------------
+ do_stage2 = bool(args.run_stage2) and not bool(args.skip_stage2)
+ if do_stage2:
+ preds_path = build_stage2_census_predictors(
+ stage2_dir=paths.stage2_dir,
state_fips=str(args.stage2_state_fips),
acs_year=int(args.stage2_acs_year),
- min_obs_per_bg=int(args.stage2_min_obs_per_bg),
- alpha=float(args.stage2_alpha),
- standardize=bool(args.stage2_standardize),
- no_ols=bool(args.stage2_no_ols),
+ decennial_year=int(args.stage2_decennial_year),
fetch_census=bool(args.stage2_fetch_census),
- baseline_cluster=(str(args.stage2_baseline_cluster) if args.stage2_baseline_cluster is not None else None),
- predictors_from=(Path(args.stage2_predictors_from) if args.stage2_predictors_from is not None else None),
census_cache_path=(Path(args.stage2_census_cache) if args.stage2_census_cache is not None else None),
+ census_predictors_path=(
+ Path(args.stage2_census_predictors) if args.stage2_census_predictors is not None else None
+ ),
)
- write_run_manifest(run_dir=run_dir, args=args)
+ run_stage2_multinom_r(
+ run_dir=paths.run_dir,
+ crosswalk_path=Path(args.stage2_crosswalk),
+ census_predictors_path=preds_path,
+ out_dir=paths.stage2_r_dir,
+ baseline_cluster=(str(args.stage2_baseline_cluster) if args.stage2_baseline_cluster is not None else None),
+ alpha=(float(args.stage2_alpha) if args.stage2_alpha is not None else None),
+ min_obs_per_bg=int(args.stage2_min_obs_per_bg),
+ allow_missing_predictors=bool(args.stage2_allow_missing_predictors),
+ standardize=int(args.stage2_standardize),
+ use_vgam=int(args.stage2_use_vgam),
+ verbose=bool(args.stage2_verbose),
+ no_emoji=bool(args.stage2_no_emoji),
+ )
+ _update_run_state(
+ paths.run_state,
+ {"status": "running", "stage2": {"complete": True, "engine": "R", "out_dir": str(paths.stage2_r_dir)}},
+ )
+ else:
+ _update_run_state(paths.run_state, {"status": "running", "stage2": {"complete": False, "skipped": True}})
+ # ----------------------------
+ # Finalize
+ # ----------------------------
+ write_run_manifest(
+ run_dir=paths.run_dir,
+ args=args,
+ selected_k=selected_k,
+ k_search_enabled=bool(args.k_search),
+ k_search_output_dir=k_search_out_dir,
+ log_path=log_path,
+ )
+ _update_run_state(paths.run_state, {"status": "complete"})
logger.info("Pipeline completed successfully.")
def build_parser() -> argparse.ArgumentParser:
- p = argparse.ArgumentParser(description="ComEd Smart Meter Analysis Pipeline Orchestrator (Phase 3)")
+ p = argparse.ArgumentParser(description="ComEd Smart Meter Analysis Pipeline Orchestrator (Production-Grade)")
+ # Core
p.add_argument("--run-name", required=True, help="Run name (directory under data/runs/)")
p.add_argument("--year-month", required=True, help="Target month in YYYYMM format (e.g., 202307)")
p.add_argument("--from-s3", action="store_true", help="Download CSVs from S3 into run_dir/raw/")
- p.add_argument("--num-files", type=int, default=10, help="Number of S3 files to download (default: 10)")
+ # Skip flags
+ p.add_argument("--skip-download", action="store_true", help="Skip S3 download step.")
+ p.add_argument("--skip-ingest", action="store_true", help="Skip CSV->parquet ingest step (use existing outputs).")
+ p.add_argument(
+ "--skip-prepare", action="store_true", help="Skip clustering preparation (use existing sampled_profiles)."
+ )
+ p.add_argument("--skip-k-search", action="store_true", help="Skip k-search run (use existing outputs if present).")
+ p.add_argument("--skip-cluster", action="store_true", help="Skip clustering (use existing cluster_assignments).")
+ p.add_argument("--skip-stage2", action="store_true", help="Skip Stage 2 even if --run-stage2 is set.")
+
+ # Download knobs
+ p.add_argument("--num-files", type=int, default=10, help="Non-hybrid: number of S3 files to download (default: 10)")
+ p.add_argument("--no-fail-fast", action="store_true", help="Download: allow errors up to --max-errors")
+ p.add_argument(
+ "--max-errors", type=int, default=10, help="Download: max errors before aborting (when --no-fail-fast)"
+ )
+ p.add_argument("--retries", type=int, default=3, help="Download: retries per file (in addition to initial attempt)")
+ p.add_argument("--backoff-factor", type=float, default=2.0, help="Download: exponential backoff factor")
+ p.add_argument("--log-every", type=int, default=100, help="Download: progress logging interval")
- p.add_argument("--sample-days", type=int, default=31, help="Days to sample for clustering (default: 31)")
+ # Hybrid
+ p.add_argument(
+ "--hybrid",
+ action="store_true",
+ help="Enable disk-safe hybrid loop (download chunk -> process -> delete raw -> repeat).",
+ )
+ p.add_argument(
+ "--download-chunk-size",
+ type=int,
+ default=DEFAULT_DOWNLOAD_CHUNK_SIZE,
+ help=f"Hybrid: keys per download chunk (default: {DEFAULT_DOWNLOAD_CHUNK_SIZE})",
+ )
+ p.add_argument(
+ "--min-free-gb",
+ type=float,
+ default=DEFAULT_MIN_FREE_GB,
+ help=f"Hybrid: minimum free disk GB (default: {DEFAULT_MIN_FREE_GB})",
+ )
+ p.add_argument(
+ "--parts-dirname",
+ type=str,
+ default=DEFAULT_PARTS_DIRNAME,
+ help=f"Hybrid: processed/ parts dir name (default: {DEFAULT_PARTS_DIRNAME})",
+ )
+ p.add_argument(
+ "--delete-raw-after-process",
+ action="store_true",
+ help="Hybrid: delete raw CSVs after each successful processing cycle (recommended).",
+ )
+ p.add_argument(
+ "--hybrid-write-monolith", action="store_true", help="Hybrid: keep per-cycle monolith parquet outputs."
+ )
+
+ # Ingest
+ p.add_argument("--day-mode", choices=["calendar", "billing"], default="calendar", help="Day attribution mode")
+
+ # Prepare
+ p.add_argument(
+ "--sample-days", type=int, default=DEFAULT_SAMPLE_DAYS, help=f"Days to sample (default: {DEFAULT_SAMPLE_DAYS})"
+ )
p.add_argument(
"--sample-households",
type=int,
default=None,
help="Households to sample (default: all). Provide an integer to limit.",
)
+ p.add_argument("--seed", type=int, default=DEFAULT_SEED, help=f"Random seed for sampling (default: {DEFAULT_SEED})")
+ p.add_argument(
+ "--clustering-seed",
+ type=int,
+ default=DEFAULT_SEED,
+ help=f"Random seed for clustering (default: {DEFAULT_SEED})",
+ )
- p.add_argument("--seed", type=int, default=42, help="Random seed for sampling (default: 42)")
- p.add_argument("--clustering-seed", type=int, default=42, help="Random seed for clustering (default: 42)")
-
- p.add_argument("--k", type=int, required=True, help="Number of clusters (single k per run)")
-
- p.add_argument("--day-mode", choices=["calendar", "billing"], default="calendar", help="Day attribution mode")
+ # K-search
+ p.add_argument("--k-search", action="store_true", help="Enable k-selection search (runs before clustering).")
+ p.add_argument("--k-min", type=int, default=DEFAULT_K_MIN, help=f"k-search: minimum k (default: {DEFAULT_K_MIN})")
+ p.add_argument("--k-max", type=int, default=DEFAULT_K_MAX, help=f"k-search: maximum k (default: {DEFAULT_K_MAX})")
+ p.add_argument(
+ "--k-values", default=None, help="k-search: explicit comma-separated k values (overrides --k-min/--k-max)."
+ )
+ p.add_argument(
+ "--k-search-sample-households", type=int, default=None, help="k-search: limit households used by k-search."
+ )
+ p.add_argument(
+ "--k-search-max-profiles",
+ type=int,
+ default=DEFAULT_K_SEARCH_MAX_PROFILES,
+ help=f"k-search: cap profiles for k-search (default: {DEFAULT_K_SEARCH_MAX_PROFILES})",
+ )
+ p.add_argument(
+ "--k-search-output-dirname",
+ type=str,
+ default=DEFAULT_K_SEARCH_DIRNAME,
+ help=f"k-search: output dir under clustering/ (default: {DEFAULT_K_SEARCH_DIRNAME})",
+ )
- p.add_argument("--run-stage2", action="store_true", help="Run Stage 2 log-ratio regression (optional)")
+ # Clustering k
+ p.add_argument(
+ "--k",
+ default=None,
+ help="Number of clusters for single-k clustering. If --k-search is enabled, use --k auto or omit --k.",
+ )
- # ----------------------------
- # Stage 2 options (user-facing)
- # ----------------------------
+ # Stage 2 (R multinomial)
+ p.add_argument("--run-stage2", action="store_true", help="Run Stage 2 regression (R multinomial blockgroup model)")
+ p.add_argument(
+ "--stage2-build-census-only",
+ action="store_true",
+ help="Stage 2: build census cache + predictors parquet only, then exit.",
+ )
p.add_argument(
"--stage2-crosswalk",
default=str(DEFAULT_CROSSWALK_PATH),
- help=f"ZIP+4 → block-group crosswalk path (default: {DEFAULT_CROSSWALK_PATH})",
+ help=f"ZIP+4 → BG crosswalk path (default: {DEFAULT_CROSSWALK_PATH})",
)
p.add_argument(
- "--stage2-state-fips",
- default=DEFAULT_STATE_FIPS,
- help=f"State FIPS (default: {DEFAULT_STATE_FIPS})",
+ "--stage2-state-fips", default=DEFAULT_STATE_FIPS, help=f"State FIPS (default: {DEFAULT_STATE_FIPS})"
)
p.add_argument(
- "--stage2-acs-year",
- type=int,
- default=DEFAULT_ACS_YEAR,
- help=f"ACS year (default: {DEFAULT_ACS_YEAR})",
+ "--stage2-acs-year", type=int, default=DEFAULT_ACS_YEAR, help=f"ACS year (default: {DEFAULT_ACS_YEAR})"
)
p.add_argument(
- "--stage2-fetch-census",
- action="store_true",
- help="Stage 2: force re-fetch Census data (ignore cache)",
+ "--stage2-decennial-year", type=int, default=2020, help="Decennial year for urban_percent (default: 2020)"
)
p.add_argument(
- "--stage2-census-cache",
+ "--stage2-fetch-census", action="store_true", help="Stage 2: force re-fetch Census data (ignore cache)"
+ )
+ p.add_argument("--stage2-census-cache", default=None, help="Optional census cache parquet path.")
+ p.add_argument(
+ "--stage2-census-predictors",
default=None,
- help=(
- "Optional census cache parquet path. If omitted, defaults to "
- "data/runs/{run_name}/stage2/census_cache_{state_fips}_{acs_year}.parquet (Option 1)."
- ),
+ help="Optional census predictors parquet path (default: stage2/census_predictors_blockgroups.parquet).",
+ )
+ p.add_argument(
+ "--stage2-min-obs-per-bg", type=int, default=50, help="Stage 2: min household-days per BG (default: 50)"
+ )
+ p.add_argument(
+ "--stage2-alpha",
+ type=float,
+ default=None,
+ help="Optional: pass alpha to Stage2 R script. If omitted, no --alpha flag is passed.",
)
-
p.add_argument(
- "--stage2-min-obs-per-bg",
+ "--stage2-allow-missing-predictors",
type=int,
- default=50,
- help="Stage 2 minimum household-day observations per block group (default: 50)",
+ choices=[0, 1],
+ default=0,
+ help="If 1, drop block groups with any predictor NA (complete-case); if 0, abort if any BG would be dropped (default: 0).",
)
- p.add_argument("--stage2-alpha", type=float, default=0.5, help="Stage 2 Laplace smoothing alpha (default: 0.5)")
- p.add_argument("--stage2-standardize", action="store_true", help="Stage 2: standardize predictors")
- p.add_argument("--stage2-no-ols", action="store_true", help="Stage 2: skip OLS robustness check")
p.add_argument(
"--stage2-baseline-cluster",
default=None,
help="Stage 2: optional baseline cluster label (default: most frequent cluster)",
)
+ p.add_argument("--stage2-verbose", action="store_true", help="Stage 2 (R): verbose logs.")
+ p.add_argument("--stage2-no-emoji", action="store_true", help="Stage 2 (R): disable unicode icons.")
+
p.add_argument(
- "--stage2-predictors-from",
- default=None,
- help="Stage 2: optional path to predictors list (one per line) to force exact predictors",
+ "--stage2-standardize",
+ type=int,
+ choices=[0, 1],
+ default=0,
+ help="Stage 2 (R): divide predictors by SD only (NO mean-centering). 0=off, 1=on (default: 0).",
+ )
+ p.add_argument(
+ "--stage2-use-vgam",
+ type=int,
+ choices=[0, 1],
+ default=1,
+ help="Stage 2 (R): 1=use VGAM::vglm (IRLS), 0=use nnet::multinom (default: 1).",
)
return p
@@ -470,11 +1570,32 @@ def build_parser() -> argparse.ArgumentParser:
def main() -> None:
args = build_parser().parse_args()
- if int(args.num_files) <= 0 and args.from_s3:
- raise ValueError("--num-files must be > 0 when using --from-s3")
-
- if int(args.k) <= 1:
- raise ValueError("--k must be >= 2")
+ # Validation
+ if args.hybrid:
+ if int(args.download_chunk_size) <= 0:
+ raise ValueError("--download-chunk-size must be > 0")
+ if float(args.min_free_gb) <= 0:
+ raise ValueError("--min-free-gb must be > 0")
+ else:
+ if args.from_s3 and not args.skip_download and int(args.num_files) <= 0:
+ raise ValueError("--num-files must be > 0 when using --from-s3 (non-hybrid mode)")
+
+ if args.k_search_max_profiles is not None and int(args.k_search_max_profiles) <= 0:
+ raise ValueError("--k-search-max-profiles must be > 0")
+
+ # k validation: only required if we will actually run clustering
+ if not bool(args.skip_cluster):
+ if not bool(args.k_search):
+ if args.k is None or str(args.k).lower() == "auto":
+ raise ValueError("--k is required when --k-search is not enabled (unless --skip-cluster is set)")
+ if int(args.k) < 2:
+ raise ValueError("--k must be >= 2")
+ else:
+ if args.k is not None and str(args.k).lower() != "auto" and int(args.k) < 2:
+ raise ValueError("--k must be >= 2")
+ else:
+ if args.k is not None and str(args.k).lower() != "auto" and int(args.k) < 2:
+ raise ValueError("--k must be >= 2")
run_pipeline(args)
diff --git a/smart_meter_analysis/census.py b/smart_meter_analysis/census.py
index 47f1ff2..2849e90 100644
--- a/smart_meter_analysis/census.py
+++ b/smart_meter_analysis/census.py
@@ -296,7 +296,7 @@ def fetch_acs_data(
logger.info(f"Retrieved {df.height} block groups from ACS raw pull")
# Debug: Check if required variables are present
- required_vars = ["B17001_001E", "B17001_002E", "B22001_001E", "B22001_002E", "B03001_001E", "B03001_003E"]
+ required_vars: list[str] = []
missing_vars = [v for v in required_vars if v not in df.columns]
if missing_vars:
logger.warning(f"Required variables not found in API response: {missing_vars}")
diff --git a/smart_meter_analysis/census_specs.py b/smart_meter_analysis/census_specs.py
index c617df2..3f81a31 100644
--- a/smart_meter_analysis/census_specs.py
+++ b/smart_meter_analysis/census_specs.py
@@ -3,7 +3,9 @@
- VARIABLE_SPECS is the spec-driven registry used by census.py to decide which ACS/Decennial
variables to request and how to engineer features.
-- STAGE2_PREDICTORS_47 is the *stable* final predictor list used by Stage 2 regression.
+- STAGE2_PREDICTORS_47 is the stable, documented predictor inventory.
+- STAGE2_PREDICTORS_MODEL_42 is the deterministic model input list after omitting one
+ reference variable per compositional share group (to ensure identifiability).
"""
from __future__ import annotations
@@ -63,6 +65,80 @@
"urban_percent",
]
+EXPECTED_STAGE2_PREDICTOR_COUNT = len(STAGE2_PREDICTORS_47)
+
+# -----------------------------------------------------------------------------
+# Deterministic compositional reference specification
+# -----------------------------------------------------------------------------
+# These groups contain "shares" that sum to (approximately) 100 within a block group.
+# To avoid a non-unique regression (perfect multicollinearity), we omit one reference
+# variable per group from the model input list. The omitted share is documented and
+# can be derived from the included shares:
+# reference_share = 100 - sum(other_shares_in_group)
+STAGE2_COMPOSITION_GROUPS: list[dict[str, Any]] = [
+ {
+ "group_id": "income_3bin",
+ "group_name": "Household income distribution (3 bins)",
+ "variables": ["pct_income_under_25k", "pct_income_25k_to_75k", "pct_income_75k_plus"],
+ "reference": "pct_income_75k_plus",
+ "rationale": (
+ "Shares sum to 100; omit one bin to ensure identifiability. "
+ "Highest-income share provides a stable reference for interpreting lower/middle shares."
+ ),
+ },
+ {
+ "group_id": "tenure_2bin",
+ "group_name": "Housing tenure (occupied units)",
+ "variables": ["pct_owner_occupied", "pct_renter_occupied"],
+ "reference": "pct_owner_occupied",
+ "rationale": "Owner + renter shares sum to 100; omit owners so renter share is interpreted relative to owners.",
+ },
+ {
+ "group_id": "vintage_3bin",
+ "group_name": "Housing vintage distribution (3 bins)",
+ "variables": ["pct_housing_built_2000_plus", "pct_housing_built_1980_1999", "old_building_pct"],
+ "reference": "pct_housing_built_2000_plus",
+ "rationale": "Shares sum to 100; omit newest stock so older stock shares are interpreted relative to newer construction.",
+ },
+ {
+ "group_id": "home_value_3bin",
+ "group_name": "Owner-occupied home value distribution (3 bins)",
+ "variables": ["pct_home_value_under_150k", "pct_home_value_150k_to_299k", "pct_home_value_300k_plus"],
+ "reference": "pct_home_value_300k_plus",
+ "rationale": "Shares sum to 100; omit highest-value share so lower/middle shares are interpreted relative to highest-value homes.",
+ },
+ {
+ "group_id": "age_6bin",
+ "group_name": "Population age distribution (6 bins)",
+ "variables": [
+ "pct_population_under_5",
+ "pct_population_5_to_17",
+ "pct_population_18_to_24",
+ "pct_population_25_to_44",
+ "pct_population_45_to_64",
+ "pct_population_65_plus",
+ ],
+ "reference": "pct_population_25_to_44",
+ "rationale": "Shares sum to 100; omit 25-44 as a central working-age reference group.",
+ },
+]
+
+STAGE2_REFERENCE_PREDICTORS: set[str] = {g["reference"] for g in STAGE2_COMPOSITION_GROUPS}
+
+# Validate references exist in the 47-list (fail fast if someone edits lists incorrectly)
+_missing_refs = sorted([r for r in STAGE2_REFERENCE_PREDICTORS if r not in STAGE2_PREDICTORS_47])
+if _missing_refs:
+ raise ValueError(f"Reference predictors not found in STAGE2_PREDICTORS_47: {_missing_refs}")
+
+# Model input list: keep original order but omit references deterministically
+STAGE2_PREDICTORS_MODEL_42: list[str] = [p for p in STAGE2_PREDICTORS_47 if p not in STAGE2_REFERENCE_PREDICTORS]
+
+EXPECTED_STAGE2_MODEL_PREDICTOR_COUNT = len(STAGE2_PREDICTORS_MODEL_42)
+if EXPECTED_STAGE2_PREDICTOR_COUNT != 47:
+ raise ValueError(f"Expected 47 predictors, found {EXPECTED_STAGE2_PREDICTOR_COUNT}")
+if EXPECTED_STAGE2_MODEL_PREDICTOR_COUNT != 42:
+ raise ValueError(f"Expected 42 model predictors, found {EXPECTED_STAGE2_MODEL_PREDICTOR_COUNT}")
+
# -----------------------------------------------------------------------------
# Complete spec registry with ACS variable codes
# -----------------------------------------------------------------------------
@@ -172,7 +248,7 @@
"name": "pct_housing_built_2000_plus",
"category": "housing",
"source": "acs",
- "numerator": "B25034_002E + B25034_003E + B25034_004E + B25034_005E", # 2000-2009, 2010-2013, 2014-2017, 2018-2020+
+ "numerator": "B25034_002E + B25034_003E + B25034_004E + B25034_005E",
"denominator": "B25034_001E",
"transformation": "none",
},
@@ -180,7 +256,7 @@
"name": "pct_housing_built_1980_1999",
"category": "housing",
"source": "acs",
- "numerator": "B25034_006E + B25034_007E", # 1990-1999, 1980-1989
+ "numerator": "B25034_006E + B25034_007E",
"denominator": "B25034_001E",
"transformation": "none",
},
@@ -188,7 +264,7 @@
"name": "old_building_pct",
"category": "housing",
"source": "acs",
- "numerator": "B25034_008E + B25034_009E + B25034_010E + B25034_011E", # Pre-1980
+ "numerator": "B25034_008E + B25034_009E + B25034_010E + B25034_011E",
"denominator": "B25034_001E",
"transformation": "none",
},
@@ -196,7 +272,7 @@
"name": "pct_structure_single_family_detached",
"category": "housing",
"source": "acs",
- "numerator": "B25024_002E", # 1-unit detached
+ "numerator": "B25024_002E",
"denominator": "B25024_001E",
"transformation": "none",
},
@@ -204,7 +280,7 @@
"name": "pct_structure_single_family_attached",
"category": "housing",
"source": "acs",
- "numerator": "B25024_003E", # 1-unit attached
+ "numerator": "B25024_003E",
"denominator": "B25024_001E",
"transformation": "none",
},
@@ -212,7 +288,7 @@
"name": "pct_structure_multifamily_2_to_4",
"category": "housing",
"source": "acs",
- "numerator": "B25024_004E + B25024_005E", # 2 units + 3-4 units
+ "numerator": "B25024_004E + B25024_005E",
"denominator": "B25024_001E",
"transformation": "none",
},
@@ -220,7 +296,7 @@
"name": "pct_structure_multifamily_5_to_19",
"category": "housing",
"source": "acs",
- "numerator": "B25024_006E + B25024_007E", # 5-9 units + 10-19 units
+ "numerator": "B25024_006E + B25024_007E",
"denominator": "B25024_001E",
"transformation": "none",
},
@@ -228,7 +304,7 @@
"name": "pct_structure_multifamily_20_plus",
"category": "housing",
"source": "acs",
- "numerator": "B25024_008E + B25024_009E", # 20-49 units + 50+ units
+ "numerator": "B25024_008E + B25024_009E",
"denominator": "B25024_001E",
"transformation": "none",
},
@@ -236,7 +312,7 @@
"name": "pct_structure_multifamily_10_plus",
"category": "housing",
"source": "acs",
- "numerator": "B25024_007E + B25024_008E + B25024_009E", # 10-19, 20-49, 50+
+ "numerator": "B25024_007E + B25024_008E + B25024_009E",
"denominator": "B25024_001E",
"transformation": "none",
},
@@ -244,7 +320,7 @@
"name": "pct_structure_mobile_home",
"category": "housing",
"source": "acs",
- "numerator": "B25024_010E", # Mobile home
+ "numerator": "B25024_010E",
"denominator": "B25024_001E",
"transformation": "none",
},
@@ -276,7 +352,7 @@
"name": "pct_rent_burden_30_plus",
"category": "housing",
"source": "acs",
- "numerator": "B25070_007E + B25070_008E + B25070_009E + B25070_010E", # 30-34.9%, 35-39.9%, 40-49.9%, 50%+
+ "numerator": "B25070_007E + B25070_008E + B25070_009E + B25070_010E",
"denominator": "B25070_001E",
"transformation": "none",
},
@@ -284,7 +360,7 @@
"name": "pct_rent_burden_50_plus",
"category": "housing",
"source": "acs",
- "numerator": "B25070_010E", # 50%+
+ "numerator": "B25070_010E",
"denominator": "B25070_001E",
"transformation": "none",
},
@@ -292,15 +368,15 @@
"name": "pct_owner_cost_burden_30_plus_mortgage",
"category": "housing",
"source": "acs",
- "numerator": "B25091_008E + B25091_009E + B25091_010E + B25091_011E", # With mortgage: 30-34.9%, 35-39.9%, 40-49.9%, 50%+
- "denominator": "B25091_002E", # With mortgage total
+ "numerator": "B25091_008E + B25091_009E + B25091_010E + B25091_011E",
+ "denominator": "B25091_002E",
"transformation": "none",
},
{
"name": "pct_owner_cost_burden_50_plus_mortgage",
"category": "housing",
"source": "acs",
- "numerator": "B25091_011E", # With mortgage: 50%+
+ "numerator": "B25091_011E",
"denominator": "B25091_002E",
"transformation": "none",
},
@@ -308,16 +384,16 @@
"name": "pct_owner_overcrowded_2plus_per_room",
"category": "housing",
"source": "acs",
- "numerator": "B25014_007E", # Owner: 2+ persons per room
- "denominator": "B25014_002E", # Owner total
+ "numerator": "B25014_007E",
+ "denominator": "B25014_002E",
"transformation": "none",
},
{
"name": "pct_renter_overcrowded_2plus_per_room",
"category": "housing",
"source": "acs",
- "numerator": "B25014_013E", # Renter: 2+ persons per room
- "denominator": "B25014_008E", # Renter total
+ "numerator": "B25014_013E",
+ "denominator": "B25014_008E",
"transformation": "none",
},
# -------------------------------------------------------------------------
@@ -341,8 +417,8 @@
"name": "pct_single_parent_households",
"category": "household",
"source": "acs",
- "numerator": "B11001_006E + B11001_007E", # Male householder + Female householder, no spouse
- "denominator": "B11001_001E", # Total households
+ "numerator": "B11001_006E + B11001_007E",
+ "denominator": "B11001_001E",
"transformation": "none",
},
# -------------------------------------------------------------------------
@@ -359,15 +435,15 @@
"name": "pct_female",
"category": "demographic",
"source": "acs",
- "numerator": "B01001_026E", # Female total
- "denominator": "B01001_001E", # Total population
+ "numerator": "B01001_026E",
+ "denominator": "B01001_001E",
"transformation": "none",
},
{
"name": "pct_white_alone",
"category": "demographic",
"source": "acs",
- "numerator": "B03002_003E", # White alone, not Hispanic/Latino
+ "numerator": "B03002_003E",
"denominator": "B03002_001E",
"transformation": "none",
},
@@ -375,7 +451,7 @@
"name": "pct_black_alone",
"category": "demographic",
"source": "acs",
- "numerator": "B03002_004E", # Black alone, not Hispanic/Latino
+ "numerator": "B03002_004E",
"denominator": "B03002_001E",
"transformation": "none",
},
@@ -383,7 +459,7 @@
"name": "pct_asian_alone",
"category": "demographic",
"source": "acs",
- "numerator": "B03002_006E", # Asian alone, not Hispanic/Latino
+ "numerator": "B03002_006E",
"denominator": "B03002_001E",
"transformation": "none",
},
@@ -391,7 +467,7 @@
"name": "pct_two_or_more_races",
"category": "demographic",
"source": "acs",
- "numerator": "B03002_009E", # Two or more races, not Hispanic/Latino
+ "numerator": "B03002_009E",
"denominator": "B03002_001E",
"transformation": "none",
},
@@ -399,7 +475,7 @@
"name": "pct_population_under_5",
"category": "demographic",
"source": "acs",
- "numerator": "B01001_003E + B01001_027E", # Male + Female under 5
+ "numerator": "B01001_003E + B01001_027E",
"denominator": "B01001_001E",
"transformation": "none",
},
@@ -450,11 +526,8 @@
"name": "urban_percent",
"category": "spatial",
"source": "decennial",
- "numerator": "H2_002N", # Urban housing units
- "denominator": "H2_002N + H2_003N", # Total (urban + rural)
+ "numerator": "H2_002N",
+ "denominator": "H2_002N + H2_003N",
"transformation": "none",
},
]
-
-# Optional convenience: expose the expected count
-EXPECTED_STAGE2_PREDICTOR_COUNT = len(STAGE2_PREDICTORS_47)
diff --git a/test b/test
new file mode 100644
index 0000000..e69de29
diff --git a/uv.lock b/uv.lock
index f299c53..67dd1a3 100644
--- a/uv.lock
+++ b/uv.lock
@@ -756,6 +756,15 @@ wheels = [
{ url = "https://files.pythonhosted.org/packages/f7/ec/67fbef5d497f86283db54c22eec6f6140243aae73265799baaaa19cd17fb/ghp_import-2.1.0-py3-none-any.whl", hash = "sha256:8337dd7b50877f163d4c0289bc1f1c7f127550241988d568c1db512c4324a619", size = 11034, upload-time = "2022-05-02T15:47:14.552Z" },
]
+[[package]]
+name = "git-filter-repo"
+version = "2.47.0"
+source = { registry = "https://pypi.org/simple" }
+sdist = { url = "https://files.pythonhosted.org/packages/b9/bb/7a283f568af6b0528ade65e8ace84bd6ba46003e429101bcd62c232d01a5/git_filter_repo-2.47.0.tar.gz", hash = "sha256:411b27e68a080c07a69c233cb526dbc2d848b09a72f10477f4444dd0822cf290", size = 275743, upload-time = "2024-12-04T03:10:48.2Z" }
+wheels = [
+ { url = "https://files.pythonhosted.org/packages/60/60/d3943f0880ebcb7e0bdf79254d10dddd39c7b656eeecae32b8806ff66dec/git_filter_repo-2.47.0-py3-none-any.whl", hash = "sha256:2cd04929b9024e83e65db571cbe36aec65ead0cb5f9ec5abe42158654af5ad83", size = 76282, upload-time = "2024-12-04T03:10:46.064Z" },
+]
+
[[package]]
name = "griffe"
version = "1.15.0"
@@ -2898,6 +2907,7 @@ dependencies = [
{ name = "boto3" },
{ name = "botocore" },
{ name = "cenpy" },
+ { name = "git-filter-repo" },
{ name = "ipykernel" },
{ name = "matplotlib" },
{ name = "memory-profiler" },
@@ -2938,6 +2948,7 @@ requires-dist = [
{ name = "boto3", specifier = ">=1.40.46" },
{ name = "botocore", specifier = ">=1.40.47" },
{ name = "cenpy", specifier = ">=1.0.1" },
+ { name = "git-filter-repo", specifier = ">=2.47.0" },
{ name = "ipykernel", specifier = ">=6.30.1" },
{ name = "matplotlib", specifier = ">=3.9.4" },
{ name = "memory-profiler", specifier = ">=0.61.0" },