{"spec_id":"contour-density","library":"highcharts","language":"javascript","code":"// anyplot.ai\n// contour-density: Density Contour Plot\n// Library: highcharts 12.6.0 | JavaScript 22.23.2\n// Quality: 90/100 | Created: 2026-09-04\n//# anyplot-orientation: landscape\n\nconst THEME = window.ANYPLOT_THEME;\nconst t = window.ANYPLOT_TOKENS;\nconst INK_MUTED = THEME === \"light\" ? \"#6B6A63\" : \"#A8A79F\";\n\n// --- Deterministic RNG (LCG) + Box-Muller normal deviates -------------------\nfunction makeRng(seed) {\n  let state = seed >>> 0;\n  return function () {\n    state = (state * 1664525 + 1013904223) >>> 0;\n    return state / 4294967296;\n  };\n}\nfunction randomNormal(rng) {\n  const u1 = Math.max(rng(), 1e-12);\n  const u2 = rng();\n  return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);\n}\nfunction generateCluster(rng, n, meanX, meanY, sdX, sdY, rho) {\n  const points = [];\n  for (let k = 0; k < n; k++) {\n    const z1 = randomNormal(rng);\n    const z2 = randomNormal(rng);\n    points.push({\n      x: meanX + sdX * z1,\n      y: meanY + sdY * (rho * z1 + Math.sqrt(1 - rho * rho) * z2),\n    });\n  }\n  return points;\n}\n\n// --- Data: machined-part QC measurements, main batch + a tool-wear drift ---\nconst rng = makeRng(42);\nconst mainBatch = generateCluster(rng, 550, 25.0, 48.0, 0.14, 0.55, 0.55);\nconst driftBatch = generateCluster(rng, 150, 25.55, 49.1, 0.12, 0.45, 0.45);\nconst parts = mainBatch.concat(driftBatch);\n\n// --- 2D kernel density estimate on a grid -----------------------------------\nfunction std(values) {\n  const mean = values.reduce((a, b) => a + b, 0) / values.length;\n  const variance = values.reduce((a, b) => a + (b - mean) ** 2, 0) / values.length;\n  return Math.sqrt(variance);\n}\nfunction computeDensityGrid(pts, gx, gy, bwX, bwY) {\n  const ny = gy.length;\n  const nx = gx.length;\n  const grid = Array.from({ length: ny }, () => new Array(nx).fill(0));\n  pts.forEach((p) => {\n    for (let j = 0; j < ny; j++) {\n      const dy = (gy[j] - p.y) / bwY;\n      const wy = Math.exp(-0.5 * dy * dy);\n      if (wy < 1e-6) continue;\n      const row = grid[j];\n      for (let i = 0; i < nx; i++) {\n        const dx = (gx[i] - p.x) / bwX;\n        row[i] += wy * Math.exp(-0.5 * dx * dx);\n      }\n    }\n  });\n  return grid;\n}\n\nconst diameters = parts.map((p) => p.x);\nconst weights = parts.map((p) => p.y);\nconst n = parts.length;\nconst bwX = std(diameters) * Math.pow(n, -1 / 6);\nconst bwY = std(weights) * Math.pow(n, -1 / 6);\nconst xMin = Math.min(...diameters) - 3 * bwX;\nconst xMax = Math.max(...diameters) + 3 * bwX;\nconst yMin = Math.min(...weights) - 3 * bwY;\nconst yMax = Math.max(...weights) + 3 * bwY;\n\nconst GRID_N = 70;\nconst gridX = Array.from({ length: GRID_N }, (_, i) => xMin + (i * (xMax - xMin)) / (GRID_N - 1));\nconst gridY = Array.from({ length: GRID_N }, (_, j) => yMin + (j * (yMax - yMin)) / (GRID_N - 1));\nconst densityGrid = computeDensityGrid(parts, gridX, gridY, bwX, bwY);\n\nlet maxDensity = 0;\ndensityGrid.forEach((row) => row.forEach((v) => { if (v > maxDensity) maxDensity = v; }));\n\n// --- Marching squares: extract iso-density contour lines from the grid -----\nfunction buildPaths(segments) {\n  const pointSegs = new Map();\n  segments.forEach((seg) => {\n    seg.forEach((p) => {\n      if (!pointSegs.has(p)) pointSegs.set(p, []);\n      pointSegs.get(p).push(seg);\n    });\n  });\n  const visited = new Set();\n  const paths = [];\n  segments.forEach((seg) => {\n    if (visited.has(seg)) return;\n    visited.add(seg);\n    const path = [seg[0], seg[1]];\n    let extended = true;\n    while (extended) {\n      extended = false;\n      const last = path[path.length - 1];\n      const next = pointSegs.get(last).find((s) => !visited.has(s));\n      if (next) {\n        visited.add(next);\n        const nextPoint = next[0] === last ? next[1] : next[0];\n        path.push(nextPoint);\n        extended = nextPoint !== path[0];\n      }\n    }\n    extended = true;\n    while (extended) {\n      extended = false;\n      const first = path[0];\n      const next = pointSegs.get(first).find((s) => !visited.has(s));\n      if (next) {\n        visited.add(next);\n        path.unshift(next[0] === first ? next[1] : next[0]);\n        extended = true;\n      }\n    }\n    paths.push(path);\n  });\n  return paths;\n}\n\nfunction marchingSquaresPaths(values, gx, gy, level) {\n  const ny = values.length;\n  const nx = values[0].length;\n\n  // Shared edge crossings, computed once so neighboring cells reference the\n  // same point object (exact chaining without floating-point epsilon compares).\n  const hCross = [];\n  for (let j = 0; j < ny; j++) {\n    const row = [];\n    for (let i = 0; i < nx - 1; i++) {\n      const v0 = values[j][i];\n      const v1 = values[j][i + 1];\n      if ((v0 >= level) !== (v1 >= level)) {\n        const frac = (level - v0) / (v1 - v0);\n        row.push({ x: gx[i] + frac * (gx[i + 1] - gx[i]), y: gy[j] });\n      } else {\n        row.push(null);\n      }\n    }\n    hCross.push(row);\n  }\n  const vCross = [];\n  for (let i = 0; i < nx; i++) {\n    const col = [];\n    for (let j = 0; j < ny - 1; j++) {\n      const v0 = values[j][i];\n      const v1 = values[j + 1][i];\n      if ((v0 >= level) !== (v1 >= level)) {\n        const frac = (level - v0) / (v1 - v0);\n        col.push({ x: gx[i], y: gy[j] + frac * (gy[j + 1] - gy[j]) });\n      } else {\n        col.push(null);\n      }\n    }\n    vCross.push(col);\n  }\n\n  const segments = [];\n  for (let j = 0; j < ny - 1; j++) {\n    for (let i = 0; i < nx - 1; i++) {\n      const v00 = values[j][i];\n      const v10 = values[j][i + 1];\n      const v11 = values[j + 1][i + 1];\n      const v01 = values[j + 1][i];\n      const c00 = v00 >= level;\n      const c10 = v10 >= level;\n      const c11 = v11 >= level;\n      const c01 = v01 >= level;\n      const caseIndex = (c00 ? 1 : 0) | (c10 ? 2 : 0) | (c11 ? 4 : 0) | (c01 ? 8 : 0);\n      if (caseIndex === 0 || caseIndex === 15) continue;\n      const edgeB = hCross[j][i];\n      const edgeT = hCross[j + 1][i];\n      const edgeL = vCross[i][j];\n      const edgeR = vCross[i + 1][j];\n      const center = (v00 + v10 + v11 + v01) / 4;\n      let pairs;\n      if (caseIndex === 5) pairs = center >= level ? [[edgeL, edgeT], [edgeB, edgeR]] : [[edgeL, edgeB], [edgeT, edgeR]];\n      else if (caseIndex === 10) pairs = center >= level ? [[edgeL, edgeB], [edgeT, edgeR]] : [[edgeL, edgeT], [edgeB, edgeR]];\n      else if (caseIndex === 1 || caseIndex === 14) pairs = [[edgeL, edgeB]];\n      else if (caseIndex === 2 || caseIndex === 13) pairs = [[edgeB, edgeR]];\n      else if (caseIndex === 3 || caseIndex === 12) pairs = [[edgeL, edgeR]];\n      else if (caseIndex === 4 || caseIndex === 11) pairs = [[edgeR, edgeT]];\n      else if (caseIndex === 6 || caseIndex === 9) pairs = [[edgeB, edgeT]];\n      else pairs = [[edgeL, edgeT]]; // caseIndex 7 or 8\n      pairs.forEach((pair) => segments.push(pair));\n    }\n  }\n  return buildPaths(segments);\n}\n\nfunction pathsToSeriesData(paths) {\n  // A bare `null` array element would get an auto-incremented index x (0, 1, 2…)\n  // from Highcharts' tuple parser, corrupting the xAxis extremes. Give the gap\n  // an explicit, in-range x instead.\n  const data = [];\n  paths.forEach((path, idx) => {\n    if (idx > 0) data.push({ x: path[0].x, y: null });\n    path.forEach((p) => data.push([p.x, p.y]));\n  });\n  return data;\n}\n\n// --- Sequential Imprint colormap for the density levels ---------------------\nfunction hexToRgb(hex) {\n  const value = parseInt(hex.slice(1), 16);\n  return [(value >> 16) & 255, (value >> 8) & 255, value & 255];\n}\nfunction lerpColor(hex1, hex2, frac) {\n  const a = hexToRgb(hex1);\n  const b = hexToRgb(hex2);\n  const r = Math.round(a[0] + (b[0] - a[0]) * frac);\n  const g = Math.round(a[1] + (b[1] - a[1]) * frac);\n  const bl = Math.round(a[2] + (b[2] - a[2]) * frac);\n  return `rgb(${r}, ${g}, ${bl})`;\n}\n\nconst levelFractions = [0.12, 0.28, 0.46, 0.64, 0.82];\nconst contourSeries = levelFractions.map((frac, idx) => {\n  const paths = marchingSquaresPaths(densityGrid, gridX, gridY, frac * maxDensity);\n  return {\n    type: \"spline\",\n    name: `Density ${Math.round(frac * 100)}%`,\n    color: lerpColor(t.seq[0], t.seq[1], idx / (levelFractions.length - 1)),\n    data: pathsToSeriesData(paths),\n    lineWidth: 2.5,\n    marker: { enabled: false },\n    enableMouseTracking: false,\n    showInLegend: true,\n  };\n});\n\nconst scatterRgb = hexToRgb(INK_MUTED);\nconst scatterSeries = {\n  type: \"scatter\",\n  name: \"Measured parts\",\n  color: `rgba(${scatterRgb[0]}, ${scatterRgb[1]}, ${scatterRgb[2]}, 0.35)`,\n  data: parts.map((p) => [p.x, p.y]),\n  marker: { radius: 2.5, symbol: \"circle\", lineWidth: 0 },\n  showInLegend: false,\n  tooltip: { pointFormat: \"Diameter: {point.x:.2f} mm<br/>Weight: {point.y:.2f} g\" },\n};\n\n// --- Chart -------------------------------------------------------------------\nconst title = \"Process QC · contour-density · javascript · highcharts · anyplot.ai\";\nconst titleFontSize = Math.round(22 * Math.min(1, 67 / title.length));\n\nHighcharts.chart(\"container\", {\n  chart: {\n    type: \"spline\",\n    backgroundColor: \"transparent\",\n    animation: false,\n    style: { fontFamily: \"inherit\" },\n  },\n  credits: { enabled: false },\n  title: {\n    text: title,\n    style: { color: t.ink, fontSize: `${titleFontSize}px`, fontWeight: \"600\" },\n  },\n  xAxis: {\n    title: { text: \"Part Diameter (mm)\", style: { color: t.inkSoft, fontSize: \"16px\" } },\n    lineColor: t.inkSoft,\n    tickColor: t.inkSoft,\n    gridLineColor: t.grid,\n    labels: { style: { color: t.inkSoft, fontSize: \"14px\" } },\n  },\n  yAxis: {\n    title: { text: \"Part Weight (g)\", style: { color: t.inkSoft, fontSize: \"16px\" } },\n    gridLineColor: t.grid,\n    labels: { style: { color: t.inkSoft, fontSize: \"14px\" } },\n  },\n  legend: {\n    itemStyle: { color: t.inkSoft, fontSize: \"14px\" },\n    itemHoverStyle: { color: t.ink },\n  },\n  plotOptions: {\n    series: { animation: false },\n  },\n  series: [scatterSeries, ...contourSeries],\n});\n"}