Skip to content

Build analysis notebook for benchmark results #13

Description

@jathavaan

Task description

Context

The benchmarking framework (doppa) runs containerised experiments on Azure Container Instances, with related benchmarks executing concurrently in the same wall-clock window (one ACI container per experiment within a batch). Each single-machine experiment uses a sequential stopping rule: the iteration loop runs until the bootstrapped 95% CI on the mean elapsed time converges within 5% relative half-width, bounded below by BENCHMARK_MIN_ITERATIONS = 10 iterations and a 60-second timed-window floor, and above by the per-query-type iteration ceiling from BenchmarkIteration and a hard timeout of BENCHMARK_MAX_TIMED_WINDOW_SECONDS = 3600s. National-scale spatial joins opt out of sequential stopping (use_sequential_stopping=False) and run a fixed 5-iteration count with warmup_iterations=1.

The full experiment consists of 3 passes (benchmark_run = 1, 2, 3), each executed as a separate orchestrator invocation. Per-iteration samples are persisted to Azure Blob Storage as Parquet files in a Hive-partitioned layout:

benchmarks/
  query_id=<id>/
    run_id=<run_id>/
      benchmark_run=<pass>/
        iteration=<n>/
          data.parquet

Run-level metadata (achieved iterations, stop reason, CI statistics) is appended to a single benchmark_metadata.parquet in the metadata container. Cost analytics are stored per-query as separate Parquet files (aci_cost.parquet, blob_cost.parquet, postgres_cost.parquet, databricks_cost.parquet).

This notebook implements the analysis pipeline that turns those raw samples into the tables, plots, and statistical claims required by Chapters 6 (Results) and 7 (Discussion), following the methodology defined in thesis Section 4.4.

Goals

  • Load benchmark results from blob storage into a tidy long-format DataFrame.
  • Compute per-pass descriptive statistics per (query_id, benchmark_run) cell, using median as the primary estimator (Section 4.4.2).
  • Run within-pass paired statistical comparisons with effect sizes (Section 4.4.4).
  • Aggregate across passes for absolute-performance claims with consistency checks (Section 4.4.3).
  • Structure results around the three research questions (Sections 4.4.5–4.4.7).
  • Generate publication-ready tables and plots.

Data shape

Per-iteration sample schema (SchemaVersion V4)

Each data.parquet file contains one row per iteration with these columns:

Column Type Description
status str "success" or "failed"
failure_reason str | None Exception message when status is failed
elapsed_time float | None Wall-clock elapsed time in seconds (None for failed iterations)
network_bytes_sent int Bytes sent during the iteration (from psutil.net_io_counters)
network_bytes_received int Bytes received during the iteration
started_at str (ISO 8601) Iteration start timestamp (UTC)
ended_at str (ISO 8601) Iteration end timestamp (UTC)
cpu_time_user_seconds float User-mode CPU time for the process
cpu_time_system_seconds float System-mode CPU time for the process
result_cardinality int | None Number of rows in the result set (-1 if result is None)
schema_version str Always "v4" for current experiments

Databricks-only columns (all None for ACI-based single-machine benchmarks):

Column Type Description
executor_input_bytes_read int | None Bytes read by Spark executors from storage
executor_run_time_ms int | None Total executor computation time
shuffle_read_bytes int | None Intra-cluster shuffle read volume
shuffle_write_bytes int | None Intra-cluster shuffle write volume
driver_collection_time_ms int | None Time spent collecting results to driver
stage_durations_ms str | None JSON-serialised list of per-stage durations

Hive partition columns (derived from the blob path, not stored in the Parquet file body):

Column Type Description
query_id str Full benchmark identifier, e.g. point-in-polygon-lookup-duckdb-small
run_id str Date-suffixed run identifier, e.g. 2026-05-16-A1B2C3
benchmark_run int Pass number (1, 2, or 3)
iteration int Global iteration number across all passes (iteration + ceiling × (benchmark_run - 1)); pass 2 starts at ceiling + 1, not 1

Derived columns (computed during load)

The query_id encodes three dimensions that must be extracted:

