Metadata-Version: 2.4
Name: dpv-calibration
Version: 0.1.2
Summary: Univariate DPV calibration: OLS fit, inverse-regression prediction, LOD/LOQ.
Project-URL: Repository, https://github.com/danvu772/dpv-calibration
License: MIT License
        
        Copyright (c) 2026 Dan Vu
        
        Permission is hereby granted, free of charge, to any person obtaining a copy
        of this software and associated documentation files (the "Software"), to deal
        in the Software without restriction, including without limitation the rights
        to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
        copies of the Software, and to permit persons to whom the Software is
        furnished to do so, subject to the following conditions:
        
        The above copyright notice and this permission notice shall be included in all
        copies or substantial portions of the Software.
        
        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
        IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
        FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
        AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
        LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
        OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
        SOFTWARE.
License-File: LICENSE
Requires-Python: >=3.10
Requires-Dist: dpv-monopeak>=0.1.1
Requires-Dist: numpy>=1.24
Requires-Dist: openpyxl>=3.1
Requires-Dist: pandas>=2.0
Requires-Dist: scikit-learn>=1.3
Requires-Dist: scipy>=1.10
Provides-Extra: plot
Requires-Dist: matplotlib>=3.7; extra == 'plot'
Description-Content-Type: text/markdown

# dpv-calibration

Univariate OLS calibration for Differential Pulse Voltammetry biosensor data: fit a calibration line from standard scans at known concentrations, predict analyte concentration from new scans via inverse regression, and receive 95% confidence intervals with LOD/LOQ classification.

