#include <psemek/app/application_base.hpp>
#include <psemek/app/default_application_factory.hpp>
#include <psemek/gfx/painter.hpp>
#include <psemek/gfx/gl.hpp>
#include <psemek/math/camera.hpp>
#include <psemek/math/gauss.hpp>
#include <psemek/math/intersection.hpp>
#include <psemek/random/generator.hpp>
#include <psemek/random/device.hpp>
#include <psemek/random/uniform_ball.hpp>
#include <psemek/random/uniform.hpp>
#include <psemek/pcg/perlin.hpp>
#include <psemek/pcg/fractal.hpp>
#include <psemek/util/ndarray.hpp>
#include <psemek/util/clock.hpp>
#include <psemek/log/log.hpp>
#include <psemek/io/file_stream.hpp>

using namespace psemek;

auto make_perlin(random::generator & rng, int min_octave, int max_octave, float power)
{
	std::vector<pcg::perlin<float, 2>> octaves;
	std::vector<float> weights;

	random::uniform_sphere_vector_distribution<float, 2> random_vector{};

	for (int octave = min_octave; octave < max_octave; ++octave)
	{
		int size = 1 << octave;
		util::ndarray<math::vector<float, 2>, 2> gradients({size + 1, size + 1});
		for (auto & g : gradients)
			g = random_vector(rng);
		octaves.emplace_back(std::move(gradients));
		weights.push_back(std::pow(power, - octave));
	}

	float weight_sum = 0.f;
	for (auto w : weights)
		weight_sum += w;
	for (auto & w : weights)
		w /= weight_sum;

	return pcg::fractal<pcg::perlin<float, 2>>(std::move(octaves), std::move(weights));
}

void make_random_vector_field(random::generator & rng, util::ndarray<math::vector<float, 2>, 2> & result, int min_octave, int max_octave, float scale)
{
	auto noise_1 = make_perlin(rng, min_octave, max_octave, 2.f);
	auto noise_2 = make_perlin(rng, min_octave, max_octave, 2.f);

	for (int y = 0; y < result.height(); ++y)
	{
		for (int x = 0; x < result.width(); ++x)
		{
			math::point p{(x + 0.5f) / result.height(), (y + 0.5f) / result.width()};
			result(x, y) = math::normalized(math::vector{2.f * noise_1(p) - 1.f, 2.f * noise_2(p) - 1.f}) * scale;
		}
	}
}

template <typename T>
void add(util::ndarray<T, 2> & dst, util::ndarray<T, 2> const & src)
{
	auto src_begin = src.begin();
	auto src_end = src.end();
	auto dst_begin = dst.begin();
	for (; src_begin != src_end;)
		*dst_begin++ += *src_begin++;
}

template <typename T>
void scale(util::ndarray<T, 2> & array, float factor)
{
	for (auto & v : array)
		v *= factor;
}

float do_lerp(float x, float y, float t)
{
	return math::lerp(x, y, t);
}

math::vector<float, 2> do_lerp(math::vector<float, 2> const & x, math::vector<float, 2> const & y, float t)
{
	return math::lerp(length(x), length(y), t) * math::slerp(math::normalized(x), math::normalized(y), t);
}

template <typename T>
void lerp(util::ndarray<T, 2> & dst, util::ndarray<T, 2> const & src0, util::ndarray<T, 2> const & src1, float t)
{
	auto src0_begin = src0.begin();
	auto src0_end = src0.end();
	auto src1_begin = src1.begin();
	auto dst_begin = dst.begin();

	for (; src0_begin != src0_end;)
		*dst_begin++ = do_lerp(*src0_begin++, *src1_begin++, t);
}

template <typename T>
T sample_bilinear(util::ndarray<T, 2> const & array, math::point<float, 2> p)
{
	p[0] = math::clamp(p[0] - 0.5f, {0.f, array.width() - 1.f});
	p[1] = math::clamp(p[1] - 0.5f, {0.f, array.height() - 1.f});

	int ix = std::min<int>(array.width() - 2, std::floor(p[0]));
	int iy = std::min<int>(array.height() - 2, std::floor(p[1]));

	float tx = p[0] - ix;
	float ty = p[1] - iy;

	return math::lerp(
		math::lerp(array(ix + 0, iy + 0), array(ix + 1, iy + 0), tx),
		math::lerp(array(ix + 0, iy + 1), array(ix + 1, iy + 1), tx),
		ty
	);
}

