conformalSource.js

import { niceStep } from './grid.js';
import { MAX_LINES, SAMPLES } from './gl/conformalGridPass.js';

const PROBE_RINGS = 48;
const PROBE_ANGLES = 48;
const PROBE_MIN_RADIUS = 1e-4;
const PROBE_MAX_RADIUS = 1e4;

const REFINE_STEPS = 24;
const REFINE_GROWTH = 1.6;

const NEWTON_HIT_THRESHOLD = 24;
const NEWTON_SEEDS = 6;
const NEWTON_GRID = 20;

const VIEW_MARGIN = 0.35;

const TARGET_PX_SPACING = 110;
const REGION_PAD = 0.15;
const MAX_CLUSTERS = 4;

const VIEW_CELLS_X = 32;
const VIEW_CELLS_Y = 20;

const MIN_CLUSTER_COVERAGE = 0.02;

const MIN_CLUSTER_NOVELTY = 0.25;

const TARGET_SEGMENT_PX = 6;

const REGION_TRIM = 0.03;

const STEP_HYSTERESIS = 1.7;

const VIEW_EPSILON = 0.02;

function medianOf(values) {
  if (values.length === 0) return 0;
  const sorted = Float64Array.from(values).sort();
  return sorted[sorted.length >> 1];
}

function quantile(values, q) {
  const sorted = Float64Array.from(values).sort();
  const idx = Math.min(sorted.length - 1, Math.max(0, Math.round(q * (sorted.length - 1))));
  return sorted[idx];
}

function derivativeMagnitude(complexFn, re, im) {
  const h = Math.max(1e-7, Math.hypot(re, im) * 1e-5);
  const a = complexFn.evaluate(re, im);
  const b = complexFn.evaluate(re + h, im);
  if (!a || !b) return null;
  const d = Math.hypot(b.re - a.re, b.im - a.im) / h;
  return Number.isFinite(d) && d > 0 ? d : null;
}

function newtonStep(complexFn, re, im, target) {
  const a = complexFn.evaluate(re, im);
  if (!a) return null;
  const h = Math.max(1e-9, Math.hypot(re, im) * 1e-6);
  const b = complexFn.evaluate(re + h, im);
  if (!b) return null;
  const dRe = (b.re - a.re) / h;
  const dIm = (b.im - a.im) / h;
  const den = dRe * dRe + dIm * dIm;
  if (!(den > 0) || !Number.isFinite(den)) return null;
  const nRe = a.re - target.re;
  const nIm = a.im - target.im;
  const stepRe = (nRe * dRe + nIm * dIm) / den;
  const stepIm = (nIm * dRe - nRe * dIm) / den;
  const outRe = re - stepRe;
  const outIm = im - stepIm;
  if (!Number.isFinite(outRe) || !Number.isFinite(outIm)) return null;
  return { re: outRe, im: outIm };
}

function solveToward(complexFn, seed, target, iterations = 12) {
  let cur = seed;
  for (let i = 0; i < iterations; i++) {
    const next = newtonStep(complexFn, cur.re, cur.im, target);
    if (!next) return cur;
    const moved = Math.hypot(next.re - cur.re, next.im - cur.im);
    cur = next;
    if (moved <= Math.max(1e-12, Math.hypot(cur.re, cur.im) * 1e-12)) break;
  }
  return cur;
}

