API Reference

def lh(Q, beta=0.925, return_exceed=False):
    """Lyne-Hollick digital filter, 2-pass (Lyne & Hollick, 1979).

    Forward pass:  b[t] = beta*b[t-1] + (1-beta)/2 * (Q[t] + Q[t-1])
    Backward pass: applied to the forward result for smoothing.

    The default beta=0.925 follows the recommendation of Nathan & McMahon
    (1990). Using lh_multi(Q, beta=0.925, num_pass=3) reproduces the
    Nathan-McMahon recommended 3-pass protocol exactly.

    Lyne, V. and Hollick, M. (1979) Stochastic time-variable rainfall-runoff
    modelling. Inst. Eng. Aust. Natl. Conf., 89-93.

    Nathan, R.J. and McMahon, T.A. (1990) Evaluation of automated techniques
    for base flow and recession analyses. Water Resour. Res. 26(7), 1465-1473.

    Args:
        Q (numpy.ndarray): Streamflow.
        beta (float): Filter parameter. Default 0.925.
        return_exceed (bool): If True, append exceed count.

    Returns:
        numpy.ndarray: Baseflow.
    """
    if return_exceed:
        b = np.zeros(Q.shape[0] + 1)
    else:
        b = np.zeros(Q.shape[0])

    # first pass (forward)
    b[0] = Q[0]
    for i in range(Q.shape[0] - 1):
        b[i + 1] = beta * b[i] + (1 - beta) / 2 * (Q[i] + Q[i + 1])
        if b[i + 1] > Q[i + 1]:
            b[i + 1] = Q[i + 1]
            if return_exceed:
                b[-1] += 1

    # second pass (backward)
    b1 = np.copy(b)
    for i in range(Q.shape[0] - 2, -1, -1):
        b[i] = beta * b[i + 1] + (1 - beta) / 2 * (b1[i + 1] + b1[i])
        if b[i] > b1[i]:
            b[i] = b1[i]
            if return_exceed:
                b[-1] += 1
    return b

def lh_multi(Q, beta=0.925, num_pass=2, return_exceed=False):
    """Lyne-Hollick digital filter, n-pass (Spongberg, 2000).

    Applies the LH filter with alternating forward/backward passes.

    Spongberg, M.E. (2000) Spectral analysis of base flow separation with
    digital filters. Water Resour. Res. 36(3), 745-752.

    Args:
        Q (numpy.ndarray): Streamflow.
        beta (float): Filter parameter. Default 0.925.
        num_pass (int): Number of passes. Default 2.
        return_exceed (bool): If True, append exceed count.

    Returns:
        numpy.ndarray: Baseflow.
    """
    if return_exceed:
        b = np.zeros(Q.shape[0] + 1)
    else:
        b = np.zeros(Q.shape[0])

    b[0] = Q[0]

    for n in range(num_pass):
        if n != 0:
            b = np.flip(b, axis=0)
            Q = b.copy()

        for i in range(Q.shape[0] - 1):
            b[i + 1] = beta * b[i] + (1 - beta) / 2 * (Q[i] + Q[i + 1])
            if b[i + 1] > Q[i + 1]:
                b[i + 1] = Q[i + 1]
                if return_exceed:
                    b[-1] += 1

    if num_pass % 2 == 0:
        b = np.flip(b, axis=0)

    return b

def eckhardt(Q, a, BFImax, initial_method='Q0', return_exceed=False):
    """Eckhardt two-parameter filter (Eckhardt, 2005).

    b[t] = ((1-BFImax)*a*b[t-1] + (1-a)*BFImax*Q[t]) / (1 - a*BFImax)

    Eckhardt, K. (2005) How to construct recursive digital filters for
    baseflow separation. Hydrol. Process. 19, 507-515.

    Args:
        Q (numpy.ndarray): Streamflow.
        a (float): Recession coefficient.
        BFImax (float): Maximum baseflow index.
        initial_method (str or float): Method to set b[0]. Default 'Q0'.
        return_exceed (bool): If True, append exceed count.

    Returns:
        numpy.ndarray: Baseflow.
    """
    denom = 1 - a * BFImax
    return _recursive_digital_filter(
        Q, alpha=(1 - BFImax) * a / denom,
        beta=(1 - a) * BFImax / denom, gamma=0.0,
        initial_method=initial_method, return_exceed=return_exceed)

def chapman_maxwell(Q, a, initial_method='Q0', return_exceed=False):
    """Chapman-Maxwell filter (Chapman & Maxwell, 1996).

    b[t] = a/(2-a) * b[t-1] + (1-a)/(2-a) * Q[t]

    Equivalent to Eckhardt with BFImax = 0.5.

    Chapman, T.G., Maxwell, A.I. (1996) - Baseflow separation - comparison of
    numerical methods with tracer experiments. Hydrol. Water Resour. Symp.,
    Inst. Eng. Australia, Hobart, 539-545.

    Args:
        Q (numpy.ndarray): Streamflow.
        a (float): Recession coefficient.
        initial_method (str or float): Method to set b[0]. Default 'Q0'.
        return_exceed (bool): If True, append exceed count.

    Returns:
        numpy.ndarray: Baseflow.
    """
    return _recursive_digital_filter(
        Q, alpha=a / (2 - a), beta=(1 - a) / (2 - a), gamma=0.0,
        initial_method=initial_method, return_exceed=return_exceed)