std::string season_str(int season)
{
	switch (season)
	{
	case 0: return "Winter";
	case 1: return "Spring";
	case 2: return "Summer";
	case 3: return "Autumn";
	}

	return "???";
}

std::string day_night_str(int day_night)
{
	switch (day_night)
	{
	case 0: return "Night";
	case 1: return "Day";
	}

	return "???";
}

struct climate_snapshot
{
	util::ndarray<math::vector<float, 2>, 2> wind_velocity;
	util::ndarray<float, 2> temperature;
	util::ndarray<float, 2> humidity;
};

static constexpr int N = 128;

struct solver
{
	util::ndarray<float, 2> const & terrain;
	random::generator & rng;
	int season;
	int day_night;

	climate_snapshot & snapshot;

	int solve_velocity_iterations = 0;
	int max_solve_velocity_iterations = N / 2;

	int solve_temperature_iterations = 0;
	int max_solve_temperature_iterations = N / 4;

	int solve_humidity_iterations = 0;
	int max_solve_humidity_iterations = 3 * N;

	util::ndarray<math::vector<float, 2>, 2> initial_random_wind_velocity;
	util::ndarray<float, 2> pressure;
	util::ndarray<float, 2> new_temperature;
	util::ndarray<float, 2> new_humidity;

	int stage = 0;
	bool finished = false;

	solver(util::ndarray<float, 2> const & terrain,	random::generator & rng, int season, int day_night, climate_snapshot & snapshot)
		: terrain(terrain)
		, rng(rng)
		, season(season)
		, day_night(day_night)
		, snapshot(snapshot)
	{
		initial_random_wind_velocity.resize({N, N});
		snapshot.wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
		pressure.resize({N, N}, 0.f);
		snapshot.temperature.resize({N, N}, 0.f);
		new_temperature.resize({N, N}, 0.f);
		snapshot.humidity.resize({N, N}, 0.f);
		new_humidity.resize({N, N}, 0.f);

		make_random_vector_field(rng, initial_random_wind_velocity, 3, 6, 0.005f);

		for (int y = 0; y < N; ++y)
			for (int x = 0; x < N; ++x)
				if (terrain(x, y) < 0.f)
					initial_random_wind_velocity(x, y) *= 5.f;
				else
					initial_random_wind_velocity(x, y) *= std::exp(- 2.f * terrain(x, y));

		snapshot.wind_velocity = initial_random_wind_velocity.copy();

		for (int y = 0; y < N; ++y)
			for (int x = 0; x < N; ++x)
				snapshot.temperature(x, y) = expected_temperature(y);
	}

	float expected_temperature(int y)
	{
		return math::lerp(40.f, -20.f, (y + 0.5f) / N) // Base latitude temperature in spring/autumn
			- 15.f * std::cos(season * math::pi * 0.5f) // Seasonal offset of +/-15C
			- 5.f * std::cos(day_night * math::pi) // Day-night fluctuation of +/-5C
			;
	}

	float expected_humidity(int x, int y)
	{
		if (terrain(x, y) < 0.f)
			return 1.f;
		else
			return 0.f;
	}