function clusterHits(hits, gridR, gridA) {
  const occupied = new Int32Array(gridR * gridA).fill(-1);
  const buckets = new Map();
  for (let i = 0; i < hits.length; i++) {
    const c = hits[i].ring * gridA + hits[i].angle;
    occupied[c] = 0;
    let list = buckets.get(c);
    if (!list) buckets.set(c, (list = []));
    list.push(i);
  }

  const label = new Int32Array(gridR * gridA).fill(-1);
  let labels = 0;
  const stack = [];
  for (let c = 0; c < gridR * gridA; c++) {
    if (occupied[c] !== 0 || label[c] !== -1) continue;
    const id = labels++;
    label[c] = id;
    stack.push(c);
    while (stack.length) {
      const cur = stack.pop();
      const ring = Math.floor(cur / gridA);
      const angle = cur - ring * gridA;
      for (let dr = -1; dr <= 1; dr++) {
        for (let da = -1; da <= 1; da++) {
          const nr = ring + dr;
          if (nr < 0 || nr >= gridR) continue;
          const na = (angle + da + gridA) % gridA;
          const next = nr * gridA + na;
          if (occupied[next] === 0 && label[next] === -1) {
            label[next] = id;
            stack.push(next);
          }
        }
      }
    }
  }

  if (labels === 0) return [];
  const clusters = Array.from({ length: labels }, () => []);
  for (const [cell, list] of buckets) {
    const id = label[cell];
    for (const i of list) clusters[id].push(hits[i]);
  }
  return clusters;
}

function clusterRegion(allMembers) {

  const first = new Map();
  const last = new Map();
  for (const m of allMembers) {
    if (!first.has(m.cell)) first.set(m.cell, m);
    last.set(m.cell, m);
  }
  const members = [...first.values(), ...last.values()];

  const res = members.map((m) => m.re);
  const ims = members.map((m) => m.im);
  if (members.length < 24) {
    return {
      minRe: Math.min(...res), maxRe: Math.max(...res),
      minIm: Math.min(...ims), maxIm: Math.max(...ims),
      midRe: medianOf(res), midIm: medianOf(ims),
    };
  }
  return {
    minRe: quantile(res, REGION_TRIM),
    maxRe: quantile(res, 1 - REGION_TRIM),
    minIm: quantile(ims, REGION_TRIM),
    maxIm: quantile(ims, 1 - REGION_TRIM),
    midRe: medianOf(res),
    midIm: medianOf(ims),
  };
}

function clipRegion(box, step, maxLines) {
  const maxSpan = maxLines * step;
  const out = { ...box };
  if (box.maxRe - box.minRe > maxSpan) {
    const c = Number.isFinite(box.midRe) ? box.midRe : (box.minRe + box.maxRe) / 2;
    out.minRe = c - maxSpan / 2;
    out.maxRe = c + maxSpan / 2;
  }
  if (box.maxIm - box.minIm > maxSpan) {
    const c = Number.isFinite(box.midIm) ? box.midIm : (box.minIm + box.maxIm) / 2;
    out.minIm = c - maxSpan / 2;
    out.maxIm = c + maxSpan / 2;
  }
  return out;
}

function padRegion(box, step) {

  const padRe = (box.maxRe - box.minRe) * REGION_PAD + step * 0.5;
  const padIm = (box.maxIm - box.minIm) * REGION_PAD + step * 0.5;
  return {
    minRe: box.minRe - padRe,
    maxRe: box.maxRe + padRe,
    minIm: box.minIm - padIm,
    maxIm: box.maxIm + padIm,
  };
}

const LOG_MIN_RADIUS = Math.log(PROBE_MIN_RADIUS);
const LOG_RADIUS_SPAN = Math.log(PROBE_MAX_RADIUS) - LOG_MIN_RADIUS;

function latticeOf(re, im) {
  const r = Math.max(Math.hypot(re, im), PROBE_MIN_RADIUS);
  const ring = Math.min(PROBE_RINGS - 1, Math.max(0, Math.round(
    ((Math.log(r) - LOG_MIN_RADIUS) / LOG_RADIUS_SPAN) * (PROBE_RINGS - 1))));
  let theta = Math.atan2(im, re);
  if (theta < 0) theta += 2 * Math.PI;
  const angle = Math.min(PROBE_ANGLES - 1, Math.floor((theta / (2 * Math.PI)) * PROBE_ANGLES));
  return { ring, angle };
}

export class ConformalSourcePlanner {
  constructor() {
    this.plan = null;
    this._expr = null;
    this._scale = 0;
    this._centerRe = 0;
    this._centerIm = 0;
  }

