#include <psemek/app/application_base.hpp>
#include <psemek/app/default_application_factory.hpp>
#include <psemek/phys/engine_2d.hpp>
#include <psemek/gfx/painter.hpp>
#include <psemek/geom/simplex.hpp>
#include <psemek/geom/math.hpp>
#include <psemek/geom/camera.hpp>
#include <psemek/geom/scale.hpp>
#include <psemek/util/hash_table.hpp>
#include <psemek/util/clock.hpp>
#include <psemek/random/generator.hpp>
#include <psemek/random/uniform.hpp>
#include <psemek/random/normal.hpp>
#include <psemek/random/weighted.hpp>
#include <psemek/log/log.hpp>
#include <psemek/util/thread.hpp>
#include <psemek/util/to_string.hpp>

#include <deque>

using namespace psemek;

static float const sim_dt = 0.005f;
static float const gravity = 25.f;
static float const air_friction = 0.1f;
static float const ground_friction = 100.f;
static float const hills_friction = 100.f;
static float const spring_damping = 10.f;
static float const max_spring_force = 25000.f;

static geom::interval<float> const cell_size_range{1.f, 1.f};
static geom::interval<float> const cell_period_range{0.25f, 4.f};
static float const muscle_expand_extent = 0.5f;

static int const population_size = 4 * 1024;

static float const simulation_time = 30.f;
static float const display_simulation_time = 30.f;

enum class cell_type
{
	hard_tissue,
	soft_tissue,
	muscle,
	dead,
};

struct cell_info
{
	cell_type type;
	float size = 1.f;
	float phase = 0.f;
	float period = 1.f;
};

float spring_force(cell_type type)
{
	switch (type)
	{
	case cell_type::soft_tissue:
		return 2000.f;
	case cell_type::hard_tissue:
		return 5000.f;
	case cell_type::muscle:
		return 2000.f;
	case cell_type::dead:
		return 100.f;
	}

	return 0.f;
}

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

struct map
{
	geom::vector<float, 2> ground_normal{0.f, 1.f};

	std::vector<geom::point<float, 2>> ground_points;
};

using creature_blueprint = std::unordered_map<geom::point<int, 2>, cell_info>;

struct creature
{
	struct cell
	{
		std::uint32_t indices[4];
		cell_info info;
	};

	creature_blueprint blueprint;
	int generation = 0;
	std::optional<float> score = std::nullopt;

	std::vector<geom::point<float, 2>> positions;
	std::vector<geom::vector<float, 2>> velocities;
	std::vector<cell> cells;

	std::vector<geom::segment<std::uint32_t>> outer_edges;

	geom::point<float, 2> center() const;
	geom::vector<float, 2> center_velocity() const;
	void translate(geom::vector<float, 2> const & delta);
	void push(geom::vector<float, 2> const & delta_velocity);

	bool dead() const;

	void update(float dt);
	void collide(map const & map, float dt);

	void compute_outer_edges();
};

geom::point<float, 2> creature::center() const
{
	geom::vector<float, 2> sum{0.f, 0.f};
	for (int i = 1; i < positions.size(); ++i)
	{
		sum += positions[i] - positions[0];
	}
	return positions[0] + sum / (1.f * positions.size());
}

geom::vector<float, 2> creature::center_velocity() const
{
	geom::vector<float, 2> sum{0.f, 0.f};
	for (int i = 1; i < velocities.size(); ++i)
	{
		sum += velocities[i];
	}
	return sum / (1.f * velocities.size());
}

void creature::translate(geom::vector<float, 2> const & delta)
{
	for (auto & pos : positions)
		pos += delta;
}

void creature::push(geom::vector<float, 2> const & delta_velocity)
{
	for (auto & vel : velocities)
		vel += delta_velocity;
}

bool creature::dead() const
{
	for (auto const & cell : cells)
		if (cell.info.type == cell_type::dead)
			return true;
	return false;
}

