# Python guide for the Difference Significance Profile (DSP) dataset

## Purpose

This document contains the code-oriented instructions for loading, checking, visualising, and reproducing the supplied DSP outputs. For the dataset overview, experimental devices, analysis settings, naming convention, folder structure, and detailed CSV descriptions, see [`README.md`](README.md).

## Recommended software environment

The examples require Python 3.10 or later and the following packages:

```bash
python -m pip install numpy pandas matplotlib
```

## Important filename convention

Processed files end with an alpha code:

| Suffix | Significance level | Simultaneous level |
|---:|---:|---:|
| `001` | α = 0.01 | 99% |
| `005` | α = 0.05 | 95% |
| `010` | α = 0.10 | 90% |

The `005` files represent the principal 95% simultaneous analysis.

---

## How to use the data

## Minimal Python workflow

The following example reads the main 95% analysis for the different-temperature comparison and prints the significant bins.

```python
from pathlib import Path
import pandas as pd

root = Path("data")
comparison = root / "T1_T2"

summary = pd.read_csv(
    comparison / "experimental_T120_1_T145_1_summary_stats_005.csv"
)

dsp = pd.read_csv(
    comparison / "experimental_T120_1_T145_1_dsp_per_bin_005.csv"
)

significant = dsp.loc[
    dsp["sig_flag"].eq(1),
    ["bin_left", "bin_right", "bin_mid", "Delta", "SE", "Z"]
]

print(summary.to_string(index=False))
print(significant.to_string(index=False))
```

## Reconstructing the PSD overlay

```python
import pandas as pd
import matplotlib.pyplot as plt

psd = pd.read_csv(
    "data/T1_T2/experimental_T120_1_T145_1_psd_bins_005.csv"
)

width = psd["bin_right"] - psd["bin_left"]

plt.bar(
    psd["bin_mid"], psd["p1hat"], width=width,
    alpha=0.5, label="120 °C"
)
plt.bar(
    psd["bin_mid"], psd["p2hat"], width=width,
    alpha=0.5, label="145 °C"
)
plt.xlabel("Particle size [µm]")
plt.ylabel("Probability")
plt.legend()
plt.tight_layout()
plt.show()
```

## Reconstructing the DSP difference plot

```python
import pandas as pd
import matplotlib.pyplot as plt

x = pd.read_csv(
    "data/T1_T2/experimental_T120_1_T145_1_dsp_per_bin_005.csv"
)

width = x["bin_right"] - x["bin_left"]
limit = x["crit_boot"] * x["SE"]

plt.fill_between(
    x["bin_mid"], -limit, limit,
    alpha=0.3, label="Simultaneous reference band"
)

nonsig = x["sig_flag"].eq(0)
sig = x["sig_flag"].eq(1)

plt.bar(
    x.loc[nonsig, "bin_mid"], x.loc[nonsig, "Delta"],
    width=width.loc[nonsig], label="Not significant"
)
plt.bar(
    x.loc[sig, "bin_mid"], x.loc[sig, "Delta"],
    width=width.loc[sig], label="Significant"
)
plt.axhline(0, linewidth=1)
plt.xlabel("Particle size [µm]")
plt.ylabel("Probability difference, Δ")
plt.legend()
plt.tight_layout()
plt.show()
```

## Loading the covariance matrix

```python
import pandas as pd

cov = pd.read_csv(
    "data/T1_T2/experimental_T120_1_T145_1_cov_matrix_005.csv"
)

assert cov.shape == (75, 75)
```
---

## Reproducing the analysis conceptually

The supplied outputs can be reconstructed from the raw particle files using the following workflow:

1. Read the two positive particle-size vectors.
2. Pool the two samples and define one common binning.
3. Count particles from each sample in every bin.
4. Calculate the empirical probabilities `p1hat` and `p2hat`.
5. Calculate the signed difference `Delta = p1hat - p2hat`.
6. Under the null hypothesis of equal PSDs, calculate the pooled probability in each bin.
7. Calculate the standard error and local statistic `Z`.
8. Generate paired multinomial bootstrap count vectors under the pooled null distribution.
9. For every bootstrap iteration, record the maximum absolute `Z` value across all bins.
10. Use the `(1 - alpha)` quantile of these maxima as the simultaneous threshold.
11. Mark bins with `abs(Z) > threshold`.
12. Calculate the global Pearson statistic and its bootstrap p-value.

The supplied analysis uses:

- 75 common bins for both comparisons;
- 40,000 parametric-bootstrap iterations;
- simultaneous levels of 99%, 95%, and 90%; and
- a consistent random-number sequence within each comparison, as evidenced by identical `M_star` files across alpha levels.

For exact reproduction, use the same bin edges, random-number seed, bootstrap algorithm, and quantile definition as the original analysis code.

---

## Data integrity checks

The following checks are useful after downloading or moving the repository:

```python
from pathlib import Path
import numpy as np
import pandas as pd

root = Path("data")

# Raw sample sizes
assert len(pd.read_csv(root / "SOPAT/experimental_T120_1.csv")) == 2881
assert len(pd.read_csv(root / "SOPAT/experimental_T120_11.csv")) == 3015
assert len(pd.read_csv(root / "SOPAT/experimental_T145_1.csv")) == 2658

# Main per-bin output
x = pd.read_csv(
    root / "T1_T2/experimental_T120_1_T145_1_dsp_per_bin_005.csv"
)
assert len(x) == 75
assert x["n1i"].sum() == 2881
assert x["n2i"].sum() == 2658
assert np.isclose(x["p1hat"].sum(), 1.0)
assert np.isclose(x["p2hat"].sum(), 1.0)
assert np.isclose(x["Delta"].sum(), 0.0, atol=1e-12)
assert (x["sig_flag"] == (x["Z"].abs() > x["crit_boot"])).all()

# Bootstrap length
m = pd.read_csv(
    root / "T1_T2/experimental_T120_1_T145_1_bootstrap_M_005.csv"
)
assert len(m) == 40000

# Covariance dimensions
cov = pd.read_csv(
    root / "T1_T2/experimental_T120_1_T145_1_cov_matrix_005.csv"
)
assert cov.shape == (75, 75)
```

---

## Notes for porosimetry files

The porosimetry exports use decimal commas and may require explicit import handling. Inspect the first rows before automated processing. The porosimeter model is not encoded in the supplied CSV files.

```python
from pathlib import Path

for filename in [
    "data/T1_T2/experimental_T120_porosity.csv",
    "data/T1_T2/experimental_T145_porosity.csv",
]:
    path = Path(filename)
    print(f"\n--- {path.name} ---")
    with path.open("r", encoding="utf-8-sig", errors="replace") as handle:
        for _ in range(5):
            print(handle.readline().rstrip())
```