#!/usr/bin/env node
// Checks the MATLAB layer against the numbers the chapters claim, through the
// numbl CLI, with no browser involved.
//
//   npm run test:matlab            (NUMBL_DIR defaults to ~/src/numbl)
//
// Every assertion here is a statement made in the prose, so if one of these
// fails either the code or the text is wrong.
import { join, dirname } from 'node:path'
import { fileURLToPath } from 'node:url'
import { runLab } from './runlab.mjs'

const root = join(dirname(fileURLToPath(import.meta.url)), '..')
const src = (p) => join(root, 'src', p)

let failures = 0
const check = (name, ok, detail = '') => {
  console.log(`${ok ? '  ok  ' : ' FAIL '} ${name}${detail ? `   ${detail}` : ''}`)
  if (!ok) failures++
}
const near = (a, b, tol) => Math.abs(a - b) <= tol
const last = (row) => {
  for (let i = row.length - 1; i >= 0; i--) if (row[i] != null) return row[i]
  return null
}
const best = (row) => Math.min(...row.filter((v) => v != null))

const ONE = src('chapters/onedim/driver.m')
const TWO = src('chapters/twodim/driver.m')

console.log('\nChapter 1: three-point differences')
{
  const { out } = runLab(ONE, src('chapters/ch01/fd2.m'), {
    kind: 'smooth',
    mode: 'sweep',
    ns: [16, 32, 64, 128, 256, 512],
    nmaxs: [1e9],
    timed: 0,
  })
  const q = last(out.orders[0])
  check('observed order is 2', near(q, 2, 0.02), `q = ${q.toFixed(3)}`)
  const { out: low } = runLab(ONE, src('chapters/ch01/fd2.m'), {
    kind: 'lowreg',
    mode: 'sweep',
    ns: [64, 128, 256, 512],
    nmaxs: [1e9],
    timed: 0,
  })
  const ql = last(low.orders[0])
  check('the low-regularity solution converges more slowly', ql < 1.7, `q = ${ql.toFixed(3)}`)
}

console.log('\nChapter 2: the roundoff floor and Richardson')
{
  const { out } = runLab(ONE, src('chapters/ch01/fd2.m'), {
    kind: 'smooth',
    mode: 'sweep',
    ns: [1024, 4096, 16384, 65536, 262144],
    nmaxs: [1e9],
    timed: 0,
  })
  const e = out.E[0]
  const bi = e.indexOf(best(e))
  check('the error stops improving before N = 262144', bi < e.length - 1, `best at N = ${out.ns[bi]}`)
  check('the best error is a few times 1e-9', best(e) < 3e-8 && best(e) > 1e-10, `${best(e).toExponential(2)}`)

  const { out: r } = runLab(ONE, src('chapters/ch02/richardson.m'), {
    kind: 'smooth',
    mode: 'sweep',
    ns: [32, 64, 128, 256, 512],
    nmaxs: [1e9, 1e9],
    timed: 0,
  })
  check('Richardson gives order 4', near(last(r.orders[1]), 4, 0.1), `q = ${last(r.orders[1]).toFixed(3)}`)
  check('and it beats the plain scheme by >1e4', last(r.E[0]) / last(r.E[1]) > 1e4)

  const { out: rl } = runLab(ONE, src('chapters/ch02/richardson.m'), {
    kind: 'lowreg',
    mode: 'sweep',
    ns: [64, 128, 256, 512],
    nmaxs: [1e9, 1e9],
    timed: 0,
  })
  check(
    'on lowreg the extrapolation buys no extra order',
    near(last(rl.orders[1]), last(rl.orders[0]), 0.1) && last(rl.orders[1]) < 1.7,
    `q = ${last(rl.orders[1]).toFixed(2)} against ${last(rl.orders[0]).toFixed(2)}`,
  )
}

console.log('\nChapter 3: the compact fourth-order scheme')
{
  const { out } = runLab(ONE, src('chapters/ch03/highorder.m'), {
    kind: 'smooth',
    mode: 'sweep',
    ns: [64, 256, 1024],
    nmaxs: [1e9, 1e9],
    timed: 0,
  })
  check('Numerov is fourth order', near(last(out.orders[1]), 4, 0.05), `q = ${last(out.orders[1]).toFixed(3)}`)
  check('at N = 1024 it is below 1e-9', last(out.E[1]) < 1e-9, last(out.E[1]).toExponential(2))
}