void creature::update(float dt)
{
	for (auto & vel : velocities)
		vel[1] -= gravity * dt;

	for (auto & vel : velocities)
	{
		auto v = geom::length(vel);
		vel -= air_friction * v * vel * dt;
	}

	for (int i = 0; i < positions.size(); ++i)
		positions[i] += velocities[i] * dt;

	static geom::vector<float, 2> const deltas[4]
	{
		{-0.5f, -0.5f},
		{ 0.5f, -0.5f},
		{ 0.5f,  0.5f},
		{-0.5f,  0.5f},
	};

	for (auto & cell : cells)
	{
		float size = cell.info.size;

		if (cell.info.type == cell_type::muscle)
		{
			cell.info.phase += dt / cell.info.period;
			while (cell.info.phase >= 1.f)
				cell.info.phase -= 1.f;
			size *= 1.f + muscle_expand_extent * std::sin(cell.info.phase * 2.f * geom::pi);
		}

		geom::point<float, 2> center{0.f, 0.f};
		for (auto idx : cell.indices)
		{
			center += (positions[idx] - geom::point{0.f, 0.f}) * 0.25f;
		}

		float A = 0.f;
		float B = 0.f;

		for (int i = 0; i < 4; ++i)
		{
			auto const p = positions[cell.indices[i]] - center;
			auto const q = deltas[i] * size;

			A += geom::dot(p, q);
			B += geom::det(p, q);
		}

		float const angle = - std::atan2(B, A);

		for (int i = 0; i < 4; ++i)
		{
			auto const target = center + geom::rotate(deltas[i] * size, angle);

			auto force = spring_force(cell.info.type) * (target - positions[cell.indices[i]]);
			if (auto f = geom::length(force); f > max_spring_force)
				cell.info.type = cell_type::dead;
			velocities[cell.indices[i]] += force * dt;
		}

		auto cmvel = geom::vector{0.f, 0.f};
		auto rotation = 0.f;
		for (auto idx : cell.indices)
		{
			cmvel += velocities[idx] * 0.25f;
			rotation += geom::det(positions[idx] - center, velocities[idx]) * (0.25f / geom::length_sqr(positions[idx] - center));
		}

		for (auto idx : cell.indices)
		{
			auto target_vel = cmvel + geom::ort(positions[idx] - center) * rotation;

			velocities[idx] += (target_vel - velocities[idx]) * (1.f - std::exp(- spring_damping * dt));
		}
	}
}

void creature::collide(map const & map, float dt)
{
	for (int i = 0; i < positions.size(); ++i)
	{
		auto delta = positions[i] - geom::point<float, 2>{0.f, 0.f};
		auto dist = geom::dot(delta, map.ground_normal);
		if (dist < 0.f)
		{
			positions[i] -= dist * map.ground_normal;

			auto tangent = geom::ort(map.ground_normal);
			auto vn = geom::dot(velocities[i], map.ground_normal);
			auto vt = geom::dot(velocities[i], tangent);

			vn = std::max(0.f, vn);
			vt *= std::exp(- dt * ground_friction);

			velocities[i] = map.ground_normal * vn + tangent * vt;
		}

		auto it = std::upper_bound(map.ground_points.begin(), map.ground_points.end(), positions[i][0], [](float v, auto const & p){ return v < p[0]; });
		if (it != map.ground_points.begin() && it != map.ground_points.end())
		{
			auto j = (it - map.ground_points.begin()) - 1;

			auto n = geom::normalized(geom::ort(map.ground_points[j + 1] - map.ground_points[j]));
			float dist = geom::dot(positions[i] - map.ground_points[j], n);

			if (dist < 0.f)
			{
				positions[i] -= dist * n;

				auto tangent = geom::ort(n);
				auto vn = geom::dot(velocities[i], n);
				auto vt = geom::dot(velocities[i], tangent);

				vn = std::max(0.f, vn);
				vt *= std::exp(- dt * hills_friction);

				velocities[i] = n * vn + tangent * vt;
			}
		}
	}
}

void creature::compute_outer_edges()
{
	util::hash_set<geom::segment<std::uint32_t>> edge_set;

	for (auto const & cell : cells)
	{
		for (int i = 0; i < 4; ++i)
		{
			geom::segment<std::uint32_t> segment;
			segment[0] = cell.indices[i];
			segment[1] = cell.indices[(i + 1) % 4];
			edge_set.insert(segment);
		}
	}

	outer_edges.clear();
	for (auto const & edge : edge_set)
	{
		auto dual = edge;
		std::swap(dual[0], dual[1]);
		if (!edge_set.contains(dual))
			outer_edges.push_back(edge);
	}
}