# Example: "point-in-polygon-lookup-duckdb-small"
#   workload_type  = "point-in-polygon-lookup"
#   configuration  = "duckdb"
#   dataset_size   = "small"
#
# Example: "national-scale-spatial-join-databricks-broadcast-8-nodes-large"
#   workload_type  = "national-scale-spatial-join"
#   configuration  = "databricks-broadcast-8-nodes"
#   dataset_size   = "large"
#
# Example: "knn-search-local-small"
#   workload_type  = "knn-search"
#   configuration  = "local"   (= Shapefile)
#   dataset_size   = "small"

A parsing function must handle the variable-length configuration segment. The canonical mapping is defined by benchmarks.yml: each entry has id and dataset_size. Load the YAML at notebook start and use it to look up dataset_size per query_id, then derive workload_type and configuration by stripping the known suffix.

Metadata schema (benchmark_metadata.parquet)

Column Type Description
id str (UUID) Metadata entry identifier
timestamp datetime When the run completed
query_id str Same as in samples
run_id str Same as in samples
achieved_iterations int Successful timed iterations (excludes warmup and failed)
failed_iterations int Iterations that raised an exception
stop_reason str One of: precision, timeout, ceiling, fixed, partial, failed
ci_half_width_seconds float | None Bootstrapped CI half-width on mean elapsed (from stopping rule)
ci_half_width_relative float | None Same as fraction of mean
mean_elapsed_seconds float | None Mean elapsed across timed iterations
median_elapsed_seconds float | None Median elapsed across timed iterations

Cost analytics schema (per *_cost.parquet)

Column Type Description
compute_cost float Compute component (ACI vCPU+RAM, Postgres, or Databricks DBU+VM)
storage_cost float Storage component
network_cost float Network egress/ingress component
operations_cost float Blob operations component
total_cost float Sum of all components

Iteration counts per query type (from BenchmarkIteration)

Query type Ceiling Stopping Warmup
point-in-polygon-lookup 2500 Sequential 5
knn-search 4000 Sequential 5
bbox-filtering 900 Sequential 5
national-scale-spatial-join (ACI) 5 Fixed 1
national-scale-spatial-join (Databricks) 5 Fixed 1

Sequential stopping parameters: min_iterations=10, min_window=60s, max_window=3600s, ci_target=0.05, bootstrap_resamples=1000, confidence=0.95.

Batch structure (from benchmarks.yml)

RQ1 batches (single-machine, 3 query types × 2 sizes):

Workload Small (3-way) Large (2-way)
point-in-polygon-lookup duckdb, postgis, local duckdb, postgis
knn-search duckdb, postgis, local duckdb, postgis
bbox-filtering duckdb, postgis, local duckdb, postgis

RQ2 batches (national-scale spatial join):

Single-node pairs: (duckdb, postgis) at small/medium/large.
Cross-architecture pairs (at 8 nodes, small/medium): (duckdb, postgis, broadcast-8, partitioned-8).
Sedona strategy pairs: (broadcast-N, partitioned-N) at matching node counts.
Sedona scaling pairs: matched node-count pairs within same strategy (e.g. broadcast-4 ↔ broadcast-12 at large).

Analysis pipeline

1. Load and validate

import pandas as pd
import yaml
from azure.storage.blob import ContainerClient

def load_samples(container_client: ContainerClient) -> pd.DataFrame:
    """
    Walk the Hive-partitioned blob tree under the benchmarks container.
    Parse query_id, run_id, benchmark_run, iteration from the blob path.
    Read each data.parquet, attach partition columns, and concatenate.
    """
    ...

def load_metadata(container_client: ContainerClient) -> pd.DataFrame:
    """Download and read benchmark_metadata.parquet from the metadata container."""
    ...

def load_benchmarks_yml(path: str = "benchmarks.yml") -> dict:
    """Load benchmarks.yml and build lookup tables for query_id parsing and batch membership."""
    ...

def parse_query_id(query_id: str, benchmarks: dict) -> tuple[str, str, str]:
    """
    Look up query_id in benchmarks.yml to get dataset_size,
    then derive workload_type and configuration.
    Returns (workload_type, configuration, dataset_size).
    """
    entry = benchmarks[query_id]
    dataset_size = entry["dataset_size"]
    suffix = f"-{dataset_size}"
    stem = query_id.removesuffix(suffix)
    # workload_type is the longest known prefix; configuration is the remainder
    ...

