package plates

import (
	"math"
	"sort"
	"testing"

	"salty/terrain/internal/world"
)

// testPlanet is a small cylinder with the same shape of arithmetic as a real one: a whole number of columns
// and a noise period that divides the circumference.
func testPlanet(t *testing.T) world.Planet {
	t.Helper()
	p, err := world.New(40000, 8, 100, 50, 0, 40000)
	if err != nil {
		t.Fatalf("planet: %v", err)
	}
	return p
}

// handModel builds a Model with a raster set by the caller, so a test can ask about a boundary in a known
// place. Nothing here goes through Build: the point is to control the partition rather than the seed.
func handModel(t *testing.T, p world.Planet, gw, gh int, cell []int16, ps []Plate) *Model {
	t.Helper()
	return &Model{
		P: p, Cfg: Default(), Plates: ps,
		GW: gw, GH: gh, GCellM: p.CircumferenceM() / float64(gw),
		Cell: cell,
	}
}

// stripes paints two vertical bands: plate 1 from column lo up to hi, plate 0 everywhere else.
//
// That is **two** contacts, not one, and it is worth saying why every test here is written in pairs. A
// cylinder cut into two strips has a margin at each end of each strip, and under a pure translation the
// plates are closing at one of them and opening at the other by exactly the same amount. There is no way to
// arrange two plates on a cylinder that only collide. The invariant is the test.
func stripes(gw, gh, lo, hi int) []int16 {
	cell := make([]int16, gw*gh)
	for gy := range gh {
		for gx := range gw {
			if gx >= lo && gx < hi {
				cell[gy*gw+gx] = 1
			}
		}
	}
	return cell
}

// rollX moves every column east by n, wrapping. Rolling the map is the whole seam test: the cylinder has no
// preferred meridian, so a partition and the same partition rolled must produce the same tectonics.
func rollX(cell []int16, gw, gh, n int) []int16 {
	out := make([]int16, len(cell))
	for gy := range gh {
		for gx := range gw {
			out[gy*gw+((gx+n)%gw)] = cell[gy*gw+gx]
		}
	}
	return out
}

// closingPair drives plate 0 east and plate 1 west, with no spin.
func closingPair() []Plate {
	return []Plate{
		{ID: 0, TransXM: 0.02, Continental: true, AreaCells: 1},
		{ID: 1, TransXM: -0.02, Continental: true, AreaCells: 1},
	}
}

// meanClosing is a boundary's average closing rate, in metres a year.
func meanClosing(b Boundary) float64 {
	if len(b.V) == 0 {
		return 0
	}
	total := 0.0
	for _, v := range b.V {
		total += v.ClosingMYr
	}
	return total / float64(len(b.V))
}

// sortedMeans is every boundary's mean closing rate, in order: the signature of a whole planet's tectonics,
// independent of which order the boundaries happened to be found in.
func sortedMeans(bs []Boundary) []float64 {
	out := make([]float64, len(bs))
	for i, b := range bs {
		out[i] = meanClosing(b)
	}
	sort.Float64s(out)
	return out
}

func kinds(bs []Boundary) []string {
	out := make([]string, len(bs))
	for i, b := range bs {
		out[i] = b.Dominant().String()
	}
	sort.Strings(out)
	return out
}