	void step()
	{
		if (stage == 0)
		{
			for (int y = 1; y + 1 < N; ++y)
			{
				for (int x = 1; x + 1 < N; ++x)
				{
					float divergence = 0.f
						+ (initial_random_wind_velocity(x + 1, y)[0] - initial_random_wind_velocity(x - 1, y)[0]) * 0.5f
						+ (initial_random_wind_velocity(x, y + 1)[1] - initial_random_wind_velocity(x, y + 1)[1]) * 0.5f
						;

					// (sum - 4.f * p(x,y)) = divergence
					// sum - 4.f * p = divergence
					// p = (sum - divergence) / 4.f
					pressure(x, y) = (pressure(x - 1, y) + pressure(x + 1, y) + pressure(x, y - 1) + pressure(x, y + 1) - divergence) / 4.f;
				}
			}

			for (int y = 1; y + 1 < N; ++y)
			{
				for (int x = 1; x + 1 < N; ++x)
				{
					math::vector pressure_gradient
					{
						(pressure(x + 1, y) - pressure(x - 1, y)) * 0.5f,
						(pressure(x, y + 1) - pressure(x, y - 1)) * 0.5f,
					};
					snapshot.wind_velocity(x, y) = initial_random_wind_velocity(x, y) - pressure_gradient;

					if (terrain(x, y) > 0.f)
						snapshot.wind_velocity(x, y) *= std::exp(- 2.f * terrain(x, y));
				}
			}

			++solve_velocity_iterations;
			if (solve_velocity_iterations == max_solve_velocity_iterations)
				++stage;
			return;
		}

		if (stage == 1)
		{
			for (int y = 0; y < N; ++y)
			{
				for (int x = 0; x < N; ++x)
				{
					float step_time = 100.f;

					math::point p{x + 0.5f, y + 0.5f};
					float advection = sample_bilinear(snapshot.temperature, p - snapshot.wind_velocity(x, y) * step_time);

					float heating_speed = (terrain(x, y) < 0.f) ? 0.001f : 0.0001f;
					float heating_delta = (expected_temperature(y) - snapshot.temperature(x, y)) * (- std::expm1(- heating_speed * step_time));

					new_temperature(x, y) = heating_delta + advection;
				}
			}
			std::swap(snapshot.temperature, new_temperature);

			for (int y = 1; y + 1 < N; ++y)
			{
				for (int x = 1; x + 1 < N; ++x)
				{
					float average = (0.f
						+ snapshot.temperature(x - 1, y)
						+ snapshot.temperature(x + 1, y)
						+ snapshot.temperature(x, y - 1)
						+ snapshot.temperature(x, y + 1)
					) * 0.25f;

					new_temperature(x, y) = math::lerp(snapshot.temperature(x, y), average, 0.5f);
				}
			}
			std::swap(snapshot.temperature, new_temperature);

			++solve_temperature_iterations;
			if (solve_temperature_iterations >= max_solve_temperature_iterations)
				++stage;
			return;
		}

		if (stage == 2)
		{
			for (int y = 0; y < N; ++y)
			{
				for (int x = 0; x < N; ++x)
				{
					float step_time = 100.f;

					math::point p{x + 0.5f, y + 0.5f};
					float advection = sample_bilinear(snapshot.humidity, p - snapshot.wind_velocity(x, y) * step_time);

					float evaporation_speed = 0.0001f;
					float evaporation_delta = (expected_humidity(x, y) - snapshot.humidity(x, y)) * (- std::expm1(- evaporation_speed * step_time));

					new_humidity(x, y) = evaporation_delta + advection;
				}
			}
			std::swap(snapshot.humidity, new_humidity);

			for (int y = 1; y + 1 < N; ++y)
			{
				for (int x = 1; x + 1 < N; ++x)
				{
					float average = (0.f
						+ snapshot.humidity(x - 1, y)
						+ snapshot.humidity(x + 1, y)
						+ snapshot.humidity(x, y - 1)
						+ snapshot.humidity(x, y + 1)
					) * 0.25f;

					new_humidity(x, y) = math::lerp(snapshot.humidity(x, y), average, 0.5f);
				}
			}
			std::swap(snapshot.humidity, new_humidity);

			++solve_humidity_iterations;
			if (solve_humidity_iterations >= max_solve_humidity_iterations)
			{
				++stage;
				finished = true;
			}
			return;
		}
	}
};

struct river_vertex
{
	// Empty for river deltas
	std::optional<math::point<int, 2>> parent;
	std::vector<math::point<int, 2>> children = {};
	float flow = 0.f;
};

struct river_net
{
	util::ndarray<std::optional<river_vertex>, 2> vertices;
	util::ndarray<bool, 2> water;

	std::vector<math::point<int, 2>> queue;
};