def chapman(Q, a=0.925, initial_method='Q0', return_exceed=False):
    """Chapman filter (Chapman, 1991).

    b[t] = (3a-1)/(3-a) * b[t-1] + (1-a)/(3-a) * (Q[t] + Q[t-1])

    Chapman, T.G. (1991) Comment on 'Evaluation of automated techniques for
    base flow and recession analyses' by R.J. Nathan and T.A. McMahon.
    Water Resour. Res. 27(7), 1783-1784.

    Args:
        Q (numpy.ndarray): Streamflow.
        a (float): Recession coefficient. Default 0.925.
        initial_method (str or float): Method to set b[0]. Default 'Q0'.
        return_exceed (bool): If True, append exceed count.

    Returns:
        numpy.ndarray: Baseflow.
    """
    return _recursive_digital_filter(
        Q, alpha=(3 * a - 1) / (3 - a), beta=(1 - a) / (3 - a), gamma=1.0,
        initial_method=initial_method, return_exceed=return_exceed)

def boughton(Q, a, C, initial_method='Q0', return_exceed=False):
    """Boughton two-parameter filter (Boughton, 1993).

    b[t] = a/(1+C) * b[t-1] + C/(1+C) * Q[t]

    Boughton W.C. (1993) - A hydrograph-based model for estimating water yield
    of ungauged catchments. Inst. Eng. Aust. Natl. Conf. Publ. 93/14, 317-324.

    Args:
        Q (numpy.ndarray): Streamflow.
        a (float): Recession coefficient.
        C (float): Calibrated parameter (see param_calibrate).
        initial_method (str or float): Method to set b[0]. Default 'Q0'.
        return_exceed (bool): If True, append exceed count.

    Returns:
        numpy.ndarray: Baseflow.
    """
    return _recursive_digital_filter(
        Q, alpha=a / (1 + C), beta=C / (1 + C), gamma=0.0,
        initial_method=initial_method, return_exceed=return_exceed)

def ewma(Q, e, initial_method='Q0', return_exceed=False):
    """Exponential weighted moving average filter (Tularam & Ilahee, 2008).

    b[t] = (1-e) * b[t-1] + e * Q[t]

    Tularam, G.A. and Ilahee, M. (2008) Exponential smoothing method of base
    flow separation and its impact on continuous loss estimates. Am. J. Environ.
    Sci. 4(2), 136-144.

    Args:
        Q (numpy.ndarray): Streamflow.
        e (float): Smoothing parameter.
        initial_method (str or float): Method to set b[0]. Default 'Q0'.
        return_exceed (bool): If True, append exceed count.

    Returns:
        numpy.ndarray: Baseflow.
    """
    return _recursive_digital_filter(
        Q, alpha=(1 - e), beta=e, gamma=0.0,
        initial_method=initial_method, return_exceed=return_exceed)

def furey(Q, a, A, initial_method='Q0', return_exceed=False):
    """Furey-Gupta physically-based filter (Furey & Gupta, 2001).

    b[t] = (a - A*(1-a)) * b[t-1] + A*(1-a) * Q[t-1]

    Note: this filter uses Q[t-1] (not Q[t]) in the second term, so gamma
    acts on the *previous* timestep's streamflow via the general form with
    beta applied to Q[t] + gamma*Q[t-1]. Here we set beta on Q[t-1] by
    mapping: alpha = a - A*(1-a), beta = A*(1-a), gamma = 0, and shifting
    the Q index.

    Furey, P.R. and Gupta, V.K. (2001) A physically based filter for
    separating base flow from streamflow time series. Water Resour. Res.
    37(11), 2709-2722.

    Args:
        Q (numpy.ndarray): Streamflow.
        a (float): Recession coefficient.
        A (float): Calibrated parameter.
        initial_method (str or float): Method to set b[0]. Default 'Q0'.
        return_exceed (bool): If True, append exceed count.

    Returns:
        numpy.ndarray: Baseflow.
    """
    # Furey uses Q[i] (previous timestep) rather than Q[i+1], so it doesn't
    # fit the standard general form directly. Implement the loop explicitly.
    n = Q.shape[0]
    if return_exceed:
        b = np.zeros(n + 1)
    else:
        b = np.zeros(n)

    b[0] = _init_baseflow(Q, initial_method)

    for i in range(n - 1):
        b[i + 1] = (a - A * (1 - a)) * b[i] + A * (1 - a) * Q[i]
        if b[i + 1] > Q[i + 1]:
            b[i + 1] = Q[i + 1]
            if return_exceed:
                b[-1] += 1
    return b

def willems(Q, a, w, initial_method='Q0', return_exceed=False):
    """Willems digital filter (Willems, 2009).

    v = (1-w)*(1-a) / (2*w)
    b[t] = (a-v)/(1+v) * b[t-1] + v/(1+v) * (Q[t] + Q[t-1])

    Willems, P. (2009) A time series tool to support the multi-criteria
    performance evaluation of rainfall-runoff models. Environ. Model. Softw.
    24(3), 311-321.

    Args:
        Q (numpy.ndarray): Streamflow.
        a (float): Recession coefficient.
        w (float): Average proportion of quick flow in streamflow.
        initial_method (str or float): Method to set b[0]. Default 'Q0'.
        return_exceed (bool): If True, append exceed count.

    Returns:
        numpy.ndarray: Baseflow.
    """
    v = (1 - w) * (1 - a) / (2 * w)
    return _recursive_digital_filter(
        Q, alpha=(a - v) / (1 + v), beta=v / (1 + v), gamma=1.0,
        initial_method=initial_method, return_exceed=return_exceed)