func TestABoundaryAcrossTheSeamIsOneBoundary(t *testing.T) {
	p := testPlanet(t)
	const gw, gh = 200, 100
	m := handModel(t, p, gw, gh, stripes(gw, gh, 50, 150), closingPair())
	bs := m.buildBoundaries()

	if len(bs) != 2 {
		t.Fatalf("two vertical contacts on a cylinder, got %d boundaries", len(bs))
	}
	// Roll the partition so one contact sits exactly on the meridian. If the seam were special the boundary
	// through it would come back cut in half - two chains of half the length - or with a whole circumference
	// of jump in the middle of it.
	rolled := handModel(t, p, gw, gh, rollX(m.Cell, gw, gh, 50), closingPair())
	rbs := rolled.buildBoundaries()
	if len(rbs) != 2 {
		t.Fatalf("after rolling the map onto the seam: %d boundaries, want 2", len(rbs))
	}

	for _, b := range rbs {
		for i := 0; i+1 < len(b.V); i++ {
			d := math.Hypot(b.V[i+1].XM-b.V[i].XM, b.V[i+1].YM-b.V[i].YM)
			if d > maxGapCells*m.GCellM*1.01 {
				t.Fatalf("a %.0f m step between neighbouring vertices: X was wrapped, not unwrapped", d)
			}
		}
		if got, want := b.LengthM(), bs[0].LengthM(); math.Abs(got-want) > m.GCellM {
			t.Errorf("rolled boundary is %.0f m, unrolled %.0f m", got, want)
		}
	}

	// The tectonics have to be the same planet, not just the same shape: the rolled map's margins close at
	// the same rates as the unrolled map's.
	before, after := sortedMeans(bs), sortedMeans(rbs)
	for i := range before {
		if math.Abs(before[i]-after[i]) > 1e-9 {
			t.Errorf("closing rate %d is %.6g before the roll and %.6g after", i, before[i], after[i])
		}
	}

	circ := p.CircumferenceM()
	crossed := false
	for _, b := range rbs {
		for _, v := range b.V {
			if v.XM < 0 || v.XM >= circ {
				crossed = true
			}
		}
	}
	if !crossed {
		t.Error("no vertex outside 0..circumference, so nothing was unwrapped and the roll tested nothing")
	}
}

func TestClosingIsTheSameWhicheverPlateIsCalledA(t *testing.T) {
	p := testPlanet(t)
	const gw, gh = 200, 100

	forward := handModel(t, p, gw, gh, stripes(gw, gh, 50, 150), closingPair())
	// Swap which stripe belongs to which plate, and swap the motions with it. Physically nothing has moved:
	// the same two materials are being driven together at the same margin. Every closing rate must therefore
	// come back identical, because the normal flips and the relative velocity flips with it.
	swapped := stripes(gw, gh, 50, 150)
	for i := range swapped {
		swapped[i] = 1 - swapped[i]
	}
	reverse := handModel(t, p, gw, gh, swapped, []Plate{
		{ID: 0, TransXM: -0.02, Continental: true, AreaCells: 1},
		{ID: 1, TransXM: 0.02, Continental: true, AreaCells: 1},
	})

	fm, rm := sortedMeans(forward.buildBoundaries()), sortedMeans(reverse.buildBoundaries())
	if len(fm) != len(rm) {
		t.Fatalf("%d boundaries one way round and %d the other", len(fm), len(rm))
	}
	for i := range fm {
		if math.Abs(fm[i]-rm[i]) > 1e-9 {
			t.Errorf("closing rate %d is %.6g one way round and %.6g the other; the sign convention is not "+
				"symmetric", i, fm[i], rm[i])
		}
	}
}

func TestOneMarginClosesAndTheOtherOpens(t *testing.T) {
	p := testPlanet(t)
	const gw, gh = 200, 100
	bs := handModel(t, p, gw, gh, stripes(gw, gh, 50, 150), closingPair()).buildBoundaries()
	if len(bs) != 2 {
		t.Fatalf("got %d boundaries, want 2", len(bs))
	}

	means := sortedMeans(bs)
	if means[0] >= 0 || means[1] <= 0 {
		t.Fatalf("closing rates %.4g and %.4g; two strips on a cylinder give one of each", means[0], means[1])
	}
	// Equal and opposite, because a pure translation is the same relative velocity at both margins and the
	// only thing that differs is which way the normal points.
	if math.Abs(means[0]+means[1]) > 1e-9 {
		t.Errorf("closing rates %.6g and %.6g are not equal and opposite", means[0], means[1])
	}
	if got := math.Abs(means[1]); math.Abs(got-0.04) > 1e-9 {
		t.Errorf("plates at 2 cm/yr each close at %.4g m/yr; 0.04 is the sum of the two speeds", got)
	}
}

