islaApocalypse-v2/Tools/Scripts/LandHistogram.cs

using System;

namespace IslaApocalypse.Tools
{
	/// <summary>
	/// ⭐ THE LAND-HEIGHT DISTRIBUTION — one instrument doing two jobs.
	///
	/// ═══ JOB 1: THE CALIBRATION ENGINE ═══
	///
	/// The redistribution curve's knots ARE percentiles of this distribution (P60/73/83/88/96/99).
	/// This class measures them. That is what makes the band shares 60/13/10/5/8/3/1 exact by
	/// construction rather than approximately right.
	///
	/// > ⚠ The reference SHIPPED the six resulting literals and threw the instrument away. No
	/// > sampler, no histogram, no percentile helper survives at the tag — so its knots could never
	/// > be re-derived, only trusted. Rebuilding the measuring device is the point of this file:
	/// > a calibration you cannot re-run is a magic number with a good story.
	///
	/// ═══ JOB 2: THE DIAGNOSTIC ═══
	///
	/// It is also what separates the three candidate causes of "flat, undramatic upper terrain",
	/// which look identical in a render and completely different here:
	///
	///   1. CURVE SQUASHING      — the SHAPED histogram piles mass at the bench/plateau output
	///                             heights. The curve is flattening ground that had relief.
	///   2. NOISE EMPTY UP HIGH  — the RAW histogram's top (above ~P96) is a thin sliver spread over
	///                             a wide range. There is nothing up there to shape.
	///   3. SHARE ALLOCATION     — little land is TARGETED into plateau+spike (3 % + 1 %) by
	///                             construction. The curve is doing exactly what it was told.
	///
	/// These have different fixes — reshape the segments, change the noise, or reallocate the
	/// shares — so guessing which one it is costs a whole task.
	///
	/// ═══ WHY A HISTOGRAM AND NOT A SORTED SAMPLE ARRAY ═══
	///
	/// Exact quantiles want every sample sorted; the reference pooled 340 M of them. At 4 bytes each
	/// that is 1.3 GB and a sort to match. A fine fixed-width histogram with IN-BIN LINEAR
	/// INTERPOLATION gives quantiles accurate to well under a bin, in flat memory, at any pool size,
	/// and streams across seeds without holding a single sample. At the default
	/// <see cref="CalibrationBinWidth"/> the resolution is 1e-4 raw ≈ 2.5 cm of world height —
	/// four decimal places on a knot, against reference literals quoted to six.
	///
	/// ⚠ Engine-free (System only). It sits in <c>Tools/</c> rather than <c>Core/</c> because it is a
	/// MEASURING INSTRUMENT for the generator, not a contract about the world — the same reasoning
	/// that keeps <c>IslandFalloff</c> here. Core carries what the world IS; Tools carries what we
	/// point at it.
	/// </summary>
	public sealed class LandHistogram
	{
		/// <summary>
		/// Bin width for CALIBRATION, in raw height units. 1e-4 raw ≈ 2.5 cm — finer than any knot
		/// distinction that could matter, and 38,500 bins over the working range is 300 KB.
		/// </summary>
		public const float CalibrationBinWidth = 1e-4f;

		/// <summary>
		/// Top of the binned range, raw. Generous: pass-1 land maxes near 1.45 and the curve's tail
		/// can exceed the 420 m cap. Anything above lands in the overflow bin and is REPORTED, never
		/// silently dropped.
		/// </summary>
		public const float DefaultTop = 4.0f;

		/// <summary>Heights at or below this are not land and are not counted. The curve's identity threshold.</summary>
		public readonly float SeaLevel;

		/// <summary>Bin width in raw height units.</summary>
		public readonly float BinWidth;

		/// <summary>Top of the binned range; samples above it go to <see cref="OverflowCount"/>.</summary>
		public readonly float Top;

		private readonly long[] _counts;

		/// <summary>Land samples at or above <see cref="Top"/>. ⚠ Reported, not hidden.</summary>
		public long OverflowCount { get; private set; }

		/// <summary>Total land samples accumulated, overflow included.</summary>
		public long TotalLand { get; private set; }

		/// <summary>Lowest and highest land sample seen, exactly (not bin-quantized).</summary>
		public float MinLand { get; private set; } = float.MaxValue;
		public float MaxLand { get; private set; } = float.MinValue;

		/// <summary>How many fields have been pooled in. The calibration pool's size.</summary>
		public int FieldsPooled { get; private set; }

		public LandHistogram(float seaLevel, float binWidth = CalibrationBinWidth, float top = DefaultTop)
		{
			if (binWidth <= 0f) throw new ArgumentOutOfRangeException(nameof(binWidth), binWidth, "Bin width must be positive.");
			if (top <= seaLevel) throw new ArgumentOutOfRangeException(nameof(top), top, "Top must exceed sea level.");

			SeaLevel = seaLevel;
			BinWidth = binWidth;
			Top = top;
			_counts = new long[(int)MathF.Ceiling((top - seaLevel) / binWidth)];
		}

		/// <summary>Number of bins (excluding overflow).</summary>
		public int BinCount => _counts.Length;

		/// <summary>Raw height at the low edge of bin <paramref name="i"/>.</summary>
		public float BinLow(int i) => SeaLevel + i * BinWidth;