Validation checks:

  • Assert schema_version == "v4" for all rows.
  • For each (query_id, benchmark_run) cell, count successful iterations and compare against achieved_iterations from metadata. Flag any discrepancies.
  • Filter rows where status == "failed". Report count per (query_id, benchmark_run). If any cell loses > 5% of its iterations to failures, label that cell as exploratory.
  • Confirm that within each batch (determined by related_script_ids in benchmarks.yml), all query_ids share the same run_id per pass — this validates concurrent execution and enables paired analysis.
  • Parse started_at/ended_at to datetime. Verify temporal overlap between paired configs within a batch.
  • Validate result_cardinality consistency across configurations within a batch — differing cardinality signals a correctness issue.

2. Descriptive statistics per cell (Section 4.4.2)

For each (query_id, benchmark_run) cell, compute over successful iterations:

Primary estimator: median — median is robust to right-skewed timing distributions and tolerant of one-sided contamination from noisy-neighbour effects (Section 4.4.2, Figure 4.4.2). Mean is reported alongside for completeness but does not anchor claims.

import numpy as np
from scipy.stats import bootstrap as scipy_bootstrap

def descriptive_stats(values: np.ndarray, n_bootstrap: int = 10_000, confidence: float = 0.95) -> dict:
    n = len(values)
    result = scipy_bootstrap(
        data=(values,),
        statistic=np.median,
        n_resamples=n_bootstrap,
        confidence_level=confidence,
        method="BCa",
    )
    return {
        "n": n,
        "median": float(np.median(values)),
        "mean": float(np.mean(values)),
        "std": float(np.std(values, ddof=1)) if n > 1 else np.nan,
        "cv": float(np.std(values, ddof=1) / np.mean(values)) if np.mean(values) > 0 and n > 1 else np.nan,
        "iqr": float(np.percentile(values, 75) - np.percentile(values, 25)),
        "min": float(np.min(values)),
        "p25": float(np.percentile(values, 25)),
        "p75": float(np.percentile(values, 75)),
        "p95": float(np.percentile(values, 95)),
        "max": float(np.max(values)),
        "median_ci_lower": float(result.confidence_interval.low),
        "median_ci_upper": float(result.confidence_interval.high),
    }

Compute independently for each metric:

  • Primary: elapsed_time, network_bytes_received, network_bytes_sent
  • Auxiliary: cpu_time_user_seconds, cpu_time_system_seconds, result_cardinality
  • Distributed phase (Databricks only): executor_input_bytes_read, executor_run_time_ms, shuffle_read_bytes, shuffle_write_bytes, driver_collection_time_ms

Output: one tidy DataFrame per metric, indexed by (query_id, benchmark_run).

3. Cross-pass aggregation for absolute-performance claims (Section 4.4.3)

For each query_id, collect the 3 per-pass medians and compute:

def cross_pass_aggregation(pass_medians: np.ndarray) -> dict:
    """pass_medians: array of shape (3,), one per benchmark_run."""
    headline = float(np.median(pass_medians))
    spread = float(np.ptp(pass_medians))
    return {
        "headline_median": headline,
        "pass_1_median": float(pass_medians[0]),
        "pass_2_median": float(pass_medians[1]),
        "pass_3_median": float(pass_medians[2]),
        "pass_range": spread,
        "pass_range_relative": spread / headline if headline > 0 else np.nan,
        "consistent": (spread / headline <= 0.10) if headline > 0 else False,
    }

Consistency flag (Georges, Buytaert & Eeckhout, 2007): the 3 pass-medians must fall within 10% of the headline median. Cells where consistent == False are labelled exploratory — the cross-pass spread is too wide to support a firm absolute-performance claim.

With only N=3 passes, the honest uncertainty bound is the range itself; no parametric CI is computed on the pass-medians.

4. Within-pass paired comparisons (Section 4.4.4)

For each (workload_type, dataset_size, benchmark_run):

Step 1 — Identify paired configurations from benchmarks.yml via related_script_ids. Pairing is valid because batch members execute concurrently under the same wall-clock window, so iteration i of config A and iteration i of config B experience approximately the same cloud conditions.