func TestContinentsCollideAndOceansSubduct(t *testing.T) {
	p := testPlanet(t)
	const gw, gh = 200, 100
	cell := stripes(gw, gh, 50, 150)

	both := handModel(t, p, gw, gh, cell, closingPair()).buildBoundaries()
	if got, want := kinds(both), []string{"collision", "rift"}; !sameStrings(got, want) {
		t.Errorf("two continental plates give %v, want %v", got, want)
	}

	// The same geometry with one side oceanic: the closing margin is a subduction zone and the arc belongs
	// to the continent, because continental crust is too buoyant to go down.
	ps := closingPair()
	ps[1].Continental = false
	oceanic := handModel(t, p, gw, gh, cell, ps).buildBoundaries()
	if got, want := kinds(oceanic), []string{"ridge", "subduction"}; !sameStrings(got, want) {
		t.Errorf("continent against ocean gives %v, want %v", got, want)
	}
	for _, b := range oceanic {
		for _, v := range b.V {
			if v.Kind == Subduction && v.Over != 0 {
				t.Fatalf("the overriding plate is %d, but plate 1 is the oceanic one", v.Over)
			}
			if v.Kind != Subduction && v.Over != -1 {
				t.Fatalf("a %q vertex carries an overriding plate", v.Kind)
			}
		}
	}
}

func sameStrings(a, b []string) bool {
	if len(a) != len(b) {
		return false
	}
	for i := range a {
		if a[i] != b[i] {
			return false
		}
	}
	return true
}

func TestMotionAlongTheLineIsATransform(t *testing.T) {
	p := testPlanet(t)
	const gw, gh = 200, 100
	// Vertical contacts, both plates sliding north and south: the relative motion is entirely along the line,
	// so neither margin closes or opens.
	m := handModel(t, p, gw, gh, stripes(gw, gh, 50, 150), []Plate{
		{ID: 0, TransYM: 0.02, Continental: true, AreaCells: 1},
		{ID: 1, TransYM: -0.02, Continental: true, AreaCells: 1},
	})
	for _, b := range m.buildBoundaries() {
		if got := b.Dominant(); got != Transform {
			t.Errorf("plates sliding past each other give %q, want %q", got, Transform)
		}
		if got := math.Abs(meanClosing(b)); got > 1e-9 {
			t.Errorf("a transform margin closes at %.3g m/yr", got)
		}
	}
}

func TestKindOf(t *testing.T) {
	oblique := 60 * math.Pi / 180
	cases := []struct {
		name            string
		closing, slip   float64
		bothContinental bool
		want            Kind
	}{
		{"head-on continental", 1, 0, true, Collision},
		{"head-on with an ocean", 1, 0, false, Subduction},
		{"opening continental", -1, 0, true, Rift},
		{"opening with an ocean", -1, 0, false, Ridge},
		{"pure slip", 0, 1, true, Transform},
		{"oblique but still closing", 1, 1.5, true, Collision},
		{"slip has taken over", 1, 2, true, Transform},
	}
	for _, c := range cases {
		if got := kindOf(c.closing, c.slip, oblique, c.bothContinental); got != c.want {
			t.Errorf("%s: got %q, want %q", c.name, got, c.want)
		}
	}
}