console.log('\nChapter 4: Chebyshev collocation')
{
  const { out } = runLab(ONE, src('chapters/ch04/spectral.m'), {
    kind: 'smooth',
    mode: 'sweep',
    ns: [16, 32, 64, 128],
    nmaxs: [1e9, 1e9, 1e9],
    timed: 0,
  })
  check('spectral reaches 1e-13 by N = 64', out.E[2][2] < 1e-12, out.E[2][2].toExponential(2))
  const { out: low } = runLab(ONE, src('chapters/ch04/spectral.m'), {
    kind: 'lowreg',
    mode: 'sweep',
    ns: [64, 128, 256],
    nmaxs: [1e9, 1e9, 1e9],
    timed: 0,
  })
  check('on lowreg it is no better than the three-point scheme', low.E[2][2] > low.E[0][2])
}

console.log('\nChapters 5 and 6: the five-point Laplacian, two solvers')
{
  const { out } = runLab(TWO, src('chapters/ch06/fastpoisson.m'), {
    kind: 'expsin',
    par: 1,
    mode: 'sweep',
    ns: [8, 16, 32, 128, 512],
    nmaxs: [32, 512],
    timed: 0,
  })
  check('the five-point scheme is second order', near(last(out.orders[1]), 2, 0.02), `q = ${last(out.orders[1]).toFixed(3)}`)
  for (let j = 0; j < 3; j++) {
    check(
      `dense and transform agree at N = ${out.ns[j]}`,
      Math.abs(out.E[0][j] - out.E[1][j]) < 1e-9 * out.E[0][j],
      `${out.E[0][j].toExponential(4)} vs ${out.E[1][j].toExponential(4)}`,
    )
  }
  check('at N = 512 the error is 4e-8', near(last(out.E[1]) / 3.99e-8, 1, 0.05), last(out.E[1]).toExponential(3))
}

console.log('\nChapter 7: the compact nine-point stencil and 2-D Chebyshev')
{
  const { out } = runLab(TWO, src('chapters/ch07/highorder2d.m'), {
    kind: 'expsin',
    par: 1,
    mode: 'sweep',
    ns: [8, 16, 32, 64],
    nmaxs: [512, 512, 64],
    timed: 0,
  })
  check('nine-point is sixth order on Laplace', near(out.orders[1][1], 6, 0.1), `q = ${out.orders[1][1].toFixed(2)}`)
  check('nine-point reaches 1e-13 by N = 32', out.E[1][2] < 1e-13, out.E[1][2].toExponential(2))
  check('Chebyshev reaches 1e-13 by N = 16', out.E[2][1] < 1e-13, out.E[2][1].toExponential(2))
}

console.log('\nChapter 8: iterative solvers')
{
  const { out } = runLab(src('chapters/ch08/driver.m'), src('chapters/ch08/iterative.m'), {
    kind: 'expsin',
    par: 1,
    N: 64,
    nit: 60,
    tol: 1e-11,
  })
  const mg = out.R[4].filter((v) => v != null)
  check('multigrid converges in fewer than 15 cycles', mg.length < 16, `${mg.length - 1} cycles`)
  const mgErr = out.E[4].filter((v) => v != null)
  check(
    'and it lands on the discretisation error',
    near(mgErr[mgErr.length - 1] / out.discErr, 1, 0.02),
    `${mgErr[mgErr.length - 1].toExponential(3)} vs ${out.discErr.toExponential(3)}`,
  )
  const jac = out.R[0].filter((v) => v != null)
  check('Jacobi is nowhere near after 60 sweeps', jac[jac.length - 1] > 1e-2, jac[jac.length - 1].toExponential(2))
}

console.log('\nChapter 9: the L-shaped domain')
{
  const NSM = [
    [16, 32, 64],
    [8, 16, 32],
  ]
  const { out: b } = runLab(src('chapters/ch09/driver.m'), src('chapters/ch09/lshape.m'), {
    kind: 'both',
    mode: 'sweep',
    nsm: NSM,
  })
  check('the difference scheme loses second order', near(last(b.orders[0]), 1.37, 0.1), `q = ${last(b.orders[0]).toFixed(3)}`)
  check('the harmonic basis reaches 1e-11', last(b.E[1]) < 1e-11, last(b.E[1]).toExponential(2))

  const { out: s } = runLab(src('chapters/ch09/driver.m'), src('chapters/ch09/lshape.m'), {
    kind: 'smooth',
    mode: 'sweep',
    nsm: NSM,
  })
  check('with smooth data the corner costs nothing', near(last(s.orders[0]), 2, 0.02), `q = ${last(s.orders[0]).toFixed(3)}`)

  const { out: g } = runLab(src('chapters/ch09/driver.m'), src('chapters/ch09/lshape.m'), {
    kind: 'sing',
    mode: 'sweep',
    nsm: NSM,
  })
  check('the singular function is in the basis exactly', g.E[1][0] < 1e-14, g.E[1][0].toExponential(2))
}