Step 2 — Align iterations. Configs in a batch may have different achieved iteration counts due to the sequential stopping rule. Truncate all configs to the minimum successful iteration count within the batch.

Step 3 — Run the appropriate test:

from scipy.stats import wilcoxon, friedmanchisquare
from statsmodels.stats.multitest import multipletests
from itertools import combinations

def vargha_delaney_a12(a: np.ndarray, b: np.ndarray) -> float:
    """
    Vargha & Delaney (2000) A_12: probability that a randomly chosen
    observation from A exceeds one from B.
    """
    m, n = len(a), len(b)
    total = 0.0
    for ai in a:
        total += np.sum(ai > b) + 0.5 * np.sum(ai == b)
    return total / (m * n)

def classify_a12(a12: float) -> str:
    """Arcuri & Briand (2014) thresholds."""
    a = max(a12, 1 - a12)
    if a >= 0.71: return "large"
    if a >= 0.64: return "medium"
    if a >= 0.56: return "small"
    return "negligible"

def pairwise_comparison(
    a: np.ndarray, b: np.ndarray,
    n_bootstrap: int = 10_000, confidence: float = 0.95,
) -> dict:
    """
    Wilcoxon signed-rank on paired samples + A_12 + ratio/difference CIs.
    a and b must be aligned by iteration index and equal length.
    """
    stat, p_value = wilcoxon(a - b, alternative="two-sided")

    a12 = vargha_delaney_a12(a, b)

    ratios = a / np.where(b == 0, np.nan, b)
    ratios = ratios[~np.isnan(ratios)]
    ratio_ci = scipy_bootstrap(
        (ratios,), statistic=np.median,
        n_resamples=n_bootstrap, confidence_level=confidence, method="BCa",
    )

    diffs = a - b
    diff_ci = scipy_bootstrap(
        (diffs,), statistic=np.median,
        n_resamples=n_bootstrap, confidence_level=confidence, method="BCa",
    )

    return {
        "n_paired": len(a),
        "wilcoxon_stat": float(stat),
        "p_value": float(p_value),
        "a12": float(a12),
        "a12_category": classify_a12(a12),
        "ratio_median": float(np.median(ratios)),
        "ratio_ci_lower": float(ratio_ci.confidence_interval.low),
        "ratio_ci_upper": float(ratio_ci.confidence_interval.high),
        "diff_median": float(np.median(diffs)),
        "diff_ci_lower": float(diff_ci.confidence_interval.low),
        "diff_ci_upper": float(diff_ci.confidence_interval.high),
    }

For ≥3 configurations in a batch, run Friedman as the omnibus, then pairwise Wilcoxon signed-rank post-hoc with Holm–Bonferroni correction:

def omnibus_and_posthoc(
    groups: dict[str, np.ndarray],
    benchmark_run: int,
    metric: str,
) -> pd.DataFrame:
    names = sorted(groups.keys())
    n = min(len(v) for v in groups.values())
    aligned = {name: groups[name][:n] for name in names}

    # Friedman omnibus (informational; pairwise tests always run)
    if len(names) >= 3:
        stat, p_omnibus = friedmanchisquare(*[aligned[name] for name in names])
    else:
        stat, p_omnibus = np.nan, np.nan

    # Pairwise Wilcoxon + effect size
    rows = []
    for name_a, name_b in combinations(names, 2):
        result = pairwise_comparison(aligned[name_a], aligned[name_b])
        result["config_a"] = name_a
        result["config_b"] = name_b
        result["benchmark_run"] = benchmark_run
        result["metric"] = metric
        result["friedman_stat"] = float(stat)
        result["friedman_p"] = float(p_omnibus)
        rows.append(result)

    df = pd.DataFrame(rows)

    # Holm–Bonferroni correction within this family
    reject, corrected_p, _, _ = multipletests(df["p_value"], method="holm")
    df["p_value_holm"] = corrected_p
    df["significant"] = reject

    return df

Repeat the entire pipeline for all three primary metrics: elapsed_time, network_bytes_received, network_bytes_sent.

5. Cross-pass consistency for relative claims (Section 4.4.5)

For each pairwise comparison, check whether all 3 passes agree:

def cross_pass_consistency(pairwise_df: pd.DataFrame) -> pd.DataFrame:
    """
    pairwise_df has rows for all (config_a, config_b, benchmark_run, metric).
    Groups by (config_a, config_b, metric) and checks 3-pass agreement.
    """
    def check(group: pd.DataFrame) -> pd.Series:
        ratios = group["ratio_median"].values
        sig = group["significant"].values
        a12s = group["a12"].values

        direction_agrees = all(r > 1 for r in ratios) or all(r < 1 for r in ratios)
        all_significant = bool(np.all(sig))
        all_nontrivial = all(classify_a12(a) != "negligible" for a in a12s)

        return pd.Series({
            "direction_consistent": direction_agrees,
            "all_significant": all_significant,
            "effect_size_consistent": all_nontrivial,
            "fully_consistent": direction_agrees and all_significant and all_nontrivial,
        })

    return pairwise_df.groupby(["config_a", "config_b", "metric"]).apply(check).reset_index()

Reporting rule: A configuration is reported as "faster" / "transferring fewer bytes" only when all 3 passes agree on direction, statistical significance (after Holm–Bonferroni), and at least a small effect size ($\hat{A}_{12} \geq 0.56$ or $\leq 0.44$). Pairs where passes disagree are themselves a finding — they indicate cloud variability dominates the configuration effect for that query — and are reported as such.

6. Outlier handling (Chen & Revels, 2016)

Default: keep all measurements. Report robust statistics (median, IQR). Do not apply parametric outlier detection — the 3-sigma rule is invalid for non-i.i.d. timing measurements (Chen & Revels, 2016).

If an outlier sensitivity check is desired:

  • Flag values above the 99th percentile within a cell (non-parametric criterion).
  • Run the full pipeline twice: once with all data, once excluding flagged values.
  • If conclusions differ, the divergence is the headline finding.
  • Report count and percentage removed per cell.

7. RQ-specific analysis

RQ1: Cloud-native vs. traditional (Section 4.4.5)

Scope: point-in-polygon-lookup, knn-search, bbox-filtering across small and large.

For each (workload_type, dataset_size):

  • 3-way Friedman + post-hoc at small (duckdb, postgis, local/shapefile).
  • 2-way Wilcoxon at large (duckdb, postgis).
  • Report: median ratio, $\hat{A}_{12}$, Holm-corrected p-value, cross-pass consistency.
  • Ranking table: which configuration is fastest / transfers fewest bytes, with consistency flag.

RQ2: Scaling against single-node (Section 4.4.6)

Scope: national-scale-spatial-join.

  • At each (dataset_size, worker_count): compare Sedona strategy vs DuckDB, vs PostGIS.
  • Scaling curve: median elapsed_time vs worker count per strategy, with cross-pass range as error bars.
  • Speedup ratio: median(single_node) / median(distributed) with bootstrap CI.
  • Within-Sedona: Friedman + post-hoc (broadcast vs partitioned vs default) at matched (dataset_size, worker_count).

RQ3: Consistency across patterns (Section 4.4.7)

  • For each configuration pair, check whether the ranking (faster/slower) is consistent across all workload types tested in RQ1.
  • Cross-workload agreement table: one row per configuration pair, one column per workload type, cells show direction + significance.

8. Visualisation

For each (workload_type, dataset_size):

  1. Box-and-whisker plots of elapsed_time per (configuration, benchmark_run) — passes side by side, configurations grouped. Makes per-pass consistency visually obvious.
  2. Same for network_bytes_received and network_bytes_sent.
  3. Scatter: elapsed_time vs network_bytes_received per configuration (all iterations, all passes) — shows CPU-bound vs network-bound separation.
  4. Time-series diagnostic: iteration index vs elapsed_time within a pass, coloured by configuration — exposes drift, warm-up artefacts, or noisy-neighbour episodes. Appendix material.

For RQ2:
5. Scaling curve: median elapsed_time vs worker count, one line per Sedona strategy, with cross-pass range as error bars. DuckDB and PostGIS plotted as horizontal reference lines.

Style: matplotlib with consistent fonts/sizes for thesis figures. Export as PDF for LaTeX \includegraphics.

9. Reporting tables

