// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

package core

import (
	"math"
	"testing"
)

// TestInterpolateMonotoneShape checks the two promises of the
// Fritsch-Carlson tangents: the curve passes through every knot, and
// interpolating monotone data never overshoots the data's range.
func TestInterpolateMonotoneShape(t *testing.T) {
	xv := []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}
	yv := []float64{0, 0.1, 1.9, 2.0, 2.1, 2.2, 3.9, 4.0, 4.1, 8}
	xs := mustFloats(t, xv, len(xv))
	ys := mustFloats(t, yv, len(yv))
	dense := make([]float64, 0, 801)
	for i := range 801 {
		dense = append(dense, -1+float64(i)*10/800)
	}
	q := mustFloats(t, dense, len(dense))
	out, err := InterpolateMonotone(xs, ys, q)
	if err != nil {
		t.Fatalf("InterpolateMonotone: %v", err)
	}
	for k := range xv {
		if got := out.FloatAt((k + 1) * 80); math.Abs(got-yv[k]) > 1e-12 {
			t.Fatalf("knot %d: %.12g, want %.12g", k, got, yv[k])
		}
	}
	for i := range out.Len() {
		if v := out.FloatAt(i); v < -1e-12 || v > 8+1e-12 {
			t.Fatalf("overshoot at %g: %.12g outside the data range", dense[i], v)
		}
	}
}

// TestInterpolateMonotoneBasics checks the two-knot limit (a
// straight line), the flat clamp outside, and the refusal of
// repeated knots.
func TestInterpolateMonotoneBasics(t *testing.T) {
	xs := mustFloats(t, []float64{0, 1}, 2)
	ys := mustFloats(t, []float64{10, 20}, 2)
	q := mustFloats(t, []float64{0.25, -3, 4}, 3)
	out, err := InterpolateMonotone(xs, ys, q)
	if err != nil {
		t.Fatalf("InterpolateMonotone: %v", err)
	}
	if got := out.FloatAt(0); math.Abs(got-12.5) > 1e-12 {
		t.Fatalf("linear read = %.12g, want 12.5", got)
	}
	if got := out.FloatAt(1); got != 10 {
		t.Fatalf("clamp below = %.12g, want 10", got)
	}
	if got := out.FloatAt(2); got != 20 {
		t.Fatalf("clamp above = %.12g, want 20", got)
	}
	bad := mustFloats(t, []float64{1, 1, 2}, 3)
	if _, err := InterpolateMonotone(bad, ys, q); err == nil {
		t.Fatal("repeated knots accepted")
	}
}

// TestInterpolateGrid checks multilinear reads: exact at the nodes,
// the bilinear average at a cell centre, clamped outside, and the
// argument guards.
func TestInterpolateGrid(t *testing.T) {
	grid := mustFloats(t, []float64{
		0, 10,
		20, 30,
	}, 2, 2)
	origins := []float64{0, 0}
	steps := []float64{1, 0.5}
	// Node reads land exactly; the cell centre averages the corners.
	queries := mustFloats(t, []float64{
		0, 0,
		1, 0.5,
		0.5, 0.25,
		-9, 9,
	}, 4, 2)
	out, err := InterpolateGrid(grid, origins, steps, queries)
	if err != nil {
		t.Fatalf("InterpolateGrid: %v", err)
	}
	if got := out.FloatAt(0); got != 0 {
		t.Fatalf("node (0,0) = %.12g, want 0", got)
	}
	if got := out.FloatAt(1); got != 30 {
		t.Fatalf("node (1,0.5) = %.12g, want 30", got)
	}
	if got := out.FloatAt(2); math.Abs(got-15) > 1e-12 {
		t.Fatalf("cell centre = %.12g, want 15", got)
	}
	if got := out.FloatAt(3); got != 10 {
		t.Fatalf("clamp = %.12g, want 10", got)
	}
	if _, err := InterpolateGrid(grid, origins, []float64{1, -1}, queries); err == nil {
		t.Fatal("negative step accepted")
	}
	badQueries := mustFloats(t, []float64{0, 0, 0}, 3, 1)
	if _, err := InterpolateGrid(grid, origins, steps, badQueries); err == nil {
		t.Fatal("wrong query columns accepted")
	}
}