Depends on [dpv-monopeak](https://github.com/danvu772/dpv-monopeak) for peak feature extraction from raw voltage/current arrays.

Developed by the Ray Research Lab at UH Manoa.

---

## Installation

```bash
pip install dpv-calibration
```

For `plot_calibration` (requires matplotlib):

```bash
pip install "dpv-calibration[plot]"
```

---

## Quick start

```python
import numpy as np
from dpv_calibration import DPVScan, CalibrationModel, plot_calibration

# Construct calibration scans (v in volts, i in µA, concentration_M in mol/L)
scans = [
    DPVScan(v=v_arr, i=i_arr, concentration_M=1e-9),
    DPVScan(v=v_arr, i=i_arr, concentration_M=1e-8),
    # one DPVScan per replicate per concentration level
]

# Fit (log scale is auto selected when dynamic range exceeds 1000x)
model = CalibrationModel(feature="Height")
model.fit(scans)
print(model)
# CalibrationModel(feature='Height', log_x=True, R²=0.9821, RMSECV=0.3104,
#                  LOD=3.14e-13 M, LOQ=1.05e-12 M)

# Predict concentration from a new raw scan
result = model.predict(v=new_v, i=new_i)
print(result.concentration_M)  # float, mol/L
print(result.status)            # "below_lod" | "detected" | "quantifiable"
print(result.ci)                # (lower_M, upper_M), 95% CI in mol/L

# Plot the calibration (requires matplotlib)
ax = plot_calibration(model)
ax.figure.savefig("calibration.png", dpi=150, bbox_inches="tight")
```

---

## API reference

### `DPVScan`

```python
@dataclass
class DPVScan:
    v: np.ndarray           # voltage array, shape (N,), volts
    i: np.ndarray           # current array, shape (N,), µA
    concentration_M: float  # analyte concentration, mol/L, must be > 0
```

Plain data container. No processing occurs at construction. `v` and `i` must be the same length. Current units must be consistent across every scan in a calibration set.

---

### `PredictionResult`

```python
@dataclass
class PredictionResult:
    concentration_M: float  # predicted analyte concentration, mol/L
    status: str             # detection classification, see table below
    ci: tuple               # (lower_M, upper_M), 95% CI in mol/L
```

**`status` values:**

| Value | Condition | Meaning |
|---|---|---|
| `"below_lod"` | `concentration_M < model.lod_M` | Signal is indistinguishable from blank noise. Do not report a numeric concentration. |
| `"detected"` | `model.lod_M <= concentration_M < model.loq_M` | Analyte presence is confirmed but concentration cannot be reliably quantified. Use for presence or absence decisions only. |
| `"quantifiable"` | `concentration_M >= model.loq_M` | Result is reliable for quantitative reporting. |

`ci` is a `(lower_M, upper_M)` tuple in mol/L. Both bounds are in linear molar units regardless of whether the model was fitted in log space. The interval widens when the predicted concentration is far from the center of the calibration range or when `n` is small.

---

### `CalibrationModel`

```python
CalibrationModel(feature: str = "Height", log_x: bool | None = None)
```

**Constructor parameters:**

| Parameter | Type | Default | Description |
|---|---|---|---|
| `feature` | `str` | `"Height"` | The dpv_monopeak feature used as the calibration signal `y`. Valid values: `"Height"`, `"Area"`, `"Width"`, `"X"`, `"YOffset"`, `"MaxSlope"`, `"MinSlope"`, `"SumSlope"`. |
| `log_x` | `bool or None` | `None` | Whether to fit in log10 concentration space. `None` triggers auto detection: `True` when `max(concentration) / min(concentration) > 1000`, otherwise `False`. Override by passing `True` or `False` explicitly. |

**Public attributes available after `.fit()`:**

| Attribute | Type | Unit | Description |
|---|---|---|---|
| `feature` | `str` | | Feature name set at construction. |
| `log_x` | `bool` | | Whether the model uses log10 concentration space. Resolved from `None` after fit. |
| `r2` | `float` | | Coefficient of determination of the OLS fit on training data. |
| `rmsecv` | `float` | log10(M) when `log_x=True`, M when `log_x=False` | Root mean square error of leave one out cross validation in the same units as the fitted `x`. |
| `lod_M` | `float` | mol/L | Limit of detection. Smallest concentration distinguishable from noise (3 standard deviation rule). |
| `loq_M` | `float` | mol/L | Limit of quantification. Smallest concentration reliably quantifiable (10 standard deviation rule). Always greater than `lod_M`. |

**`repr` examples:**

```
# Fitted
CalibrationModel(feature='Height', log_x=True, R²=0.9821, RMSECV=0.3104,
                 LOD=3.14e-13 M, LOQ=1.05e-12 M)

# Not yet fitted
CalibrationModel(feature='Height', not fitted)
```

#### `.fit(scans: list[DPVScan]) -> CalibrationModel`

Fit the calibration model. Returns `self` to support chaining: `model.fit(scans).predict(v, i)`.

Internally:
1. Calls `dpv_monopeak.extract_dpv_features(scan.v, scan.i)` for each scan. Scans where no valid Faradaic peak is found (None return) are silently skipped.
2. Resolves `log_x` if it was `None`.
3. Fits OLS in concentration vs feature space (optionally in log scale).
4. Computes and stores `r2`, `rmsecv`, `lod_M`, `loq_M`.

Multiple replicates at the same concentration are supported. Each `DPVScan` is an independent observation.

#### `.predict(v, i, m=1, alpha=0.05) -> PredictionResult`

```python
model.predict(
    v: np.ndarray,        # voltage array of the unknown scan, volts
    i: np.ndarray,        # current array of the unknown scan, µA
    m: int = 1,           # number of replicates averaged to produce this scan
    alpha: float = 0.05,  # CI significance level; 0.05 gives 95% CI, 0.01 gives 99%
) -> PredictionResult
```

Predicts concentration from a raw DPV scan using inverse regression.

Raises `RuntimeError` if called before `.fit()`. Raises `ValueError` if no valid peak is found in the scan.

**`m` (replicates):** if the scan is the mean of `m` replicate measurements, pass `m` accordingly. CI width scales as `1 / sqrt(m)`. Pass `m=1` for a single scan.

**How inverse regression works:** standard OLS predicts feature from concentration (`y = b0 + b1 * x`). Prediction reverses this by solving for `x`: `x0_hat = (y0 - b0) / b1`. This is not a reverse regression of x on y. It correctly propagates calibration uncertainty into the CI. When `log_x=True`, `x0_hat` is in log10(M) and is back transformed to mol/L before being stored in `PredictionResult.concentration_M`.

---

### `plot_calibration`

```python
from dpv_calibration import plot_calibration  # requires [plot] extra

plot_calibration(
    model: CalibrationModel,
    ax=None,              # matplotlib Axes to draw into; creates a new figure if None
    alpha: float = 0.05,  # significance level for the confidence band
) -> matplotlib.axes.Axes
```

Raises `ImportError` if matplotlib is not installed. Raises `RuntimeError` if the model is not fitted.

**Elements drawn:**

| Element | Description |
|---|---|
| Scatter | Training points in the fitted x space |
| OLS line | Fitted line over the calibration range plus 10% padding on each side |
| Confidence band | Shaded region for the mean response at each x value (Lavagnini & Magno 2007, Eq. 9b) |
| Orange dashed line | LOD position in fitted x space |
| Red dashed line | LOQ position in fitted x space |
| Text annotation | R² and RMSECV in the lower right corner |

**X axis:** `"log₁₀(Concentration / M)"` when `model.log_x=True`; `"Concentration (M)"` when `False`.

**Y axis:** value of `model.feature`.

---

## Loading your data

The following is a complete, ready to paste function for **PSTrace multi-curve Excel exports** from PalmSens instruments. Adapt it to your own file format. The only requirement for the rest of the API is a list of `DPVScan` objects with correct `v`, `i`, and `concentration_M` values.

```python
import re
import numpy as np
import pandas as pd
from dpv_calibration import DPVScan


def load_palmsens_xlsx(
    path: str,
    keyword: str = "cort",
    conc_pattern: str = r"(\d+(?:\.\d+)?)\s*(pM|nM|uM|mM)",
) -> list:
    """
    Load DPV scans from a PSTrace multi-curve Excel export.

    Expected file layout:
        Row 0:  scan name  | (blank)  | scan name  | (blank) | ...
        Row 1:  "V"        | "uA"     | "V"        | "uA"    | ...
        Row 2+: voltage    | current  | voltage    | current | ...

    Columns come in pairs (voltage, current). Scans must have exactly 81
    data rows after dropping NaN values (standard PSTrace DPV export length).

    Parameters
    ----------
    path         : path to the .xlsx file
    keyword      : case insensitive substring filter on scan names;
                   use "" to load all scans
    conc_pattern : regex with two capture groups (numeric value, unit);
                   applied to each scan name to extract concentration;
                   scans with no match are skipped
    """
    _units = {"pm": 1e-12, "nm": 1e-9, "um": 1e-6, "mm": 1e-3}

    def _parse_conc(label):
        m = re.search(conc_pattern, label, re.IGNORECASE)
        if not m:
            return None
        value, unit = float(m.group(1)), m.group(2).lower()
        return value * _units[unit]

    raw = pd.read_excel(path, header=None)
    scans = []
    for col in range(0, raw.shape[1], 2):
        name = str(raw.iloc[0, col])
        if keyword.lower() not in name.lower():
            continue
        data = raw.iloc[2:, col:col + 2].dropna()
        data.columns = ["V", "I"]
        data = data.apply(pd.to_numeric, errors="coerce").dropna()
        if len(data) != 81:
            continue
        conc = _parse_conc(name)
        if conc is None:
            continue
        scans.append(DPVScan(
            v=data["V"].values,
            i=data["I"].values,
            concentration_M=conc,
        ))
    return scans
```

**Concentration labels recognised by the default regex** (case insensitive, first match in the label name):

| Label | Parsed concentration |
|---|---|
| `"1pM Cort 5mM redox det"` | `1e-12` M |
| `"10nM cort"` | `1e-8` M |
| `"0.5uM Cortisol"` | `5e-7` M |
| `"1mM Cort PBS det"` | `1e-3` M |

---

## Mathematical reference

> Full derivations with LaTeX: [MATH.md](https://github.com/danvu772/dpv-calibration/blob/main/MATH.md)

All equations follow Lavagnini & Magno (2007).

**Variables:** `n` = number of calibration points after skipping invalid scans. `xi` = concentration of point i, or `log10(concentration)` when `log_x=True`. `yi` = feature value of point i. `x_mean` = mean of all `xi`. `Sxx = sum((xi - x_mean)^2)`.

### OLS fit

```
b1   = sum((xi - x_mean) * yi) / Sxx
b0   = y_mean - b1 * x_mean
sy_x = sqrt(sum((yi - (b0 + b1 * xi))^2) / (n - 2))
```

`sy_x` is the residual standard deviation. It measures scatter around the fitted line and is the basis for LOD, LOQ, and all confidence intervals. The denominator `n - 2` reflects two degrees of freedom used to estimate `b0` and `b1`.

### Inverse regression

Given a new feature value `y0` extracted from an unknown scan:

```
x0_hat = (y0 - b0) / b1
```

Standard error of `x0_hat`, accounting for calibration fit uncertainty and the distance of the prediction from the calibration center:

```
se = (sy_x / abs(b1)) * sqrt(1/m + 1/n + (x0_hat - x_mean)^2 / Sxx)
```

where `m` is the number of replicates averaged in the unknown scan.

The `(1 - alpha) * 100%` confidence interval:

```
t     = t_distribution.ppf(1 - alpha/2, df=n - 2)
ci_lo = x0_hat - t * se
ci_hi = x0_hat + t * se
```

When `log_x=True`, back transform all results to mol/L:

```
concentration_M = 10^x0_hat
ci = (10^ci_lo, 10^ci_hi)
```

### Confidence band in `plot_calibration`

The shaded region shows the 95% confidence interval for the mean response at each x value (Lavagnini & Magno 2007, Eq. 9b):

```
half_band(x) = t * sy_x * sqrt(1/n + (x - x_mean)^2 / Sxx)
upper(x) = b0 + b1*x + half_band(x)
lower(x) = b0 + b1*x - half_band(x)
```

### LOD and LOQ

Computed in fitted x space, back transformed when `log_x=True`:

```
x_lod = 3  * sy_x / b1    # limit of detection in fitted x space
x_loq = 10 * sy_x / b1    # limit of quantification in fitted x space
```

```
# when log_x=False:
lod_M = x_lod
loq_M = x_loq

# when log_x=True:
lod_M = 10^x_lod
loq_M = 10^x_loq
```

LOD corresponds to a net signal of `3 * sy_x` above blank. Below this threshold, signal cannot be distinguished from noise. LOQ uses factor 10 and marks the threshold for reliable quantitative accuracy. These map to `PredictionResult.status` as follows: `concentration_M < lod_M` gives `"below_lod"`, `lod_M <= concentration_M < loq_M` gives `"detected"`, `concentration_M >= loq_M` gives `"quantifiable"`.

### RMSECV

Leave one out cross validation in fitted x space:

```
for each point i in 1..n:
    fit b1_loo, b0_loo on all points except i  (same OLS formulas)
    x_pred_i = (yi - b0_loo) / b1_loo
    error_i  = (xi - x_pred_i)^2

RMSECV = sqrt(mean(error_i for i in 1..n))
```

When `log_x=True`, RMSECV is in units of log10(M). An RMSECV of 0.5 means predictions are off by half a decade on average. RMSECV is not back transformed because it is a symmetric error metric defined in fitting space.

### R squared

```
R2 = 1 - sum((yi - y_hat_i)^2) / sum((yi - y_mean)^2)
```

Reported as `model.r2` for convention. Lavagnini & Magno (2007) caution that R² near 1.0 does not guarantee a valid linear model. The primary diagnostics are `sy_x` and RMSECV.

---

## Usage with AI assistants

This README is designed to be pasted directly into an LLM coding assistant as context. All types, units, parameter names, return values, and mathematical definitions are documented explicitly so that an LLM can generate correct, complete code against this API without access to the source.

---

## Citation

Lavagnini, I. & Magno, F. (2007). A statistical overview on univariate calibration, inverse regression, and detection limits: application to gas chromatography/mass spectrometry technique. *Mass Spectrometry Reviews*, 26, 1-18. https://doi.org/10.1002/mas.20100

---

## License

MIT © 2026 Dan Vu