		/// <summary>Sample count in bin <paramref name="i"/>.</summary>
		public long BinCountAt(int i) => _counts[i];

		/// <summary>
		/// Pool one field's LAND samples in. Call repeatedly to build a multi-seed pool — the
		/// reference calibrated across a pooled batch, and one seed's distribution is not the
		/// island's.
		///
		/// ⚠ Strictly <c>&gt; SeaLevel</c>, matching the curve's own <c>h &lt;= Sea → identity</c>
		/// test. A pixel exactly at sea is not land, and counting it would put a spike in bin 0 that
		/// drags every low percentile down.
		/// </summary>
		public void Accumulate(float[,] field, int mapSize)
		{
			for (int x = 0; x < mapSize; x++)
			{
				for (int y = 0; y < mapSize; y++)
				{
					float h = field[x, y];
					if (h <= SeaLevel) continue;

					TotalLand++;
					if (h < MinLand) MinLand = h;
					if (h > MaxLand) MaxLand = h;

					int bin = (int)((h - SeaLevel) / BinWidth);
					if (bin >= _counts.Length) OverflowCount++;
					else _counts[bin]++;
				}
			}
			FieldsPooled++;
		}

		/// <summary>
		/// The quantile at <paramref name="percent"/> (0..100) — a raw height, interpolated inside
		/// its bin so the answer is not quantized to <see cref="BinWidth"/>.
		///
		/// ⚠ Throws on an empty pool rather than returning sea level. An all-ocean seed silently
		/// calibrating every knot to 0.15 is exactly the kind of quiet nonsense that ships.
		/// </summary>
		public float Quantile(double percent)
		{
			if (TotalLand == 0)
				throw new InvalidOperationException(
					"[LandHistogram] No land samples pooled — cannot take a quantile. " +
					"Check the sea level and that pass 1 actually produced an island.");
			if (percent < 0.0 || percent > 100.0)
				throw new ArgumentOutOfRangeException(nameof(percent), percent, "A percentile is 0..100.");

			double target = percent / 100.0 * TotalLand;
			long cum = 0;

			for (int i = 0; i < _counts.Length; i++)
			{
				long c = _counts[i];
				if (c == 0) continue;
				if (cum + c >= target)
				{
					// Linear position inside the bin: the samples in it are assumed uniform, which
					// is the standard histogram-quantile assumption and is a sub-bin error.
					double within = (target - cum) / c;
					return BinLow(i) + (float)(within * BinWidth);
				}
				cum += c;
			}

			// Only reachable when the quantile falls in the overflow — a real answer we cannot
			// resolve, so say so rather than returning Top as if it were measured.
			throw new InvalidOperationException(
				$"[LandHistogram] P{percent} falls above the binned range (top {Top}); " +
				$"{OverflowCount} of {TotalLand} samples overflowed. Raise `top` and re-measure.");
		}

		/// <summary>Fraction of land strictly below <paramref name="h"/>, interpolated within the bin.</summary>
		public double FractionBelow(float h)
		{
			if (TotalLand == 0) return 0.0;
			if (h <= SeaLevel) return 0.0;

			int bin = (int)((h - SeaLevel) / BinWidth);
			if (bin >= _counts.Length) return 1.0;

			long cum = 0;
			for (int i = 0; i < bin; i++) cum += _counts[i];

			double within = (h - BinLow(bin)) / BinWidth;
			return (cum + within * _counts[bin]) / TotalLand;
		}

		/// <summary>Fraction of land in <c>[lo, hi)</c>.</summary>
		public double FractionBetween(float lo, float hi) => Math.Max(0.0, FractionBelow(hi) - FractionBelow(lo));

		/// <summary>The tallest bin's count — the y-axis a plot needs.</summary>
		public long PeakBinCount()
		{
			long peak = 0;
			foreach (long c in _counts) if (c > peak) peak = c;
			return peak;
		}

		/// <summary>
		/// Re-bin into a coarser histogram for DISPLAY. The calibration histogram has 38,500 bins;
		/// a plot has room for a few hundred, and drawing one bin per pixel column of a 1,200 px plot
		/// would alias the distribution into noise.
		/// </summary>
		public LandHistogram Rebin(float displayBinWidth)
		{
			var outH = new LandHistogram(SeaLevel, displayBinWidth, Top)
			{
				TotalLand = TotalLand,
				OverflowCount = OverflowCount,
				MinLand = MinLand,
				MaxLand = MaxLand,
				FieldsPooled = FieldsPooled,
			};
			for (int i = 0; i < _counts.Length; i++)
			{
				if (_counts[i] == 0) continue;
				// Bin centre, so a sample does not systematically bias to the low edge.
				int j = (int)((BinLow(i) + BinWidth * 0.5f - SeaLevel) / displayBinWidth);
				if (j >= outH._counts.Length) outH.OverflowCount += _counts[i];
				else outH._counts[j] += _counts[i];
			}
			return outH;
		}

		public override string ToString() =>
			$"LandHistogram({TotalLand:N0} land samples from {FieldsPooled} field(s), " +
			$"[{MinLand:F4} .. {MaxLand:F4}] raw, bin {BinWidth:G3}, overflow {OverflowCount})";
	}
}