/*
This Source Code Form is subject to the
terms of the Mozilla Public License, v.
2.0. If a copy of the MPL was not
distributed with this file, You can
obtain one at
http://mozilla.org/MPL/2.0/.
*/

#include "stdafx.h"
#include "editor/editorProfile.hpp"
#include "editor/editorFormat.hpp"
#include "utilities/translation.h"

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <limits>
#include <optional>
#include <tuple>

namespace profile
{

namespace
{

// a switch, and the timbers it shares, may lie on a vertical curve: a concave one, or a convex one no sharper than this
double constexpr switch_crest{2000.0}; // m
// per mille: the change of the grade on a switch taken for a break of it, 3 mm over 30 m of it; below it the switch lies in one plane
double constexpr switch_break{0.1};

double polyline(line const &Line, double const Chainage)
{
	auto const &points{Line.points};
	if (points.empty())
		return 0.0;
	if (points.size() == 1 || Chainage <= points.front().chainage)
		return points.front().elevation + (points.size() > 1 ? grade_after(Line, 0) * (Chainage - points.front().chainage) : 0.0);
	if (Chainage >= points.back().chainage)
		return points.back().elevation + grade_after(Line, points.size() - 2) * (Chainage - points.back().chainage);
	auto const next{std::upper_bound(points.begin(), points.end(), Chainage, [](double const Value, pvi const &Point) { return Value < Point.chainage; })};
	auto const index{static_cast<std::size_t>(std::distance(points.begin(), next)) - 1};
	return points[index].elevation + grade_after(Line, index) * (Chainage - points[index].chainage);
}

bool overlaps(double const From, double const To, context::zone const &Zone)
{
	return From < Zone.to - 1e-6 && To > Zone.from + 1e-6;
}

} // namespace

double grade_after(line const &Line, std::size_t const Index)
{
	auto const &points{Line.points};
	if (Index + 1 >= points.size())
		return 0.0;
	auto const run{points[Index + 1].chainage - points[Index].chainage};
	return run > 1e-9 ? (points[Index + 1].elevation - points[Index].elevation) / run : 0.0;
}

double required_radius(line const &Line)
{
	return std::max({10.0, Line.norms.radius_min, Line.norms.radius_factor * Line.speed * Line.speed});
}

void assign_radii(line &Line)
{
	auto const required{std::ceil(required_radius(Line) / 10.0) * 10.0};
	for (std::size_t k = 1; k + 1 < Line.points.size(); ++k)
	{
		auto &point{Line.points[k]};
		if (false == point.automatic || point.joint)
			continue;
		auto const change{std::abs(grade_after(Line, k) - grade_after(Line, k - 1)) * 1000.0};
		point.radius = change > Line.norms.curve_threshold + 1e-6 ? required : 0.0;
	}
}

void complete(line &Line, context const &Context)
{
	auto &points{Line.points};
	points.erase(std::remove_if(points.begin(), points.end(), [](pvi const &Point) { return Point.joint; }), points.end());
	if (points.size() < 2)
		return;
	bool const start{Line.join_start && Context.start_joined};
	bool const end{Line.join_end && Context.end_joined};
	if (start)
		points.front().elevation = Context.start_elevation;
	if (end)
		points.back().elevation = Context.end_elevation;

	// the vertical curve of the joint starts at the end of the route, off the adjoining track, or past the switches and the shared
	// timbers the end lies on: over them the grade of the adjoining track goes on unchanged
	auto const joint = [&](pvi const &End, pvi const &Neighbour, double const Grade, double const Direction) -> std::optional<pvi> {
		auto lead{0.0};
		for (bool grown{true}; grown;)
		{
			grown = false;
			for (auto const &zone : Context.switches)
			{
				auto const to_start{Direction > 0.0 ? zone.from - End.chainage : End.chainage - zone.to};
				auto const to_end{Direction > 0.0 ? zone.to - End.chainage : End.chainage - zone.from};
				if (to_start <= lead + 0.5 && to_end > lead)
				{
					lead = to_end;
					grown = true;
				}
			}
		}
		auto const run{(Neighbour.chainage - End.chainage) * Direction - lead};
		if (run < 1.0)
			return std::nullopt;
		auto const offset{std::abs(Neighbour.elevation - (End.elevation + Grade * (lead + run) * Direction))};
		if (offset / run * 1000.0 <= Line.norms.curve_threshold + 1e-6)
		{
			// no curve needed, but the grade breaks only past the timbers the end lies on
			if (lead < 0.5)
				return std::nullopt;
			pvi result;
			result.chainage = End.chainage + lead * Direction;
			result.elevation = End.elevation + Grade * lead * Direction;
			result.automatic = false;
			result.joint = true;
			return result;
		}
		auto radius{std::ceil(required_radius(Line) / 10.0) * 10.0};
		double tangent;
		if (run * run >= 2.0 * radius * offset)
			tangent = (run - std::sqrt(run * run - 2.0 * radius * offset)) * 0.5;
		else
		{
			tangent = run * 0.5;
			radius = run * run / (2.0 * offset);
		}
		pvi result;
		result.chainage = End.chainage + (lead + tangent) * Direction;
		result.elevation = End.elevation + Grade * (lead + tangent) * Direction;
		result.radius = radius;
		result.automatic = false;
		result.joint = true;
		return result;
	};
	if (start)
		if (auto const point{joint(points.front(), points[1], Context.start_grade, 1.0)})
			points.insert(points.begin() + 1, *point);
	if (end)
		if (auto const point{joint(points.back(), points[points.size() - 2], Context.end_grade, -1.0)})
			points.insert(points.end() - 1, *point);
	assign_radii(Line);
}

std::vector<curve> curves(line const &Line)
{
	std::vector<curve> result;
	auto const &points{Line.points};
	for (std::size_t k = 1; k + 1 < points.size(); ++k)
	{
		auto const radius{points[k].radius};
		if (radius <= 0.0)
			continue;
		curve c;
		c.point = k;
		c.change = grade_after(Line, k) - grade_after(Line, k - 1);
		c.tangent = radius * std::abs(c.change) * 0.5;
		if (c.tangent < 1e-6)
			continue;
		c.sagitta = c.tangent * c.tangent / (2.0 * radius);
		c.start = points[k].chainage - c.tangent;
		c.end = points[k].chainage + c.tangent;
		result.push_back(c);
	}
	return result;
}

double elevation(line const &Line, double const Chainage)
{
	auto height{polyline(Line, Chainage)};
	for (auto const &c : curves(Line))
	{
		auto const distance{std::abs(Chainage - Line.points[c.point].chainage)};
		if (distance >= c.tangent)
			continue;
		height += c.change * (c.tangent - distance) * (c.tangent - distance) / (4.0 * c.tangent);
	}
	return height;
}

double grade(line const &Line, double const Chainage)
{
	auto const &points{Line.points};
	if (points.size() < 2)
		return 0.0;
	std::size_t index{0};
	while (index + 2 < points.size() && Chainage >= points[index + 1].chainage)
		++index;
	auto slope{grade_after(Line, index)};
	for (auto const &c : curves(Line))
	{
		auto const offset{Chainage - points[c.point].chainage};
		if (std::abs(offset) >= c.tangent)
			continue;
		slope -= c.change * (c.tangent - std::abs(offset)) / (2.0 * c.tangent) * (offset < 0.0 ? -1.0 : 1.0);
	}
	return slope;
}

double radius_for_sagitta(line const &Line, std::size_t const Point, double const Sagitta)
{
	if (Point == 0 || Point + 1 >= Line.points.size())
		return 0.0;
	auto const change{std::abs(grade_after(Line, Point) - grade_after(Line, Point - 1))};
	if (change < 1e-9)
		return Line.points[Point].radius;
	return std::max(0.0, 8.0 * std::abs(Sagitta) / (change * change));
}

std::vector<issue> check(line const &Line, context const &Context)
{
	std::vector<issue> result;
	auto const &points{Line.points};
	auto const &norms{Line.norms};
	auto const add = [&](bool const Error, int const Point, double const Chainage, std::string Text) {
		result.push_back({Error, Point, Chainage, std::move(Text)});
	};
	if (points.size() < 2)
	{
		add(true, -1, 0.0, STR_C("The grade line needs at least the two ends"));
		return result;
	}
	for (std::size_t k = 0; k + 1 < points.size(); ++k)
	{
		if (points[k + 1].chainage - points[k].chainage < 0.5)
			add(true, static_cast<int>(k + 1), points[k + 1].chainage, format(STR_C("Points %zu and %zu are closer than 0.5 m"), k + 1, k + 2));
		auto const slope{grade_after(Line, k) * 1000.0};
		if (std::abs(slope) > norms.grade_max + 1e-6)
			add(false, static_cast<int>(k), points[k].chainage, format(STR_C("Grade %.2f per mille between points %zu and %zu exceeds %.2f"), slope, k + 1, k + 2, norms.grade_max));
	}

	auto const list{curves(Line)};
	auto const required{required_radius(Line)};
	for (std::size_t k = 1; k + 1 < points.size(); ++k)
	{
		auto const change{std::abs(grade_after(Line, k) - grade_after(Line, k - 1)) * 1000.0};
		auto const radius{points[k].radius};
		if (radius <= 0.0 && change > norms.curve_threshold + 1e-6)
			add(false, static_cast<int>(k), points[k].chainage, format(STR_C("Point %zu: change of grade %.2f per mille needs a vertical curve"), k + 1, change));
		else if (radius > 0.0 && radius < required - 1e-6)
			add(false, static_cast<int>(k), points[k].chainage, format(STR_C("Point %zu: R %.0f m is below the required %.0f m for %.0f km/h"), k + 1, radius, required, Line.speed));
		// a change as small as the grade lines laid one after another meet with, a few millimetres over the length of a switch, is no
		// break of it: the tracks of a slip met so lie in one plane all the same (pila.scn: 0.05 per mille inside shared timbers)
		if (radius <= 0.0 && change > switch_break)
			for (auto const &zone : Context.switches)
				// a break on the bound of a switch is let be, within what the chainages of routes laid over each other differ by
				if (points[k].chainage > zone.from + 0.05 && points[k].chainage < zone.to - 0.05)
					add(true, static_cast<int>(k), points[k].chainage, format(STR_C("Point %zu: the grade changes inside switch %s"), k + 1, zone.name.c_str()));
	}
	for (std::size_t i = 0; i < list.size(); ++i)
	{
		auto const &c{list[i]};
		auto const point{static_cast<int>(c.point)};
		if (c.start < points.front().chainage - 1e-6 || c.end > points.back().chainage + 1e-6)
			add(true, point, points[c.point].chainage, format(STR_C("Point %zu: the vertical curve reaches past the end of the route"), c.point + 1));
		if (i + 1 < list.size() && c.end > list[i + 1].start + 1e-6)
			add(true, point, c.end, format(STR_C("Vertical curves of points %zu and %zu overlap by %.2f m"), c.point + 1, list[i + 1].point + 1, c.end - list[i + 1].start));
		// a switch may lie on a concave curve, and on a convex one no sharper than switch_crest, though better off them; within
		// what the fit rounds the radius to its room by
		for (auto const &zone : Context.switches)
			if (overlaps(c.start, c.end, zone))
			{
				if (c.change < 0.0 && points[c.point].radius < switch_crest * 0.99)
					add(true, point, points[c.point].chainage, format(STR_C("Point %zu: the convex vertical curve on switch %s is sharper than %.0f m"), c.point + 1, zone.name.c_str(), switch_crest));
				else
					add(false, point, points[c.point].chainage, format(STR_C("Point %zu: the vertical curve lies on switch %s"), c.point + 1, zone.name.c_str()));
				break;
			}
	}
	if (norms.element_min > 0.0)
	{
		for (std::size_t k = 1; k + 2 < points.size(); ++k)
		{
			auto const tangent = [&](std::size_t const Point) {
				for (auto const &c : list)
					if (c.point == Point)
						return c.tangent;
				return 0.0;
			};
			auto const free{(points[k + 1].chainage - tangent(k + 1)) - (points[k].chainage + tangent(k))};
			if (free < norms.element_min - 1e-6)
				add(false, static_cast<int>(k), points[k].chainage, format(STR_C("Constant grade between points %zu and %zu is %.1f m long, less than %.0f m"), k + 1, k + 2, std::max(0.0, free), norms.element_min));
		}
	}
	for (auto const &zone : Context.switches)
	{
		auto const slope{std::abs(grade(Line, (zone.from + zone.to) * 0.5)) * 1000.0};
		if (norms.grade_switch_max > 0.0 && slope > norms.grade_switch_max + 1e-6)
			add(false, -1, zone.from, format(STR_C("Grade %.2f per mille over switch %s exceeds %.2f"), slope, zone.name.c_str(), norms.grade_switch_max));
	}
	for (auto const &fixed : Context.fixed)
	{
		auto const departure{elevation(Line, fixed.first) - fixed.second};
		if (std::abs(departure) > norms.fixed_tolerance)
		{
			add(true, -1, fixed.first, format(STR_C("At %s the grade line departs %.3f m from a fixed point"), format_chainage(fixed.first).c_str(), departure));
			break;
		}
	}
	auto const ends = [&](bool const Atend) {
		auto const joined{Atend ? Context.end_joined : Context.start_joined};
		if (false == joined)
			return;
		auto const chainage{Atend ? points.back().chainage : points.front().chainage};
		auto const departure{elevation(Line, chainage) - (Atend ? Context.end_elevation : Context.start_elevation)};
		if (std::abs(departure) > norms.fixed_tolerance)
			add(false, Atend ? static_cast<int>(points.size()) - 1 : 0, chainage, format(Atend ? STR_C("The end of the grade line is %.3f m off the adjoining track") : STR_C("The start of the grade line is %.3f m off the adjoining track"), departure));
		auto const slope{grade(Line, chainage)};
		auto const change{std::abs(slope - (Atend ? Context.end_grade : Context.start_grade)) * 1000.0};
		if (change > norms.curve_threshold + 1e-6)
			add(false, Atend ? static_cast<int>(points.size()) - 1 : 0, chainage, format(Atend ? STR_C("Change of grade %.2f per mille against the adjoining track at the end") : STR_C("Change of grade %.2f per mille against the adjoining track at the start"), change));
	};
	ends(false);
	ends(true);
	std::stable_sort(result.begin(), result.end(), [](issue const &A, issue const &B) { return A.error && false == B.error; });
	return result;
}

line recognize(std::vector<double> const &Chainages, std::vector<double> const &Elevations, std::vector<double> const &Grades, recognition const &Options)
{
	line result;
	auto const count{std::min({Chainages.size(), Elevations.size(), Grades.size()})};
	if (count < 2)
		return result;

	struct part
	{
		std::size_t first, last;
		double slope{0.0}, offset{0.0};
	};
	auto const fit = [&](part &Part) {
		double sx{0.0}, sy{0.0}, sxx{0.0}, sxy{0.0};
		auto const origin{Chainages[Part.first]};
		auto const n{static_cast<double>(Part.last - Part.first + 1)};
		for (auto i = Part.first; i <= Part.last; ++i)
		{
			auto const x{Chainages[i] - origin};
			sx += x;
			sy += Elevations[i];
			sxx += x * x;
			sxy += x * Elevations[i];
		}
		auto const denominator{n * sxx - sx * sx};
		Part.slope = std::abs(denominator) > 1e-12 ? (n * sxy - sx * sy) / denominator : 0.0;
		Part.offset = (sy - Part.slope * sx) / n - Part.slope * origin;
	};

	std::vector<part> parts;
	for (std::size_t i = 0; i < count;)
	{
		auto j{i};
		auto low{Grades[i]}, high{Grades[i]};
		while (j + 1 < count && std::max(high, Grades[j + 1]) - std::min(low, Grades[j + 1]) <= Options.grade_tolerance)
		{
			++j;
			low = std::min(low, Grades[j]);
			high = std::max(high, Grades[j]);
		}
		if (Chainages[j] - Chainages[i] >= Options.minimum_length)
		{
			part p{i, j};
			fit(p);
			parts.push_back(p);
			i = j + 1;
		}
		else
			++i;
	}
	auto const apart = [&](part const &A, part const &B, double const Chainage) { return std::abs((A.slope - B.slope) * Chainage + A.offset - B.offset); };
	for (std::size_t i = 0; i + 1 < parts.size();)
	{
		if (std::abs(parts[i].slope - parts[i + 1].slope) < Options.grade_tolerance * 0.5 && apart(parts[i], parts[i + 1], (Chainages[parts[i].last] + Chainages[parts[i + 1].first]) * 0.5) <= Options.step)
		{
			parts[i].last = parts[i + 1].last;
			fit(parts[i]);
			parts.erase(parts.begin() + i + 1);
		}
		else
			++i;
	}

	result.points.push_back({Chainages.front(), Elevations.front(), 0.0});
	for (std::size_t i = 0; i + 1 < parts.size(); ++i)
	{
		auto const &a{parts[i]};
		auto const &b{parts[i + 1]};
		auto const from{Chainages[a.last]};
		auto const to{Chainages[b.first]};
		auto const change{b.slope - a.slope};
		auto chainage{std::abs(change) > 1e-9 ? (a.offset - b.offset) / change : (from + to) * 0.5};
		chainage = std::clamp(chainage, from, std::max(from, to));
		// the grades don't meet between the parts, as at a stretch lifted with the grade it had: a break where each ends
		if (apart(a, b, chainage) > Options.step)
		{
			for (auto const &[at, slope, offset] : {std::tuple{from, a.slope, a.offset}, std::tuple{to, b.slope, b.offset}})
				if (at - result.points.back().chainage >= 0.5)
				{
					pvi point;
					point.automatic = false;
					point.chainage = at;
					point.elevation = slope * at + offset;
					result.points.push_back(point);
				}
			continue;
		}
		pvi point;
		point.automatic = false;
		point.chainage = chainage;
		point.elevation = a.slope * chainage + a.offset;
		auto const gap{to - from};
		if (gap >= Options.sharp_gap && std::abs(change) > 1e-9)
			point.radius = std::round(gap / std::abs(change));
		if (point.chainage - result.points.back().chainage >= 0.5)
			result.points.push_back(point);
	}
	if (Chainages[count - 1] - result.points.back().chainage < 0.5 && result.points.size() > 1)
		result.points.pop_back();
	result.points.push_back({Chainages[count - 1], Elevations[count - 1], 0.0});

	// the gap between constant grades is shorter than the curve, the sagitta gives the radius
	for (std::size_t i = 1; i + 1 < result.points.size(); ++i)
	{
		auto &point{result.points[i]};
		if (point.radius <= 0.0)
			continue;
		auto const next{std::lower_bound(Chainages.begin(), Chainages.begin() + count, point.chainage)};
		if (next == Chainages.begin() || next == Chainages.begin() + count)
			continue;
		auto const index{static_cast<std::size_t>(std::distance(Chainages.begin(), next))};
		auto const run{Chainages[index] - Chainages[index - 1]};
		auto const measured{run > 1e-9 ? Elevations[index - 1] + (Elevations[index] - Elevations[index - 1]) * (point.chainage - Chainages[index - 1]) / run : Elevations[index]};
		auto const radius{radius_for_sagitta(result, i, measured - point.elevation)};
		if (radius > 0.0)
			point.radius = std::round(radius);
	}
	if (Options.detail <= 0.0)
		return result;
	std::vector<char> off(count, 0);
	for (std::size_t i = 0; i < count; ++i)
		off[i] = std::abs(elevation(result, Chainages[i]) - Elevations[i]) > Options.detail ? 1 : 0;
	std::vector<pvi> bends;
	for (std::size_t i = 0; i < count;)
	{
		if (off[i] == 0)
		{
			++i;
			continue;
		}
		auto const first{i > 0 ? i - 1 : i};
		while (i < count && off[i] != 0)
			++i;
		auto const last{std::min(i, count - 1)};
		std::vector<char> keep(count, 0);
		keep[first] = keep[last] = 1;
		std::vector<std::pair<std::size_t, std::size_t>> stack{{first, last}};
		while (false == stack.empty())
		{
			auto const [a, b]{stack.back()};
			stack.pop_back();
			if (b <= a + 1)
				continue;
			auto const run{std::max(1e-9, Chainages[b] - Chainages[a])};
			double worst{0.0};
			std::size_t index{a};
			for (auto k = a + 1; k < b; ++k)
			{
				auto const departure{std::abs(Elevations[k] - (Elevations[a] + (Elevations[b] - Elevations[a]) * (Chainages[k] - Chainages[a]) / run))};
				if (departure > worst)
				{
					worst = departure;
					index = k;
				}
			}
			if (worst > Options.detail * 0.5)
			{
				keep[index] = 1;
				stack.push_back({a, index});
				stack.push_back({index, b});
			}
		}
		auto const low{Chainages[first]}, high{Chainages[last]};
		result.points.erase(std::remove_if(result.points.begin() + 1, result.points.end() - 1, [&](pvi const &Point) { return Point.chainage >= low - 0.5 && Point.chainage <= high + 0.5; }), result.points.end() - 1);
		for (auto k = first; k <= last; ++k)
			if (keep[k] != 0 && k > 0 && k + 1 < count)
			{
				pvi point;
				point.automatic = false;
				point.chainage = Chainages[k];
				point.elevation = Elevations[k];
				bends.push_back(point);
			}
	}
	if (bends.empty())
		return result;
	result.points.insert(result.points.end() - 1, bends.begin(), bends.end());
	std::sort(result.points.begin(), result.points.end(), [](pvi const &Left, pvi const &Right) { return Left.chainage < Right.chainage; });
	std::vector<pvi> points;
	for (auto const &point : result.points)
		if (points.empty() || point.chainage - points.back().chainage >= 0.5 || &point == &result.points.back())
		{
			if (false == points.empty() && point.chainage - points.back().chainage < 0.5)
				points.pop_back();
			points.push_back(point);
		}
	result.points = std::move(points);
	return result;
}

std::vector<hump> humps(std::vector<double> const &Chainages, std::vector<double> const &Ground, context const &Context, double const Tolerance, std::vector<double> const &Tolerances)
{
	std::vector<hump> result;
	auto const count{std::min(Chainages.size(), Ground.size())};
	if (count < 2)
		return result;
	// the ground every few metres, the mean of the model about each, past its noise
	double constexpr spacing{5.0};
	auto const origin{Chainages.front()};
	auto const end{Chainages[count - 1]};
	auto const cells{static_cast<std::ptrdiff_t>(std::lround((end - origin) / spacing)) + 1};
	std::vector<double> level(cells, 0.0);
	std::vector<int> taken(cells, 0);
	for (std::size_t i = 0; i < count; ++i)
	{
		if (false == std::isfinite(Ground[i]) || (i < Tolerances.size() && Tolerances[i] > 0.0 && Tolerances[i] < Tolerance))
			continue;
		auto const k{std::clamp<std::ptrdiff_t>(std::lround((Chainages[i] - origin) / spacing), 0, cells - 1)};
		level[k] += Ground[i];
		++taken[k];
	}
	for (std::ptrdiff_t k = 0; k < cells; ++k)
		level[k] = taken[k] > 0 ? level[k] / taken[k] : std::numeric_limits<double>::quiet_NaN();
	auto const at = [&](std::ptrdiff_t const K) { return K >= 0 && K < cells ? level[K] : std::numeric_limits<double>::quiet_NaN(); };
	auto const cells_of = [](double const Distance) { return static_cast<std::ptrdiff_t>(std::lround(Distance / spacing)); };
	// the crest is the highest of the ground about it, the feet where the ground has fallen nearly as far as it falls within
	// the reach: some metres lower than the crest, not as many as of a hill or a bridge
	double constexpr around{60.0}; // m
	double constexpr farthest{300.0}; // m
	double constexpr lowest{1.0}; // m
	double constexpr highest{6.0}; // m
	// the grades just past the rounding of the crest fall away from it both ways, one of them more steeply than lines and
	// stations are laid, the two together a change of grade over some 50 m which no curve of 1250 m or more makes
	double constexpr inner{10.0}; // m
	double constexpr outer{40.0}; // m
	double constexpr falling{5.0}; // per mille
	double constexpr steep{25.0}; // per mille
	double constexpr sharp{40.0}; // per mille
	// nowhere about the crest as steep as the banks of a bridge, the walls of a cutting
	double constexpr wall{80.0}; // per mille
	// no switch on the crest; past it the fan, switches of their own length, not the shared timbers of the ones beside
	double constexpr clear{15.0}; // m
	double constexpr fan{150.0}; // m
	double constexpr body{15.0}; // m
	std::size_t constexpr fanning{2};
	auto const crest_of = [&](std::ptrdiff_t const C) -> std::optional<hump> {
		auto const top{at(C)};
		if (false == std::isfinite(top))
			return std::nullopt;
		for (auto k = C - cells_of(around); k <= C + cells_of(around); ++k)
			if (auto const other{at(k)}; std::isfinite(other) && (other > top || (other == top && k < C)))
				return std::nullopt;
		hump found;
		found.crest = origin + spacing * static_cast<double>(C);
		for (auto const &zone : Context.switches)
			if (zone.from < found.crest + clear && zone.to > found.crest - clear)
				return std::nullopt;
		double grades{0.0};
		double steepest{0.0};
		bool fanned{false};
		for (int const side : {-1, 1})
		{
			auto low{top};
			std::ptrdiff_t seen{0};
			auto const reach{cells_of(farthest)};
			for (std::ptrdiff_t d = 1; d <= reach; ++d)
				if (auto const height{at(C + side * d)}; std::isfinite(height))
				{
					low = std::min(low, height);
					++seen;
				}
			if (seen < cells_of(around) || top - low < lowest || top - low > highest)
				return std::nullopt;
			found.height = std::max(found.height, top - low);
			auto foot{reach};
			for (std::ptrdiff_t d = 1; d <= reach; ++d)
				if (auto const height{at(C + side * d)}; std::isfinite(height) && top - height >= 0.9 * (top - low))
				{
					foot = d;
					break;
				}
			// and on past the foot as far as the curve into the grade beyond
			(side < 0 ? found.from : found.to) = std::clamp(found.crest + side * spacing * static_cast<double>(foot + cells_of(25.0)), origin, end);
			auto const off{at(C + side * cells_of(inner))};
			auto const on{at(C + side * cells_of(outer))};
			if (false == std::isfinite(off) || false == std::isfinite(on))
				return std::nullopt;
			auto const grade{(off - on) / (outer - inner) * 1000.0};
			if (grade < falling)
				return std::nullopt;
			grades += grade;
			steepest = std::max(steepest, grade);
			for (std::ptrdiff_t d = 0; d + 2 <= cells_of(around); ++d)
				if (auto const a{at(C + side * d)}, b{at(C + side * (d + 2))}; std::isfinite(a) && std::isfinite(b) && std::abs(a - b) / (2.0 * spacing) * 1000.0 > wall)
					return std::nullopt;
			std::size_t past{0};
			for (auto const &zone : Context.switches)
				if (auto const distance{side > 0 ? zone.from - found.crest : found.crest - zone.to}; zone.to - zone.from >= body && distance >= clear && distance <= fan)
					++past;
			fanned = fanned || past >= fanning;
		}
		if (steepest < steep || grades < sharp || false == fanned)
			return std::nullopt;
		return found;
	};
	for (std::ptrdiff_t c = 0; c < cells; ++c)
		if (auto const found{crest_of(c)}; found && (result.empty() || found->crest > result.back().to))
			result.push_back(*found);
	return result;
}

line fit_ground(std::vector<double> const &Chainages, std::vector<double> const &Ground, context const &Context, line const &Settings, ground_fit const &Options, std::vector<double> const &Tolerances,
                std::vector<double> const &Reaches)
{
	line result;
	result.speed = Settings.speed;
	result.norms = Settings.norms;
	result.join_start = Settings.join_start;
	result.join_end = Settings.join_end;
	auto const count{std::min(Chainages.size(), Ground.size())};
	if (count < 2 || std::none_of(Ground.begin(), Ground.begin() + count, [](double const Height) { return std::isfinite(Height); }))
		return result;
	auto const origin{Chainages.front()};
	auto const length{Chainages[count - 1] - origin};
	if (length < 1.0)
		return result;

	// the line is sought on nodes a few metres apart, straight between them: the grade can change only at a node
	auto const nodes{std::max<std::size_t>(3, static_cast<std::size_t>(std::ceil(length / 5.0)) + 1)};
	auto const step{length / static_cast<double>(nodes - 1)};
	std::vector<std::size_t> cell(count);
	std::vector<double> part(count);
	for (std::size_t i = 0; i < count; ++i)
	{
		auto const at{std::clamp((Chainages[i] - origin) / step, 0.0, static_cast<double>(nodes - 1))};
		cell[i] = std::min(nodes - 2, static_cast<std::size_t>(at));
		part[i] = at - static_cast<double>(cell[i]);
	}
	auto const node = [&](std::size_t const K) { return origin + step * static_cast<double>(K); };
	// a change of grade at a node is the curve of the required radius around it. On a hump the least radius is the one of the
	// hump, though a convex curve over a switch is no sharper than switch_crest
	auto const required{std::ceil(required_radius(result) / 10.0) * 10.0};
	auto const on_hump = [&](double const Chainage) { return std::any_of(Options.humps.begin(), Options.humps.end(), [&](hump const &Hump) { return Chainage >= Hump.from && Chainage <= Hump.to; }); };
	auto const least = [&](double const Chainage) { return Options.hump_radius > 0.0 && on_hump(Chainage) ? std::min(required, std::ceil(Options.hump_radius / 10.0) * 10.0) : required; };
	auto const on_switch = [&](double const From, double const To) { return std::any_of(Context.switches.begin(), Context.switches.end(), [&](auto const &Zone) { return Zone.to > From && Zone.from < To; }); };
	auto const sharpest{step / (required * 1.001)};
	std::vector<double> sag_max(nodes), crest_max(nodes);
	for (std::size_t k = 0; k < nodes; ++k)
	{
		auto const radius{least(node(k))};
		sag_max[k] = crest_max[k] = step / (radius * 1.001);
		// a sharper one keeps its curve off the switches
		if (radius < switch_crest)
			for (auto const &zone : Context.switches)
				crest_max[k] = std::min(crest_max[k], std::max(step / (switch_crest * 1.001), 2.0 * std::max(zone.from - node(k), node(k) - zone.to) / (radius * 1.001)));
	}
	struct held
	{
		std::size_t cell;
		double part;
		double low;
		double high;
		bool yielding;
	};
	std::vector<held> bounds;
	for (auto const &bound : Options.bounds)
	{
		auto const at{std::clamp((bound.chainage - origin) / step, 0.0, static_cast<double>(nodes - 1))};
		auto const cell{std::min(nodes - 2, static_cast<std::size_t>(at))};
		bounds.push_back({cell, at - static_cast<double>(cell), std::min(bound.low, bound.high), std::max(bound.low, bound.high), bound.yielding});
	}
	// the bounds which contradict each other about one place, as of two crossovers into lines fitted apart in height, hold as far
	// as both can: the line keeps between them. Held both, it would tear itself apart over them
	std::sort(bounds.begin(), bounds.end(), [](held const &A, held const &B) { return A.cell + A.part < B.cell + B.part; });
	for (std::size_t first = 0; first < bounds.size();)
	{
		auto last{first + 1};
		while (last < bounds.size() && bounds[last].cell + bounds[last].part - bounds[first].cell - bounds[first].part <= 1.0)
			++last;
		auto low{-std::numeric_limits<double>::max()};
		auto high{std::numeric_limits<double>::max()};
		for (auto i = first; i < last; ++i)
		{
			low = std::max(low, bounds[i].low);
			high = std::min(high, bounds[i].high);
		}
		if (low > high)
			for (auto i = first; i < last; ++i)
			{
				bounds[i].low = high;
				bounds[i].high = low;
			}
		first = last;
	}
	bool const start_fixed{Settings.join_start && Context.start_joined};
	bool const end_fixed{Settings.join_end && Context.end_joined};

	// a change of grade costs the price times its size, so the line breaks where the ground keeps away from it long enough to
	// be worth it. The ground costs the square of its distance over its tolerance; past its reach ever less, up to a bound:
	// what parts from the line, a bridge, a viaduct, a wall, a track on another level, is gone past rather than followed. The
	// bends are no sharper than the radius of the norms and the grades no steeper than the limit
	auto const tolerance{std::max(0.01, Options.tolerance)};
	auto const told = [&](std::size_t const I) { return I < Tolerances.size() && Tolerances[I] > 0.0 ? std::max(0.005, Tolerances[I]) : tolerance; };
	// the ground the line follows however far it parts, and as closely as asked, where it parted longer than a structure would
	std::vector<double> followed(count, 0.0);
	auto const close_to{Options.follow_tolerance > 0.0 ? std::max(0.005, Options.follow_tolerance) : tolerance};
	auto const own = [&](std::size_t const I) { return followed[I] > 0.0 ? close_to : told(I); };
	auto const reach = [&](std::size_t const I) { return followed[I] > 0.0 ? followed[I] : I < Reaches.size() && Reaches[I] > 0.0 ? std::max(0.005, Reaches[I]) : told(I); };
	for (auto const &hump : Options.humps)
		for (std::size_t i = 0; i < count; ++i)
			if (std::isfinite(Ground[i]) && told(i) >= tolerance && Chainages[i] >= hump.from && Chainages[i] <= hump.to)
				followed[i] = hump.height + tolerance;
	auto const price{24000.0 * tolerance * tolerance};
	auto const grade_max{(Options.grade_max > 0.0 ? Options.grade_max : result.norms.grade_max) / 1000.0};
	// what has to hold, the heights of the joined ends, holds before the grades of the ends, and these before the limits of the
	// norms. A yielding bound counts as much as a kilometre of the ground off as far, and gives
	// way to the norms
	double constexpr stiff{1e11};
	double constexpr firm{1e9};
	double constexpr leaning{1e3};

	std::vector<double> x(nodes);
	{
		std::size_t i{0};
		for (std::size_t k = 0; k < nodes; ++k)
		{
			while (i + 1 < count && Chainages[i] < node(k))
				++i;
			auto j{i};
			while (j > 0 && false == std::isfinite(Ground[j]))
				--j;
			while (j + 1 < count && false == std::isfinite(Ground[j]))
				++j;
			x[k] = Ground[j];
		}
	}
	// the cost of a line, with its slope and the band of its curvature when asked. The distance of the ground counts up to
	// its tolerance times the root of the spread, the size of a change of grade is rounded off within the flat
	auto spread{256.0};
	auto flat{sharpest};
	std::vector<double> gradient(nodes), diagonal(nodes), first(nodes), second(nodes);
	auto const cost = [&](std::vector<double> const &Line, bool const Derivatives) {
		if (Derivatives)
		{
			std::fill(gradient.begin(), gradient.end(), 0.0);
			std::fill(diagonal.begin(), diagonal.end(), 1e-9);
			std::fill(first.begin(), first.end(), 0.0);
			std::fill(second.begin(), second.end(), 0.0);
		}
		auto total{0.0};
		auto const add = [&](std::initializer_list<std::pair<std::size_t, double>> Terms, double const Value, double const Slope, double const Curvature) {
			total += Value;
			if (false == Derivatives)
				return;
			for (auto const &[a, ca] : Terms)
			{
				gradient[a] += Slope * ca;
				for (auto const &[b, cb] : Terms)
				{
					if (b == a)
						diagonal[a] += Curvature * ca * cb;
					else if (b == a + 1)
						first[a] += Curvature * ca * cb;
					else if (b == a + 2)
						second[a] += Curvature * ca * cb;
				}
			}
		};
		auto const square = [&](std::initializer_list<std::pair<std::size_t, double>> Terms, double const Off, double const Weight) { add(Terms, Weight * Off * Off, 2.0 * Weight * Off, 2.0 * Weight); };
		for (std::size_t i = 0; i < count; ++i)
		{
			if (false == std::isfinite(Ground[i]))
				continue;
			auto const off{Line[cell[i]] * (1.0 - part[i]) + Line[cell[i] + 1] * part[i] - Ground[i]};
			auto const agreeing{spread * reach(i) * reach(i)};
			auto const share{agreeing / (agreeing + off * off)};
			auto const weight{tolerance * tolerance / (own(i) * own(i))};
			add({{cell[i], 1.0 - part[i]}, {cell[i] + 1, part[i]}}, weight * off * off * share, weight * 2.0 * off * share * share, weight * 2.0 * share * share);
		}
		for (std::size_t k = 1; k + 1 < nodes; ++k)
		{
			auto const terms = {std::pair{k - 1, 1.0 / step}, std::pair{k, -2.0 / step}, std::pair{k + 1, 1.0 / step}};
			auto const bend{(Line[k - 1] - 2.0 * Line[k] + Line[k + 1]) / step};
			auto const root{std::sqrt(bend * bend + flat * flat)};
			add(terms, price * root, price * bend / root, price * flat * flat / (root * root * root));
			if (auto const most{bend > 0.0 ? sag_max[k] : crest_max[k]}; std::abs(bend) > most)
				square(terms, bend - std::copysign(most, bend), firm);
		}
		for (std::size_t k = 0; k + 1 < nodes; ++k)
		{
			auto const grade{(Line[k + 1] - Line[k]) / step};
			if (std::abs(grade) > grade_max)
				square({{k, -1.0 / step}, {k + 1, 1.0 / step}}, grade - std::copysign(grade_max, grade), firm);
		}
		if (start_fixed)
		{
			square({{0, 1.0}}, Line[0] - Context.start_elevation, stiff);
			if (Options.hold_grades)
				square({{0, -1.0 / step}, {1, 1.0 / step}}, (Line[1] - Line[0]) / step - Context.start_grade, firm * 10.0);
		}
		if (end_fixed)
		{
			square({{nodes - 1, 1.0}}, Line[nodes - 1] - Context.end_elevation, stiff);
			if (Options.hold_grades)
				square({{nodes - 2, -1.0 / step}, {nodes - 1, 1.0 / step}}, (Line[nodes - 1] - Line[nodes - 2]) / step - Context.end_grade, firm * 10.0);
		}
		for (auto const &bound : bounds)
		{
			auto const height{Line[bound.cell] * (1.0 - bound.part) + Line[bound.cell + 1] * bound.part};
			if (height < bound.low || height > bound.high)
				square({{bound.cell, 1.0 - bound.part}, {bound.cell + 1, bound.part}}, height - std::clamp(height, bound.low, bound.high), bound.yielding ? leaning : firm * 10.0);
		}
		return total;
	};
	// the step of Newton, the band factored as L D L^T with two diagonals below the main one
	auto const newton = [&]() {
		std::vector<double> d(nodes), l1(nodes, 0.0), l2(nodes, 0.0), z(nodes);
		for (std::size_t i = 0; i < nodes; ++i)
		{
			d[i] = diagonal[i] - (i >= 1 ? l1[i - 1] * l1[i - 1] * d[i - 1] : 0.0) - (i >= 2 ? l2[i - 2] * l2[i - 2] * d[i - 2] : 0.0);
			if (i + 1 < nodes)
				l1[i] = (first[i] - (i >= 1 ? l2[i - 1] * l1[i - 1] * d[i - 1] : 0.0)) / d[i];
			if (i + 2 < nodes)
				l2[i] = second[i] / d[i];
		}
		for (std::size_t i = 0; i < nodes; ++i)
			z[i] = -gradient[i] - (i >= 1 ? l1[i - 1] * z[i - 1] : 0.0) - (i >= 2 ? l2[i - 2] * z[i - 2] : 0.0);
		for (std::size_t i = 0; i < nodes; ++i)
			z[i] /= d[i];
		for (auto i = nodes - 1; i-- > 0;)
			z[i] -= l1[i] * z[i + 1] + (i + 2 < nodes ? l2[i] * z[i + 2] : 0.0);
		return z;
	};
	// the spread and the flat narrow down level by level, each level going on from the line of the one before: from nearly
	// all the ground agreeing and changes of grade costing their square, to the tolerance and their size
	auto const solve = [&]() {
		spread = 256.0;
		flat = sharpest;
		for (;;)
		{
			auto current{cost(x, true)};
			for (int pass = 0; pass < 100; ++pass)
			{
				auto const direction{newton()};
				auto descent{0.0};
				for (std::size_t k = 0; k < nodes; ++k)
					descent += gradient[k] * direction[k];
				std::vector<double> trial(nodes);
				auto length{1.0};
				for (;; length *= 0.5)
				{
					for (std::size_t k = 0; k < nodes; ++k)
						trial[k] = x[k] + length * direction[k];
					if (cost(trial, false) <= current + 1e-4 * length * descent || length < 1e-6)
						break;
				}
				auto moved{0.0};
				for (std::size_t k = 0; k < nodes; ++k)
					moved = std::max(moved, std::abs(length * direction[k]));
				x = std::move(trial);
				current = cost(x, true);
				if (moved < 1e-5)
					break;
			}
			if (spread <= 1.0 && flat <= 1e-6)
				break;
			spread = std::max(1.0, spread * 0.5);
			flat = std::max(1e-6, flat * 0.5);
		}
	};
	solve();
	// the ground parted from the line farther along than a structure is long, not as far as a viaduct or a hill over a tunnel
	// keep it, and nowhere as steep as the banks and the walls under a bridge, is none: a hump, a rise, a hollow the line goes
	// over. It agrees with the line however far it is, and the line is sought again. The ground which tells the height no more
	// closely than the tolerance counts, the model, not the rail tops beside; its steepness over a few metres, past the noise
	// of the model
	double constexpr chord{5.0};
	for (int round = 0; round < 3 && Options.follow_length > 0.0; ++round)
	{
		bool more{false};
		std::vector<std::size_t> run;
		std::size_t back{0};
		std::size_t previous{count};
		int side{0};
		double farthest{0.0};
		double steepest{0.0};
		auto const close = [&]() {
			if (side != 0 && Chainages[run.back()] - Chainages[run.front()] > Options.follow_length && farthest <= Options.follow_height && steepest <= Options.follow_steepest / 1000.0)
			{
				for (auto const i : run)
					if (followed[i] <= 0.0)
						followed[i] = farthest + tolerance;
				more = true;
			}
			side = 0;
			run.clear();
		};
		auto const extend = [&](std::size_t const I) {
			while (back + 1 < run.size() && Chainages[I] - Chainages[run[back + 1]] >= chord)
				++back;
			if (false == run.empty() && Chainages[I] - Chainages[run[back]] >= chord)
				steepest = std::max(steepest, std::abs(Ground[I] - Ground[run[back]]) / (Chainages[I] - Chainages[run[back]]));
			run.push_back(I);
		};
		for (std::size_t i = 0; i < count; ++i)
		{
			if (false == std::isfinite(Ground[i]) || told(i) < tolerance)
				continue;
			auto const off{x[cell[i]] * (1.0 - part[i]) + x[cell[i] + 1] * part[i] - Ground[i]};
			auto const now{followed[i] > 0.0 ? 0 : off > reach(i) ? 1 : off < -reach(i) ? -1 : 0};
			if (now != side)
			{
				close();
				if (now != 0)
				{
					side = now;
					back = 0;
					farthest = 0.0;
					steepest = 0.0;
					if (previous < count)
						run.push_back(previous);
				}
			}
			if (now != 0)
			{
				farthest = std::max(farthest, std::abs(off));
				extend(i);
			}
			previous = i;
		}
		close();
		if (false == more)
			break;
		solve();
	}

	// the design: the changes of grade along a bend, one after another the same way, are one point of intersection where
	// the grades before and after them meet. No change at a node is sharper than the required radius, so the curve of that
	// radius stays within the nodes of its bend; where the ends held the line sharper, the curve is as sharp as the line
	auto const change = [&](std::size_t const K) { return (x[K - 1] - 2.0 * x[K] + x[K + 1]) / step; };
	auto const significant{sharpest * 0.01};
	auto const bends = [&](std::size_t const K) { return std::abs(change(K)) >= significant; };
	// an end whose grade the line couldn't take over keeps only its height, the grade breaks at the joint
	if (start_fixed && std::abs((x[1] - x[0]) / step - Context.start_grade) > significant)
		result.join_start = false;
	if (end_fixed && std::abs((x[nodes - 1] - x[nodes - 2]) / step - Context.end_grade) > significant)
		result.join_end = false;
	result.points.push_back({origin, x[0]});
	std::vector<double> rooms{0.0};
	std::vector<double> lengths{0.0};
	for (std::size_t k = 1; k + 1 < nodes;)
	{
		if (false == bends(k))
		{
			++k;
			continue;
		}
		auto const first{k};
		auto const sign{change(k) > 0.0};
		double sum{0.0}, moment{0.0};
		for (; k + 1 < nodes && bends(k) && (change(k) > 0.0) == sign; ++k)
		{
			sum += change(k);
			moment += change(k) * node(k);
		}
		auto const at{moment / sum};
		pvi point{at, x[first] + (x[first] - x[first - 1]) / step * (at - node(first))};
		auto const room{std::min(at - node(first) + step * 0.5, node(k - 1) + step * 0.5 - at)};
		result.points.push_back(point);
		rooms.push_back(room);
		lengths.push_back(node(k - 1) - node(first) + step);
	}
	result.points.push_back({node(nodes - 1), x[nodes - 1]});
	rooms.push_back(0.0);
	lengths.push_back(0.0);
	// the curve of the required radius where it has room, off the curves beside; where it hasn't, as sharp as the room lets it
	// be. On a hump as long as the bend of the line, no sharper than the radius of the hump, nor than a switch lets a convex one
	for (std::size_t i = 1; i + 1 < result.points.size(); ++i)
	{
		auto const turn{std::abs(grade_after(result, i) - grade_after(result, i - 1))};
		if (turn * 1000.0 <= result.norms.curve_threshold + 1e-6)
		{
			// a change no curve needs breaks bare, but not on a switch
			if (auto const at{result.points[i].chainage}; turn > 1e-9 && on_switch(at, at))
			{
				result.points[i].radius = std::min(required, 2.0 * std::max(0.0, rooms[i] - 1e-4) / turn);
				result.points[i].automatic = false;
			}
			continue;
		}
		auto const radius{least(result.points[i].chainage)};
		auto bent{radius < required ? std::clamp(std::ceil(lengths[i] / turn / 10.0) * 10.0, radius, required) : required};
		if (auto const at{result.points[i].chainage}; bent < switch_crest && grade_after(result, i) < grade_after(result, i - 1) && on_switch(at - bent * turn * 0.5, at + bent * turn * 0.5))
			bent = switch_crest;
		if (bent * turn * 0.5 > rooms[i] - 1e-4)
		{
			result.points[i].radius = 2.0 * std::max(0.0, rooms[i] - 1e-4) / turn;
			result.points[i].automatic = false;
		}
		else if (bent < required)
		{
			result.points[i].radius = bent;
			result.points[i].automatic = false;
		}
	}
	assign_radii(result);

	// then the simplest design which keeps within half the tolerance of the line found, the one of the ground agreeing with
	// it: a point is dropped, or two bending the same way become one, as long as that departs the least and adds no fault
	std::vector<double> limit(nodes, tolerance);
	for (std::size_t i = 0; i < count; ++i)
		if (std::isfinite(Ground[i]) && std::abs(Ground[i] - (x[cell[i]] * (1.0 - part[i]) + x[cell[i] + 1] * part[i])) <= reach(i))
			for (auto const k : {cell[i], cell[i] + 1})
				limit[k] = std::min(limit[k], own(i));
	for (auto const &bound : bounds)
		for (auto const k : {bound.cell, bound.cell + 1})
			limit[k] = std::min(limit[k], std::max(0.001, (bound.high - bound.low) * 0.5));
	auto const errors = [&](line const &Line) {
		auto completed{Line};
		complete(completed, Context);
		auto const found{check(completed, Context)};
		return std::count_if(found.begin(), found.end(), [](issue const &Issue) { return Issue.error; });
	};
	auto const departure = [&](line const &Line, double const From, double const To) {
		auto worst{0.0};
		for (std::size_t k = 0; k < nodes; ++k)
			if (node(k) >= From && node(k) <= To)
				worst = std::max(worst, std::abs(elevation(Line, node(k)) - x[k]) / limit[k]);
		return worst;
	};
	auto faults{errors(result)};
	// a convex curve over a switch which the room of its bend cut sharper than a switch may lie on takes the radius it may have,
	// as far as the curves beside leave it room
	for (std::size_t i = 1; faults > 0 && i + 1 < result.points.size(); ++i)
	{
		auto const turn{grade_after(result, i) - grade_after(result, i - 1)};
		auto const radius{result.points[i].radius};
		auto const at{result.points[i].chainage};
		if (radius <= 0.0 || radius >= switch_crest || turn >= 0.0 || false == on_switch(at + radius * turn * 0.5, at - radius * turn * 0.5))
			continue;
		auto wider{result};
		wider.points[i].radius = switch_crest;
		wider.points[i].automatic = false;
		if (auto const fewer{errors(wider)}; fewer < faults)
		{
			result = std::move(wider);
			faults = fewer;
		}
	}
	// an end whose grade the curve of the joint can't take over without a fault breaks at the joint
	for (auto *join : {&result.join_start, &result.join_end})
		if (faults > 0 && *join)
		{
			*join = false;
			if (auto const fewer{errors(result)}; fewer < faults)
				faults = fewer;
			else
				*join = true;
		}
	// where the line can't bend from the grades of its ends into the ground within the norms, as a short piece between two
	// lines, the grades break at the joints: the line keeps only the heights of the ends
	if (faults > 0 && Options.hold_grades && (start_fixed || end_fixed))
	{
		auto loose{Options};
		loose.hold_grades = false;
		if (auto again{fit_ground(Chainages, Ground, Context, Settings, loose, Tolerances, Reaches)}; errors(again) < faults)
			return again;
	}
	for (;;)
	{
		auto const &points{result.points};
		auto const last{points.size() - 1};
		std::vector<std::pair<double, line>> candidates;
		for (std::size_t i = 1; i < last; ++i)
		{
			auto const from{points[i >= 2 ? i - 2 : 0].chainage};
			auto dropped{result};
			dropped.points.erase(dropped.points.begin() + i);
			assign_radii(dropped);
			candidates.emplace_back(departure(dropped, from, points[std::min(last, i + 2)].chainage), std::move(dropped));
			if (i + 1 >= last)
				continue;
			auto const before{grade_after(result, i - 1)};
			auto const between{grade_after(result, i)};
			auto const after{grade_after(result, i + 1)};
			if ((between - before) * (after - between) <= 0.0)
				continue;
			auto const at{(points[i + 1].elevation - points[i].elevation + before * points[i].chainage - after * points[i + 1].chainage) / (before - after)};
			auto merged{result};
			merged.points[i] = {at, points[i].elevation + before * (at - points[i].chainage)};
			merged.points.erase(merged.points.begin() + i + 1);
			assign_radii(merged);
			candidates.emplace_back(departure(merged, from, points[std::min(last, i + 3)].chainage), std::move(merged));
		}
		std::sort(candidates.begin(), candidates.end(), [](auto const &A, auto const &B) { return A.first < B.first; });
		auto const chosen{std::find_if(candidates.begin(), candidates.end(), [&](auto const &Candidate) { return Candidate.first <= 0.5 && errors(Candidate.second) <= faults; })};
		if (chosen == candidates.end())
			break;
		result = std::move(chosen->second);
		faults = errors(result);
	}
	return result;
}

std::vector<departure> departures(line const &Line, std::vector<double> const &Chainages, std::vector<double> const &Ground, double const Tolerance)
{
	std::vector<departure> result;
	auto const count{std::min(Chainages.size(), Ground.size())};
	auto const off = [&](std::size_t const I) { return std::isfinite(Ground[I]) ? Ground[I] - elevation(Line, Chainages[I]) : 0.0; };
	for (std::size_t i = 0; i < count;)
	{
		if (std::abs(off(i)) <= Tolerance)
		{
			++i;
			continue;
		}
		departure stretch;
		stretch.from = Chainages[i];
		auto const below{off(i) < 0.0};
		for (; i < count && std::abs(off(i)) > Tolerance && (off(i) < 0.0) == below; ++i)
		{
			if (std::abs(off(i)) > std::abs(stretch.depth))
				stretch.depth = off(i);
			stretch.to = Chainages[i];
		}
		result.push_back(stretch);
	}
	return result;
}

std::vector<departure> bridged(std::vector<double> const &Chainages, std::vector<double> const &Heights, std::vector<double> const &Ground, double const Depth, double const Steepest)
{
	std::vector<departure> result;
	auto const count{std::min({Chainages.size(), Heights.size(), Ground.size()})};
	auto const under = [&](std::size_t const I) { return std::isfinite(Ground[I]) && Heights[I] - Ground[I] > Depth; };
	// the steepness over a few metres, past the noise of the model, in and about the stretch
	double constexpr chord{5.0};
	double constexpr about{10.0};
	auto const steepest = [&](double const From, double const To) {
		double most{0.0};
		for (std::size_t i = 0, j = 0; i < count; ++i)
		{
			if (Chainages[i] < From - about || Chainages[i] > To + about || false == std::isfinite(Ground[i]))
				continue;
			j = std::max(j, i);
			while (j + 1 < count && Chainages[j] - Chainages[i] < chord)
				++j;
			if (auto const run{Chainages[j] - Chainages[i]}; run >= chord && run <= 2.0 * chord && std::isfinite(Ground[j]))
				most = std::max(most, std::abs(Ground[j] - Ground[i]) / run);
		}
		return most * 1000.0;
	};
	for (std::size_t i = 0; i < count;)
	{
		if (false == under(i))
		{
			++i;
			continue;
		}
		departure stretch;
		stretch.from = Chainages[i];
		for (; i < count && under(i); ++i)
		{
			stretch.depth = std::min(stretch.depth, Ground[i] - Heights[i]);
			stretch.to = Chainages[i];
		}
		if (steepest(stretch.from, stretch.to) > Steepest)
			result.push_back(stretch);
	}
	return result;
}

std::vector<double> breaks(line const &Line)
{
	std::vector<double> result;
	auto const &points{Line.points};
	for (std::size_t k = 1; k + 1 < points.size(); ++k)
	{
		if (points[k].radius > 0.0)
			continue;
		if (std::abs(grade_after(Line, k) - grade_after(Line, k - 1)) > 1e-7)
			result.push_back(points[k].chainage);
	}
	for (auto const &c : curves(Line))
	{
		result.push_back(c.start);
		result.push_back(c.end);
	}
	std::sort(result.begin(), result.end());
	return result;
}

std::vector<double> fewest_breaks(line const &Line, std::vector<double> Joints, double const Tolerance)
{
	std::sort(Joints.begin(), Joints.end());
	Joints.erase(std::unique(Joints.begin(), Joints.end(), [](double const A, double const B) { return std::abs(A - B) < 1e-6; }), Joints.end());
	auto const all{breaks(Line)};
	auto const inward = [&](double const Chainage, double const Other) { return grade(Line, Chainage + std::min(1e-3, std::abs(Other - Chainage) * 0.5) * (Other > Chainage ? 1.0 : -1.0)); };
	auto const follows = [&](double const From, double const To) {
		auto const length{To - From};
		auto const h0{elevation(Line, From)};
		auto const h1{elevation(Line, To)};
		auto const g0{inward(From, To) * length};
		auto const g1{inward(To, From) * length};
		auto const count{std::max(4, static_cast<int>(std::ceil(length)))};
		for (int i = 1; i < count; ++i)
		{
			auto const u{static_cast<double>(i) / count};
			auto const u2{u * u};
			auto const u3{u2 * u};
			auto const height{(2.0 * u3 - 3.0 * u2 + 1.0) * h0 + (u3 - 2.0 * u2 + u) * g0 + (-2.0 * u3 + 3.0 * u2) * h1 + (u3 - u2) * g1};
			if (std::abs(height - elevation(Line, From + u * length)) > Tolerance)
				return false;
		}
		return true;
	};
	std::vector<double> result;
	for (std::size_t k = 0; k + 1 < Joints.size(); ++k)
	{
		std::vector<double> stops;
		for (auto const chainage : all)
			if (chainage > Joints[k] + 1e-6 && chainage < Joints[k + 1] - 1e-6)
				stops.push_back(chainage);
		stops.push_back(Joints[k + 1]);
		auto from{Joints[k]};
		std::size_t next{0};
		while (next + 1 < stops.size())
		{
			auto reach{next};
			for (auto i{next + 1}; i < stops.size() && follows(from, stops[i]); ++i)
				reach = i;
			if (reach + 1 == stops.size())
				break;
			result.push_back(stops[reach]);
			from = stops[reach];
			next = reach + 1;
		}
	}
	return result;
}

std::string format_chainage(double const Chainage)
{
	auto const value{std::abs(Chainage)};
	auto const kilometres{static_cast<long>(value / 1000.0)};
	return format("%s%ld+%06.2f", Chainage < 0.0 ? "-" : "", kilometres, value - kilometres * 1000.0);
}

}