230 lines
10 KiB
Go
230 lines
10 KiB
Go
// 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
|
|
}
|