struct weather_app
	: app::application_base
{
	random::generator rng{random::device{}};
	// random::generator rng{0, 0};

	gfx::painter painter;

	math::box<float, 2> simulation_box;

	bool paused = true;
	bool show_velocity = false;
	bool show_temperature = false;
	bool show_humidity = false;
	bool show_biomes = true;
	bool show_rivers = true;

	gfx::pixmap_rgba biomes_map;
	util::ndarray<float, 2> terrain;

	// 0 - Winter
	// 1 - Spring
	// 2 - Summer
	// 3 - Autumn
	int season = 0;

	// 0 - Night
	// 1 - Day
	int day_night = 0;

	climate_snapshot snapshots[4][2];
	climate_snapshot average;
	climate_snapshot display_snapshot;

	bool need_update_display_snapshot = false;

	std::optional<struct solver> solver;

	river_net rivers;

	float display_season = 0.f;

	util::clock<> frame_clock;

	weather_app(options const &, context const & context)
	{
		simulation_box = {{{0.f, N}, {0.f, N}}};

		terrain.resize({N, N}, 0.f);
		auto heightmap = gfx::read_image<std::uint8_t>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "heightmap_seed_3.png"});
		for (int y = 0; y < N; ++y)
			for (int x = 0; x < N; ++x)
				terrain(x, y) = ((heightmap(x, y) / 255.f) * 2048.f - 512.f) / 1024.f;

		biomes_map = gfx::read_image<gfx::color_rgba>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "biomes.png"});

		context.vsync(false);
	}

	void on_event(app::key_event const & event) override
	{
		application_base::on_event(event);

		if (event.down && event.key == app::keycode::SPACE)
			paused ^= true;

		if (event.down && event.key == app::keycode::V)
			show_velocity ^= true;

		if (event.down && event.key == app::keycode::T)
			show_temperature ^= true;

		if (event.down && event.key == app::keycode::W)
			show_humidity ^= true;

		if (event.down && event.key == app::keycode::B)
			show_biomes ^= true;

		if (event.down && event.key == app::keycode::R)
			show_rivers ^= true;
	}

	void update() override
	{
		float const dt = frame_clock.restart().count();

		if (solver && solver->finished)
		{
			solver.reset();
			++day_night;
			if (day_night == 2)
			{
				day_night = 0;
				++season;
			}

			if (season == 4)
			{
				compute_average();
				init_river_deltas();
				need_update_display_snapshot = true;
			}
		}

		if (!solver && season < 4)
			solver.emplace(terrain, rng, season, day_night, snapshots[season][day_night]);

		if (!paused)
		{
			if (solver)
			{
				for (int i = 0; i < 16; ++i)
				{
					solver->step();
					if (solver->finished)
						break;
				}
			}
			else if (season == 4 && !rivers.queue.empty())
			{
				for (int i = 0; i < 16; ++i)
				{
					propagate_rivers();
					if (rivers.queue.empty())
					{
						compute_river_flow();
						break;
					}
				}
			}
		}

		if (state().key_down.contains(app::keycode::LEFT))
		{
			display_season -= 0.5f * dt;
			need_update_display_snapshot = true;
		}

		if (state().key_down.contains(app::keycode::RIGHT))
		{
			display_season += 0.5f * dt;
			need_update_display_snapshot = true;
		}

		while (display_season < 0.f)
			display_season += 1.f;
		while (display_season >= 1.f)
			display_season -= 1.f;
	}

	void compute_average()
	{
		average.wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
		average.temperature.resize({N, N}, 0.f);
		average.humidity.resize({N, N}, 0.f);

		for (auto const & season : snapshots)
			for (auto const & s : season)
			{
				add(average.wind_velocity, s.wind_velocity);
				add(average.temperature, s.temperature);
				add(average.humidity, s.humidity);
			}

		scale(average.wind_velocity, 1.f / 8.f);
		scale(average.temperature, 1.f / 8.f);
		scale(average.humidity, 1.f / 8.f);
	}

	void init_river_deltas()
	{
		rivers.vertices.resize({N, N});
		rivers.water.resize({N, N}, false);

		for (int y = 1; y + 1 < N; ++y)
		{
			for (int x = 1; x + 1 < N; ++x)
			{
				int land = 0;
				land += (terrain(x - 1, y - 1) >= 0.f ? 1 : 0);
				land += (terrain(x + 0, y - 1) >= 0.f ? 1 : 0);
				land += (terrain(x - 1, y + 0) >= 0.f ? 1 : 0);
				land += (terrain(x + 0, y + 0) >= 0.f ? 1 : 0);

				if (land > 0 && land < 4)
				{
					rivers.queue.push_back({x, y});
					rivers.vertices(x, y).emplace();
				}

				if (land < 4)
					rivers.water(x, y) = true;
			}
		}
	}

	void propagate_rivers()
	{
		auto i = random::uniform<int>(rng, 0, rivers.queue.size() - 1);
		if (i + 1 < rivers.queue.size())
			std::swap(rivers.queue[i], rivers.queue.back());

		auto p = rivers.queue.back();

		math::vector<int, 2> const neighbours[4]
		{
			{-1,  0},
			{ 1,  0},
			{ 0, -1},
			{ 0,  1},
		};

		float height = sample_bilinear(terrain, math::cast<float>(p));

		std::vector<math::point<int, 2>> candidates;
		for (auto n : neighbours)
		{
			auto q = p + n;
			if (q[0] == 0 || q[0] == N || q[1] == 0 || q[1] == N)
				continue;

			if (rivers.water(q[0], q[1]))
				continue;

			if (rivers.vertices(q[0], q[1]))
				continue;

			float nheight = sample_bilinear(terrain, math::cast<float>(q));

			if (nheight + 0.02f < height)
				continue;

			candidates.push_back(q);
		}

		if (candidates.size() <= 1)
			rivers.queue.pop_back();

		if (candidates.empty())
		{
			// TODO: remove if no parent and not going anywhere?
			return;
		}

		auto next = random::uniform_from(rng, candidates);

		rivers.vertices(p[0], p[1])->children.push_back(next);
		rivers.vertices(next[0], next[1]) = river_vertex{.parent = p};
		rivers.queue.push_back(next);
	}

	void compute_river_flow()
	{
		for (int y = 0; y < N; ++y)
		{
			for (int x = 0; x < N; ++x)
			{
				if (auto & v = rivers.vertices(x, y); v && !v->parent)
				{
					compute_river_flow_at({x, y});
				}
			}
		}
	}

	float compute_river_flow_at(math::point<int, 2> const & p)
	{
		auto & v = rivers.vertices(p[0], p[1]);
		if (!v) return 0.f;

		if (v->children.empty())
		{
			v->flow = std::max(0.f, sample_bilinear(average.humidity, math::cast<float>(p)) - 0.05f);
		}
		else
		{
			for (auto c : v->children)
				v->flow += compute_river_flow_at(c);
		}

		return v->flow;
	}

	void update_display_snapshot()
	{
		display_snapshot.wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
		display_snapshot.temperature.resize({N, N}, 0.f);
		display_snapshot.humidity.resize({N, N}, 0.f);

		int i = std::min<int>(3, std::floor(4.f * display_season));
		int j = (i + 1) % 4;
		float t = (4.f * display_season) - i;

		auto & s0 = snapshots[i][0];
		auto & s1 = snapshots[j][0];

		lerp(display_snapshot.wind_velocity, s0.wind_velocity, s1.wind_velocity, t);
		lerp(display_snapshot.temperature, s0.temperature, s1.temperature, t);
		lerp(display_snapshot.humidity, s0.humidity, s1.humidity, t);
	}

	void present() override
	{
		if (need_update_display_snapshot)
		{
			update_display_snapshot();
			need_update_display_snapshot = false;
		}

		gl::ClearColor(0.f, 0.f, 0.f, 0.f);
		gl::Clear(gl::COLOR_BUFFER_BIT);

		float const aspect_ratio = state().size[0] * 1.f / state().size[1];
		math::box<float, 2> view_box = math::expand(simulation_box, 1.f);
		if (view_box[0].length() / view_box[1].length() > aspect_ratio)
			view_box[1] = math::expand(view_box[1], (view_box[0].length() / aspect_ratio - view_box[1].length()) / 2.f);
		else
			view_box[0] = math::expand(view_box[0], (view_box[1].length() * aspect_ratio - view_box[0].length()) / 2.f);

		std::optional<math::vector<int, 2>> mouseover_cell;
		{
			auto mouse = math::lerp(view_box, math::vector{state().mouse[0] * 1.f / state().size[0], 1.f - state().mouse[1] * 1.f / state().size[1]});
			int x = std::floor(mouse[0]);
			int y = std::floor(mouse[1]);

			if (x >= 0 && x < N && y >= 0 && y < N)
				mouseover_cell = {x, y};
		}

		[[maybe_unused]] float const pixel_size = view_box[0].length() / state().size[0];

		auto map_color = [](float value, gfx::color_4f const & negative, gfx::color_4f const & positive){
			return math::lerp(negative, positive, 1.f / (1.f + std::exp(- value)));
		};

		auto map_color_positive = [](float value, gfx::color_4f const & zero, gfx::color_4f const & positive){
			return math::lerp(zero, positive, 2.f / (1.f + std::exp(- value)) - 1.f);
		};

		auto map_temperature = [&](float temperature)
		{
			return map_color(0.5f * std::floor(temperature / 10.f), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f});
		};

		auto map_humidity = [&](float humidity)
		{
			return map_color_positive(std::round(10.f * humidity), {0.f, 1.f, 1.f, 0.f}, {0.f, 1.f, 1.f, 0.75f});
		};

		auto map_biome = [this](float temperature, float precipitation)
		{
			auto x = math::clamp<int>(math::unlerp({0.f, 0.5f}, precipitation) * biomes_map.width() , {0, biomes_map.width()  - 1});
			auto y = math::clamp<int>(math::unlerp({-10.f, 30.f}, temperature  ) * biomes_map.height(), {0, biomes_map.height() - 1});

			return gfx::to_colorf(biomes_map(x, y));
		};

		auto & avg_snapshot = (season < 4) ? snapshots[season][day_night] : average;
		auto & snapshot = (season < 4) ? snapshots[season][day_night] : display_snapshot;

		for (int y = 0; y < N; ++y)
		{
			for (int x = 0; x < N; ++x)
			{
				gfx::color_4f color = gfx::color_4f::zero();

				bool const is_water = terrain(x, y) < 0.f;
				bool const is_land = !is_water;

				if (is_water)
					color = map_color(8.f * terrain(x, y), {0.f, 0.f, 0.125f, 1.f}, {0.f, 1.f, 1.5f, 1.f});
				else if (show_biomes)
				{
					float temperature = avg_snapshot.temperature(x, y) - std::max(0.f, terrain(x, y)) * 30.f;
					color = map_biome(temperature, avg_snapshot.humidity(x, y));
				}
				else
					color = gfx::color_4f{0.3f, 0.5f, 0.1f, 1.f};

				if (is_land && x > 0 && x + 1 < N && y > 0 && y + 1 < N)
				{
					math::vector terrain_gradient
					{
						(terrain(x + 1, y) - terrain(x - 1, y)) / 2.f,
						(terrain(x, y + 1) - terrain(x, y - 1)) / 2.f,
					};

					auto terrain_normal = math::normalized(math::vector{-terrain_gradient[0], -terrain_gradient[1], 0.125f});

					float lightness = 0.5f + 0.5f * math::dot(terrain_normal, math::normalized(math::vector{1.f, 2.f, 3.f}));

					color = gfx::dark(color, 1.f - lightness);
				}

				if (show_temperature)
				{
					color = gfx::blend(color, map_temperature(snapshot.temperature(x, y)));
				}

				if (show_humidity)
				{
					color = gfx::blend(color, map_humidity(snapshot.humidity(x, y)));
				}

				painter.rect({{{x, x + 1.f}, {y, y + 1.f}}}, gfx::to_coloru8(color));
			}
		}

		if (show_velocity)
		{
			for (int y = 0; y < N; ++y)
			{
				for (int x = 0; x < N; ++x)
				{
					math::point center{x + 0.5f, y + 0.5f};
					auto v = snapshot.wind_velocity(x, y);
					if (auto l = math::length(v); l > 0.f)
					{
						float const magnification = 200.f;
						float const max_length = 1.5f;
						v *= 0.5f * max_length * (1.f - std::exp(- magnification * l)) / l;
					}
					auto n = math::ort(v) * 0.3f;

					auto color = gfx::color_4f{1.f, 1.f, 1.f, 1.f};
					if (!show_temperature)
					{
						float t = 0.2f * snapshot.temperature(x, y);
						if (t >= 0.f)
							color = map_color_positive( t, color, {1.f, 0.5f, 0.125f, 1.f});
						else
							color = map_color_positive(-t, color, {0.125f, 0.5f, 1.f, 1.f});
					}

					painter.triangle(center - v - n, center - v + n, center + v, gfx::to_coloru8(color));
				}
			}
		}

		if (show_rivers)
		{
			for (auto const & p : rivers.queue)
			{
				painter.circle(math::cast<float>(p), 3.f * pixel_size, {0, 255, 255, 255}, 12);
			}

			if (!rivers.vertices.empty())
				for (int y = 0; y < N; ++y)
				{
					for (int x = 0; x < N; ++x)
					{
						auto & v = rivers.vertices(x, y);
						if (!v || !v->parent) continue;

						float width = 2.f;
						if (rivers.queue.empty())
							width = 2.f * std::sqrt(v->flow);

						painter.line({x, y}, math::cast<float>(*v->parent), width * pixel_size, {0, 255, 255, 255}, false);
					}
				}
		}

		int row = 0;
		auto push_text = [&](std::string const & text) mutable
		{
			painter.text(view_box.corner(0.f, 1.f) - math::vector{0.f, row * pixel_size * 2.f * 12.f}, text, {.scale = {2.f * pixel_size, - 2.f * pixel_size}, .x = gfx::painter::x_align::left, .y = gfx::painter::y_align::top, .c = {255, 255, 255, 255}});
			++row;
		};

		if (solver)
		{
			push_text(std::format("{} {}", season_str(season), day_night_str(day_night)));

			if (solver->stage >= 0)
				push_text(std::format("Velocity {}/{}", solver->solve_velocity_iterations, solver->max_solve_velocity_iterations));
			else
				push_text("");

			if (solver->stage >= 1)
				push_text(std::format("Temperature {}/{}", solver->solve_temperature_iterations, solver->max_solve_temperature_iterations));
			else
				push_text("");

			if (solver->stage >= 2)
				push_text(std::format("Humidity {}/{}", solver->solve_humidity_iterations, solver->max_solve_humidity_iterations));
			else
				push_text("");

			push_text("");
		}

		if (mouseover_cell)
		{
			int x = (*mouseover_cell)[0];
			int y = (*mouseover_cell)[1];
			painter.rect({{{x, x + 1.f}, {y, y + 1.f}}}, {255, 255, 255, 127});

			push_text(std::format("{} {}", x, y));
			push_text(std::format("V = {:.3f}", length(snapshot.wind_velocity(x, y)) * 1000.f));
			push_text(std::format("H = {:.3f}", terrain(x, y) * 1024.f));
			push_text(std::format("T = {:.3f}", snapshot.temperature(x, y)));
			push_text(std::format("W = {:.3f}", snapshot.humidity(x, y)));
			push_text("");
		}

		if (season == 4)
		{
			++row;

			float axis_width = 300.f;
			float axis_height = 24.f;

			float axis_left = view_box[0].min + 30.f * pixel_size;
			float axis_bottom = view_box[1].max - row * pixel_size * 2.f * 12.f - axis_height * pixel_size;

			auto line = [&](float x0, float y0, float x1, float y1, float w, gfx::color_rgba const & c)
			{
				painter.line(
					{axis_left + x0 * axis_width * pixel_size, axis_bottom + y0 * axis_height * pixel_size},
					{axis_left + x1 * axis_width * pixel_size, axis_bottom + y1 * axis_height * pixel_size},
					pixel_size * w, c, false
				);
			};

			line(0.f, 0.5f, 1.f, 0.5f, 2.f, {255, 255, 255, 127});
			line(display_season, 0.f, display_season, 1.f, 4.f, {255, 255, 255, 255});

			++row;
			++row;
		}

		if (mouseover_cell && season == 4)
		{
			int x = (*mouseover_cell)[0];
			int y = (*mouseover_cell)[1];

			float plot_width = 300.f;
			float plot_height = 150.f;

			float plot_left = view_box[0].min + 30.f * pixel_size;
			float plot_bottom = view_box[1].max - row * pixel_size * 2.f * 12.f - plot_height * 1.5f * pixel_size;

			auto line = [&](float x0, float y0, float x1, float y1, gfx::color_rgba const & c)
			{
				painter.line(
					{plot_left + x0 * plot_width * pixel_size, plot_bottom + y0 * plot_height * pixel_size},
					{plot_left + x1 * plot_width * pixel_size, plot_bottom + y1 * plot_height * pixel_size},
					pixel_size * 2.f, c, false
				);
			};

			auto begin_plot = [&]{
				painter.rect({{{plot_left, plot_left + plot_width * pixel_size}, {plot_bottom, plot_bottom + plot_height * pixel_size}}}, {255, 255, 255, 31});

				for (int x = 0; x <= 4; ++x)
					line(x / 4.f, 0.f, x / 4.f, 1.f, {255, 255, 255, (x == 0 || x == 4) ? 127 : 31});
			};

			auto draw_plot = [&](std::vector<float> const & ys, gfx::color_rgba const & c)
			{
				std::vector<float> xs(5);
				for (int i = 0; i <= 4; ++i)
					xs[i] = i / 4.f;

				auto system = math::matrix<float, 16, 16>::zero();
				auto coeffs = math::vector<float, 16>::zero();

				for (int i = 0; i < 4; ++i)
				{
					float x0 = xs[i];
					system[2 * i + 0][4 * i + 0] = 1.f;
					system[2 * i + 0][4 * i + 1] = x0;
					system[2 * i + 0][4 * i + 2] = x0 * x0;
					system[2 * i + 0][4 * i + 3] = x0 * x0 * x0;
					coeffs[2 * i + 0] = ys[i];

					float x1 = xs[i + 1];
					system[2 * i + 1][4 * i + 0] = 1.f;
					system[2 * i + 1][4 * i + 1] = x1;
					system[2 * i + 1][4 * i + 2] = x1 * x1;
					system[2 * i + 1][4 * i + 3] = x1 * x1 * x1;
					coeffs[2 * i + 1] = ys[i + 1];
				}

				for (int i = 0; i < 4; ++i)
				{
					int j = (i + 1) % 4;
					float x0 = xs[i + 1];
					float x1 = xs[j];

					system[8 + 2 * i + 0][4 * i + 1] = 1.f;
					system[8 + 2 * i + 0][4 * i + 2] = 2.f * x0;
					system[8 + 2 * i + 0][4 * i + 3] = 3.f * x0 * x0;
					system[8 + 2 * i + 0][4 * j + 1] = - 1.f;
					system[8 + 2 * i + 0][4 * j + 2] = - 2.f * x1;
					system[8 + 2 * i + 0][4 * j + 3] = - 3.f * x1 * x1;

					system[8 + 2 * i + 1][4 * i + 2] = 2.f;
					system[8 + 2 * i + 1][4 * i + 3] = 6.f * x0;
					system[8 + 2 * i + 1][4 * j + 2] = - 2.f;
					system[8 + 2 * i + 1][4 * j + 3] = - 6.f * x1;
				}

				math::gauss(system, coeffs);

				auto eval = [&](int i, float x)
				{
					return 0.f
						+ coeffs[4 * i + 0]
						+ coeffs[4 * i + 1] * x
						+ coeffs[4 * i + 2] * x * x
						+ coeffs[4 * i + 3] * x * x * x
						;
				};

				for (int i = 0; i < 4; ++i)
				{
					int D = 8;
					for (int j = 0; j < D; ++j)
					{
						float x0 = (i + (j + 0.f) / D) / 4.f;
						float x1 = (i + (j + 1.f) / D) / 4.f;
						line(x0, eval(i, x0), x1, eval(i, x1), c);
					}

					if (i == std::min<int>(3, std::floor(4.f * display_season)))
					{
						float x = display_season;
						float y = eval(i, x);
						painter.circle({plot_left + x * plot_width * pixel_size, plot_bottom + y * plot_height * pixel_size}, 4.f * pixel_size, c);
					}
				}
			};

			// Temperature
			{
				begin_plot();

				for (int y = 0; y <= 10; ++y)
				{
					line(0.f, y / 10.f, 1.f, y / 10.f, {255, 255, 255, (y == 5) ? 127 : 31});
				}

				std::vector<float> ys[2];

				for (int d = 0; d < 2; ++d)
				{
					for (int s = 0; s <= 4; ++s)
					{
						float value = 0.5f + snapshots[s % 4][d].temperature(x, y) / 100.f;
						ys[d].push_back(value);
					}
				}

				draw_plot(ys[0], gfx::color_rgba{0, 127, 255, 255});
				draw_plot(ys[1], gfx::color_rgba{255, 127, 0, 255});
			}

			plot_bottom -= (plot_height + 50.f) * pixel_size;

			// Humidity
			{
				begin_plot();

				for (int y = 0; y <= 5; ++y)
				{
					line(0.f, y / 5.f, 1.f, y / 5.f, {255, 255, 255, (y == 0) ? 127 : 31});
				}

				std::vector<float> ys;

				for (int s = 0; s <= 4; ++s)
				{
					float value = (snapshots[s % 4][0].humidity(x, y) + snapshots[s % 4][1].humidity(x, y)) / 2.f;
					ys.push_back(value);
				}

				draw_plot(ys, {0, 255, 255, 255});
			}
		}

		painter.render(math::orthographic_camera{view_box}.transform());
	}
};

namespace psemek::app
{

	std::unique_ptr<application::factory> make_application_factory()
	{
		return default_application_factory<weather_app>({.name = "Weather simulation test"});
	}

}