def ihacres(Q, a, C, alpha_s, initial_method='Q0', return_exceed=False):
    """Jakeman-Hornberger / IHACRES baseflow filter (Jakeman & Hornberger, 1993).

    b[t] = a/(1+C) * b[t-1] + C/(1+C) * (Q[t] + alpha_s * Q[t-1])

    Three-parameter extension of the Boughton filter. When alpha_s = 0,
    reduces exactly to boughton(). The alpha_s term adds dependence on
    the previous timestep's total flow, accounting for delayed baseflow
    response.

    In the generalized filter taxonomy: alpha = a/(1+C), beta = C/(1+C),
    gamma = alpha_s. This filter occupies the continuum between the
    gamma=0 and gamma=1 families.

    Jakeman, A.J. and Hornberger, G.M. (1993) How much complexity is
    warranted in a rainfall-runoff model? Water Resour. Res. 29(8),
    2637-2649.

    Args:
        Q (numpy.ndarray): Streamflow.
        a (float): Recession coefficient for slow flow (0.9-0.999).
        C (float): Partitioning parameter (0.1-1.0).
        alpha_s (float): Previous-timestep dependence (-0.99 to 0).
            Set to 0 to recover the Boughton filter.
        initial_method (str or float): Method to set b[0]. Default 'Q0'.
        return_exceed (bool): If True, append exceed count.

    Returns:
        numpy.ndarray: Baseflow.
    """
    return _recursive_digital_filter(
        Q, alpha=a / (1 + C), beta=C / (1 + C), gamma=alpha_s,
        initial_method=initial_method, return_exceed=return_exceed)

def what(Q, BFImax, a):
    """WHAT method (Lim et al., 2005). Alias for eckhardt().

    Mathematically identical to the Eckhardt filter.

    Lim, K.J., et al. (2005) Automated Web GIS based hydrograph analysis tool,
    WHAT. JAWRA 41(6), 1407-1416.

    Args:
        Q (numpy.ndarray): Streamflow.
        BFImax (float): Maximum baseflow index.
        a (float): Recession coefficient.

    Returns:
        numpy.ndarray: Baseflow.
    """
    return eckhardt(Q, a, BFImax)

def fixed(Q, area=None):
    """Fixed interval graphical method from HYSEP (Sloto & Crouse, 1996).

    Sloto, R.A. & Crouse, M.Y. (1996) HYSEP: A computer program for
    streamflow hydrograph separation and analysis. USGS WRI 96-4040.

    Args:
        Q (numpy.ndarray): Streamflow.
        area (float): Basin area in km^2.

    Returns:
        numpy.ndarray: Baseflow.
    """
    inN = _hysep_interval(area)
    return _fixed_interpolation(Q, inN)

def slide(Q, area):
    """Sliding interval graphical method from HYSEP (Sloto & Crouse, 1996).

    Sloto, R.A. & Crouse, M.Y. (1996) HYSEP: A computer program for
    streamflow hydrograph separation and analysis. USGS WRI 96-4040.

    Args:
        Q (numpy.ndarray): Streamflow.
        area (float): Basin area in km^2.

    Returns:
        numpy.ndarray: Baseflow.
    """
    inN = _hysep_interval(area)
    return _slide_interpolation(Q, inN)

def local(Q, b_LH, area=None, return_exceed=False):
    """Local minimum graphical method from HYSEP (Sloto & Crouse, 1996).

    Sloto, R.A. & Crouse, M.Y. (1996) HYSEP: A computer program for
    streamflow hydrograph separation and analysis. USGS WRI 96-4040.

    Args:
        Q (numpy.ndarray): Streamflow.
        b_LH (numpy.ndarray): Lyne-Hollick baseflow (used for edge filling).
        area (float): Basin area in km^2.
        return_exceed (bool): If True, append exceed count.

    Returns:
        numpy.ndarray: Baseflow.
    """
    idx_turn = _local_turn(Q, _hysep_interval(area))
    if idx_turn.shape[0] < 3:
        raise IndexError('Less than 3 turning points found')
    b = _linear_interpolation(Q, idx_turn, return_exceed=return_exceed)
    b[:idx_turn[0]] = b_LH[:idx_turn[0]]
    b[idx_turn[-1] + 1:] = b_LH[idx_turn[-1] + 1:]
    return b

def ukih(Q, b_LH, return_exceed=False):
    """UK Institute of Hydrology smoothed minima method (UKIH, 1980).

    Aksoy, H., Kurt, I. and Eris, E. (2009) Filtered smoothed minima baseflow
    separation method. J. Hydrol. 372(1), 94-101.

    Args:
        Q (numpy.ndarray): Streamflow.
        b_LH (numpy.ndarray): Lyne-Hollick baseflow (used for edge filling).
        return_exceed (bool): If True, append exceed count.

    Returns:
        numpy.ndarray: Baseflow.
    """
    N = 5
    block_end = Q.shape[0] // N * N
    idx_min = np.argmin(Q[:block_end].reshape(-1, N), axis=1)
    idx_min = idx_min + np.arange(0, block_end, N)
    idx_turn = _ukih_turn(Q, idx_min)
    if idx_turn.shape[0] < 3:
        raise IndexError('Less than 3 turning points found')
    b = _linear_interpolation(Q, idx_turn, return_exceed=return_exceed)
    b[:idx_turn[0]] = b_LH[:idx_turn[0]]
    b[idx_turn[-1] + 1:] = b_LH[idx_turn[-1] + 1:]
    return b

