// Package stats is how a run is judged. "It reads as real geology" is not a screenshot; it is a straight
// slope-area plot and an S-shaped hypsometric curve, and this package produces both.
//
// The project's habit already is to measure rather than eyeball — a slope histogram settled the noise tuning
// and attributed rill damage to a specific pass — and these are the two standard checks the incoming spec
// added on top. The slope-area exponent in particular is the direct test of whether the fluvial pass did the
// thing it exists to do, so it is the proof that closes build-order step 4.
package stats

import (
	"fmt"
	"math"
)

type Bin struct {
	LogA float64 `json:"log_area"`
	LogS float64 `json:"log_slope"`
	N    int     `json:"n"`
}

// SlopeArea is the stream-power signature. At steady state S = (U/K)^(1/n) * A^(-m/n), so log S against
// log A is a straight line of gradient -m/n: -0.5 at the defaults. A curved or scattered plot means K, m, n
// or the run length is wrong, and no amount of detail noise will hide it.
//
// The catch, and it took a bad R2 to notice: that relation has the same gradient but a *different intercept*
// for every uplift rate. This map's uplift spans 0.2 to 5 mm/yr, so regressing every channel together stacks
// twenty-five-fold-separated parallel lines into a cloud and fits nonsense to it. Slope is therefore
// normalised by (U/K)^(1/n) first, which collapses every regime onto one line through the origin and tests
// the exponent rather than the uplift field's heterogeneity. RawExponent keeps the unnormalised fit, which is
// what a single-uplift map would report and is worth seeing next to it.
type SlopeArea struct {
	Bins        []Bin   `json:"bins"`
	Exponent    float64 `json:"exponent"`
	RawExponent float64 `json:"raw_exponent"`
	Expected    float64 `json:"expected"`
	R2          float64 `json:"r2"`
	RawR2       float64 `json:"raw_r2"`
	Channels    int     `json:"channel_cells"`
	ThreshKm2   float64 `json:"threshold_km2"`
}

// Hypsometry is the second check: cumulative area against normalised elevation should be S-shaped. The
// integral is the single number - convex and high means too young or too much uplift, concave and low means
// over-eroded. Mature landscapes sit near 0.4 to 0.6.
type Hypsometry struct {
	Integral float64   `json:"integral"`
	Curve    []float64 `json:"curve"` // area fraction at 11 elevation fractions, 0.0 to 1.0
}

type Slopes struct {
	Under15Deg float64 `json:"under_15_deg"`
	Under30Deg float64 `json:"under_30_deg"`
	Over50Deg  float64 `json:"over_50_deg"`
	MedianDeg  float64 `json:"median_deg"`
}

type Report struct {
	LandFraction float64 `json:"land_fraction"`
	ClipFraction float64 `json:"clip_fraction"`
	// Min and Max span the whole field, sea floor included, because that is what the 16-bit encoding has to
	// fit. Land relief is the number that says anything about the terrain, and they are not the same: a
	// -180 m sea floor flatters the relief by 180 m for free.
	ReliefM         float64    `json:"relief_m"`
	MinM            float64    `json:"min_m"`
	MaxM            float64    `json:"max_m"`
	LandMinM        float64    `json:"land_min_m"`
	LandMaxM        float64    `json:"land_max_m"`
	LandReliefM     float64    `json:"land_relief_m"`
	Slopes          Slopes     `json:"slopes"`
	SlopeArea       SlopeArea  `json:"slope_area"`
	Hypsometry      Hypsometry `json:"hypsometry"`
	DrainageDensity float64    `json:"drainage_density_per_km"`

	// LeafFraction is the share of land cells that drain nothing but themselves. See accumulate.go: it is
	// the one number that separates a drainage network from a comb of parallel non-converging flow lines.
	LeafFraction float64 `json:"leaf_fraction"`

	// LandCells is how much land the world has and MeasuredLandCells how much of it the land statistics
	// below actually walked. They differ only on a partial run - `bake --only` leaves most of a planet at sea
	// level - and when they do, every distribution here describes the part that was solved while the extent
	// above describes the whole cylinder. Summary says so rather than leaving the two to be compared.
	LandCells         int64 `json:"land_cells"`
	MeasuredLandCells int64 `json:"measured_land_cells"`

	// Buckets is the whole-map aggregates split by the uplift class that caused them.
	// The map-wide median above cannot tell a mountain belt from a plain, and that is the question.
	Buckets []UpliftBucket `json:"uplift_buckets"`
}

// expectedGradient is the -m/n the theory predicts, kept in one place so the verdict compares the fit against
// the exponents the run was actually configured with rather than against the defaults.
var expectedGradient = -0.5

// SetExpected is called once from the command before any report is computed.
func SetExpected(m, n float64) {
	if n != 0 {
		expectedGradient = -m / n
	}
}

// fitLine is an ordinary least-squares fit returning the gradient and R².
func fitLine(x, y []float64) (float64, float64) {
	n := float64(len(x))
	if n < 3 {
		return 0, 0
	}
	var sx, sy, sxx, sxy float64
	for i := range x {
		sx += x[i]
		sy += y[i]
		sxx += x[i] * x[i]
		sxy += x[i] * y[i]
	}
	den := n*sxx - sx*sx
	if math.Abs(den) < 1e-12 {
		return 0, 0
	}
	grad := (n*sxy - sx*sy) / den
	intercept := (sy - grad*sx) / n
	mean := sy / n
	var ssRes, ssTot float64
	for i := range x {
		pred := grad*x[i] + intercept
		ssRes += (y[i] - pred) * (y[i] - pred)
		ssTot += (y[i] - mean) * (y[i] - mean)
	}
	if ssTot < 1e-12 {
		return grad, 0
	}
	return grad, 1 - ssRes/ssTot
}