  _stale(viewport, expr) {
    if (!this.plan || expr !== this._expr) return true;
    if (Math.abs(Math.log(viewport.scale / this._scale)) > VIEW_EPSILON) return true;
    const spanRe = viewport.width / viewport.scale;
    const spanIm = viewport.height / viewport.scale;
    if (Math.abs(viewport.centerRe - this._centerRe) > spanRe * VIEW_EPSILON) return true;
    if (Math.abs(viewport.centerIm - this._centerIm) > spanIm * VIEW_EPSILON) return true;
    return false;
  }

  update(viewport, complexFn) {
    if (!complexFn.compiled) return null;
    if (!this._stale(viewport, complexFn.expr)) return this.plan;

    const next = this._build(viewport, complexFn);
    this._expr = complexFn.expr;
    this._scale = viewport.scale;
    this._centerRe = viewport.centerRe;
    this._centerIm = viewport.centerIm;
    this.plan = next;
    return next;
  }

  _build(viewport, complexFn) {
    const halfRe = viewport.width / viewport.scale / 2;
    const halfIm = viewport.height / viewport.scale / 2;
    const loRe = viewport.centerRe - halfRe * (1 + VIEW_MARGIN);
    const hiRe = viewport.centerRe + halfRe * (1 + VIEW_MARGIN);
    const loIm = viewport.centerIm - halfIm * (1 + VIEW_MARGIN);
    const hiIm = viewport.centerIm + halfIm * (1 + VIEW_MARGIN);

    const hits = [];

    const cellRep = new Map();

    const target = { re: viewport.centerRe, im: viewport.centerIm };

    const nearMisses = [];

    const probe = (re, im, ring, angle) => {
      const w = complexFn.evaluate(re, im);
      if (!w) return;
      if (w.re < loRe || w.re > hiRe || w.im < loIm || w.im > hiIm) {
        if (ring >= 0) {
          const d = Math.hypot(w.re - target.re, w.im - target.im);
          if (Number.isFinite(d)) nearMisses.push({ re, im, d });
        }
        return;
      }
      const cx = Math.min(VIEW_CELLS_X - 1, Math.max(0,
        Math.floor(((w.re - loRe) / (hiRe - loRe)) * VIEW_CELLS_X)));
      const cy = Math.min(VIEW_CELLS_Y - 1, Math.max(0,
        Math.floor(((w.im - loIm) / (hiIm - loIm)) * VIEW_CELLS_Y)));
      const cell = cy * VIEW_CELLS_X + cx;
      let r = ring;
      let a = angle;
      if (r < 0) ({ ring: r, angle: a } = latticeOf(re, im));
      const index = hits.length;
      hits.push({ re, im, ring: r, angle: a, cell });
      if (!cellRep.has(cell)) cellRep.set(cell, index);
    };

    const ratio = Math.pow(PROBE_MAX_RADIUS / PROBE_MIN_RADIUS, 1 / (PROBE_RINGS - 1));
    let radius = PROBE_MIN_RADIUS;
    for (let ring = 0; ring < PROBE_RINGS; ring++, radius *= ratio) {
      for (let a = 0; a < PROBE_ANGLES; a++) {
        const theta = (2 * Math.PI * a) / PROBE_ANGLES;
        probe(radius * Math.cos(theta), radius * Math.sin(theta), ring, a);
      }
    }

    if (cellRep.size < NEWTON_HIT_THRESHOLD && nearMisses.length > 0) {
      nearMisses.sort((a, b) => a.d - b.d);
      const seeds = [];
      for (const candidate of nearMisses) {
        if (seeds.length >= NEWTON_SEEDS) break;

        const apart = seeds.every(
          (sd) => Math.hypot(sd.re - candidate.re, sd.im - candidate.im) >
                  Math.hypot(candidate.re, candidate.im) * 0.25
        );
        if (apart) seeds.push(candidate);
      }
      for (const seed of seeds) {
        const root = solveToward(complexFn, seed, target);
        const d = derivativeMagnitude(complexFn, root.re, root.im);
        if (!d) continue;
        const spanRe = (halfRe * (1 + VIEW_MARGIN)) / d;
        const spanIm = (halfIm * (1 + VIEW_MARGIN)) / d;
        if (!Number.isFinite(spanRe) || !Number.isFinite(spanIm)) continue;
        for (let j = 0; j < NEWTON_GRID; j++) {
          const im = root.im - spanIm + (2 * spanIm * j) / (NEWTON_GRID - 1);
          for (let i = 0; i < NEWTON_GRID; i++) {
            const re = root.re - spanRe + (2 * spanRe * i) / (NEWTON_GRID - 1);
            probe(re, im, -1, 0);
          }
        }
      }
    }

    const prior = this.plan && this.plan.bounds;
    if (prior) {
      const cRe = (prior.minRe + prior.maxRe) / 2;
      const cIm = (prior.minIm + prior.maxIm) / 2;
      const hRe = ((prior.maxRe - prior.minRe) / 2) * REFINE_GROWTH || 1;
      const hIm = ((prior.maxIm - prior.minIm) / 2) * REFINE_GROWTH || 1;
      for (let j = 0; j < REFINE_STEPS; j++) {
        const im = cIm - hIm + (2 * hIm * j) / (REFINE_STEPS - 1);
        for (let i = 0; i < REFINE_STEPS; i++) {
          probe(cRe - hRe + (2 * hRe * i) / (REFINE_STEPS - 1), im, -1, 0);
        }
      }
    }

    if (cellRep.size === 0) return null;

    const derivatives = [];
    const reps = [...cellRep.values()];
    const stride = Math.max(1, Math.floor(reps.length / 96));
    for (let i = 0; i < reps.length; i += stride) {
      const h = hits[reps[i]];
      const d = derivativeMagnitude(complexFn, h.re, h.im);
      if (d !== null) derivatives.push(d);
    }
    const scaleFactor = medianOf(derivatives) || 1;

    let step = niceStep(TARGET_PX_SPACING / (viewport.scale * scaleFactor));
    if (!Number.isFinite(step) || step <= 0) return null;

    if (this.plan && this.plan.majorStep && this._expr === complexFn.expr) {
      const drift = step / this.plan.majorStep;
      if (drift < STEP_HYSTERESIS && drift > 1 / STEP_HYSTERESIS) {
        step = this.plan.majorStep;
      }
    }

    const ranked = clusterHits(hits, PROBE_RINGS, PROBE_ANGLES)
      .map((members) => {
        const cells = new Set();
        for (const m of members) cells.add(m.cell);
        return { members, coverage: cells.size };
      })
      .filter((c) => c.coverage >= Math.max(2, cellRep.size * MIN_CLUSTER_COVERAGE))
      .sort((a, b) => b.coverage - a.coverage);

    const claimed = new Set();
    const clusters = [];
    for (const c of ranked) {
      if (clusters.length >= MAX_CLUSTERS) break;
      const cells = new Set();
      for (const m of c.members) cells.add(m.cell);
      let fresh = 0;
      for (const cell of cells) if (!claimed.has(cell)) fresh++;
      if (clusters.length > 0 && fresh / cells.size < MIN_CLUSTER_NOVELTY) continue;
      for (const cell of cells) claimed.add(cell);
      clusters.push(c);
    }

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

    const boxes = clusters.map((c) => clusterRegion(c.members));

    const sampleLimitedSpan = (SAMPLES * TARGET_SEGMENT_PX) / (viewport.scale * scaleFactor);
    const spanLimit = Math.min(MAX_LINES - 4, Math.max(4, sampleLimitedSpan / step));
    const regions = boxes.map((b) => padRegion(clipRegion(b, step, spanLimit), step));

    const bounds = regions.reduce(
      (acc, r) => ({
        minRe: Math.min(acc.minRe, r.minRe),
        maxRe: Math.max(acc.maxRe, r.maxRe),
        minIm: Math.min(acc.minIm, r.minIm),
        maxIm: Math.max(acc.maxIm, r.maxIm),
      }),
      { minRe: Infinity, maxRe: -Infinity, minIm: Infinity, maxIm: -Infinity }
    );

    return { majorStep: step, regions, bounds };
  }
}