package stats

import (
	"math"
	"math/rand/v2"
	"sort"
	"testing"

	"salty/terrain/internal/field"
)

// A small world with real structure in it: a coast, a range, a plain, and a sea the statistics have to leave
// out. Deterministic, so both halves of every comparison see the same ground.
func testWorld(t *testing.T, w, h int, cellM float64) (*field.Field, []bool, []float32) {
	t.Helper()
	r := rand.New(rand.NewPCG(11, 13))
	f := field.New(w, h, cellM)
	land := make([]bool, w*h)
	up := make([]float32, w*h)
	for y := 0; y < h; y++ {
		for x := 0; x < w; x++ {
			i := y*w + x
			if x < w/8 {
				f.Data[i] = -40 // the sea, which must not appear in any land statistic
				continue
			}
			land[i] = true
			t := float64(x) / float64(w)
			// A range towards the east, a plain in the middle, and enough noise to give the slopes a spread.
			// The 40 m base keeps every land cell above sea level, so "the land minimum is positive" is a
			// statement about the mask rather than about this formula.
			f.Data[i] = float32(40 + 300*t*t + 18*math.Sin(float64(x)/9)*math.Cos(float64(y)/7) +
				r.NormFloat64()*3)
			up[i] = float32((0.02 + 0.9*t*t*t) / 1000)
		}
	}
	return f, land, up
}

func testOptions() Options {
	return Options{ElevMin: -1024, ElevMax: 2048, TalusDeg: 35, ReliefWindowM: 500,
		ChannelM2: 1e6, K: 5e-5, M: 0.5, N: 1}
}

// The histograms replaced sorts, and the whole point is that nothing a run is judged by moved. This is the
// same data through both: the old implementation is reproduced here as the reference, so that a future change
// to the fast path has something to be wrong against.
func TestTheHistogramsAgreeWithSorting(t *testing.T) {
	const w, h, cellM = 220, 160, 8.0
	f, land, up := testWorld(t, w, h, cellM)

	acc := New(testOptions())
	acc.Add(Input{H: f, Land: land, UpliftMYr: up})
	acc.AddExtent(f.Data, land, 0)
	got := acc.Report(cellM)

	// --- the reference, by sorting, exactly as the package used to do it -----------------------------
	slope := f.Slope()
	var degs, elevs []float64
	for i := range f.Data {
		if !land[i] {
			continue
		}
		degs = append(degs, math.Atan(float64(slope.Data[i]))*180/math.Pi)
		elevs = append(elevs, float64(f.Data[i]))
	}
	sort.Float64s(degs)
	sort.Float64s(elevs)
	frac := func(v []float64, limit float64) float64 {
		return float64(sort.SearchFloat64s(v, limit)) / float64(len(v))
	}

	const slopeTol = 90.0 / slopeBins // one bin: the whole error budget of a histogram quantile
	if d := math.Abs(got.Slopes.MedianDeg - degs[len(degs)/2]); d > slopeTol {
		t.Errorf("median slope %.4f against %.4f", got.Slopes.MedianDeg, degs[len(degs)/2])
	}
	for _, c := range []struct {
		name string
		got  float64
		want float64
	}{
		{"under 15", got.Slopes.Under15Deg, frac(degs, 15)},
		{"under 30", got.Slopes.Under30Deg, frac(degs, 30)},
		{"over 50", got.Slopes.Over50Deg, 1 - frac(degs, 50)},
	} {
		if math.Abs(c.got-c.want) > 0.002 {
			t.Errorf("slopes %s: %.4f against %.4f", c.name, c.got, c.want)
		}
	}

	// The hypsometric integral is a mean and is carried exactly, so it has to match to the bit of a float sum.
	lo, hi := elevs[0], elevs[len(elevs)-1]
	var sum float64
	for _, v := range elevs {
		sum += (v - lo) / (hi - lo)
	}
	if d := math.Abs(got.Hypsometry.Integral - sum/float64(len(elevs))); d > 1e-9 {
		t.Errorf("hypsometric integral %.6f against %.6f", got.Hypsometry.Integral, sum/float64(len(elevs)))
	}
	if got.LandMinM != lo || got.LandMaxM != hi {
		t.Errorf("land range %.3f..%.3f against %.3f..%.3f", got.LandMinM, got.LandMaxM, lo, hi)
	}

	// And the per-class breakdown, which is the block that matters most.
	for _, b := range got.Buckets {
		var bdeg []float64
		for i := range f.Data {
			if !land[i] {
				continue
			}
			mm := float64(up[i]) * 1000
			if mm < b.LoMmYr || mm >= b.HiMmYr {
				continue
			}
			bdeg = append(bdeg, math.Atan(float64(slope.Data[i]))*180/math.Pi)
		}
		if len(bdeg) != b.Cells {
			t.Errorf("bucket %s holds %d cells, the reference found %d", b.Name, b.Cells, len(bdeg))
		}
		sort.Float64s(bdeg)
		if d := math.Abs(b.MedianDeg - bdeg[len(bdeg)/2]); d > slopeTol {
			t.Errorf("bucket %s median %.4f against %.4f", b.Name, b.MedianDeg, bdeg[len(bdeg)/2])
		}
		p90 := bdeg[min(len(bdeg)*9/10, len(bdeg)-1)]
		if d := math.Abs(b.P90Deg - p90); d > slopeTol {
			t.Errorf("bucket %s P90 %.4f against %.4f", b.Name, b.P90Deg, p90)
		}
	}
	if len(got.Buckets) < 2 {
		t.Fatalf("only %d buckets came out; this test measured almost nothing", len(got.Buckets))
	}
}

