package fluvial

import (
	"math"
	"testing"
)

// TestClampToReposeCutsACone is the smallest possible check on the constraint: a cone far steeper than the
// repose angle must come back at or under it.
func TestClampToReposeCutsACone(t *testing.T) {
	const (
		w, h  = 81, 81
		cellM = 10.0
		talus = 0.4 // about 22 degrees
	)
	base := make([]bool, w*h)
	field := make([]float32, w*h)
	for y := 0; y < h; y++ {
		for x := 0; x < w; x++ {
			d := math.Hypot(float64(x-w/2), float64(y-h/2)) * cellM
			field[y*w+x] = float32(math.Max(0, 800-2.0*d)) // slope 2.0, five times repose
		}
	}
	g := NewGrid(w, h, cellM, base)
	g.SetElevationRange(-100, 2000)
	removed := g.ClampToRepose(field, talus)

	worst := 0.0
	for y := 0; y < h; y++ {
		for x := 0; x < w; x++ {
			for k := 0; k < 8; k++ {
				nx, ny := x+dx8[k], y+dy8[k]
				if nx < 0 || ny < 0 || nx >= w || ny >= h {
					continue
				}
				d := cellM
				if dx8[k] != 0 && dy8[k] != 0 {
					d = cellM * math.Sqrt2
				}
				if s := float64(field[y*w+x]-field[ny*w+nx]) / d; s > worst {
					worst = s
				}
			}
		}
	}
	t.Logf("removed %.2f m mean; steepest slope now %.3f (repose %.3f)", removed, worst, talus)
	if worst > talus*1.02 {
		t.Errorf("steepest slope %.3f exceeds repose %.3f", worst, talus)
	}
}

// TestClampToReposeIsIsotropic asks what shape is left, which the test above cannot: a four-sided pyramid
// satisfies "no slope exceeds repose" exactly, so the constraint check says nothing about whether the clamp
// cut a cone or cut a pyramid.
//
// Clamp a cone far above repose, then for each of 360 azimuths find by bisection the radius at which the
// surface falls through a fixed height. A cone gives a constant radius; the amplitudes of the four-fold and
// eight-fold Fourier components of that radius, as a fraction of its mean, say how far from one it is.
//
// What the numbers turn out to be, and what they are not. Measured 0.97 % four-fold and 2.39 % eight-fold -
// and *identical* with the pop-order jitter, with the allowance jitter, with both and with neither. On a cone
// no two cells share a bucket, because the surface falls twenty metres a cell against a one-centimetre
// bucket, so the ordering bias this file's jitter removes has nothing to bite on here. The residual is
// geometry: a path to a point at 22.5 degrees has to be built of cardinal and diagonal steps, and the octile
// distance it accumulates exceeds the straight line by up to 8 %, so an eight-connected clamp cuts an
// octagon out of a cone whatever order it works in. That is irreducible without a wider neighbourhood, and
// the thresholds below sit above it: this test guards against a regression to something far worse, and
// TestBucketPQDoesNotPreferRasterOrder is what actually holds the ordering honest.
func TestClampToReposeIsIsotropic(t *testing.T) {
	const (
		w, h  = 201, 201
		cellM = 10.0
		talus = 0.4
		level = 300.0
	)
	field := make([]float32, w*h)
	for y := 0; y < h; y++ {
		for x := 0; x < w; x++ {
			d := math.Hypot(float64(x-w/2), float64(y-h/2)) * cellM
			field[y*w+x] = float32(math.Max(0, 2000-2.0*d))
		}
	}
	g := NewGrid(w, h, cellM, make([]bool, w*h))
	g.SetSeed(37125)
	g.SetElevationRange(-100, 4000)
	g.ClampToRepose(field, talus)

	at := func(fx, fy float64) float64 { // bilinear, in cells
		x0, y0 := int(fx), int(fy)
		if x0 < 0 || y0 < 0 || x0 >= w-1 || y0 >= h-1 {
			return 0
		}
		tx, ty := fx-float64(x0), fy-float64(y0)
		return (1-ty)*((1-tx)*float64(field[y0*w+x0])+tx*float64(field[y0*w+x0+1])) +
			ty*((1-tx)*float64(field[(y0+1)*w+x0])+tx*float64(field[(y0+1)*w+x0+1]))
	}

	const rays = 360
	var sum, c4r, c4i, c8r, c8i float64
	for i := 0; i < rays; i++ {
		th := 2 * math.Pi * float64(i) / rays
		cs, sn := math.Cos(th), math.Sin(th)
		lo, hi := 0.0, float64(w/2-2)
		for n := 0; n < 40; n++ { // bisect on the radius where the surface crosses `level`
			mid := (lo + hi) / 2
			if at(float64(w/2)+mid*cs, float64(h/2)+mid*sn) > level {
				lo = mid
			} else {
				hi = mid
			}
		}
		r := (lo + hi) / 2
		sum += r
		c4r += r * math.Cos(4*th)
		c4i += r * math.Sin(4*th)
		c8r += r * math.Cos(8*th)
		c8i += r * math.Sin(8*th)
	}
	a4 := 2 * math.Hypot(c4r, c4i) / sum
	a8 := 2 * math.Hypot(c8r, c8i) / sum
	t.Logf("clamped cone: mean radius %.2f cells, four-fold %.2f%%, eight-fold %.2f%%",
		sum/rays, a4*100, a8*100)
	if a4 > 0.02 || a8 > 0.04 {
		t.Errorf("the clamped cone is %.2f%% four-fold and %.2f%% eight-fold against 0.97 and 2.39 measured: "+
			"it is a pyramid, not an octagon", a4*100, a8*100)
	}
}

// TestBucketPQDoesNotPreferRasterOrder is the unit underneath it. Pushed plain, cells at one elevation come
// back in exactly reverse insertion order, which is a Spearman correlation of -1.
func TestBucketPQDoesNotPreferRasterOrder(t *testing.T) {
	const n = 4096
	order := func(jitter bool) float64 {
		q := newBucketPQ(0, 100)
		for i := 0; i < n; i++ {
			if jitter {
				q.pushJittered(50, int32(i), (hashXY(1, int32(i%64), int32(i/64), jitterReposeOrder)-0.5)*2*reposeOrderBuckets)
			} else {
				q.push(50, int32(i))
			}
		}
		var sum float64
		for pos := 0; pos < n; pos++ {
			idx := float64(q.pop())
			sum += (float64(pos) - float64(n-1)/2) * (idx - float64(n-1)/2)
		}
		var varr float64
		for i := 0; i < n; i++ {
			d := float64(i) - float64(n-1)/2
			varr += d * d
		}
		return sum / varr
	}
	plain, jittered := order(false), order(true)
	t.Logf("pop order against flat index: plain %.3f, jittered %.3f", plain, jittered)
	if plain > -0.99 {
		t.Errorf("plain push no longer pops in reverse insertion order (%.3f); this test's premise is gone", plain)
	}
	if math.Abs(jittered) > 0.05 {
		t.Errorf("jittered push still correlates with flat index at %.3f", jittered)
	}
}