Table 1 — Descriptive statistics per cell.
One row per (query_id, benchmark_run, metric).
Columns: n, median, mean, std, cv, iqr, min, p25, p75, p95, max, median_ci_lower, median_ci_upper.

Table 2 — Cross-pass medians and consistency.
One row per (query_id, metric).
Columns: headline_median, pass_1_median, pass_2_median, pass_3_median, pass_range, pass_range_relative, consistent.

Table 3 — Pairwise comparisons.
One row per (config_a, config_b, workload_type, dataset_size, benchmark_run, metric).
Columns: n_paired, ratio_median, ratio_ci_lower, ratio_ci_upper, diff_median, diff_ci_lower, diff_ci_upper, wilcoxon_stat, p_value, p_value_holm, significant, a12, a12_category, friedman_stat, friedman_p.

Table 4 — Cross-pass consistency for relative claims.
One row per (config_a, config_b, workload_type, dataset_size, metric).
Columns: direction_consistent, all_significant, effect_size_consistent, fully_consistent.

All tables exportable via pandas.DataFrame.to_latex(escape=True).

Methodology references

Cite in the notebook header:

  • Arcuri & Briand (2014) — Wilcoxon signed-rank with $\hat{A}_{12}$; test choice and effect-size thresholds.
  • Chen & Revels (2016) — non-i.i.d. timing measurements; rationale for non-parametric methods and against 3-sigma outlier removal.
  • Georges, Buytaert & Eeckhout (2007) — CI-based reporting; precision-based "exploratory" labelling.
  • Hoefler & Belli (2015) — twelve rules for reporting parallel-system performance; report distributions, not just means.
  • Iosup et al. (2011), Schad et al. (2010), Leitner & Cito (2016), Laaber et al. (2019) — cloud variability; rationale for per-pass analysis and paired comparisons.
  • JCGM 200:2012 (VIM) — measurement vocabulary.
  • Kalibera & Jones (2013) — multi-level experiment design; framework for the 3-pass × N-iteration structure.
  • Raasveldt et al. (2018) — fair benchmarking pitfalls; rationale for warm-up and result-set materialisation control.
  • Vargha & Delaney (2000) — original $\hat{A}_{12}$ definition.

Acceptance criteria

  • Raw results from blob storage load into a single tidy DataFrame with derived workload_type, configuration, dataset_size columns, validated for schema V4 and completeness against metadata.
  • Metadata from benchmark_metadata.parquet loaded, cross-referenced with per-cell sample counts, and stop-reason distribution reported.
  • Tables 1–4 produced for primary metrics (elapsed_time, network_bytes_received, network_bytes_sent).
  • Tables 1–2 also produced for auxiliary metrics (cpu_time_user_seconds, cpu_time_system_seconds) and Databricks phase metrics where applicable.
  • Box plots, scatter plots, and scaling curves produced for every (workload_type, dataset_size) combination.
  • Pairwise comparisons report Wilcoxon signed-rank (or Friedman + post-hoc), Holm-corrected p-values, and $\hat{A}_{12}$ with effect-size categories.
  • Cross-pass consistency flags computed for every absolute and relative claim.
  • RQ-specific analysis sections (RQ1, RQ2, RQ3) produce the results needed for thesis Chapters 6 and 7.
  • LaTeX tables exportable directly into the Results chapter via to_latex.
  • Notebook header documents the methodology and cites all 9 sources above.
  • Re-running the notebook from scratch on the raw results reproduces every figure and table without manual intervention.

Open questions

  • Decide whether to include a 4th or 5th pass for benchmark groups where cross-pass spread exceeds 10%. Cost: each extra pass adds ~24 hours of ACI time + Databricks cluster costs.
  • Decide whether to report cold-start performance separately. Currently excluded by warmup. If ACI container provisioning latency is interesting, capture from ACI API metadata rather than iteration measurements.
  • Confirm whether result_cardinality consistency should be asserted within each batch (differing cardinality = potential correctness bug).

Out of scope

  • Running the actual benchmark experiment (handled by the framework).
  • Methodology-chapter prose (separate issue).
  • Discussion-chapter interpretation (separate issue).
  • Cross-provider comparison (Azure-only is a stated limitation of the thesis).

Metadata

Metadata

Assignees

Labels

enhancementNew feature or request

Type

No type

Projects

No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions