sejeteralo/frontend/app/utils/bezier-math.ts

/**
 * TypeScript port of the Bézier tariff math engine.
 *
 * Mirrors backend/app/engine/integrals.py and pricing.py.
 * Uses Cardano's formula + Newton-Raphson polish for cubic solving.
 */

// ── Cubic solver ──

/**
 * Solve ax³ + bx² + cx + d = 0 for real roots in [0, 1].
 * Uses Cardano's method with Newton-Raphson refinement.
 */
function solveCubicInUnit(a: number, b: number, c: number, d: number): number | null {
  if (Math.abs(a) < 1e-12) {
    // Degenerate: quadratic
    if (Math.abs(b) < 1e-12) {
      // Linear
      if (Math.abs(c) < 1e-12) return null
      const t = -d / c
      return t >= -1e-10 && t <= 1 + 1e-10 ? clamp01(t) : null
    }
    const disc = c * c - 4 * b * d
    if (disc < 0) return null
    const sqrtDisc = Math.sqrt(disc)
    const t1 = (-c + sqrtDisc) / (2 * b)
    const t2 = (-c - sqrtDisc) / (2 * b)
    if (t1 >= -1e-10 && t1 <= 1 + 1e-10) return clamp01(t1)
    if (t2 >= -1e-10 && t2 <= 1 + 1e-10) return clamp01(t2)
    return null
  }

  // Normalize: t³ + pt² + qt + r = 0
  const p = b / a
  const q = c / a
  const r = d / a

  // Depressed cubic: u³ + pu + q = 0 via substitution t = u - p/3
  const p1 = q - p * p / 3
  const q1 = r - p * q / 3 + 2 * p * p * p / 27

  const discriminant = q1 * q1 / 4 + p1 * p1 * p1 / 27

  const roots: number[] = []

  if (discriminant > 1e-12) {
    // One real root
    const sqrtD = Math.sqrt(discriminant)
    const u = cbrt(-q1 / 2 + sqrtD)
    const v = cbrt(-q1 / 2 - sqrtD)
    roots.push(u + v - p / 3)
  } else if (discriminant < -1e-12) {
    // Three real roots (casus irreducibilis)
    const m = Math.sqrt(-p1 / 3)
    const theta = Math.acos((-q1 / 2) / (m * m * m)) / 3
    roots.push(
      2 * m * Math.cos(theta) - p / 3,
      2 * m * Math.cos(theta - 2 * Math.PI / 3) - p / 3,
      2 * m * Math.cos(theta - 4 * Math.PI / 3) - p / 3,
    )
  } else {
    // Double or triple root
    const u = cbrt(-q1 / 2)
    roots.push(2 * u - p / 3, -u - p / 3)
  }

  // Find root in [0,1] and refine with Newton-Raphson
  for (const root of roots) {
    if (root >= -0.1 && root <= 1.1) {
      let t = clamp01(root)
      // Newton-Raphson polish (3 iterations)
      for (let i = 0; i < 3; i++) {
        const f = ((a * t + b) * t + c) * t + d
        const fp = (3 * a * t + 2 * b) * t + c
        if (Math.abs(fp) < 1e-14) break
        t = clamp01(t - f / fp)
      }
      return t
    }
  }

  return null
}

function cbrt(x: number): number {
  return x < 0 ? -Math.pow(-x, 1 / 3) : Math.pow(x, 1 / 3)
}

function clamp01(t: number): number {
  return Math.max(0, Math.min(1, t))
}

// ── Integral computation ──

export interface IntegralResult {
  alpha1: number
  alpha2: number
  beta2: number
}

export function computeIntegrals(
  volume: number,
  vinf: number,
  vmax: number,
  pmax: number,
  a: number,
  b: number,
  c: number,
  d: number,
  e: number,
): IntegralResult {
  if (volume <= vinf) {
    const T = solveTier1T(volume, vinf, b)
    const alpha1 = computeAlpha1(T, vinf, a, b)
    return { alpha1, alpha2: 0, beta2: 0 }
  } else {
    const alpha1 = computeAlpha1(1.0, vinf, a, b)
    const wmax = vmax - vinf
    const T = solveTier2T(volume - vinf, wmax, c, d)
    const uu = computeUU(T, c, d, e)
    const alpha2 = (volume - vinf) - 3 * uu * wmax
    const beta2 = 3 * pmax * wmax * uu
    return { alpha1, alpha2, beta2 }
  }
}

function solveTier1T(volume: number, vinf: number, b: number): number {
  if (volume <= 0) return 0
  if (volume >= vinf) return 1
  const ratio = volume / vinf
  const t = solveCubicInUnit(1 - 3 * b, 3 * b, 0, -ratio)
  return t ?? 0
}

function solveTier2T(w: number, wmax: number, c: number, d: number): number {
  if (w <= 0) return 0
  if (w >= wmax) return 1
  const ratio = w / wmax
  const t = solveCubicInUnit(
    3 * (c + d - c * d) - 2,
    3 * (1 - 2 * c - d + c * d),
    3 * c,
    -ratio,
  )
  return t ?? 0
}

function computeAlpha1(T: number, vinf: number, a: number, b: number): number {
  return 3 * vinf * (
    Math.pow(T, 6) / 6 * (-9 * a * b + 3 * a + 6 * b - 2) +
    Math.pow(T, 5) / 5 * (24 * a * b - 6 * a - 13 * b + 3) +
    3 * Math.pow(T, 4) / 4 * (-7 * a * b + a + 2 * b) +
    Math.pow(T, 3) / 3 * 6 * a * b
  )
}

function computeUU(T: number, c: number, d: number, e: number): number {
  return (
    (-3 * c * d + 9 * e * c * d + 3 * c - 9 * e * c + 3 * d - 9 * e * d + 6 * e - 2) * Math.pow(T, 6) / 6 +
    (2 * c * d - 15 * e * c * d - 4 * c + 21 * e * c - 2 * d + 15 * e * d - 12 * e + 2) * Math.pow(T, 5) / 5 +
    (6 * e * c * d + c - 15 * e * c - 6 * e * d + 6 * e) * Math.pow(T, 4) / 4 +
    (3 * e * c) * Math.pow(T, 3) / 3
  )
}

// ── Pricing computation ──

export interface HouseholdData {
  volume_m3: number
  status: string
}

export interface PricingResult {
  p0: number
  curveVolumes: number[]
  curvePricesM3: number[]
}

export interface ImpactRow {
  volume: number
  oldPrice: number
  newPriceRP: number
  newPriceRS: number
}

export function computeP0(
  households: HouseholdData[],
  recettes: number,
  abop: number,
  abos: number,
  vinf: number,
  vmax: number,
  pmax: number,
  a: number,
  b: number,
  c: number,
  d: number,
  e: number,
): number {
  let totalAbo = 0
  let totalAlpha = 0
  let totalBeta = 0

  for (const h of households) {
    const abo = h.status === 'RS' ? abos : abop
    totalAbo += abo
    const vol = Math.max(h.volume_m3, 1e-5)
    const { alpha1, alpha2, beta2 } = computeIntegrals(vol, vinf, vmax, pmax, a, b, c, d, e)
    totalAlpha += alpha1 + alpha2
    totalBeta += beta2
  }

  if (totalAbo >= recettes) return 0
  if (totalAlpha === 0) return 0

  return (recettes - totalAbo - totalBeta) / totalAlpha
}

/**
 * Generate price curve points (price per m³ vs volume).
 */
export function generateCurve(
  vinf: number,
  vmax: number,
  pmax: number,
  p0: number,
  a: number,
  b: number,
  c: number,
  d: number,
  e: number,
  nbpts: number = 200,
): PricingResult {
  const curveVolumes: number[] = []
  const curvePricesM3: number[] = []

  const dt = 1 / (nbpts - 1)

  // Tier 1
  for (let i = 0; i < nbpts; i++) {
    const t = Math.min(i * dt, 1 - 1e-6)
    const v = vinf * ((1 - 3 * b) * t * t * t + 3 * b * t * t)
    const p = p0 * ((3 * a - 2) * t * t * t + (-6 * a + 3) * t * t + 3 * a * t)
    curveVolumes.push(v)
    curvePricesM3.push(p)
  }

  // Tier 2
  for (let i = 0; i < nbpts; i++) {
    const t = Math.min(i * dt, 1 - 1e-6)
    const v = vinf + (vmax - vinf) * (
      (3 * (c + d - c * d) - 2) * t * t * t +
      3 * (1 - 2 * c - d + c * d) * t * t +
      3 * c * t
    )
    const p = p0 + (pmax - p0) * ((1 - 3 * e) * t * t * t + 3 * e * t * t)
    curveVolumes.push(v)
    curvePricesM3.push(p)
  }

  return { p0, curveVolumes, curvePricesM3 }
}

/**
 * Compute price impacts at reference volume levels.
 */
export function computeImpacts(
  households: HouseholdData[],
  recettes: number,
  abop: number,
  abos: number,
  vinf: number,
  vmax: number,
  pmax: number,
  a: number,
  b: number,
  c: number,
  d: number,
  e: number,
  referenceVolumes: number[] = [30, 60, 90, 150, 300, 600],
): { p0: number; impacts: ImpactRow[] } {
  const p0 = computeP0(households, recettes, abop, abos, vinf, vmax, pmax, a, b, c, d, e)

  // Linear baseline
  const totalVol = households.reduce((s, h) => s + Math.max(h.volume_m3, 1e-5), 0)
  const totalAbo = households.reduce((s, h) => s + (h.status === 'RS' ? abos : abop), 0)
  const oldPM3 = totalVol > 0 ? (recettes - totalAbo) / totalVol : 0

  const impacts: ImpactRow[] = referenceVolumes.map((vol) => {
    const { alpha1, alpha2, beta2 } = computeIntegrals(vol, vinf, vmax, pmax, a, b, c, d, e)
    return {
      volume: vol,
      oldPrice: abop + oldPM3 * vol,
      newPriceRP: abop + (alpha1 + alpha2) * p0 + beta2,
      newPriceRS: abos + (alpha1 + alpha2) * p0 + beta2,
    }
  })

  return { p0, impacts }
}

// ── Control point mapping ──

export interface ControlPoints {
  // Tier 1: P1(fixed), P2, P3, P4
  p1: { x: number; y: number }
  p2: { x: number; y: number }
  p3: { x: number; y: number }
  p4: { x: number; y: number }
  // Tier 2: P4(shared), P5, P6, P7(fixed)
  p5: { x: number; y: number }
  p6: { x: number; y: number }
  p7: { x: number; y: number }
}

export function paramsToControlPoints(
  vinf: number,
  vmax: number,
  pmax: number,
  p0: number,
  a: number,
  b: number,
  c: number,
  d: number,
  e: number,
): ControlPoints {
  return {
    p1: { x: 0, y: 0 },
    p2: { x: 0, y: a * p0 },
    p3: { x: b * vinf, y: p0 },
    p4: { x: vinf, y: p0 },
    p5: { x: vinf + c * (vmax - vinf), y: p0 },
    p6: {
      x: vinf + (vmax - vinf) * (1 - d + c * d),
      y: p0 + e * (pmax - p0),
    },
    p7: { x: vmax, y: pmax },
  }
}

export function controlPointsToParams(
  cp: ControlPoints,
  vmax: number,
  pmax: number,
  p0: number,
): { vinf: number; a: number; b: number; c: number; d: number; e: number } {
  const vinf = cp.p4.x

  const a = p0 > 0 ? clamp01(cp.p2.y / p0) : 0.5
  const b = vinf > 0 ? clamp01(cp.p3.x / vinf) : 0.5

  const wmax = vmax - vinf
  const c = wmax > 0 ? clamp01((cp.p5.x - vinf) / wmax) : 0.5

  const qmax = pmax - p0
  const e = qmax > 0 ? clamp01((cp.p6.y - p0) / qmax) : 0.5

  // d from p6.x: x6 = vinf + wmax * (1 - d + c*d) => d = (1 - (x6-vinf)/wmax) / (1-c)
  let d_val = 0.5
  if (wmax > 0) {
    const ratio = (cp.p6.x - vinf) / wmax
    if (Math.abs(1 - c) > 1e-10) {
      d_val = clamp01((1 - ratio) / (1 - c))
    }
  }

  return { vinf, a, b, c, d: d_val, e }
}