gfx::color_rgba cell_color(cell_type type)
{
	switch (type)
	{
	case cell_type::soft_tissue:
		return {255, 255, 255, 255};
	case cell_type::hard_tissue:
		return {127, 127, 127, 255};
	case cell_type::muscle:
		return {255, 127, 127, 255};
	case cell_type::dead:
		return {0, 0, 0, 0};
	}

	return {255, 0, 255, 255};
}

void draw(gfx::painter & painter, creature const & creature)
{
	for (auto const & cell : creature.cells)
	{
		geom::point<float, 2> center{0.f, 0.f};
		for (auto idx : cell.indices)
			center += 0.25f * (creature.positions[idx] - geom::point{0.f, 0.f});

		gfx::color_rgba const color = cell_color(cell.info.type);
		gfx::color_rgba const bgcolor = gfx::dark(color, 0.25f);

		geom::point<float, 2> ps[4];
		for (int i = 0; i < 4; ++i)
			ps[i] = creature.positions[cell.indices[i]];

		painter.triangle(ps[0], ps[1], ps[2], bgcolor);
		painter.triangle(ps[2], ps[0], ps[3], bgcolor);

		for (int i = 0; i < 4; ++i)
			ps[i] = geom::lerp(ps[i], center, 0.25f);

		painter.triangle(ps[0], ps[1], ps[2], color);
		painter.triangle(ps[2], ps[0], ps[3], color);
	}

	for (auto const & edge : creature.outer_edges)
	{
		gfx::color_rgba const color{0, 0, 0, 255};
		float const width = 0.125f;
		painter.line(creature.positions[edge.points[0]], creature.positions[edge.points[1]], width, color, false);
	}
}

struct creature_builder
{
	creature_builder() = default;

	creature_builder(creature_blueprint blueprint)
		: blueprint_(std::move(blueprint))
	{}

	void add(geom::point<int, 2> const position, cell_info const & info)
	{
		blueprint_[position] = info;
	}

	bool contains(geom::point<int, 2> const & position) const
	{
		return blueprint_.contains(position);
	}

	creature build(int generation)
	{
		static geom::vector<int, 2> const vertex_delta[4]
		{
			{0, 0},
			{1, 0},
			{1, 1},
			{0, 1},
		};

		creature result;
		result.generation = generation;

		util::hash_map<geom::point<int, 2>, std::uint32_t> vertex_id;

		for (auto const & cell : blueprint_)
		{
			auto & cell_out = result.cells.emplace_back();
			cell_out.info = cell.second;
			for (int i = 0; i < 4; ++i)
			{
				geom::point<int, 2> vertex = cell.first + vertex_delta[i];
				if (auto it = vertex_id.find(vertex); it != vertex_id.end())
				{
					cell_out.indices[i] = it->second;
				}
				else
				{
					cell_out.indices[i] = (vertex_id[vertex] = result.positions.size());
					result.positions.push_back(geom::cast<float>(vertex));
					result.velocities.push_back(geom::vector<float, 2>::zero());
				}
			}
		}

		result.blueprint = std::move(blueprint_);
		result.compute_outer_edges();
		return result;
	}

	bool connected() const
	{
		if (blueprint_.empty())
			return false;

		util::hash_set<geom::point<int, 2>> visited;

		std::deque<geom::point<int, 2>> queue;
		queue.push_back(blueprint_.begin()->first);
		visited.insert(queue.back());

		while (!queue.empty())
		{
			auto cur = queue.front();
			queue.pop_front();

			for (auto n : neighbours)
			{
				auto nn = cur + n;
				if (blueprint_.contains(nn) && !visited.contains(nn))
				{
					visited.insert(nn);
					queue.push_back(nn);
				}
			}
		}

		return visited.size() == blueprint_.size();
	}

private:
	creature_blueprint blueprint_;
};

int const display_creature_count = 4;

cell_info random_cell(random::generator & rng)
{
	cell_info result;
	result.size = random::uniform(rng, cell_size_range);
	if (auto t = random::uniform(rng, 0, 2); t == 0)
		result.type = cell_type::hard_tissue;
	else if (t == 1)
		result.type = cell_type::soft_tissue;
	else
	{
		result.type = cell_type::muscle;
		result.phase = random::uniform<float>(rng);
		result.period = random::uniform<float>(rng, cell_period_range);
	}
	return result;
}