// Summary is the one block a run prints. Written so the numbers that decide whether the run was any good are
// the ones you see without asking.
func (r Report) Summary() string {
	sa := r.SlopeArea
	verdict := "no channels: the solve did not organise"
	if sa.Channels > 0 && len(sa.Bins) >= 3 {
		switch {
		case sa.R2 >= 0.9 && math.Abs(sa.Exponent-sa.Expected) < 0.15:
			verdict = "straight and at the expected gradient: stream power is doing its job"
		case sa.R2 >= 0.9:
			verdict = "straight but off gradient: K, m or n is wrong, or the run is too short"
		default:
			verdict = "scattered: not at steady state, or pits are routing badly"
		}
	}
	hyp := "mature"
	switch {
	case r.Hypsometry.Integral > 0.6:
		hyp = "convex: too young, or too much uplift"
	case r.Hypsometry.Integral < 0.35:
		hyp = "concave: over-eroded"
	}
	// A partial run measures the whole cylinder's extent and only the solved landmasses' ground, and the two
	// sitting next to each other invite exactly the wrong comparison. Say so, rather than leave somebody to
	// work out afterwards why the drainage density looked impossible.
	partial := ""
	if r.LandCells > 0 && r.MeasuredLandCells > 0 && r.MeasuredLandCells < r.LandCells {
		partial = fmt.Sprintf(
			"  PARTIAL: the line above is the whole world; everything below is the %.0f%% of its land that\n"+
				"           was actually solved (%d of %d cells). The two are not comparable.\n",
			100*float64(r.MeasuredLandCells)/float64(r.LandCells), r.MeasuredLandCells, r.LandCells)
	}
	return fmt.Sprintf(
		"  field %.0f..%.0f m; land %.0f..%.0f m (relief %.0f m), %.0f%% land, %.2f%% clipped\n"+
			"%s"+
			"  slopes: %.0f%% under 15 deg, %.0f%% under 30, %.1f%% over 50, median %.1f deg\n"+
			"  slope-area: exponent %.3f (expect %.3f), R2 %.3f over %d bins, %d channel cells above %.2f km2\n"+
			"              unnormalised %.3f, R2 %.3f (heterogeneous uplift, so this one is expected to be worse)\n"+
			"              %s\n"+
			"  hypsometric integral %.3f (%s); drainage density %.2f /km; %.1f%% of land drains nothing\n"+
			"%s",
		r.MinM, r.MaxM, r.LandMinM, r.LandMaxM, r.LandReliefM, r.LandFraction*100, r.ClipFraction*100,
		partial,
		r.Slopes.Under15Deg*100, r.Slopes.Under30Deg*100, r.Slopes.Over50Deg*100, r.Slopes.MedianDeg,
		sa.Exponent, sa.Expected, sa.R2, len(sa.Bins), sa.Channels, sa.ThreshKm2,
		sa.RawExponent, sa.RawR2,
		verdict, r.Hypsometry.Integral, hyp, r.DrainageDensity, r.LeafFraction*100,
		BucketSummary(r.Buckets))
}

// Uplift buckets: the measurement that decides whether a plain is a plain.
//
// Every aggregate above is taken over the whole land mask, and that is exactly what hid the problem this
// bucketing was added to find. A continent whose mountains are at 35 degrees and whose plains are at 32
// reports a median of 31 and looks, from the summary, like a mountainous map — which it is, but not for the
// reason anyone assumed. Splitting by the uplift rate that *caused* the slope separates the two questions:
// "are the mountains right" and "are the plains plains".
//
// Uplift is the right axis rather than elevation. Elevation is the output of the solve, so bucketing by it
// mixes a low mountain valley in with a plain and moves the boundary every time a constant changes; uplift
// is an input, fixed before the first step, and it is the term that sets steady-state slope through
// S = U/(K*A^m). A cell's bucket therefore does not move when the run does.

// UpliftBucket is one class of the uplift field and what the landscape did with it.
type UpliftBucket struct {
	Name        string  `json:"name"`
	LoMmYr      float64 `json:"lo_mm_yr"`
	HiMmYr      float64 `json:"hi_mm_yr"`
	LandFrac    float64 `json:"land_fraction"` // share of land in this bucket
	MedianDeg   float64 `json:"median_deg"`
	P90Deg      float64 `json:"p90_deg"`
	MedianRelM  float64 `json:"median_relief_m"` // local relief, max-min over the window below
	WindowM     float64 `json:"relief_window_m"`
	NearTalus   float64 `json:"near_talus_fraction"` // within 2 degrees of the angle of repose
	MedianElevM float64 `json:"median_elev_m"`
	Cells       int     `json:"cells"`
}

// BucketSummary is the block the buckets print. Kept separate from Summary so a run that has no uplift field
// to hand still prints the rest.
func BucketSummary(bs []UpliftBucket) string {
	if len(bs) == 0 {
		return ""
	}
	s := "  by uplift class:\n"
	for _, b := range bs {
		hi := fmt.Sprintf("%.2f", b.HiMmYr)
		if b.HiMmYr >= 100 {
			hi = "  up"
		}
		s += fmt.Sprintf("    %-9s %.2f..%s mm/yr  %4.0f%% of land  slope %4.1f deg median, %4.1f P90  "+
			"relief %5.0f m/%.0f m  at talus %4.0f%%  median %.0f m\n",
			b.Name, b.LoMmYr, hi, b.LandFrac*100, b.MedianDeg, b.P90Deg, b.MedianRelM, b.WindowM,
			b.NearTalus*100, b.MedianElevM)
	}
	return s
}