venkatnsn

"""Regional Sensitivity Analysis + GLUE uncertainty in one script."""
import numpy as np
from scipy.stats import qmc

# 5,000-sample Latin Hypercube over the prior parameter bounds
N = 5000
sampler = qmc.LatinHypercube(d=5, seed=2026)
u = sampler.random(N)
bounds = np.array(BOUNDS)
samples = bounds[:, 0] + u * (bounds[:, 1] - bounds[:, 0])

venkatnsn

"""Differential Evolution against KGE, with split-sample validation."""
import numpy as np
from scipy.optimize import differential_evolution

WARMUP = 365  # days excluded from the metric

def kge(sim, obs):
    """Kling-Gupta Efficiency (Gupta et al. 2009)."""
    mask = np.isfinite(sim) & np.isfinite(obs)
    sim, obs = sim[mask], obs[mask]

venkatnsn

"""5-parameter two-bucket rainfall-runoff model, Numba-JIT'd.

Soil bucket -> ET (atm) + saturation excess -> Quickflow + Percolation
Groundwater bucket -> Baseflow

Total streamflow Q = Quickflow + Baseflow.
Five free parameters: S_max, f_quick, k_perc, tau_b, gamma.
"""
import numpy as np
from numba import njit, prange

venkatnsn

"""Real-time daily data for any U.S. catchment, free and no-auth."""
import requests, pandas as pd
from dataretrieval import nwis

CFS_TO_M3S = 0.0283168
SECONDS_PER_DAY = 86400.0

def fetch_streamflow(site_no, start, end, area_km2):
    """USGS NWIS daily mean discharge -> mm/day (depth equivalent)."""
    df, _ = nwis.get_dv(sites=site_no, parameterCd="00060", start=start, end=end)

venkatnsn

behave_idx = np.where(behave)[0]
likelihood = nse_results[behave_idx] - nse_results[behave_idx].min()
weights    = likelihood / likelihood.sum()

# Future scenario: 20% drier rainfall
P_future = P * 0.8
Q_future = np.zeros((len(behave_idx), days))
for k, i in enumerate(behave_idx):
    p = dict(S_max=S_max_vals[i], f=f_vals[i], tau=tau_vals[i])
    Q_future[k], _ = run_bucket(P_future, PET, p)

venkatnsn

from scipy.stats import qmc
import matplotlib.pyplot as plt

N = 5000
sampler = qmc.LatinHypercube(d=3, seed=1)
u = sampler.random(N)
S_max_vals = 20.0 + (400.0 - 20.0) * u[:, 0]
f_vals     = 0.0  + (1.0   - 0.0)  * u[:, 1]
tau_vals   = 1.0  + (100.0 - 1.0)  * u[:, 2]

venkatnsn

half = days // 2
P_cal,    P_val    = P[:half],    P[half:]
PET_cal,  PET_val  = PET[:half],  PET[half:]
Qobs_cal, Qobs_val = Q_obs[:half], Q_obs[half:]

def neg_nse_cal(x):
    p = dict(S_max=x[0], f=x[1], tau=x[2])
    Q_sim, _ = run_bucket(P_cal, PET_cal, p)
    return -nse(Q_sim, Qobs_cal)

venkatnsn

from scipy.optimize import differential_evolution

bounds = [(20, 400),     # S_max
          (0.0, 1.0),    # f
          (1.0, 100.0)]  # tau

result = differential_evolution(neg_nse, bounds, tol=1e-7, seed=42)
p_opt  = dict(S_max=result.x[0], f=result.x[1], tau=result.x[2])
print(f"DE NSE = {-result.fun:.3f}")
print(f"  Params: {p_opt}")

venkatnsn

from scipy.optimize import minimize

def neg_nse(x):
    p = dict(S_max=x[0], f=x[1], tau=x[2])
    Q_sim, _ = run_bucket(P, PET, p)
    return -nse(Q_sim, Q_obs)

x0     = [100, 0.5, 20]
result = minimize(neg_nse, x0, method='Nelder-Mead')
print(f"Nelder-Mead NSE = {-result.fun:.3f}")

venkatnsn

guess = dict(S_max=100, f=0.5, tau=20)
Q_guess, _ = run_bucket(P, PET, guess)
print(f"Manual guess NSE = {nse(Q_guess, Q_obs):.3f}")
# Manual guess NSE = 0.71