func TestSpinMakesAMarginChangeAlongItsLength(t *testing.T) {
	p := testPlanet(t)
	const gw, gh = 200, 100
	cell := stripes(gw, gh, 50, 150)

	// Pure translation: the relative velocity is the same everywhere, so a straight margin closes at one rate
	// from end to end. That uniformity is exactly what the in-plane pole exists to break.
	flat := handModel(t, p, gw, gh, cell, closingPair()).buildBoundaries()
	if spread := worstSpread(flat); spread > 1e-9 {
		t.Errorf("without spin a margin varies by %.3g m/yr along its own length; it should not", spread)
	}

	// The same plates with one rotating about a pole off to the side close at one end and slide at the other.
	spun := closingPair()
	spun[0].PoleXM, spun[0].PoleYM = 0, 0
	spun[0].OmegaRadYr = 2e-6
	spinning := handModel(t, p, gw, gh, cell, spun).buildBoundaries()
	if spread := worstSpread(spinning); spread < 1e-3 {
		t.Errorf("with spin a margin varies by only %.3g m/yr; the rotation is not reaching the boundary",
			spread)
	}
}

// worstSpread is the largest range of closing rates found *within* a single boundary. Within, not across:
// two margins of the same pair legitimately differ, and measuring across them would report that difference
// as variation along a line.
func worstSpread(bs []Boundary) float64 {
	worst := 0.0
	for _, b := range bs {
		lo, hi := math.Inf(1), math.Inf(-1)
		for _, v := range b.V {
			lo = math.Min(lo, v.ClosingMYr)
			hi = math.Max(hi, v.ClosingMYr)
		}
		if !math.IsInf(lo, 1) && hi-lo > worst {
			worst = hi - lo
		}
	}
	return worst
}

func TestVelocityIsContinuousAcrossTheSeam(t *testing.T) {
	p := testPlanet(t)
	circ := p.CircumferenceM()
	// A plate whose pole sits just east of the meridian. Measured without wrapping, the lever arm a metre
	// west of the seam would be a whole circumference long, and the plate would spin the wrong way there.
	pl := Plate{ID: 0, PoleXM: 10, PoleYM: 0, OmegaRadYr: 1e-6, TransXM: 0.01}

	ax, ay := pl.VelocityAt(p, circ-1, 0)
	bx, by := pl.VelocityAt(p, 1, 0)
	acrossSeam := math.Hypot(ax-bx, ay-by)

	// The same two-metre gap in open map, away from the meridian: the seam must cost nothing extra.
	cx, cy := pl.VelocityAt(p, circ/2-1, 0)
	dx, dy := pl.VelocityAt(p, circ/2+1, 0)
	elsewhere := math.Hypot(cx-dx, cy-dy)

	if math.Abs(acrossSeam-elsewhere) > 1e-12 {
		t.Errorf("velocity changes by %.3g m/yr over two metres at the meridian and %.3g m/yr over two "+
			"metres anywhere else", acrossSeam, elsewhere)
	}
	// And the failure this guards against is enormous, not subtle: an unwrapped lever arm would be a whole
	// circumference and give a jump of omega*circ.
	if acrossSeam > pl.OmegaRadYr*circ/100 {
		t.Errorf("velocity jumps by %.3g m/yr at the meridian; the lever arm was not wrapped", acrossSeam)
	}
}

func TestBuildCoversThePlanetAndReadsTheLandMask(t *testing.T) {
	p := testPlanet(t)
	cfg := Default()
	cfg.Count = 6

	allSea, err := Build(p, 7, cfg, func(xM, yM float64) bool { return false })
	if err != nil {
		t.Fatalf("build: %v", err)
	}
	total := 0
	for _, pl := range allSea.Plates {
		total += pl.AreaCells
		if pl.Continental {
			t.Errorf("plate %d is continental on a planet with no land", pl.ID)
		}
	}
	if total != allSea.GW*allSea.GH {
		t.Errorf("plates cover %d cells of %d; the partition has holes", total, allSea.GW*allSea.GH)
	}
	if len(allSea.Boundaries) == 0 {
		t.Fatal("six plates and no boundaries between them")
	}
	for _, b := range allSea.Boundaries {
		if b.A >= b.B {
			t.Errorf("boundary pair (%d, %d) is not ordered", b.A, b.B)
		}
		for _, v := range b.V {
			if v.Kind == Collision || v.Kind == Rift {
				t.Errorf("a %q on a planet with no continental plate at all", v.Kind)
			}
		}
	}

	allLand, err := Build(p, 7, cfg, func(xM, yM float64) bool { return true })
	if err != nil {
		t.Fatalf("build: %v", err)
	}
	for _, pl := range allLand.Plates {
		if !pl.Continental {
			t.Errorf("plate %d is oceanic on a planet that is all land", pl.ID)
		}
	}
	for _, b := range allLand.Boundaries {
		for _, v := range b.V {
			if v.Kind == Subduction || v.Kind == Ridge {
				t.Errorf("a %q with no oceanic plate to make it", v.Kind)
			}
		}
	}
}