def part(Q, area, log_cycle_threshold=0.1):
    """PART baseflow separation (Rutledge, 1998).

    Recession-based graphical method that identifies days where streamflow
    equals baseflow using an antecedent recession requirement, then
    interpolates in log-space between qualifying days.

    The procedure is run three times with different recession durations
    (floor(N), ceil(N), ceil(N)+1) and combined via linear interpolation
    evaluated at the exact N = A^0.2.

    Rutledge, A.T. (1998) Computer programs for describing the recession of
    ground-water discharge and for estimating mean ground-water recharge and
    discharge from streamflow records — Update. USGS WRI 98-4148.

    Args:
        Q (numpy.ndarray): Daily streamflow.
        area (float): Drainage area in km^2.
        log_cycle_threshold (float): Maximum daily log10 decline for a
            qualifying day. Default 0.1.

    Returns:
        numpy.ndarray: Daily baseflow array.
    """
    # Compute real-valued N (days of antecedent recession)
    area_sqmi = 0.3861022 * area
    N_real = area_sqmi ** 0.2

    N_floor = max(int(np.floor(N_real)), 1)
    N_ceil = max(int(np.ceil(N_real)), 2)
    if N_ceil == N_floor:
        N_ceil = N_floor + 1
    N_values = [N_floor, N_ceil, N_ceil + 1]

    results = []
    for N in N_values:
        b = _part_single(Q, N, log_cycle_threshold)
        results.append(b)

    # Linear interpolation between floor and ceil results
    if N_ceil > N_floor:
        frac = (N_real - N_floor) / (N_ceil - N_floor)
    else:
        frac = 0.0
    b_out = results[0] * (1 - frac) + results[1] * frac

    return b_out

def strict_baseflow(Q, ice=None, quantile=0.9):
    """Identify strict baseflow periods in a streamflow time series.

    Applies heuristic rules based on the derivative of the hydrograph to
    identify timesteps that are unambiguously baseflow-dominated. Used by
    recession_coefficient() for recession constant estimation.

    Args:
        Q (numpy.ndarray): Streamflow.
        ice (numpy.ndarray, optional): Boolean mask for ice-affected periods.
        quantile (float): Quantile threshold for major events. Default 0.9.

    Returns:
        numpy.ndarray: Boolean mask (True = strict baseflow).
    """
    dQ = (Q[2:] - Q[:-2]) / 2

    # 1. flow data associated with positive and zero values of dy / dt
    wet1 = np.concatenate([[True], dQ >= 0, [True]])

    # 2. previous 2 points before points with dy/dt>=0, plus the next 3 points
    idx_first = np.where(wet1[1:].astype(int) - wet1[:-1].astype(int) == 1)[0] + 1
    idx_last = np.where(wet1[1:].astype(int) - wet1[:-1].astype(int) == -1)[0]
    idx_before = np.repeat([idx_first], 2) - np.tile(range(1, 3), idx_first.shape)
    idx_next = np.repeat([idx_last], 3) + np.tile(range(1, 4), idx_last.shape)
    idx_remove = np.concatenate([idx_before, idx_next])
    wet2 = np.full(Q.shape, False)
    wet2[idx_remove.clip(min=0, max=Q.shape[0] - 1)] = True

    # 3. five data points after major events (quantile)
    growing = np.concatenate([[True], (Q[1:] - Q[:-1]) >= 0, [True]])
    idx_major = np.where((Q >= np.quantile(Q, quantile)) & growing[:-1] & ~growing[1:])[0]
    idx_after = np.repeat([idx_major], 5) + np.tile(range(1, 6), idx_major.shape)
    wet3 = np.full(Q.shape, False)
    wet3[idx_after.clip(min=0, max=Q.shape[0] - 1)] = True

    # 4. flow data followed by a data point with a larger value of -dy / dt
    wet4 = np.concatenate([[True], dQ[1:] - dQ[:-1] < 0, [True, True]])

    # dry points = strict baseflow
    dry = ~(wet1 + wet2 + wet3 + wet4)

    if ice is not None:
        dry[ice] = False

    return dry

def bn77(Q, L_min, snow_freeze_period, observational_precision, quantile=0.9):
    """Brutsaert-Nieber drought flow identification (Cheng et al., 2016).

    Identifies drought (pure baseflow) points in a discharge time series
    using recession slope analysis and multiple elimination criteria.

    Cheng, L., Zhang, L. and Brutsaert, W. (2016) Automated selection of pure
    base flows from regular daily streamflow data: objective algorithm.
    J. Hydrol. Eng. 21(11), 06016008.

    Args:
        Q (numpy.ndarray): Discharge time series.
        L_min (int): Minimum recession episode length.
        snow_freeze_period (tuple): (start_index, end_index) of snow/freeze period.
        observational_precision (float): Minimum Q threshold.
        quantile (float): Quantile for major event identification. Default 0.9.

    Returns:
        numpy.ndarray: Indices of drought flow points.
    """
    S = _estimate_recession_slope(Q)
    recession_episodes = _identify_recession_episodes(S, L_min)
    drought_flow_points = _eliminate_points(
        recession_episodes, L_min, snow_freeze_period,
        observational_precision, Q, S, quantile)
    return drought_flow_points