// The property the planet depends on: a world cut into pieces and accumulated piece by piece has to report
// what one pass over the whole thing would. Everything here is additive by construction, and this is the
// assertion that says so end to end rather than one histogram at a time.
func TestPoolingPiecesMatchesOnePass(t *testing.T) {
	const w, h, cellM = 240, 120, 8.0
	f, land, up := testWorld(t, w, h, cellM)

	whole := New(testOptions())
	whole.Add(Input{H: f, Land: land, UpliftMYr: up})
	whole.AddExtent(f.Data, land, 0)

	// The same ground in three horizontal strips. Slope and relief read neighbours, so a strip's own edge
	// rows differ from the whole - which is exactly the seam a region has, and the reason the comparison
	// below is on the *distributions* rather than cell by cell.
	pooled := New(testOptions())
	for _, band := range [][2]int{{0, 40}, {40, 80}, {80, 120}} {
		sub := field.New(w, band[1]-band[0], cellM)
		subLand := make([]bool, w*(band[1]-band[0]))
		subUp := make([]float32, len(subLand))
		copy(sub.Data, f.Data[band[0]*w:band[1]*w])
		copy(subLand, land[band[0]*w:band[1]*w])
		copy(subUp, up[band[0]*w:band[1]*w])
		pooled.Add(Input{H: sub, Land: subLand, UpliftMYr: subUp})
		pooled.AddExtent(sub.Data, subLand, 0)
	}

	a, b := whole.Report(cellM), pooled.Report(cellM)
	if a.LandFraction != b.LandFraction {
		t.Errorf("land fraction %.6f pooled against %.6f whole", b.LandFraction, a.LandFraction)
	}
	if a.LandMinM != b.LandMinM || a.LandMaxM != b.LandMaxM {
		t.Errorf("land range %.3f..%.3f pooled against %.3f..%.3f",
			b.LandMinM, b.LandMaxM, a.LandMinM, a.LandMaxM)
	}
	// Elevation does not read neighbours at all, so it has to pool to the bit.
	if math.Abs(a.Hypsometry.Integral-b.Hypsometry.Integral) > 1e-12 {
		t.Errorf("hypsometric integral %.9f pooled against %.9f", b.Hypsometry.Integral, a.Hypsometry.Integral)
	}
	// Slope reads one cell either side, so six rows of a 120-row world are clamped differently. The
	// distribution has to survive that; a tenth of a degree is far inside anything Summary turns on.
	if d := math.Abs(a.Slopes.MedianDeg - b.Slopes.MedianDeg); d > 0.1 {
		t.Errorf("median slope %.3f pooled against %.3f", b.Slopes.MedianDeg, a.Slopes.MedianDeg)
	}
	for i := range a.Buckets {
		if i >= len(b.Buckets) {
			t.Fatalf("pooling lost a bucket: %d against %d", len(b.Buckets), len(a.Buckets))
		}
		if a.Buckets[i].Cells != b.Buckets[i].Cells {
			t.Errorf("bucket %s: %d cells pooled against %d", a.Buckets[i].Name,
				b.Buckets[i].Cells, a.Buckets[i].Cells)
		}
		if d := math.Abs(a.Buckets[i].MedianDeg - b.Buckets[i].MedianDeg); d > 0.2 {
			t.Errorf("bucket %s median %.3f pooled against %.3f", a.Buckets[i].Name,
				b.Buckets[i].MedianDeg, a.Buckets[i].MedianDeg)
		}
	}
}