func TestTheSameSeedGivesTheSamePlanet(t *testing.T) {
	p := testPlanet(t)
	cfg := Default()
	land := func(xM, yM float64) bool { return yM > 5000 && yM < 12000 }

	a, err := Build(p, 9342, cfg, land)
	if err != nil {
		t.Fatalf("build: %v", err)
	}
	b, err := Build(p, 9342, cfg, land)
	if err != nil {
		t.Fatalf("build: %v", err)
	}
	if len(a.Boundaries) != len(b.Boundaries) {
		t.Fatalf("%d boundaries then %d; the set is not deterministic", len(a.Boundaries), len(b.Boundaries))
	}
	for i := range a.Boundaries {
		if a.Boundaries[i].A != b.Boundaries[i].A || a.Boundaries[i].B != b.Boundaries[i].B {
			t.Fatalf("boundary %d is a different pair on the second run", i)
		}
		if math.Abs(a.Boundaries[i].LengthM()-b.Boundaries[i].LengthM()) > 1e-9 {
			t.Fatalf("boundary %d is a different length on the second run", i)
		}
	}
}

func TestTheWarpBendsMarginsAndTooMuchOfItBreaksThem(t *testing.T) {
	p := testPlanet(t)

	// Sinuosity - the line's own length over the distance between its ends - is what "a margin is not a
	// ruled line" means as a number. A weighted Voronoi edge is a circular arc even at no warp, so the
	// baseline is a little over 1 rather than exactly 1.
	measure := func(warp float64) (sinuosity float64, boundaries int) {
		cfg := Default()
		cfg.Count = 7
		cfg.WarpFraction = warp
		m, err := Build(p, 3630, cfg, func(xM, yM float64) bool { return false })
		if err != nil {
			t.Fatalf("warp %.2f: %v", warp, err)
		}
		total, n := 0.0, 0
		for _, b := range m.Boundaries {
			if len(b.V) < 10 {
				continue
			}
			last := b.V[len(b.V)-1]
			if straight := math.Hypot(last.XM-b.V[0].XM, last.YM-b.V[0].YM); straight > 0 {
				total += b.LengthM() / straight
				n++
			}
		}
		if n == 0 {
			return 0, len(m.Boundaries)
		}
		return total / float64(n), len(m.Boundaries)
	}

	straight, straightCount := measure(0)
	warped, warpedCount := measure(Default().WarpFraction)
	if straight > 1.05 {
		t.Errorf("an unwarped partition already has a sinuosity of %.3f; it should be close to a polygon net",
			straight)
	}
	if warped <= straight*1.03 {
		t.Errorf("the warp takes sinuosity from %.3f to %.3f, which is no bend at all", straight, warped)
	}

	// Past about half a plate spacing the displacement folds back on itself and the partition grows islands
	// of one plate inside another, which the tracer faithfully chains into extra rings. The count is the
	// symptom, and this is the bound warpFineGain and the default are set under.
	_, tooMuch := measure(0.5)
	if tooMuch <= warpedCount {
		t.Skipf("no fragmentation at warp 0.5 on this seed (%d boundaries against %d); the bound still holds "+
			"but this seed does not show it", tooMuch, warpedCount)
	}
	if warpedCount != straightCount {
		t.Errorf("the default warp changed the boundary count from %d to %d; it should bend margins, not "+
			"create them", straightCount, warpedCount)
	}
}