def recession_coefficient(Q, strict):
    """
    Calculates the recession coefficient `K` from the given discharge `Q` and a boolean mask `strict` indicating which values to use.

    The recession coefficient `K` is calculated as follows:
    1. Extract the middle values of `Q` (`cQ`) and the centered finite difference of `Q` (`dQ`) using the `strict` mask.
    2. Sort `dQ / cQ` in descending order and take the value at the 5th percentile.
    3. Calculate `K` as the negative ratio of `cQ` to `dQ` at the selected index.
    4. Return the exponential of `-1 / K` as the final recession coefficient.

    Args:
        Q (numpy.ndarray): Array of discharge values.
        strict (numpy.ndarray): Boolean mask indicating which values of `Q` to use.

    Returns:
        float: The calculated recession coefficient.
    """
    cQ, dQ = Q[1:-1], (Q[2:] - Q[:-2]) / 2
    cQ, dQ = cQ[strict[1:-1]], dQ[strict[1:-1]]

    idx = np.argsort(-dQ / cQ)[np.floor(dQ.shape[0] * 0.05).astype(int)]
    K = - cQ[idx] / dQ[idx]
    return np.exp(-1 / K)

def param_calibrate(param_range, method, Q, b_LH, a):
    """
    Calibrates the parameters for a baseflow estimation method.

    Args:
        param_range (numpy.ndarray): The range of parameter values to test.
        method (callable): The baseflow estimation method to use.
        Q (numpy.ndarray): The discharge values.
        b_LH (numpy.ndarray): The low-flow baseflow values.
        a (float): The parameter for the baseflow estimation method.

    Returns:
        float: The optimal parameter value from the given range.
    """
    idx_rec = recession_period(Q)
    idx_oth = np.full(Q.shape[0], True)
    idx_oth[idx_rec] = False
    return param_calibrate_jit(param_range, method, Q, b_LH, a, idx_rec, idx_oth)

def recession_period(Q):
    """
    Identifies the recession periods in the discharge time series.

    The function takes the discharge time series `Q` as input and returns the indices of the beginning and end of the recession periods. The recession periods are identified by finding the local maxima in the 3-point moving average of the discharge time series. The function keeps only the recession periods that are at least 10 time steps long, and trims the beginning of each recession period by 60% of the recession period duration.

    Args:
        Q (numpy.ndarray): The discharge time series.

    Returns:
        numpy.ndarray: The indices of the beginning and end of the recession periods.
    """
    idx_dec = np.zeros(Q.shape[0] - 1, dtype=np.int64)
    Q_ave = moving_average(Q, 3)
    idx_dec[1:-1] = (Q_ave[:-1] - Q_ave[1:]) > 0
    idx_beg = np.where(idx_dec[:-1] - idx_dec[1:] == -1)[0] + 1
    idx_end = np.where(idx_dec[:-1] - idx_dec[1:] == 1)[0] + 1
    idx_keep = (idx_end - idx_beg) >= 10
    idx_beg = idx_beg[idx_keep]
    idx_end = idx_end[idx_keep]
    duration = idx_end - idx_beg
    idx_beg = idx_beg + np.ceil(duration * 0.6).astype(np.int64)
    return multi_arange(idx_beg, idx_end)

def maxmium_BFI(Q, b_LH, a, date=None):
    """
    Calculates the maximum baseflow index (BFI) for a given discharge time series.

    The function takes the discharge time series `Q`, the baseflow time series `b_LH`, and the recession coefficient `a` as input. It calculates the annual baseflow and discharge, and then computes the maximum BFI. If the maximum BFI is greater than 0.9, the function returns the ratio of the total baseflow to the total discharge instead.

    Args:
        Q (numpy.ndarray): The discharge time series.
        b_LH (numpy.ndarray): The baseflow time series.
        a (float): The recession coefficient.
        date (datetime.datetime, optional): The date associated with the discharge time series. If provided, the function will compute the annual BFI for each year.

    Returns:
        float: The maximum baseflow index.
    """
    b = backward(Q, b_LH, a)

    if date is None:
        idx_end = b.shape[0] // 365 * 365
        annual_b = np.mean(b[:idx_end].reshape(-1, 365), axis=1)
        annual_Q = np.mean(Q[:idx_end].reshape(-1, 365), axis=1)
        annual_BFI = annual_b / annual_Q
    else:
        idx_year = date.year - date.year.min()
        counts = np.bincount(idx_year)
        idx_valid = counts > 0
        annual_b = np.bincount(idx_year, weights=b)[idx_valid] / counts[idx_valid]
        annual_Q = np.bincount(idx_year, weights=Q)[idx_valid] / counts[idx_valid]
        annual_BFI = annual_b / annual_Q

    BFI_max = np.max(annual_BFI)
    BFI_max = BFI_max if BFI_max < 0.9 else np.sum(annual_b) / np.sum(annual_Q)
    return BFI_max