// TestDiff checks the finite differences on a known series and along
// a chosen axis of a matrix, plus the argument guards.
func TestDiff(t *testing.T) {
	a := mustFloats(t, []float64{1, 4, 9, 16}, 4)
	d1, err := Diff(a, 1, 0)
	if err != nil {
		t.Fatalf("Diff: %v", err)
	}
	want := []float64{3, 5, 7}
	for i := range want {
		if d1.FloatAt(i) != want[i] {
			t.Fatalf("first difference %d = %g, want %g", i, d1.FloatAt(i), want[i])
		}
	}
	d2, err := Diff(a, 2, 0)
	if err != nil {
		t.Fatalf("Diff order 2: %v", err)
	}
	if d2.FloatAt(0) != 2 || d2.FloatAt(1) != 2 {
		t.Fatalf("second difference = [%g, %g], want [2, 2]", d2.FloatAt(0), d2.FloatAt(1))
	}
	m := mustFloats(t, []float64{1, 4, 2, 8}, 2, 2)
	dm, err := Diff(m, 1, 0)
	if err != nil {
		t.Fatalf("Diff axis 0: %v", err)
	}
	if dm.Shape()[0] != 1 || dm.Shape()[1] != 2 {
		t.Fatalf("shape %v, want [1 2]", dm.Shape())
	}
	if dm.FloatAt(1) != 4 {
		t.Fatalf("axis-0 difference column 1 = %g, want 4", dm.FloatAt(1))
	}
	if _, err := Diff(a, 4, 0); err == nil {
		t.Fatal("order at the axis length accepted")
	}
	if _, err := Diff(m, 1, 5); err == nil {
		t.Fatal("axis past the rank accepted")
	}
}

// TestDiffNonLeadingAxis pins the run walk of Diff along an axis that is
// neither leading nor trailing: several outer positions (head) each hold
// a run of trailing elements (tail), so a run's destination is not its
// own start and the source offsets step by the line length.
func TestDiffNonLeadingAxis(t *testing.T) {
	vals := make([]int64, 24)
	for i := range vals {
		vals[i] = int64(i)
	}
	src := mustFromInts(t, vals, 2, 3, 4)

	// Axis 1: head 2, line 3, tail 4.
	got, err := Diff(src, 1, 1)
	if err != nil {
		t.Fatalf("Diff axis 1: %v", err)
	}
	if sh := got.Shape(); len(sh) != 3 || sh[0] != 2 || sh[1] != 2 || sh[2] != 4 {
		t.Fatalf("Diff axis 1 shape: %v, want [2 2 4]", sh)
	}
	for i := range 2 {
		for j := range 2 {
			for k := range 4 {
				hi, _ := IntAt(src, i, j+1, k)
				lo, _ := IntAt(src, i, j, k)
				if v, _ := IntAt(got, i, j, k); v != hi-lo {
					t.Fatalf("Diff axis 1 [%d %d %d] = %d, want %d", i, j, k, v, hi-lo)
				}
			}
		}
	}

	// Axis 2: head 6, line 4, tail 1.
	got2, err := Diff(src, 1, 2)
	if err != nil {
		t.Fatalf("Diff axis 2: %v", err)
	}
	if sh := got2.Shape(); len(sh) != 3 || sh[0] != 2 || sh[1] != 3 || sh[2] != 3 {
		t.Fatalf("Diff axis 2 shape: %v, want [2 3 3]", sh)
	}
	for i := range 2 {
		for j := range 3 {
			for k := range 3 {
				hi, _ := IntAt(src, i, j, k+1)
				lo, _ := IntAt(src, i, j, k)
				if v, _ := IntAt(got2, i, j, k); v != hi-lo {
					t.Fatalf("Diff axis 2 [%d %d %d] = %d, want %d", i, j, k, v, hi-lo)
				}
			}
		}
	}
}