'''
HYSTERESIS ANALYSIS TOOLBOX
==========================
Hysteresis index calculation and visualization for C-Q analysis of contaminants in watersheds
This toolbox provides comprehensive methods for analyzing concentration-discharge (C-Q)
relationships and hysteresis patterns in environmental data, based on established scientific
frameworks.
CORE REFERENCES:
- Lloyd et al. (2016) HESS 20:625-632 - Improved hysteresis index
- Zuecco et al. (2016) Hydrol Process 30:1449-1466 - Versatile index for event-scale analysis
- Williams (1989) J Hydrol 111:89-106 - Sediment hysteresis fundamentals
- Evans & Davies (1998) Water Resour Res 34:129-137 - C-Q hysteresis classification
- Vaughan et al. (2017) Water Resour Res 53:5345-5363 - Storm-scale modifications
- Wymore et al. (2019) Front Earth Sci 7:126 - Tropical watershed applications
- Peña et al. (2023) J Hydrol 626:130262 - HARP method for complex patterns
- Liu et al. (2021) Water Res 200:117254 - Pattern classification framework
Version 2.0 (2025-07-08)
(cc) conrad.jackisch@tbt.tu-freiberg.de (supported by Claude Opus 4)
USAGE EXAMPLES:
--------------
# Single compound analysis
fig = create_hysteresis_plot(
data, sites=['C4', 'A12'],
ccol='Cd_mgL', qcol='Q_mLs',
compound='Cd', conc_unit='mg L⁻¹', flow_unit='mL s⁻¹'
)
# Multi-compound comparison with compound-specific ranges
fig = create_multi_compound_hysteresis_plot(
data, sites=['C4', 'A12', 'B4', 'B3'],
ccols=['Cd_mgL', 'Zn_mgL', 'MPI'],
compounds=['Cd', 'Zn', 'MPI'],
conc_units=['mg L⁻¹', 'mg L⁻¹', '-'],
qcol='Q_mLs', flow_unit='mL s⁻¹'
)
# Summary statistics
stats = create_hysteresis_summary_stats(
data, sites=['C4'], ccol='Cd_mgL', qcol='Q_mLs'
)
# Timeline visualization
timeline = create_hysteresis_timeline(
data, sites=['C4', 'A12'], ccol='Cd_mgL', qcol='Q_mLs'
)
'''
import pandas as pd
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from datetime import datetime
from typing import List, Tuple, Dict, Optional, Union
# ==============================================================================
# CORE SCIENTIFIC METHODS
# ==============================================================================
def calculate_hysteresis_index(
flow_start: float,
flow_end: float,
conc_start: float,
conc_end: float,
flow_range: Optional[Tuple[float, float]] = None,
conc_range: Optional[Tuple[float, float]] = None,
method: str = 'lloyd'
) -> float:
"""
Calculate hysteresis index using established scientific methods.
Parameters
----------
flow_start, flow_end : float
Discharge values at segment start and end
conc_start, conc_end : float
Concentration values at segment start and end
flow_range : tuple of float, optional
(min, max) flow values for normalization
conc_range : tuple of float, optional
(min, max) concentration values for normalization
method : {'lloyd', 'zuecco', 'harp'}
Calculation method to use
Returns
-------
float
Hysteresis index value
- Positive: Clockwise hysteresis (concentration leads discharge)
- Negative: Counter-clockwise (concentration lags discharge)
- Near zero: Minimal hysteresis
Notes
-----
Lloyd method: HI = ΔC_norm - ΔQ_norm
Zuecco method: Angle-based approach with magnitude weighting
HARP method: Area-based calculation
"""
# Handle edge cases
if flow_start == flow_end and conc_start == conc_end:
return 0.0
# Normalize values if ranges provided
if flow_range and conc_range:
flow_min, flow_max = flow_range
conc_min, conc_max = conc_range
# Protect against zero ranges
flow_delta = flow_max - flow_min if (flow_max - flow_min) > 1e-10 else 1
conc_delta = conc_max - conc_min if (conc_max - conc_min) > 1e-10 else 1
flow_start_norm = (flow_start - flow_min) / flow_delta
flow_end_norm = (flow_end - flow_min) / flow_delta
conc_start_norm = (conc_start - conc_min) / conc_delta
conc_end_norm = (conc_end - conc_min) / conc_delta
else:
flow_start_norm = flow_start
flow_end_norm = flow_end
conc_start_norm = conc_start
conc_end_norm = conc_end
delta_flow = flow_end_norm - flow_start_norm
delta_conc = conc_end_norm - conc_start_norm
if method == 'lloyd':
# Lloyd et al. (2016): Simple normalized difference
return delta_conc - delta_flow
elif method == 'zuecco':
# Zuecco et al. (2016): Angle-based with magnitude
angle = np.arctan2(delta_conc, delta_flow)
magnitude = np.sqrt(delta_flow**2 + delta_conc**2)
return magnitude * np.sin(angle - np.pi/4)
elif method == 'harp':
# HARP method: Area-based (simplified)
return (flow_start_norm * conc_end_norm) - (flow_end_norm * conc_start_norm)
else:
raise ValueError(f"Unknown method: {method}")
def calculate_harp_metrics(
flow: np.ndarray,
conc: np.ndarray,
time: Optional[np.ndarray] = None
) -> Dict[str, float]:
"""
Calculate comprehensive HARP metrics for hysteresis characterization.
Based on Peña et al. (2023) J Hydrol 626:130262
Parameters
----------
flow : array-like
Discharge time series
conc : array-like
Concentration time series
time : array-like, optional
Time values for peak timing calculation
Returns
-------
dict
Dictionary containing:
- hysteresis_area: Normalized loop area
- residual: End-state deviation
- peak_timing: Time lag between peaks (hours)
- loop_width: Maximum loop width
"""
flow = np.asarray(flow)
conc = np.asarray(conc)
# Normalize to [0, 1]
flow_norm = (flow - flow.min()) / (flow.max() - flow.min()) if flow.max() > flow.min() else flow
conc_norm = (conc - conc.min()) / (conc.max() - conc.min()) if conc.max() > conc.min() else conc
# 1. Hysteresis area (shoelace formula)
n = len(flow_norm)
area = 0.0
for i in range(n):
j = (i + 1) % n
area += flow_norm[i] * conc_norm[j]
area -= flow_norm[j] * conc_norm[i]
hysteresis_area = abs(area) / 2.0
# 2. Residual
residual = conc_norm[-1] - conc_norm[0]
# 3. Peak timing
flow_peak_idx = np.argmax(flow)
conc_peak_idx = np.argmax(conc)
if time is not None:
time = np.asarray(time)
peak_timing = (time[conc_peak_idx] - time[flow_peak_idx]).total_seconds() / 3600
else:
peak_timing = conc_peak_idx - flow_peak_idx
# 4. Loop width
loop_width = 0.0
for c_level in np.linspace(conc_norm.min(), conc_norm.max(), 20):
tolerance = 0.05
mask = np.abs(conc_norm - c_level) < tolerance
if np.sum(mask) >= 2:
flows_at_level = flow_norm[mask]
width = flows_at_level.max() - flows_at_level.min()
loop_width = max(loop_width, width)
return {
'hysteresis_area': hysteresis_area,
'residual': residual,
'peak_timing': peak_timing,
'loop_width': loop_width
}
def classify_cq_behavior(
flow_diff: float,
conc_diff: float,
flow_range: Tuple[float, float],
conc_range: Tuple[float, float],
threshold_factor: float = 0.01
) -> str:
"""
Classify C-Q behavior based on segment changes.
Based on Williams (1989) and Evans & Davies (1998)
Parameters
----------
flow_diff : float
Change in flow between points
conc_diff : float
Change in concentration between points
flow_range : tuple
(min, max) flow values for significance testing
conc_range : tuple
(min, max) concentration values for significance testing
threshold_factor : float
Relative threshold for significant change
Returns
-------
str
Behavior classification:
- 'flushing': Q↑ C↑ (mobilization)
- 'dilution': Q↑ C↓ (dilution dominates)
- 'concentration': Q↓ C↑ (evaporation/point sources)
- 'recession': Q↓ C↓ (system recovery)
- 'chemostatic': Q changes, C stable
- 'source_variation': C changes, Q stable
- 'static': No significant changes
"""
flow_delta = flow_range[1] - flow_range[0]
conc_delta = conc_range[1] - conc_range[0]
# Determine if changes are significant
is_flow_changing = abs(flow_diff) > (threshold_factor * flow_delta) if flow_delta > 1e-10 else False
is_conc_changing = abs(conc_diff) > (threshold_factor * conc_delta) if conc_delta > 1e-10 else False
if not is_flow_changing and not is_conc_changing:
return 'static'
elif is_flow_changing and not is_conc_changing:
return 'chemostatic'
elif not is_flow_changing and is_conc_changing:
return 'source_variation'
else:
# Both changing - determine relationship
if flow_diff > 0 and conc_diff > 0:
return 'flushing'
elif flow_diff < 0 and conc_diff < 0:
return 'recession'
elif flow_diff > 0 and conc_diff < 0:
return 'dilution'
else:
return 'concentration'
def calculate_log_thickness(
hi_values: np.ndarray,
min_thickness: float = 1,
max_thickness: float = 8
) -> np.ndarray:
"""
Calculate log-normalized line thickness for visualization.
Uses log transformation to handle wide range of HI values
while maintaining visual distinction.
Parameters
----------
hi_values : array-like
Absolute hysteresis index values
min_thickness : float
Minimum line thickness
max_thickness : float
Maximum line thickness
Returns
-------
np.ndarray
Thickness values scaled between min and max
"""
hi_abs = np.abs(np.asarray(hi_values))
if len(hi_abs) == 0:
return np.array([])
# Scale factor for good visual separation
scale_factor = 100
log_transformed = np.log10(1 + hi_abs * scale_factor)
# Find range excluding near-zero values
has_variance = np.any(hi_abs > 1e-10)
if has_variance:
log_min = np.min(log_transformed[hi_abs > 1e-10])
log_max = np.max(log_transformed)
else:
return np.full(len(hi_values), min_thickness)
if log_max - log_min < 1e-10:
return np.full(len(hi_values), (min_thickness + max_thickness) / 2)
# Normalize and scale
normalized = np.zeros_like(log_transformed)
mask = hi_abs > 1e-10
normalized[mask] = (log_transformed[mask] - log_min) / (log_max - log_min)
thickness = min_thickness + (max_thickness - min_thickness) * normalized
return np.clip(thickness, min_thickness, max_thickness)
# ==============================================================================
# CORE ANALYSIS ENGINE
# ==============================================================================
def analyze_hysteresis(
data: pd.DataFrame,
sites: List[str],
ccol: str,
qcol: str,
min_dilution_threshold: float = 0.05,
method: str = 'lloyd',
include_harp: bool = False
) -> pd.DataFrame:
"""
Analyze hysteresis patterns for a specific compound across sites.
This function performs compound-specific hysteresis analysis, ensuring
that each compound is analyzed with its own concentration ranges and
characteristics.
Parameters
----------
data : pd.DataFrame
Data with columns: site_id, date, [qcol], [ccol]
sites : list of str
Site IDs to analyze
ccol : str
Concentration column for specific compound
qcol : str
Flow/discharge column
min_dilution_threshold : float
Minimum relative change for significant dilution
method : str
HI calculation method ('lloyd', 'zuecco', 'harp')
include_harp : bool
Whether to include HARP metrics
Returns
-------
pd.DataFrame
Segment-wise analysis results with:
- Behavior classification
- Hysteresis indices
- Potential pollution indicators
- Optional HARP metrics
"""
# Filter and prepare compound-specific data
analysis_data = data[data['site_id'].isin(sites)].copy()
analysis_data = analysis_data.dropna(subset=[qcol, ccol, 'date'])
if len(analysis_data) < 2:
return pd.DataFrame()
analysis_data = analysis_data.sort_values(['site_id', 'date'])
# Calculate compound-specific ranges
flow_range = (analysis_data[qcol].min(), analysis_data[qcol].max())
conc_range = (analysis_data[ccol].min(), analysis_data[ccol].max())
conc_delta = conc_range[1] - conc_range[0]
results = []
for site in sites:
site_data = analysis_data[analysis_data['site_id'] == site].reset_index(drop=True)
if len(site_data) < 2:
continue
# Calculate HARP metrics for entire event if requested
if include_harp:
harp_metrics = calculate_harp_metrics(
site_data[qcol].values,
site_data[ccol].values,
site_data['date'].values if 'date' in site_data else None
)
# Analyze segments
for i in range(len(site_data) - 1):
p1, p2 = site_data.iloc[i], site_data.iloc[i+1]
# Calculate changes
flow_diff = p2[qcol] - p1[qcol]
conc_diff = p2[ccol] - p1[ccol]
# Classify behavior
behavior = classify_cq_behavior(flow_diff, conc_diff, flow_range, conc_range)
# Calculate hysteresis indices for all methods
hi_lloyd = calculate_hysteresis_index(
p1[qcol], p2[qcol], p1[ccol], p2[ccol],
flow_range, conc_range, method='lloyd'
)
hi_zuecco = calculate_hysteresis_index(
p1[qcol], p2[qcol], p1[ccol], p2[ccol],
flow_range, conc_range, method='zuecco'
)
hi_harp = calculate_hysteresis_index(
p1[qcol], p2[qcol], p1[ccol], p2[ccol],
flow_range, conc_range, method='harp'
)
# Get primary HI based on selected method
hi_primary = {'lloyd': hi_lloyd, 'zuecco': hi_zuecco, 'harp': hi_harp}[method]
# Classify HI
if hi_primary > 0.01:
hi_class = 'positive'
elif hi_primary < -0.01:
hi_class = 'negative'
else:
hi_class = 'near-zero'
# Build result
result = {
'site_id': site,
'compound': ccol, # Track which compound
'segment_id': i,
'start_date': p1['date'],
'end_date': p2['date'],
'start_flow': p1[qcol],
'end_flow': p2[qcol],
'start_conc': p1[ccol],
'end_conc': p2[ccol],
'flow_diff': flow_diff,
'conc_diff': conc_diff,
'behavior': behavior,
'hysteresis_index_lloyd': hi_lloyd,
'hysteresis_index_zuecco': hi_zuecco,
'hysteresis_index_harp': hi_harp,
f'hysteresis_index_{method}': hi_primary,
'hi_classification': hi_class,
'hi_magnitude': abs(hi_primary)
}
# Add HydPhase if available
if 'HydPhase' in site_data.columns:
result['start_hyphase'] = p1.get('HydPhase', 'unknown')
result['end_hyphase'] = p2.get('HydPhase', 'unknown')
# Add HARP metrics if requested
if include_harp:
result.update({
'harp_area': harp_metrics['hysteresis_area'],
'harp_residual': harp_metrics['residual'],
'harp_peak_timing': harp_metrics['peak_timing'],
'harp_loop_width': harp_metrics['loop_width']
})
results.append(result)
df = pd.DataFrame(results)
# Identify pollution potential points (compound-specific)
if len(df) > 0:
# Initialize potential flags
for m in ['lloyd', 'zuecco', 'harp']:
df[f'is_potential_{m}'] = False
# Check for pollution potential patterns
for site in df['site_id'].unique():
site_df = df[df['site_id'] == site].reset_index(drop=True)
for i in range(len(site_df) - 1):
prev_seg = site_df.iloc[i]
next_seg = site_df.iloc[i+1]
# Check for concentration → dilution pattern
if (prev_seg['behavior'] == 'concentration' and
next_seg['behavior'] == 'dilution'):
# Check dilution magnitude
dilution_magnitude = abs(next_seg['conc_diff'])
is_significant = (dilution_magnitude / conc_delta > min_dilution_threshold) if conc_delta > 0 else False
# Check each method
for m in ['lloyd', 'zuecco', 'harp']:
hi_col = f'hysteresis_index_{m}'
if (prev_seg[hi_col] > 0.01 and
next_seg[hi_col] < -0.01 and
is_significant):
df.loc[site_df.index[i], f'is_potential_{m}'] = True
# Mark unanimous potential points
df['is_unanimous_potential'] = (
df['is_potential_lloyd'] &
df['is_potential_zuecco'] &
df['is_potential_harp']
)
return df
# ==============================================================================
# VISUALIZATION FUNCTIONS
# ==============================================================================
def create_hysteresis_plot(
data: pd.DataFrame,
sites: List[str],
ccol: str,
qcol: str,
compound: str,
conc_unit: str = 'mg L⁻¹',
flow_unit: str = 'L s⁻¹',
hi_method: str = 'lloyd',
show_timeline: bool = True
) -> go.Figure:
"""
Create hysteresis plot for a single compound.
Combines C-Q hysteresis loops with optional time series visualization.
"""
# Get analysis results
analysis_df = analyze_hysteresis(data, sites, ccol, qcol, method=hi_method)
if len(analysis_df) == 0:
fig = go.Figure()
fig.add_annotation(text="No data available for analysis", x=0.5, y=0.5)
return fig
# Prepare data
plot_data = data[data['site_id'].isin(sites)].copy()
plot_data = plot_data.dropna(subset=[ccol, qcol])
plot_data = plot_data.sort_values(['site_id', 'date'])
# Calculate ranges
flow_min = plot_data[qcol].min() * 0.8
flow_max = plot_data[qcol].max() * 1.2
conc_min = plot_data[ccol].min() * 0.8
conc_max = plot_data[ccol].max() * 1.2
# Color schemes
behavior_colors = {
'flushing': '#d55e00',
'dilution': '#0072b2',
'concentration': '#e69f00',
'recession': '#56b4e9',
'chemostatic': '#009e73',
'source_variation': '#cc79a7',
'static': '#999999'
}
hyphase_colors = {
'low flow': '#FEDF57',
'declining': '#51B848',
'flush': '#1F77B4'
}
# Create subplots
n_sites = len(sites)
if show_timeline:
fig = make_subplots(
rows=2, cols=n_sites,
row_heights=[0.8, 0.2],
subplot_titles=[f'Site {s}' for s in sites] + [''] * n_sites,
vertical_spacing=0.1
)
else:
fig = make_subplots(
rows=1, cols=n_sites,
subplot_titles=[f'Site {s}' for s in sites]
)
# Plot each site
for idx, site in enumerate(sites):
col = idx + 1
site_data = plot_data[plot_data['site_id'] == site]
site_analysis = analysis_df[analysis_df['site_id'] == site]
if len(site_data) < 2:
continue
# Calculate load for dot sizing
site_data['load'] = site_data[ccol] * site_data[qcol] * 86.4
load_95 = site_data['load'].quantile(0.95)
load_5 = site_data['load'].quantile(0.05)
load_range = max(load_95 - load_5, 1e-10)
# Add data points
if 'HydPhase' in site_data.columns:
for phase in site_data['HydPhase'].dropna().unique():
phase_data = site_data[site_data['HydPhase'] == phase]
sizes = 8 + 20 * np.clip((phase_data['load'] - load_5) / load_range, 0, 1)
fig.add_trace(
go.Scatter(
x=phase_data[qcol],
y=phase_data[ccol],
mode='markers',
marker=dict(
size=sizes,
color=hyphase_colors.get(phase, 'gray'),
line=dict(width=1, color='white')
),
name=phase,
showlegend=(idx == 0),
hovertemplate=f'{phase}
Q: %{{x:.2f}}
C: %{{y:.3f}}
Load: %{{customdata:.1f}}',
customdata=phase_data['load']
),
row=1, col=col
)
else:
sizes = 8 + 20 * np.clip((site_data['load'] - load_5) / load_range, 0, 1)
fig.add_trace(
go.Scatter(
x=site_data[qcol],
y=site_data[ccol],
mode='markers',
marker=dict(size=sizes, color='gray'),
name='Data',
showlegend=(idx == 0)
),
row=1, col=col
)
# Add hysteresis lines
hi_values = site_analysis[f'hysteresis_index_{hi_method}'].values
thicknesses = calculate_log_thickness(hi_values)
for i, seg in site_analysis.iterrows():
thickness = thicknesses[i - site_analysis.index[0]]
# Line style based on HI
if seg[f'hysteresis_index_{hi_method}'] > 0.01:
dash = 'solid'
elif seg[f'hysteresis_index_{hi_method}'] < -0.01:
dash = '2px 1px'
else:
dash = '1px 1px'
fig.add_trace(
go.Scatter(
x=[seg['start_flow'], seg['end_flow']],
y=[seg['start_conc'], seg['end_conc']],
mode='lines',
line=dict(
color=behavior_colors[seg['behavior']],
width=thickness,
dash=dash
),
showlegend=False,
hovertemplate=f"{seg['behavior']}
HI: {seg[f'hysteresis_index_{hi_method}']:.3f}"
),
row=1, col=col
)
# Add timeline if requested
if show_timeline:
fig.add_trace(
go.Scatter(
x=[seg['start_date'], seg['end_date']],
y=[seg['start_conc'], seg['end_conc']],
mode='lines',
line=dict(
color=behavior_colors[seg['behavior']],
width=2,
dash=dash
),
showlegend=False
),
row=2, col=col
)
# Highlight pollution potential
potential = site_analysis[site_analysis[f'is_potential_{hi_method}']]
if len(potential) > 0:
fig.add_trace(
go.Scatter(
x=potential['end_flow'],
y=potential['end_conc'],
mode='markers',
marker=dict(
size=18,
color='deeppink',
symbol='pentagon',
line=dict(width=2, color='black')
),
name='Pollution Potential',
showlegend=(idx == 0),
hovertext='Maximum pollution potential point'
),
row=1, col=col
)
# Update layout
fig.update_xaxes(type='log', title_text=f'Flow ({flow_unit})', row=1)
fig.update_yaxes(type='log', title_text=f'{compound} ({conc_unit})', row=1, col=1)
if show_timeline:
fig.update_xaxes(title_text='Date', row=2)
fig.update_yaxes(title_text=f'C ({conc_unit})', row=2, col=1)
# Add behavior legend
for behavior, color in behavior_colors.items():
fig.add_trace(
go.Scatter(
x=[None], y=[None],
mode='lines',
line=dict(color=color, width=4),
name=behavior.title(),
showlegend=True
)
)
fig.update_layout(
title=f'{compound} Hysteresis Analysis ({hi_method.title()} method)',
height=600 if show_timeline else 400,
width=400 * n_sites,
showlegend=True,
legend=dict(orientation='v', x=1.02, y=1),
template = 'none'
)
return fig
def create_multi_compound_hysteresis_plot(
data: pd.DataFrame,
sites: List[str],
ccols: List[str],
compounds: List[str],
conc_units: List[str],
qcol: str,
flow_unit: str = 'L s⁻¹',
hi_method: str = 'lloyd',
cxmin: Optional[List[float]] = None,
cxmax: Optional[List[float]] = None,
qxmin: Optional[float] = None,
qxmax: Optional[float] = None
) -> go.Figure:
"""
Create multi-compound hysteresis comparison plot.
IMPORTANT: Each compound is analyzed with its own concentration ranges
to ensure compound-specific hysteresis patterns are preserved.
Parameters
----------
data : pd.DataFrame
Input data
sites : list
Sites to plot
ccols : list
Concentration columns (one per compound)
compounds : list
Compound names for labeling
conc_units : list
Units for each compound
qcol : str
Flow column name
flow_unit : str
Flow unit
hi_method : str
Hysteresis calculation method
cxmin, cxmax : list of float, optional
Min/max concentration for each compound
qxmin, qxmax : float, optional
Min/max flow range (shared across compounds)
Returns
-------
go.Figure
Multi-panel plot with compounds as rows, sites as columns
"""
n_compounds = len(ccols)
n_sites = len(sites)
# Create subplot structure
subplot_titles = []
for comp in compounds:
for site in sites:
subplot_titles.append(f'{comp} - {site}')
for _ in sites:
subplot_titles.append('') # Empty for time series
fig = make_subplots(
rows=n_compounds * 2,
cols=n_sites,
row_heights=[0.85, 0.15] * n_compounds,
subplot_titles=subplot_titles,
vertical_spacing=0.03,
horizontal_spacing=0.05
)
# Color schemes
behavior_colors = {
'flushing': '#d55e00',
'dilution': '#0072b2',
'concentration': '#e69f00',
'recession': '#56b4e9',
'chemostatic': '#009e73',
'source_variation': '#cc79a7',
'static': '#999999'
}
hyphase_colors = {
'low flow': '#FEDF57',
'declining': '#51B848',
'flush': '#1F77B4'
}
# Determine flow range (shared across compounds)
if qxmin is None or qxmax is None:
all_flows = data[qcol].dropna()
if len(all_flows) > 0:
qxmin = all_flows.min() * 0.8 if qxmin is None else qxmin
qxmax = all_flows.max() * 1.2 if qxmax is None else qxmax
else:
qxmin, qxmax = 0.1, 100
# Process each compound separately
for comp_idx, (ccol, compound, conc_unit) in enumerate(zip(ccols, compounds, conc_units)):
# COMPOUND-SPECIFIC ANALYSIS
compound_data = data.dropna(subset=[ccol, qcol])
# Determine concentration range for THIS compound
if cxmin is None or len(cxmin) <= comp_idx:
comp_cmin = compound_data[ccol].min() * 0.8
else:
comp_cmin = cxmin[comp_idx]
if cxmax is None or len(cxmax) <= comp_idx:
comp_cmax = compound_data[ccol].max() * 1.2
else:
comp_cmax = cxmax[comp_idx]
# Analyze hysteresis for THIS compound
analysis_df = analyze_hysteresis(
compound_data, sites, ccol, qcol, method=hi_method
)
# Process each site
for site_idx, site in enumerate(sites):
h_row = comp_idx * 2 + 1
t_row = comp_idx * 2 + 2
col = site_idx + 1
site_data = compound_data[compound_data['site_id'] == site].copy()
site_analysis = analysis_df[analysis_df['site_id'] == site]
if len(site_data) < 2:
continue
# Calculate load for THIS compound
site_data['load'] = site_data[ccol] * site_data[qcol] * 86.4
load_95 = site_data['load'].quantile(0.95)
load_5 = site_data['load'].quantile(0.05)
load_range = max(load_95 - load_5, 1e-10)
# Plot data points
if 'HydPhase' in site_data.columns:
for phase in site_data['HydPhase'].dropna().unique():
phase_data = site_data[site_data['HydPhase'] == phase]
sizes = 8 + 20 * np.clip((phase_data['load'] - load_5) / load_range, 0, 1)
fig.add_trace(
go.Scatter(
x=phase_data[qcol],
y=phase_data[ccol],
mode='markers',
marker=dict(
size=sizes,
color=hyphase_colors.get(phase, 'gray'),
line=dict(width=1, color='white')
),
name=phase,
showlegend=(comp_idx == 0 and site_idx == 0),
hovertemplate=f'{compound} - {phase}
Q: %{{x:.2f}}
C: %{{y:.3f}}'
),
row=h_row, col=col
)
# Plot hysteresis lines with compound-specific thickness
hi_values = site_analysis[f'hysteresis_index_{hi_method}'].values
thicknesses = calculate_log_thickness(hi_values)
for i, seg in site_analysis.iterrows():
thickness = thicknesses[i - site_analysis.index[0]]
# Line style based on HI
if seg[f'hysteresis_index_{hi_method}'] > 0.01:
dash = 'solid'
elif seg[f'hysteresis_index_{hi_method}'] < -0.01:
dash = '2px 1px'
else:
dash = '1px 1px'
fig.add_trace(
go.Scatter(
x=[seg['start_flow'], seg['end_flow']],
y=[seg['start_conc'], seg['end_conc']],
mode='lines',
line=dict(
color=behavior_colors[seg['behavior']],
width=thickness,
dash=dash
),
showlegend=False
),
row=h_row, col=col
)
# Time series
fig.add_trace(
go.Scatter(
x=[seg['start_date'], seg['end_date']],
y=[seg['start_conc'], seg['end_conc']],
mode='lines',
line=dict(
color=behavior_colors[seg['behavior']],
width=2,
dash=dash
),
showlegend=False
),
row=t_row, col=col
)
# Highlight pollution potential
potential = site_analysis[site_analysis[f'is_potential_{hi_method}']]
unanimous = site_analysis[site_analysis['is_unanimous_potential']]
if len(potential) > 0:
# Non-unanimous potential
non_unanimous = potential[~potential['is_unanimous_potential']]
if len(non_unanimous) > 0:
fig.add_trace(
go.Scatter(
x=non_unanimous['end_flow'],
y=non_unanimous['end_conc'],
mode='markers',
marker=dict(
size=18,
color='deeppink',
symbol='pentagon'
),
name='Potential',
showlegend=(comp_idx == 0 and site_idx == 0)
),
row=h_row, col=col
)
# Unanimous potential
if len(unanimous) > 0:
fig.add_trace(
go.Scatter(
x=unanimous['end_flow'],
y=unanimous['end_conc'],
mode='markers',
marker=dict(
size=18,
color='deeppink',
symbol='pentagon',
line=dict(width=3, color='black')
),
name='Unanimous',
showlegend=(comp_idx == 0 and site_idx == 0)
),
row=h_row, col=col
)
# Update axes with compound-specific ranges
fig.update_xaxes(
type='log',
range=[np.log10(qxmin), np.log10(qxmax)],
row=h_row, col=col
)
fig.update_yaxes(
type='log',
range=[np.log10(comp_cmin), np.log10(comp_cmax)],
title_text=f'{compound} ({conc_unit})' if site_idx == 0 else '',
row=h_row, col=col
)
# Time series axes
fig.update_yaxes(
range=[comp_cmin, comp_cmax],
row=t_row, col=col
)
# Add x-axis labels only on bottom
for col in range(1, n_sites + 1):
fig.update_xaxes(
title_text=f'Flow ({flow_unit})',
row=n_compounds * 2 - 1, col=col
)
# Add legend for behaviors
for behavior, color in behavior_colors.items():
fig.add_trace(
go.Scatter(
x=[None], y=[None],
mode='lines',
line=dict(color=color, width=4),
name=behavior.title()
)
)
# Add HI legend
fig.add_trace(
go.Scatter(
x=[None], y=[None],
mode='lines',
line=dict(color='gray', width=3, dash='solid'),
name='Positive HI'
)
)
fig.add_trace(
go.Scatter(
x=[None], y=[None],
mode='lines',
line=dict(color='gray', width=3, dash='2px 1px'),
name='Negative HI'
)
)
fig.update_layout(
title=f'Multi-Compound Hysteresis Analysis ({hi_method.title()} method)
' +
'Each compound analyzed with its own concentration range',
height=450 * n_compounds,
width=350 * n_sites + 200,
showlegend=True,
legend=dict(x=1.01, y=1),
template = 'none'
)
return fig
def create_hysteresis_summary_stats(
data: pd.DataFrame,
sites: List[str],
ccol: str,
qcol: str,
hi_method: str = 'lloyd'
) -> pd.DataFrame:
"""
Generate summary statistics for hysteresis behavior.
Returns
-------
pd.DataFrame
Summary with behavior frequencies, durations, and HI statistics
"""
analysis_df = analyze_hysteresis(data, sites, ccol, qcol, method=hi_method)
if len(analysis_df) == 0:
return pd.DataFrame()
# Calculate statistics by site and behavior
summary = []
for site in sites:
site_df = analysis_df[analysis_df['site_id'] == site]
if len(site_df) == 0:
continue
# Overall statistics
total_segments = len(site_df)
date_range = site_df['end_date'].max() - site_df['start_date'].min()
# Behavior statistics
for behavior in site_df['behavior'].unique():
behavior_df = site_df[site_df['behavior'] == behavior]
# Calculate durations
durations = (behavior_df['end_date'] - behavior_df['start_date']).dt.total_seconds() / 86400
summary.append({
'site_id': site,
'compound': ccol,
'behavior': behavior,
'count': len(behavior_df),
'percentage': len(behavior_df) / total_segments * 100,
'mean_duration_days': durations.mean(),
'total_duration_days': durations.sum(),
f'mean_hi_{hi_method}': behavior_df[f'hysteresis_index_{hi_method}'].mean(),
f'std_hi_{hi_method}': behavior_df[f'hysteresis_index_{hi_method}'].std(),
'monitoring_days': date_range.days,
'total_segments': total_segments
})
return pd.DataFrame(summary)
def create_hysteresis_timeline(
data: pd.DataFrame,
sites: List[str],
ccol: str,
qcol: str,
compound: str,
hi_method: str = 'lloyd'
) -> go.Figure:
"""
Create timeline visualization of hysteresis behaviors.
"""
analysis_df = analyze_hysteresis(data, sites, ccol, qcol, method=hi_method)
if len(analysis_df) == 0:
fig = go.Figure()
fig.add_annotation(text="No data available", x=0.5, y=0.5)
return fig
# Color scheme
behavior_colors = {
'flushing': '#d55e00',
'dilution': '#0072b2',
'concentration': '#e69f00',
'recession': '#56b4e9',
'chemostatic': '#009e73',
'source_variation': '#cc79a7',
'static': '#999999'
}
fig = go.Figure()
# Plot timeline for each site
for site_idx, site in enumerate(sites):
site_df = analysis_df[analysis_df['site_id'] == site]
for _, seg in site_df.iterrows():
# Line style based on HI
if seg[f'hysteresis_index_{hi_method}'] > 0.01:
dash = 'solid'
elif seg[f'hysteresis_index_{hi_method}'] < -0.01:
dash = '2px 1px'
else:
dash = '1px 1px'
fig.add_trace(
go.Scatter(
x=[seg['start_date'], seg['end_date']],
y=[site, site],
mode='lines',
line=dict(
color=behavior_colors[seg['behavior']],
width=20,
dash=dash
),
showlegend=False,
hovertemplate=(
f"{seg['behavior']}
"
f"HI: {seg[f'hysteresis_index_{hi_method}']:.3f}
"
f"Duration: {(seg['end_date'] - seg['start_date']).days} days"
)
)
)
# Add legend
for behavior, color in behavior_colors.items():
fig.add_trace(
go.Scatter(
x=[None], y=[None],
mode='lines',
line=dict(color=color, width=10),
name=behavior.title()
)
)
fig.update_layout(
title=f'{compound} Hysteresis Behavior Timeline',
xaxis_title='Date',
yaxis_title='Site',
yaxis=dict(categoryorder='array', categoryarray=sites[::-1]),
height=100 + 60 * len(sites),
showlegend=True,
legend=dict(orientation='h', y=1.1),
template = 'none'
)
return fig