def bflow(Q, beta=0.925):
    """BFlow automated baseflow separation and recession analysis.

    Runs the Lyne-Hollick filter with 3 passes (the BFlow protocol)
    and performs recession analysis to compute the SWAT alpha factor.

    Arnold, J.G. and Allen, P.M. (1999) Automated methods for estimating
    baseflow and ground water recharge from streamflow records. JAWRA
    35(2), 411-424.

    Args:
        Q (numpy.ndarray): Daily streamflow.
        beta (float): Filter parameter. Default 0.925.

    Returns:
        dict: Dictionary with keys:
            - 'baseflow': numpy array of baseflow values (3-pass result)
            - 'alpha_factor': exponential recession constant (SWAT ALPHA_BF)
            - 'baseflow_days': days for recession to decline one log cycle
            - 'n_segments': number of recession segments used
            - 'BFI': baseflow fraction from the 3-pass filter
            - 'BFI_pass1': baseflow fraction from pass 1
            - 'BFI_pass2': baseflow fraction from pass 2
    """
    from baseflowx.separation import lh_multi

    total_Q = np.sum(Q)

    # Run 3 passes
    b1 = lh_multi(Q, beta=beta, num_pass=1)
    b2 = lh_multi(Q, beta=beta, num_pass=2)
    b3 = lh_multi(Q, beta=beta, num_pass=3)

    # Recession analysis on the 3-pass result
    recession = bflow_recession_analysis(Q, b3)

    return {
        'baseflow': b3,
        'alpha_factor': recession['alpha_factor'],
        'baseflow_days': recession['baseflow_days'],
        'n_segments': recession['n_segments'],
        'BFI': np.sum(b3) / total_Q if total_Q > 0 else np.nan,
        'BFI_pass1': np.sum(b1) / total_Q if total_Q > 0 else np.nan,
        'BFI_pass2': np.sum(b2) / total_Q if total_Q > 0 else np.nan,
    }

def bflow_recession_analysis(Q, b, ndmin=10, ndmax=300):
    """BFlow recession analysis for SWAT alpha factor estimation.

    Identifies baseflow recession segments, assembles a master recession
    curve, and fits the exponential recession constant (alpha factor).
    This maps directly to SWAT's ALPHA_BF parameter.

    Arnold, J.G. and Allen, P.M. (1999) Automated methods for estimating
    baseflow and ground water recharge from streamflow records. JAWRA
    35(2), 411-424.

    Args:
        Q (numpy.ndarray): Total streamflow.
        b (numpy.ndarray): Separated baseflow (e.g., from lh_multi).
        ndmin (int): Minimum recession segment length (days). Default 10.
        ndmax (int): Maximum recession segment length (days). Default 300.

    Returns:
        dict: Dictionary with keys:
            - 'alpha_factor': exponential recession constant
            - 'baseflow_days': days for recession to decline one log cycle
            - 'n_segments': number of recession segments used
            - 'BFI': baseflow fraction (sum(b)/sum(Q))
    """
    n = len(b)

    # Identify recession segments: periods where baseflow is declining
    declining = np.zeros(n, dtype=bool)
    declining[1:] = b[1:] < b[:-1]

    # Find contiguous declining segments
    segments = []
    seg_start = None
    for i in range(n):
        if declining[i]:
            if seg_start is None:
                seg_start = i - 1  # include the point before decline starts
        else:
            if seg_start is not None:
                seg_len = i - seg_start
                if ndmin <= seg_len <= ndmax:
                    segments.append((seg_start, i))
                seg_start = None
    # Handle segment ending at array end
    if seg_start is not None:
        seg_len = n - seg_start
        if ndmin <= seg_len <= ndmax:
            segments.append((seg_start, n))

    if len(segments) == 0:
        return {
            'alpha_factor': np.nan,
            'baseflow_days': np.nan,
            'n_segments': 0,
            'BFI': np.sum(b) / np.sum(Q) if np.sum(Q) > 0 else np.nan,
        }

    # Fit exponential recession to each segment and compute alpha
    # For each segment: b(t) = b(0) * exp(-alpha * t)
    # Taking log: ln(b(t)) = ln(b(0)) - alpha * t
    # Fit slope via least squares across all segments (master recession curve)
    all_t = []
    all_logb = []
    for start, end in segments:
        seg_b = b[start:end]
        valid = seg_b > 0
        if np.sum(valid) < 2:
            continue
        t_rel = np.arange(end - start)[valid].astype(float)
        logb = np.log(seg_b[valid])
        # Normalize: shift t so each segment starts at t=0
        all_t.append(t_rel)
        # Normalize: shift logb so each segment starts at logb=0
        all_logb.append(logb - logb[0])

    if len(all_t) == 0:
        return {
            'alpha_factor': np.nan,
            'baseflow_days': np.nan,
            'n_segments': len(segments),
            'BFI': np.sum(b) / np.sum(Q) if np.sum(Q) > 0 else np.nan,
        }

    t_all = np.concatenate(all_t)
    logb_all = np.concatenate(all_logb)

    # Least squares fit: logb = -alpha * t  (no intercept since normalized)
    # alpha = -sum(t * logb) / sum(t^2)
    alpha = -np.sum(t_all * logb_all) / (np.sum(t_all ** 2) + 1e-30)
    alpha = max(alpha, 1e-10)  # ensure positive

    # Baseflow days: days for one log cycle decline
    # One log cycle = ln(10), so alpha * BFD = ln(10)
    baseflow_days = np.log(10) / alpha

    BFI = np.sum(b) / np.sum(Q) if np.sum(Q) > 0 else np.nan

    return {
        'alpha_factor': alpha,
        'baseflow_days': baseflow_days,
        'n_segments': len(segments),
        'BFI': BFI,
    }