struct soft_creatures_2d_app
	: app::application_base
{
	soft_creatures_2d_app(options const &, context const &)
	{
		random::generator rng{random::device{}};

//		map_.ground_points.push_back({5.f, 0.f});
//		for (int x = 10; x <= 200; x += 1)
//			map_.ground_points.push_back({x / 2.f, random::uniform(rng, 0.f, std::min(x, 200) / 200.f * 4.f)});
//			map_.ground_points.push_back({x, std::min(2.f, geom::sqr(x - 10) * 0.25f)});
//			map_.ground_points.push_back({x, random::uniform(rng, 0.f, 4.f)});

		for (int species = 0; species < population_size; ++species)
		{
			while (true)
			{
				int x_size = random::uniform(rng, 2, 7);
				int y_size = random::uniform(rng, 2, 5);

				x_size = 1;
				y_size = 1;

				creature_builder builder;
				for (int x = 0; x < x_size; ++x)
				{
					for (int y = 0; y < y_size; ++y)
					{
						if (random::uniform<float>(rng) > 0.75f)
							continue;

						builder.add({x, y}, random_cell(rng));
					}
				}

				if (builder.connected())
				{
					population_.push_back(builder.build(0));
					break;
				}
			}

			if (species < display_creature_count)
				display_creatures_.push_back(population_.back());
		}

		simulator_thread_ = util::thread([this]{ run_simulation(); });
	}

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

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

	void stop() override
	{
		app::application_base::stop();
		stopping_.store(true);
		simulator_thread_.join();
	}

	void update() override
	{
		if (state().key_down.contains(app::keycode::ESCAPE))
			stop();

		if (state().key_down.contains(app::keycode::RIGHT) || simulation_time_ >= display_simulation_time)
			get_next_display_creatures();

		float const frame_dt = frame_clock_.restart().count();

		if (!paused_) physics_lag_ += frame_dt;

		while (physics_lag_ >= sim_dt)
		{
			physics_lag_ -= sim_dt;
			simulation_time_ += sim_dt;
			for (auto & creature : display_creatures_)
			{
				creature.update(sim_dt);
				creature.collide(map_, sim_dt);
			}
		}

		float max_view_x = 0.f;
		for (int c = 0; c < display_creatures_.size(); ++c)
			for (auto const & p : display_creatures_[c].positions)
				geom::make_max(max_view_x, p[0]);
		max_view_x = std::max(max_view_x + 10.f, max_view_x_);

		max_view_x_ += (max_view_x - max_view_x_) * (1.f - std::exp(- 4.f * frame_dt));
	}

	void present() override
	{
		gl::ClearColor(0.5f, 0.6f, 0.8f, 0.f);
		gl::Clear(gl::COLOR_BUFFER_BIT);

		for (int c = 0; c < display_creatures_.size(); ++c)
		{
			float aspect_ratio = state().size[0] * 1.f / state().size[1] * display_creatures_.size();
			geom::box<float, 2> view_box = {{{0.f, 0.f}, {0.f, 0.f}}};
			view_box[0] = {max_view_x_ - display_width, max_view_x_};
			view_box[1] = {0.f, view_box[0].length() / aspect_ratio};
			view_box[1] -= 2.f;

			geom::box<float, 2> ground_box;
			ground_box[0] = geom::interval<float>::full();
			ground_box[1] = {geom::limits<float>::min(), 0.f};
			ground_box = ground_box & view_box;
			gfx::color_rgba ground_color{127, 91, 65, 255};
			painter_.rect(ground_box, ground_color);

			gfx::color_rgba hills_color{91, 127, 65, 255};
			for (int i = 0; i + 1 < map_.ground_points.size(); ++i)
			{
				float x0 = map_.ground_points[i    ][0];
				float x1 = map_.ground_points[i + 1][0];

				float y0 = map_.ground_points[i    ][1];
				float y1 = map_.ground_points[i + 1][1];

				painter_.triangle({x0, 0.f}, {x1, 0.f}, {x1, y1}, hills_color);
				painter_.triangle({x0, 0.f}, {x1, y1}, {x0, y0}, hills_color);
			}

			int ruler_min = std::ceil(view_box[0].min / 5.f);
			int ruler_max = std::floor(view_box[0].max / 5.f);

			for (int r = ruler_min; r <= ruler_max; ++r)
			{
				float x = r * 5.f;
				float y = ((r % 2) == 0) ? -1.f : -0.5f;
				y = std::max(y, view_box[1].min);

				float scale = 2.f * view_box[0].length() / state().size[0];

				painter_.line({x, 0.f}, {x, y}, 0.125f, {255, 255, 255, 255}, false);
				painter_.text3d({x + 0.1f, -0.5f, 0.f}, util::to_string(r * 5), {.scale = scale, .x = gfx::painter::x_align::left, .c = {255, 255, 255, 127}}, geom::scale<float, 3>({1.f, -1.f, 1.f}).linear_matrix());
			}

			draw(painter_, display_creatures_[c]);

			float height = view_box[1].length();
			view_box[1].max += c * height;
			view_box[1].min = view_box[1].max - display_creatures_.size() * height;
			painter_.render(geom::orthographic_camera(view_box).transform());
		}

		int text_row = 0;
		auto put_line = [&](std::string const & line)
		{
			painter_.text({13.f, 10.f + text_row * 24.f}, line, {.scale = 2.f, .x = gfx::painter::x_align::left, .y = gfx::painter::y_align::top, .c = {0, 0, 0, 255}});
			painter_.text({12.f,  9.f + text_row * 24.f}, line, {.scale = 2.f, .x = gfx::painter::x_align::left, .y = gfx::painter::y_align::top, .c = {255, 255, 255, 255}});
			++text_row;
		};

		int simulated_generation = 0;
		float best_score = 0.f;
		int best_generation = 0;
		float avg_score = 0.f;
		{
			std::lock_guard lock{next_display_mutex_};
			simulated_generation = next_display_generation_ + 1;
			best_score = generation_best_score_;
			best_generation = generation_best_generation_;
			avg_score = generation_avg_score_;
		}

		put_line(util::to_string("Simulating generation: ", simulated_generation));
		put_line(util::to_string("Best score: ", best_score, " (gen ", best_generation, ")"));
		put_line(util::to_string("Average score: ", avg_score));
		put_line(util::to_string("Display generation: ", display_generation_));
		put_line(util::to_string("Time: ", simulation_time_));

		for (int c = 0; c < display_creatures_.size(); ++c)
		{
			float x = state().size[0] - 12.f;
			float y = (c * state().size[1] * 1.f) / display_creatures_.size() + 9.f;

			auto text = util::to_string("Gen ", display_creatures_[c].generation);
			painter_.text({x + 1.f, y + 1.f}, text, {.scale = 2.f, .x = gfx::painter::x_align::right, .y = gfx::painter::y_align::top, .c = {0, 0, 0, 255}});
			painter_.text({x, y}, text, {.scale = 2.f, .x = gfx::painter::x_align::right, .y = gfx::painter::y_align::top, .c = {255, 255, 255, 255}});

			if (display_creatures_[c].score)
			{
				y += 24.f;
				auto text = util::to_string("Score ", *display_creatures_[c].score);
				painter_.text({x + 1.f, y + 1.f}, text, {.scale = 2.f, .x = gfx::painter::x_align::right, .y = gfx::painter::y_align::top, .c = {0, 0, 0, 255}});
				painter_.text({x, y}, text, {.scale = 2.f, .x = gfx::painter::x_align::right, .y = gfx::painter::y_align::top, .c = {255, 255, 255, 255}});
			}
		}

		for (int i = 1; i < display_creature_count; ++i)
		{
			float y = (state().size[1] * i * 1.f) / display_creature_count;
			painter_.line({0.f, y}, {state().size[0], y}, 4.f, {0, 0, 0, 255}, false);
		}

		painter_.render(geom::window_camera(state().size[0], state().size[1]).transform());
	}

private:
	gfx::painter painter_;

	util::clock<> frame_clock_;

	map map_;

	util::thread simulator_thread_;
	std::vector<creature> population_;

	std::mutex next_display_mutex_;
	std::vector<creature> next_display_creatures_;
	bool has_next_creatures_ = false;
	int next_display_generation_ = 0;
	float generation_best_score_ = 0.f;
	float generation_avg_score_ = 0.f;
	int generation_best_generation_ = 0;

	float simulation_time_ = 0.f;
	float physics_lag_ = 0.f;
	bool paused_ = false;
	std::vector<creature> display_creatures_;
	int display_generation_ = 0;
	float const display_width = 70.f;
	float max_view_x_ = display_width - 5.f;

	std::atomic<bool> stopping_{false};

	void get_next_display_creatures()
	{
		bool new_creatures = false;

		std::lock_guard lock{next_display_mutex_};
		if (has_next_creatures_)
		{
			display_creatures_ = std::move(next_display_creatures_);
			physics_lag_ = 0.f;
			simulation_time_ = 0.f;
			has_next_creatures_ = false;
			display_generation_ = next_display_generation_;
			max_view_x_ = display_width - 5.f;
			new_creatures = true;
		}

		if (new_creatures)
		{
			for (auto & creature : display_creatures_)
			{
				auto center = creature.center();
				float radius = -std::numeric_limits<float>::infinity();
				for (auto const & pos : creature.positions)
					geom::make_max(radius, geom::distance(pos, center));

				float angle = 0.f;

				for (auto & pos : creature.positions)
					pos = geom::vector{0.f, radius} + center + geom::rotate(pos - center, angle);
			}
		}
	}

	void run_simulation()
	{
		static thread_local random::generator rng{random::device{}};
		random::normal_distribution<float> randn;

		int generation = 0;

		while (!stopping_)
		{
			int const test_iterations = std::round(simulation_time / sim_dt);

			++generation;

			std::vector<creature> evaluated_creatures(population_.size());

			#pragma omp parallel for
			for (int i = 0; i < population_.size(); ++i)
			{
				auto & creature = population_[i];
				if (!creature.score)
				{
					auto sim_creature = creature;

					float muscle_power = 0.f;
					for (auto const & cell : sim_creature.cells)
						if (cell.info.type == cell_type::muscle)
							muscle_power += std::pow(cell.info.period, -1.f);

					auto center = sim_creature.center();
					float radius = -std::numeric_limits<float>::infinity();
					for (auto const & pos : sim_creature.positions)
						geom::make_max(radius, geom::distance(pos, center));

					auto angle = random::uniform<float>(rng, -1.f, 1.f) * geom::rad(30.f);
					angle = 0.f;

					for (auto & pos : sim_creature.positions)
						pos = geom::vector{0.f, radius} + center + geom::rotate(pos - center, angle);

					auto start_pos = sim_creature.center();
					float max_y = start_pos[1];
					float sum_speed_x = 0.f;
					float sum_speed_x2 = 0.f;

					for (int iteration = 0; iteration < test_iterations; ++iteration)
					{
						sim_creature.update(sim_dt);
						sim_creature.collide(map_, sim_dt);
						geom::make_max(max_y, sim_creature.center()[1]);

						auto speed = sim_creature.center_velocity();
						sum_speed_x += speed[0];
						sum_speed_x2 += speed[0] * speed[0];
					}
					auto end_pos = sim_creature.center();

					float avg_speed = sum_speed_x / test_iterations;
					float var_speed = sum_speed_x2 / test_iterations - avg_speed * avg_speed;

					float total_speed = (end_pos[0] - start_pos[0]) / simulation_time;

					(void)total_speed;
					(void)avg_speed;
					(void)var_speed;
					(void)muscle_power;

					float score = muscle_power > 0.f ? (total_speed) : 0.f;
//					float score = muscle_power > 0.f ? (total_speed / std::sqrt(muscle_power)) : 0.f;
//					float score = muscle_power > 0.f ? (total_speed * (sim_creature.cells.size())) : 0.f;
//					float score = muscle_power > 0.f ? (total_speed / std::sqrt(muscle_power) * std::sqrt(sim_creature.cells.size())) : 0.f;
	//				float score = muscle_power > 0.f ? (total_speed / muscle_power / std::sqrt(sim_creature.cells.size())) : 0.f;

					if (sim_creature.dead())
						score = 0.f;

					creature.score = score;
				}
				evaluated_creatures[i] = std::move(creature);
			}

			if (stopping_)
				break;

			int const best = population_.size() / 4;
			int const children_count = population_.size() - best;

			population_.clear();

			std::partial_sort(evaluated_creatures.begin(), evaluated_creatures.begin() + best, evaluated_creatures.end(), [](auto const & a, auto const & b){ return *a.score > *b.score; });
			evaluated_creatures.resize(best);
			population_ = std::move(evaluated_creatures);

			float best_score = *population_.front().score;
			int best_generation = population_.front().generation;

			float avg_score = 0.f;
			for (auto const & p : population_)
				avg_score += *p.score;
			avg_score /= population_.size();

			float const mutation_boost = 1.f; //std::exp2(std::min(5.f, (generation - best_generation) / 10.f));

			float const change_type_probability = 1.f / 8.f * mutation_boost;
			float const erase_probability = 1.f / 8.f * mutation_boost;
			float const grow_probability = 1.f / 8.f * mutation_boost;
			float const sex_probability = 0.9f * mutation_boost;

			std::vector<creature> next_display;
			for (int i = 0; i < display_creature_count; ++i)
//				next_display.push_back(creatures_with_score[i * creatures_with_score.size() / display_creature_count].first);
				next_display.push_back(population_[i]);

			std::vector<float> scores(population_.size());
			for (int i = 0; i < population_.size(); ++i)
				scores[i] = *population_[i].score;

			random::weighted_distribution<float> random_parent(std::move(scores));

			for (int i = 0; i < children_count; ++i)
			{
				auto const & parent1 = population_[random_parent(rng)];
				auto const & parent2 = population_[random_parent(rng)];

				creature_blueprint blueprint = parent1.blueprint;

				if (random::uniform<float>(rng) < sex_probability)
				{
					for (auto const & cell : parent2.blueprint)
						if (random::uniform<int>(rng, 0, 1) == 0)
							blueprint[cell.first] = cell.second;
				}

				std::vector<geom::point<int, 2>> erase_cells;
				for (auto const & cell : blueprint)
				{
					if (random::uniform<float>(rng) < erase_probability)
						erase_cells.push_back(cell.first);
				}

				for (auto cell : erase_cells)
					blueprint.erase(cell);

				std::vector<std::pair<geom::point<int, 2>, cell_info>> grow_cells;
				for (auto const & cell : blueprint)
				{
					for (auto n : neighbours)
					{
						auto ncell = cell.first + n;

						if (blueprint.contains(ncell))
							continue;

						if (random::uniform<float>(rng) < grow_probability)
							grow_cells.push_back({ncell, random_cell(rng)});
					}
				}

				for (auto cell : grow_cells)
					blueprint.insert(cell);

				creature new_creature;

				{
					creature_builder builder(std::move(blueprint));
					if (builder.connected())
						new_creature = builder.build(generation);
					else
						new_creature = parent1;
				}

				for (auto & cell : new_creature.cells)
				{
					if (random::uniform<float>(rng) < change_type_probability)
					{
						cell.info = random_cell(rng);
						new_creature.score = std::nullopt;
						new_creature.generation = generation;
					}

					if (cell.info.type == cell_type::muscle)
					{
						cell.info.size += randn(rng) * 0.1f * mutation_boost;
						cell.info.phase += randn(rng) * 0.1f * mutation_boost;
						cell.info.period += randn(rng) * 0.1f * mutation_boost;

						cell.info.size = geom::clamp(cell.info.size, cell_size_range);
						cell.info.phase = std::fmod(cell.info.phase, 1.f);
						cell.info.period = geom::clamp(cell.info.period, cell_period_range);

						new_creature.score = std::nullopt;
						new_creature.generation = generation;
					}
				}

				population_.push_back(std::move(new_creature));
			}

			std::lock_guard lock{next_display_mutex_};
			next_display_creatures_ = std::move(next_display);
			next_display_generation_ = generation;
			has_next_creatures_ = true;
			generation_best_score_ = best_score;
			generation_avg_score_ = avg_score;
			generation_best_generation_ = best_generation;
		}
	}
};

namespace psemek::app
{

	std::unique_ptr<application::factory> make_application_factory()
	{
		return default_application_factory<soft_creatures_2d_app>({.name = "Soft-body creatures", .multisampling = 4});
	}

}