console.log('\nChapter 10: boundary integral equations')
{
  const { out } = runLab(src('chapters/ch10/driver.m'), src('chapters/ch10/bie.m'), {
    a: 0.3,
    nlobe: 5,
    par: 2.2,
    mode: 'sweep',
    ns: [64, 128, 256, 512],
  })
  check('Nyström reaches machine precision by n = 512', last(out.E[0]) < 1e-14, last(out.E[0]).toExponential(2))
  check('convergence is geometric, not algebraic', out.E[0][0] / out.E[0][1] > 100)
  check(
    'close evaluation is orders of magnitude worse',
    last(out.Enear[0]) / last(out.E[0]) > 1e6,
    `${last(out.Enear[0]).toExponential(2)} at 98% radius`,
  )
  check('fundamental solutions reach 1e-13', best(out.E[1]) < 1e-13, best(out.E[1]).toExponential(2))
  check(
    'and they have no close-evaluation problem',
    best(out.Enear[1]) < 1e-12,
    best(out.Enear[1]).toExponential(2),
  )
}

console.log('\nChapter 11: the Pac-Man benchmark')
{
  // zeros in a ladder skip that resolution (or the whole method)
  const { out } = runLab(src('chapters/ch11/driver.m'), src('chapters/ch11/pacman.m'), {
    kind: 'both',
    mode: 'sweep',
    nsm: [
      [16, 32, 64, 0],
      [4, 8, 16, 32],
      [4, 8, 0, 0],
      [8, 16, 32, 0],
    ],
  })
  const qfd = last(out.orders[0])
  check('Shortley-Weller is corner-limited near order 4/3', near(qfd, 1.37, 0.12), `q = ${qfd.toFixed(3)}`)
  check('the harmonic basis reaches 1e-13 with 97 unknowns', out.E[1][3] < 1e-13, out.E[1][3].toExponential(2))
  check('and converges geometrically', out.E[1][1] / out.E[1][2] > 100)
  const ratio = out.E[2][0] / out.E[2][1]
  check(
    'graded panels gain about 2^(8/3) per four levels of depth',
    ratio > 4 && ratio < 9,
    `ratio = ${ratio.toFixed(2)}`,
  )
  check('and reach 3e-5 by depth 8', out.E[2][1] < 3e-5, out.E[2][1].toExponential(2))
  check('polar Chebyshev is algebraic on the singular data', out.E[3][2] > 1e-5 && out.E[3][2] < 3e-4, out.E[3][2].toExponential(2))

  const { out: sm } = runLab(src('chapters/ch11/driver.m'), src('chapters/ch11/pacman.m'), {
    kind: 'smooth',
    mode: 'sweep',
    nsm: [
      [16, 32, 64, 0],
      [0, 0, 0, 0],
      [0, 0, 0, 0],
      [6, 16, 40, 0],
    ],
  })
  const qsw = last(sm.orders[0])
  check('with smooth data Shortley-Weller is second order', near(qsw, 2, 0.15), `q = ${qsw.toFixed(3)}`)
  check('and polar Chebyshev is geometric again, 1e-9 by N = 40', sm.E[3][2] < 1e-9, sm.E[3][2].toExponential(2))

  const { out: sg } = runLab(src('chapters/ch11/driver.m'), src('chapters/ch11/pacman.m'), {
    kind: 'sing',
    mode: 'sweep',
    nsm: [
      [0, 0, 0, 0],
      [4, 0, 0, 0],
      [0, 0, 0, 0],
      [0, 0, 0, 0],
    ],
  })
  check('the singular function is in the basis exactly', sg.E[1][0] < 1e-13, sg.E[1][0].toExponential(2))
}

console.log(failures === 0 ? '\nall checks passed' : `\n${failures} check(s) failed`)
process.exit(failures === 0 ? 0 : 1)