import numpy as np


def reference_baselines(features, rewards, current_count, bandwidth=0.5, threshold=1e-3):
    x = np.array(features, dtype=np.float64)
    r = np.array(rewards, dtype=np.float64)
    n = len(r)
    if x.shape != (n, 2) or not 0 < current_count <= n or n < 2:
        raise ValueError("invalid window")
    if not np.isfinite(x).all() or not np.isfinite(r).all() or bandwidth <= 0:
        raise ValueError("invalid numbers")
    z = np.zeros_like(x)
    for k in range(2):
        mean = sum(x[j, k] for j in range(n)) / n
        std = np.sqrt(sum((x[j, k] - mean) ** 2 for j in range(n)) / n)
        for j in range(n):
            z[j, k] = (x[j, k] - mean) / std if std > 0 else 0
    baselines, neighbors = [], []
    for i in range(n - current_count, n):
        mass = np.float64(0)
        squared = np.float64(0)
        numerator = np.float64(0)
        other_rewards = np.float64(0)
        for j in range(n):
            if j == i:
                continue
            distance = sum((z[i, k] - z[j, k]) ** 2 for k in range(2))
            weight = np.exp(-distance / (2 * bandwidth**2))
            mass += weight
            squared += weight * weight
            numerator += weight * r[j]
            other_rewards += r[j]
        baselines.append(numerator / mass if mass >= threshold else other_rewards / (n - 1))
        neighbors.append(mass**2 / squared if squared > 0 else 0)
    return np.array(baselines), np.array(neighbors)