def cmb(Q, SC, SC_BF=None, SC_RO=None, sc_bf_percentile=99,
        sc_ro_percentile=1):
    """Conductivity Mass Balance baseflow separation (Stewart et al., 2007).

    Two-component mixing model using specific conductance as a tracer:
        Q_b(t) = Q(t) * (SC(t) - SC_RO) / (SC_BF - SC_RO)

    If end-member conductivities are not provided, they are estimated from
    the SC record using percentiles.

    Stewart, M.T., Cimino, J. and Ross, M. (2007) Calibration of base flow
    separation methods with streamflow conductivity. Groundwater 45(1), 17-27.

    Args:
        Q (numpy.ndarray): Streamflow.
        SC (numpy.ndarray): Concurrent specific conductance (uS/cm).
        SC_BF (float, optional): Baseflow end-member conductivity. If None,
            estimated as the sc_bf_percentile of SC.
        SC_RO (float, optional): Surface runoff end-member conductivity.
            If None, estimated as the sc_ro_percentile of SC.
        sc_bf_percentile (float): Percentile for SC_BF estimation. Default 99.
        sc_ro_percentile (float): Percentile for SC_RO estimation. Default 1.

    Returns:
        numpy.ndarray: Daily baseflow array.
    """
    if SC_BF is None or SC_RO is None:
        est_BF, est_RO = estimate_endmembers(
            SC, bf_percentile=sc_bf_percentile,
            ro_percentile=sc_ro_percentile)
        if SC_BF is None:
            SC_BF = est_BF
        if SC_RO is None:
            SC_RO = est_RO

    if SC_BF == SC_RO:
        raise ValueError("SC_BF and SC_RO are equal; cannot separate.")

    # Core CMB equation
    b = Q * (SC - SC_RO) / (SC_BF - SC_RO)

    # Physical constraints
    b = np.clip(b, 0.0, Q)

    # Propagate NaN from SC
    b[~np.isfinite(SC)] = np.nan

    return b

def estimate_endmembers(SC, bf_percentile=99, ro_percentile=1):
    """Estimate baseflow and runoff end-member conductivities.

    Args:
        SC (numpy.ndarray): Specific conductance time series (uS/cm).
        bf_percentile (float): Percentile for baseflow end-member. Default 99.
        ro_percentile (float): Percentile for runoff end-member. Default 1.

    Returns:
        tuple: (SC_BF, SC_RO) estimated end-member conductivities.
    """
    SC_valid = SC[np.isfinite(SC)]
    SC_BF = np.percentile(SC_valid, bf_percentile)
    SC_RO = np.percentile(SC_valid, ro_percentile)
    return SC_BF, SC_RO

def calibrate_eckhardt_from_cmb(Q, SC, a=None, SC_BF=None, SC_RO=None,
                                sc_bf_percentile=99, sc_ro_percentile=1):
    """Calibrate Eckhardt BFImax using CMB as a reference.

    Runs CMB to get a reference baseflow, computes BFI, and returns a
    BFImax suitable for use with eckhardt(). Optionally estimates the
    recession coefficient if not provided.

    Args:
        Q (numpy.ndarray): Streamflow.
        SC (numpy.ndarray): Concurrent specific conductance (uS/cm).
        a (float, optional): Recession coefficient. If None, estimated
            from the streamflow using strict_baseflow + recession_coefficient.
        SC_BF (float, optional): Baseflow end-member conductivity.
        SC_RO (float, optional): Surface runoff end-member conductivity.
        sc_bf_percentile (float): Percentile for SC_BF estimation. Default 99.
        sc_ro_percentile (float): Percentile for SC_RO estimation. Default 1.

    Returns:
        dict: Dictionary with keys:
            - 'BFImax': calibrated maximum baseflow index
            - 'BFI_cmb': BFI computed from CMB
            - 'a': recession coefficient used
            - 'SC_BF': baseflow end-member conductivity used
            - 'SC_RO': runoff end-member conductivity used
    """
    from baseflowx.separation import strict_baseflow
    from baseflowx.estimate import recession_coefficient

    # Get CMB baseflow
    b_cmb = cmb(Q, SC, SC_BF=SC_BF, SC_RO=SC_RO,
                sc_bf_percentile=sc_bf_percentile,
                sc_ro_percentile=sc_ro_percentile)

    # Use only valid (non-NaN) pairs
    valid = np.isfinite(b_cmb) & np.isfinite(Q) & (Q > 0)
    BFI_cmb = np.sum(b_cmb[valid]) / np.sum(Q[valid])
    BFI_cmb = np.clip(BFI_cmb, 0.01, 0.99)

    # Estimate recession coefficient if not provided
    if a is None:
        strict = strict_baseflow(Q)
        a = recession_coefficient(Q, strict)

    # Resolve end-members used
    if SC_BF is None or SC_RO is None:
        est_BF, est_RO = estimate_endmembers(
            SC, bf_percentile=sc_bf_percentile,
            ro_percentile=sc_ro_percentile)
        SC_BF = SC_BF if SC_BF is not None else est_BF
        SC_RO = SC_RO if SC_RO is not None else est_RO

    return {
        'BFImax': BFI_cmb,
        'BFI_cmb': BFI_cmb,
        'a': a,
        'SC_BF': SC_BF,
        'SC_RO': SC_RO,
    }