// Merge is the other way pieces arrive - a planet's regions are accumulated separately and folded together -
// and it has to be the same as adding them to one accumulator.
func TestMergeMatchesAddingToOne(t *testing.T) {
	const w, h, cellM = 160, 60, 8.0
	f, land, up := testWorld(t, w, h, cellM)

	one := New(testOptions())
	one.Add(Input{H: f, Land: land, UpliftMYr: up})
	one.AddExtent(f.Data, land, 0)
	one.Add(Input{H: f, Land: land, UpliftMYr: up})
	one.AddExtent(f.Data, land, 0)

	a, b := New(testOptions()), New(testOptions())
	a.Add(Input{H: f, Land: land, UpliftMYr: up})
	a.AddExtent(f.Data, land, 0)
	b.Add(Input{H: f, Land: land, UpliftMYr: up})
	b.AddExtent(f.Data, land, 0)
	a.Merge(b)

	x, y := one.Report(cellM), a.Report(cellM)
	if x.LandFraction != y.LandFraction || x.Slopes.MedianDeg != y.Slopes.MedianDeg ||
		x.LandMinM != y.LandMinM || x.LandMaxM != y.LandMaxM {
		t.Errorf("merged report differs from one built by adding twice:\n %+v\n %+v", x.Slopes, y.Slopes)
	}
	for i := range x.Buckets {
		if x.Buckets[i].Cells != y.Buckets[i].Cells || x.Buckets[i].MedianDeg != y.Buckets[i].MedianDeg {
			t.Errorf("bucket %s differs after a merge", x.Buckets[i].Name)
		}
	}
}

// Sea cells are counted for the land fraction and for the encoding range, and are in nothing else. A single
// -40 m sea floor in a land statistic would flatter every relief number by forty metres for free.
func TestTheSeaIsNotLand(t *testing.T) {
	const w, h, cellM = 120, 80, 8.0
	f, land, up := testWorld(t, w, h, cellM)
	acc := New(testOptions())
	acc.Add(Input{H: f, Land: land, UpliftMYr: up})
	acc.AddExtent(f.Data, land, 0)
	r := acc.Report(cellM)

	if r.LandMinM < 0 {
		t.Errorf("land minimum is %.1f m; the sea got into the land statistics", r.LandMinM)
	}
	if r.MinM > -39 {
		t.Errorf("whole-field minimum is %.1f m; the sea should still bound the encoding range", r.MinM)
	}
	wantLand := 0
	for _, v := range land {
		if v {
			wantLand++
		}
	}
	if got := int(r.LandFraction*float64(w*h) + 0.5); got != wantLand {
		t.Errorf("land fraction says %d cells, the mask has %d", got, wantLand)
	}
}