How the numbers are actually computed

import {
  distanceMi,
  distanceBasedConfidence,
  idwInterpolate,
  type InterpolatedPoint,
} from "@/lib/interpolation";

export type InterpolationSample = { lat: number; lon: number; value: number };

export type VariogramModel = {
  type: "spherical" | "exponential";
  /** Semivariance at h=0+ (measurement/micro-scale error), feet². */
  nugget: number;
  /** Semivariance at h >= range — the total sill, feet². */
  sill: number;
  /** Miles. For the exponential model this is the "practical range" (h at
   *  which ~95% of the sill is reached), not the raw decay constant. */
  range: number;
};

export type KrigedPoint = InterpolatedPoint & {
  /** Ordinary kriging estimation variance, feet² — the real basis for
   *  `confidence`, kept separate in case it's useful on its own. */
  variance: number;
};

const NUM_LAG_BINS = 15;
const MIN_SAMPLES_TO_FIT = 10;
const MIN_BINS_TO_FIT = 4;
const NUM_RANGE_CANDIDATES = 25;

/**
 * Normalized variogram shape, 0 at h=0 rising to 1 at h >= range — the
 * nugget/sill scaling is applied by the caller so this can be shared
 * between fitting (candidate ranges) and evaluation (the fitted model).
 */
function shapeFn(h: number, range: number, type: VariogramModel["type"]): number {
  if (type === "spherical") {
    if (h >= range) return 1;
    const hr = h / range;
    return 1.5 * hr - 0.5 * hr ** 3;
  }
  // Exponential, "practical range" parameterization (matches the common
  // gstat/ArcGIS convention: g(range) ≈ 0.95, not the raw 1/e decay).
  return 1 - Math.exp((-3 * h) / range);
}

/**
 * The theoretical semivariogram, γ(h). Deliberately special-cases h≈0 to
 * return exactly 0 rather than the nugget: the nugget is a discontinuity
 * AT the origin (micro-scale variation between arbitrarily close points),
 * not the variogram's value at zero separation — the standard convention,
 * and what makes kriging an exact interpolator at a data point itself.
 */
function gammaAt(h: number, model: VariogramModel): number {
  if (h < 1e-6) return 0;
  return model.nugget + (model.sill - model.nugget) * shapeFn(h, model.range, model.type);
}

type VariogramBin = { lag: number; semivariance: number; count: number };

function computeEmpiricalVariogram(
  samples: InterpolationSample[]
): { bins: VariogramBin[]; maxLag: number } {
  const pairSquaredDiffs: { h: number; sv: number }[] = [];
  for (let i = 0; i < samples.length; i++) {
    for (let j = i + 1; j < samples.length; j++) {
      const h = distanceMi(samples[i].lat, samples[i].lon, samples[j].lat, samples[j].lon);
      if (h < 1e-6) continue; // coincident wells carry no lag information
      pairSquaredDiffs.push({ h, sv: 0.5 * (samples[i].value - samples[j].value) ** 2 });
    }
  }
  if (pairSquaredDiffs.length === 0) return { bins: [], maxLag: 0 };

  // Standard geostatistical practice: only trust the empirical variogram
  // out to roughly half the data's maximum extent — pairs get sparse and
  // the estimate gets noisy beyond that.
  const maxPairDist = pairSquaredDiffs.reduce((m, p) => Math.max(m, p.h), 0);
  const maxLag = maxPairDist * 0.5;
  const binWidth = maxLag / NUM_LAG_BINS;
  if (binWidth <= 0) return { bins: [], maxLag: 0 };

  const sums = new Array(NUM_LAG_BINS).fill(0);
  const counts = new Array(NUM_LAG_BINS).fill(0);
  for (const { h, sv } of pairSquaredDiffs) {
    if (h > maxLag) continue;
    const bin = Math.min(NUM_LAG_BINS - 1, Math.floor(h / binWidth));
    sums[bin] += sv;
    counts[bin] += 1;
  }

  const bins: VariogramBin[] = [];
  for (let b = 0; b < NUM_LAG_BINS; b++) {
    if (counts[b] === 0) continue;
    bins.push({ lag: (b + 0.5) * binWidth, semivariance: sums[b] / counts[b], count: counts[b] });
  }
  return { bins, maxLag };
}

type FittedCandidate = {
  type: VariogramModel["type"];
  nugget: number;
  sill: number;
  range: number;
  weightedSSE: number;
};

/**
 * Fits nugget/sill for a FIXED candidate range via weighted least squares.
 * For a fixed range, γ(h) = nugget·(1-g) + sill·g is linear in (nugget,
 * sill) — a plain 2-variable weighted regression, no iterative nonlinear
 * solver needed. Trying a spread of candidate ranges and keeping the best
 * (see fitModelForType) is a standard, auditable way to handle the one
 * genuinely nonlinear parameter without pulling in an optimization library.
 */
function fitNuggetSillForRange(
  bins: VariogramBin[],
  range: number,
  type: VariogramModel["type"]
): { nugget: number; sill: number } | null {
  let Sww11 = 0;
  let Sww12 = 0;
  let Sww22 = 0;
  let Swy1 = 0;
  let Swy2 = 0;
  for (const { lag, semivariance, count } of bins) {
    const g = shapeFn(lag, range, type);
    const x1 = 1 - g;
    const x2 = g;
    Sww11 += count * x1 * x1;
    Sww12 += count * x1 * x2;
    Sww22 += count * x2 * x2;
    Swy1 += count * x1 * semivariance;
    Swy2 += count * x2 * semivariance;
  }

  const det = Sww11 * Sww22 - Sww12 * Sww12;
  let nugget: number;
  let sill: number;
  if (Math.abs(det) < 1e-10) return null;
  nugget = (Swy1 * Sww22 - Swy2 * Sww12) / det;
  sill = (Sww11 * Swy2 - Sww12 * Swy1) / det;

  // Physically, nugget can't be negative and sill must exceed it. If the
  // unconstrained fit violates that, re-fit with nugget pinned to 0 (a
  // 1-variable regression for sill alone) rather than accept nonsense.
  if (nugget < 0 || sill <= nugget) {
    nugget = 0;
    let Swgg = 0;
    let Swgy = 0;
    for (const { lag, semivariance, count } of bins) {
      const g = shapeFn(lag, range, type);
      Swgg += count * g * g;
      Swgy += count * g * semivariance;
    }
    if (Swgg < 1e-10) return null;
    sill = Swgy / Swgg;
    if (sill <= 0) return null;
  }

  return { nugget, sill };
}

function fitModelForType(
  bins: VariogramBin[],
  type: VariogramModel["type"],
  maxLag: number
): FittedCandidate | null {
  const minRange = bins[0].lag;
  let best: FittedCandidate | null = null;

  for (let i = 0; i < NUM_RANGE_CANDIDATES; i++) {
    const t = i / (NUM_RANGE_CANDIDATES - 1);
    // Search up to 1.5x the fitted-lag window — the true range can exceed
    // the window we trusted for the empirical estimate itself.
    const range = minRange + t * (maxLag * 1.5 - minRange);
    if (range <= 0) continue;

    const fit = fitNuggetSillForRange(bins, range, type);
    if (!fit) continue;

    let weightedSSE = 0;
    for (const { lag, semivariance, count } of bins) {
      const predicted = fit.nugget + (fit.sill - fit.nugget) * shapeFn(lag, range, type);
      weightedSSE += count * (semivariance - predicted) ** 2;
    }

    if (!best || weightedSSE < best.weightedSSE) {
      best = { type, nugget: fit.nugget, sill: fit.sill, range, weightedSSE };
    }
  }

  return best;
}

/**
 * Fits an ordinary-kriging variogram model from real well data — spherical
 * and exponential are both tried, and whichever gives the lower weighted
 * SSE against the empirical (binned) semivariogram wins. Returns null when
 * there isn't enough data to fit anything meaningful (too few wells, or a
 * degenerate/flat empirical variogram) — callers should fall back to a
 * simpler method rather than krige with a fabricated model.
 */
export function fitVariogram(samples: InterpolationSample[]): VariogramModel | null {
  if (samples.length < MIN_SAMPLES_TO_FIT) return null;

  const { bins, maxLag } = computeEmpiricalVariogram(samples);
  if (bins.length < MIN_BINS_TO_FIT || maxLag <= 0) return null;

  const candidates = (["spherical", "exponential"] as const)
    .map((type) => fitModelForType(bins, type, maxLag))
    .filter((c): c is FittedCandidate => c !== null);

  if (candidates.length === 0) return null;

  const best = candidates.reduce((a, b) => (a.weightedSSE <= b.weightedSSE ? a : b));
  return { type: best.type, nugget: best.nugget, sill: best.sill, range: best.range };
}

const K_NEIGHBORS = 14;
// Ordinary kriging is well-posed with fewer, but the local system gets
// numerically flaky and the estimate loses meaning — below this, callers
// should fall back rather than trust the result.
const MIN_NEIGHBORS = 3;
const DUPLICATE_WELL_THRESHOLD_MI = 0.05;
// A spherical/exponential variogram is essentially saturated (γ ≈ sill)
// beyond its own range — a well at 3x the range carries the same "no
// correlation" signal as one at 30x, while making the local system worse
// conditioned (rows for mutually-saturated, far-apart wells become nearly
// identical). Bounding the neighbor search to this radius is standard
// local-kriging practice, not just a numerical patch: found empirically
// that without it, even well-supported points had kriging variance
// meeting or exceeding the sill everywhere, because the fixed "14 nearest
// regardless of distance" neighborhood was pulling in far, uninformative
// wells that destabilized the system for every query point.
const SEARCH_RADIUS_FACTOR = 1.5;

/** Gaussian elimination with partial pivoting on the augmented [A|b] system. */
function solveLinearSystem(A: number[][], b: number[]): number[] | null {
  const n = b.length;
  const M = A.map((row, i) => [...row, b[i]]);

  for (let col = 0; col < n; col++) {
    let pivotRow = col;
    let maxAbs = Math.abs(M[col][col]);
    for (let r = col + 1; r < n; r++) {
      if (Math.abs(M[r][col]) > maxAbs) {
        maxAbs = Math.abs(M[r][col]);
        pivotRow = r;
      }
    }
    if (maxAbs < 1e-10) return null; // singular — caller should handle null
    if (pivotRow !== col) {
      [M[col], M[pivotRow]] = [M[pivotRow], M[col]];
    }
    for (let r = col + 1; r < n; r++) {
      const factor = M[r][col] / M[col][col];
      for (let c = col; c <= n; c++) M[r][c] -= factor * M[col][c];
    }
  }

  const x = new Array(n).fill(0);
  for (let row = n - 1; row >= 0; row--) {
    let sum = M[row][n];
    for (let c = row + 1; c < n; c++) sum -= M[row][c] * x[c];
    x[row] = sum / M[row][row];
  }
  return x;
}

/**
 * Ordinary kriging at a single point, using a LOCAL neighborhood (up to
 * k=14 wells within SEARCH_RADIUS_FACTOR × the fitted range, not the full
 * statewide set) — solving an n×n system per query point is the expensive
 * part of kriging, and with up to ~2,500 query points per surface render,
 * keeping n small is what keeps this interactive. See facility-map.tsx /
 * this session's history for the O(gridPoints × wells × readingsPerWell)
 * bug this project already hit once with naive per-point IDW — this is
 * the same category of risk with a higher per-point cost, hence the small
 * bounded neighborhood.
 *
 * Falls back to plain IDW, per query point, when the local neighborhood
 * is too sparse (fewer than MIN_NEIGHBORS wells within range) or the
 * resulting system is singular — both are genuine "not enough real,
 * correlated data here" situations, not something a bigger or
 * differently-chosen neighborhood would fix (see SEARCH_RADIUS_FACTOR's
 * comment for why simply grabbing more distant wells makes this worse,
 * not better).
 *
 * Uses the classic Isaaks & Srivastava formulation directly in terms of
 * the semivariogram γ (no covariance conversion needed): solve
 *   Σ_j λ_j·γ(i,j) + μ = γ(i,0)  for each neighbor i, subject to Σ λ_j = 1
 * for the kriging weights λ and Lagrange multiplier μ, then
 *   estimate = Σ λ_i · value_i
 *   variance = Σ λ_i · γ(i,0) + μ
 */
export function krigeAt(
  targetLat: number,
  targetLon: number,
  samples: InterpolationSample[],
  model: VariogramModel
): KrigedPoint | null {
  if (samples.length === 0) return null;

  const withDist = samples
    .map((s) => ({ s, dist: distanceMi(targetLat, targetLon, s.lat, s.lon) }))
    .sort((a, b) => a.dist - b.dist);

  // A grid point essentially on top of a well reproduces that well's
  // reading exactly — kriging's exact-interpolator property, and the same
  // override IDW uses for coincident points.
  if (withDist[0].dist < 0.01) {
    return {
      depthFt: withDist[0].s.value,
      confidence: 1,
      nearestWellMi: withDist[0].dist,
      variance: 0,
    };
  }

  // Candidates within the search radius, deduped (two wells a few hundred
  // feet apart carry almost no independent information and can make the
  // local system numerically unstable), capped at K_NEIGHBORS.
  const searchRadiusMi = model.range * SEARCH_RADIUS_FACTOR;
  const neighbors: typeof withDist = [];
  for (const candidate of withDist) {
    if (candidate.dist > searchRadiusMi) break; // sorted ascending — nothing further qualifies either
    const isDuplicate = neighbors.some(
      (n) => distanceMi(n.s.lat, n.s.lon, candidate.s.lat, candidate.s.lon) < DUPLICATE_WELL_THRESHOLD_MI
    );
    if (isDuplicate) continue;
    neighbors.push(candidate);
    if (neighbors.length >= K_NEIGHBORS) break;
  }

  if (neighbors.length < MIN_NEIGHBORS) {
    // Not enough real, correlated data nearby for kriging to mean anything
    // here — fall back to plain IDW rather than force an unstable system
    // using far-away, uninformative wells. This is the same fallback
    // fitVariogram's null case uses at the dataset level, applied per
    // query point: genuinely undersampled areas get IDW's simpler, always-
    // graceful degradation instead of a numerically shaky kriging result.
    const idw = idwInterpolate(targetLat, targetLon, samples);
    return idw ? { ...idw, variance: model.sill } : null;
  }

  const k = neighbors.length;
  const A: number[][] = Array.from({ length: k + 1 }, () => new Array(k + 1).fill(0));
  const b: number[] = new Array(k + 1).fill(0);

  for (let i = 0; i < k; i++) {
    for (let j = 0; j < k; j++) {
      A[i][j] =
        i === j
          ? 0
          : gammaAt(
              distanceMi(neighbors[i].s.lat, neighbors[i].s.lon, neighbors[j].s.lat, neighbors[j].s.lon),
              model
            );
    }
    A[i][k] = 1;
    A[k][i] = 1;
    b[i] = gammaAt(neighbors[i].dist, model);
  }
  b[k] = 1;

  const solution = solveLinearSystem(A, b);
  if (!solution) {
    // Singular even after radius-bounding and dedup — same fallback as
    // the too-sparse case above, rather than leaving a hole in the surface.
    const idw = idwInterpolate(targetLat, targetLon, samples);
    return idw ? { ...idw, variance: model.sill } : null;
  }

  const weights = solution.slice(0, k);
  const mu = solution[k];

  let depthFt = 0;
  for (let i = 0; i < k; i++) depthFt += weights[i] * neighbors[i].s.value;

  let variance = mu;
  for (let i = 0; i < k; i++) variance += weights[i] * b[i];
  variance = Math.max(0, variance); // guard tiny negative numerical noise

  // Confidence uses the same distance-based heuristic as IDW rather than
  // variance/sill, by deliberate choice, not oversight. Measured against
  // this well network's actual fitted model, nugget/sill is consistently
  // ~0.7 — most of the semivariance is unexplained even at vanishingly
  // small separation, so γ(h) (and thus kriging variance) jumps to nearly
  // the sill within a fraction of a mile of any well. A variance-based
  // confidence would floor almost the entire surface, which is an honest
  // reflection of how noisy/weakly-correlated this data is, but unusable
  // as a smooth visual fade — most of the map would render at minimum
  // opacity regardless of well proximity. The real variance is still
  // returned below for anyone who wants the statistically rigorous value.
  const confidence = distanceBasedConfidence(withDist[0].dist);

  return { depthFt, confidence, nearestWellMi: withDist[0].dist, variance };
}

import { monthlyMedians, type Reading } from "@/lib/attribution/did";

/**
 * Per-well time-series projection: given a well's own reading history,
 * project forward to a future date. This is classical time-series
 * decomposition — a median-based seasonal index plus a Theil-Sen robust
 * trend line — not a trained model. It follows the same house rule as
 * kriging.ts and attribution/did.ts: return null rather than fabricate
 * a trend from too little data.
 *
 * The output is a plain depth value at a location, so it slots into
 * kriging.ts's existing InterpolationSample shape unchanged — a forecast
 * surface is just kriging fed projected values instead of measured ones.
 * This file has no dependency on, and is never imported by,
 * src/lib/attribution/* beyond reusing Reading/monthlyMedians as a type
 * and a deseasonalizing helper.
 */

const MS_PER_YEAR = 365.25 * 24 * 60 * 60 * 1000;

/**
 * Minimum distinct months of data required before fitting a trend. This
 * counts MONTHS, not raw readings — mirrors did.ts's MIN_MATCHED_MONTHS
 * reasoning: a hundred readings clustered in one week tell you nothing
 * about a multi-year trend, only how many times a well was checked.
 */
export const MIN_READINGS_TO_PROJECT = 8;
/** A "trend" fit from a record shorter than this is noise, not signal. */
export const MIN_SPAN_YEARS_TO_PROJECT = 2;
/**
 * A well observed for N years can support a projection out to roughly
 * this many multiples of N — beyond that, the well's own history simply
 * can't back the extrapolation and it's dropped rather than guessed.
 */
export const MAX_HORIZON_TO_SPAN_RATIO = 2;

function median(values: number[]): number {
  const sorted = [...values].sort((a, b) => a - b);
  const mid = Math.floor(sorted.length / 2);
  return sorted.length % 2 === 1
    ? sorted[mid]
    : (sorted[mid - 1] + sorted[mid]) / 2;
}

/** "YYYY-MM" -> fractional year, for a common numeric time axis. */
function monthKeyToYears(key: string): number {
  const year = Number(key.slice(0, 4));
  const month = Number(key.slice(5, 7));
  return year + (month - 1) / 12;
}

function yearsBetween(fromIso: string, toIso: string): number {
  return (Date.parse(toIso) - Date.parse(fromIso)) / MS_PER_YEAR;
}

export type TrendModel = {
  /** Deseasonalized trend, feet per year (positive = deepening/declining). */
  slopeFtPerYear: number;
  /** The well's own last reading — the projection's anchor point, not a
   *  regression intercept, so projecting to this exact date reproduces
   *  this exact value (the same "exact interpolator at a data point"
   *  property kriging.ts guarantees spatially). */
  anchorDate: string;
  anchorValue: number;
  /** "01".."12" -> typical deviation from this well's overall median.
   *  Months never observed are simply absent (treated as 0), not
   *  fabricated from neighboring months. */
  seasonalOffsets: Map<string, number>;
  observedSpanYears: number;
};

/**
 * Fits a well's deseasonalized trend from its own reading history.
 * Returns null when there isn't enough real data to fit anything
 * meaningful — callers should drop the well from the forecast rather
 * than krige with a fabricated trend.
 */
export function fitTrend(readings: Reading[]): TrendModel | null {
  const monthly = monthlyMedians(readings);
  if (monthly.size < MIN_READINGS_TO_PROJECT) return null;

  const entries = [...monthly.entries()].sort(([a], [b]) => (a < b ? -1 : 1));
  const times = entries.map(([key]) => monthKeyToYears(key));
  const observedSpanYears = times[times.length - 1] - times[0];
  if (observedSpanYears < MIN_SPAN_YEARS_TO_PROJECT) return null;

  const overallMedian = median(entries.map(([, v]) => v));

  // Seasonal index: each calendar month's typical deviation from this
  // well's overall median, median-based (same robustness choice did.ts
  // makes for its own per-month deltas) rather than a mean.
  const byMonthOfYear = new Map<string, number[]>();
  for (const [key, value] of entries) {
    const mm = key.slice(5, 7);
    const list = byMonthOfYear.get(mm);
    if (list) list.push(value - overallMedian);
    else byMonthOfYear.set(mm, [value - overallMedian]);
  }
  const seasonalOffsets = new Map<string, number>();
  for (const [mm, deviations] of byMonthOfYear) {
    seasonalOffsets.set(mm, median(deviations));
  }

  const deseasonalized = entries.map(
    ([key, value]) => value - (seasonalOffsets.get(key.slice(5, 7)) ?? 0)
  );

  // Theil-Sen: the median of all pairwise slopes. Same "median of
  // pairwise/bucketed quantities, no iterative solver" register as
  // did.ts's per-month deltas and kriging.ts's small fixed weighted-
  // least-squares system -- O(n^2) pairs is the same complexity class
  // as kriging's own O(wells^2) empirical variogram.
  const slopes: number[] = [];
  for (let i = 0; i < times.length; i++) {
    for (let j = i + 1; j < times.length; j++) {
      const dt = times[j] - times[i];
      if (dt < 1e-9) continue;
      slopes.push((deseasonalized[j] - deseasonalized[i]) / dt);
    }
  }
  if (slopes.length === 0) return null;
  const slopeFtPerYear = median(slopes);

  // Anchor at the well's own actual last reading (exact date/value, not
  // the last month's median) so the projection reproduces it exactly.
  const sortedReadings = [...readings].sort((a, b) => (a.date < b.date ? -1 : 1));
  const last = sortedReadings[sortedReadings.length - 1];

  return {
    slopeFtPerYear,
    anchorDate: last.date,
    anchorValue: last.value,
    seasonalOffsets,
    observedSpanYears,
  };
}

export type ForecastPoint = {
  value: number;
  horizonYears: number;
  /**
   * 0-1, purely about how much THIS WELL's own trend can be trusted at
   * this horizon -- independent of, and never folded into, kriging's
   * spatial confidence/variance (kriging.ts keeps those conceptually
   * distinct on purpose; this is a second, parallel axis of doubt, not
   * folded into either).
   */
  trendConfidence: number;
};

/**
 * Projects a fitted trend to a target date. Returns null when the
 * requested horizon exceeds what this well's own observed span can
 * support -- MAX_HORIZON_TO_SPAN_RATIO times its observed span -- rather
 * than extrapolating indefinitely from a short record.
 */
export function projectReading(
  model: TrendModel,
  targetDate: string
): ForecastPoint | null {
  const horizonYears = yearsBetween(model.anchorDate, targetDate);
  const maxHorizon = MAX_HORIZON_TO_SPAN_RATIO * model.observedSpanYears;
  if (Math.abs(horizonYears) > maxHorizon) return null;

  const targetMonth = targetDate.slice(5, 7);
  const anchorMonth = model.anchorDate.slice(5, 7);
  const seasonalDelta =
    (model.seasonalOffsets.get(targetMonth) ?? 0) -
    (model.seasonalOffsets.get(anchorMonth) ?? 0);

  const value = model.anchorValue + model.slopeFtPerYear * horizonYears + seasonalDelta;
  const trendConfidence = maxHorizon > 0 ? 1 - Math.abs(horizonYears) / maxHorizon : 0;

  return { value, horizonYears, trendConfidence };
}

/** Convenience: fitTrend + projectReading, null-propagating. */
export function projectWellSeries(
  readings: Reading[],
  targetDate: string
): ForecastPoint | null {
  const model = fitTrend(readings);
  return model ? projectReading(model, targetDate) : null;
}

/** Stable key for matching the same well between two InterpolationSample
 *  arrays (e.g. current vs. forecast values) by coordinate. */
export function wellKey(lat: number, lon: number): string {
  return `${lat},${lon}`;
}

/** A "range" from a single well can't express disagreement between wells. */
export const MIN_WELLS_FOR_SPREAD = 2;

export type ForecastSpread = {
  rangeLowFt: number;
  rangeHighFt: number;
  contributingWells: number;
};

/**
 * The spread of individual near wells' own projected change at a point --
 * NOT the single kriged/IDW point estimate, but how much the raw
 * contributing wells disagree with each other, matched to their current
 * readings by lat/lon. A real, auditable number (map-view.tsx's hover
 * readout shows it as "wells individually project X to Y ft"), not a
 * fabricated confidence interval. Returns null below MIN_WELLS_FOR_SPREAD
 * matched wells, same "nothing rather than invent a number" rule as the
 * rest of this file.
 */
export function computeForecastSpread(
  nearbyForecastWells: { lat: number; lon: number; forecastDepthFt: number }[],
  currentDepthByKey: Map<string, number>
): ForecastSpread | null {
  const changes: number[] = [];
  for (const well of nearbyForecastWells) {
    const currentDepthFt = currentDepthByKey.get(wellKey(well.lat, well.lon));
    if (currentDepthFt == null) continue;
    changes.push(well.forecastDepthFt - currentDepthFt);
  }
  if (changes.length < MIN_WELLS_FOR_SPREAD) return null;

  return {
    rangeLowFt: Math.min(...changes),
    rangeHighFt: Math.max(...changes),
    contributingWells: changes.length,
  };
}