def fetch_usgs(site_id, start_date, end_date, parameter='discharge'):
    """Fetch daily values from the USGS NWIS Water Services API.

    Args:
        site_id (str): USGS site number (e.g., '01013500').
        start_date (str): Start date as 'YYYY-MM-DD'.
        end_date (str): End date as 'YYYY-MM-DD'.
        parameter (str): Parameter to fetch. Accepts convenience names
            ('discharge', 'streamflow', 'sc', 'specific_conductance') or
            a raw NWIS parameter code (e.g., '00060'). Default 'discharge'.

    Returns:
        dict: Dictionary with keys:
            - 'dates': numpy array of datetime.date objects
            - 'values': numpy array of float64 values (NaN for missing)
            - 'units': str describing the measurement units
            - 'site_id': str echoing the requested site
            - 'parameter': str echoing the parameter code used
            - 'qualifiers': list of qualifier strings per timestep

    Raises:
        ValueError: If the parameter name is unrecognized or no data is found.
        ConnectionError: If the NWIS API request fails.

    Example:
        >>> from baseflowx.io import fetch_usgs
        >>> data = fetch_usgs('01013500', '2015-01-01', '2015-12-31')
        >>> data['values'][:5]
        array([...])
    """
    # Resolve parameter code
    param_code = _PARAM_CODES.get(parameter, parameter)
    if not param_code.isdigit():
        raise ValueError(
            f"Unrecognized parameter: {parameter!r}. Use 'discharge', 'sc', "
            f"or a numeric NWIS code like '00060'.")

    # Build request URL
    params = {
        'format': 'json',
        'sites': site_id,
        'startDT': start_date,
        'endDT': end_date,
        'parameterCd': param_code,
        'siteStatus': 'all',
    }
    url = _NWIS_DV_URL + '?' + urlencode(params)

    # Fetch data
    try:
        with urlopen(url, timeout=60) as response:
            raw = json.loads(response.read().decode('utf-8'))
    except Exception as exc:
        raise ConnectionError(f"NWIS request failed: {exc}") from exc

    # Parse JSON response
    try:
        ts = raw['value']['timeSeries']
        if not ts:
            raise ValueError(f"No data returned for site {site_id}, "
                             f"parameter {param_code}, "
                             f"{start_date} to {end_date}.")
        series = ts[0]
        values_list = series['values'][0]['value']
        unit = series['variable']['unit']['unitCode']
    except (KeyError, IndexError) as exc:
        raise ValueError(
            f"Unexpected NWIS response structure: {exc}") from exc

    # Extract dates, values, and qualifiers
    dates = []
    values = []
    qualifiers = []
    for entry in values_list:
        dates.append(datetime.strptime(
            entry['dateTime'][:10], '%Y-%m-%d').date())
        val = entry['value']
        # NWIS uses -999999 or empty string for missing
        if val == '' or val == '-999999' or val is None:
            values.append(np.nan)
        else:
            values.append(float(val))
        qualifiers.append(entry.get('qualifiers', ['']))

    return {
        'dates': np.array(dates),
        'values': np.array(values, dtype=np.float64),
        'units': unit,
        'site_id': site_id,
        'parameter': param_code,
        'qualifiers': qualifiers,
    }

def load_sample_data():
    """Load the bundled Fish River sample dataset.

    Returns daily discharge for USGS site 01013500 (Fish River near Fort Kent,
    Maine) for 2019-2020.

    Returns:
        dict: Dictionary with keys:
            - 'dates': numpy array of datetime.date objects
            - 'Q': numpy array of discharge in ft³/s
            - 'site_id': '01013500'
            - 'units': 'ft3/s'
    """
    path = _DATA_DIR / 'fish_river.csv'
    dates = []
    values = []
    with open(path, newline='') as f:
        reader = csv.DictReader(f)
        for row in reader:
            dates.append(date.fromisoformat(row['date']))
            val = row['discharge_cfs']
            values.append(float(val) if val else np.nan)
    return {
        'dates': np.array(dates),
        'Q': np.array(values, dtype=np.float64),
        'site_id': '01013500',
        'units': 'ft3/s',
    }

def clean_streamflow(series):
    """
    Cleans a streamflow time series by removing invalid values and keeping only years with at least 120 data points.

    Args:
        series (pandas.Series): The streamflow time series to be cleaned.

    Returns:
        tuple: A tuple containing the cleaned streamflow values and the corresponding dates.
    """
    date, Q = series.index, series.values.astype(float)
    has_value = np.isfinite(Q)
    date, Q = date[has_value], np.abs(Q[has_value])
    year_unique, counts = np.unique(date.year, return_counts=True)
    keep = np.isin(date.year, year_unique[counts >= 120])
    return Q[keep], date[keep]

def moving_average(x, w):
    res = np.convolve(x, np.ones(w)) / w
    return res[w - 1:-w + 1]

def backward(Q, b_LH, a):
    """
    Calculates the baseflow time series `b` from the discharge time series `Q` and the baseflow time series `b_LH` using a backward recursive approach.

    The function iterates through the discharge time series in reverse order, calculating the baseflow at each time step based on the baseflow at the next time step and the recession coefficient `a`. If the calculated baseflow exceeds the discharge at the current time step, the baseflow is set to the discharge.

    Args:
        Q (numpy.ndarray): The discharge time series.
        b_LH (numpy.ndarray): The baseflow time series.
        a (float): The recession coefficient.

    Returns:
        numpy.ndarray: The baseflow time series.
    """
    b = np.zeros(Q.shape[0])
    b[-1] = b_LH[-1]
    for i in range(Q.shape[0] - 1, 0, -1):
        b[i - 1] = b[i] / a
        if b[i] == 0:
            b[i - 1] = Q[i - 1]
        if b[i - 1] > Q[i - 1]:
            b[i - 1] = Q[i - 1]
    return b