#define KBLIB_DEF_MACROS 1
#include "kblib/convert.h"
#include "kblib/iterators.h"
#include "kblib/memory.h"
#include "kblib/stringops.h"

#include "map.h"
#include "olsennoise.hpp"
#include "turbo_colormap.hpp"

#include <algorithm>
#include <chrono>
#include <cmath>
#include <complex>
#include <cstdint>
#include <deque>
#include <fstream>
#include <iostream>
#include <regex>
#include <sstream>

using namespace std::chrono;
using std::cerr;
using std::cin;
using std::clog;
using std::cout;
using std::endl;
using std::string;
using std::vector;

static constexpr inline auto plateau(const float r) -> float {
	// $ plateau(r) = r^2 / (r^2+1) $
	return r * r / (r * r + 1);
}

template <typename T>
static constexpr inline auto square(const T x) {
	return x * x;
}

template <typename T>
constexpr static decltype(auto) do_atan2(T x, T y) {
	using std::atan2;
	return atan2(x, y);
};

template <typename T1, typename T2, typename T1_S = typename T1::Scalar,
          typename T2_S = typename T2::Scalar>
constexpr decltype(auto) atan2(const T1& x, const T2& y) {
	return Eigen::CwiseBinaryOp<
	    auto (*)(T1_S, T2_S)->std::common_type_t<T1_S, T2_S>, T1, T2>{x, y,
	                                                                  do_atan2};
}

template <typename T1, typename T2>
static constexpr inline decltype(auto) polar_t(const T1 x, const T2 y) {
	using ::atan2, std::atan2;
	// $ polar\_t(x, y) = \mathrm{atan2}(x, y) = \arg(x + iy) $
	return atan2(x, y);
}

template <typename T1, typename T2>
static constexpr inline decltype(auto) polar_r(const T1 x, const T2 y) {
	using std::sqrt;
	// $ polar\_r(x,y) = \sqrt{x^2 + y^2}\cdot\sqrt{0.1} $
	// I have no idea why I thought this scale factor was a good idea, it should
	// clearly have gone into the spiral equations instead
	constexpr auto scale_factor = 0.3162277660168379332f; // sqrt(.1)
	return sqrt(x * x + y * y) * scale_factor;
}

template <typename T>
static constexpr inline auto to_polar(const std::complex<T> p)
    -> std::complex<T> {
	// $ to\_polar(p) = |\sqrt{0.1}p| + \arg(p) $
	return {polar_r(p.real, p.imag), polar_t(p.real(), p.imag())};
}

template <typename T>
static constexpr inline auto polar(const T x, const T y) -> std::complex<T> {
	// $ polar(x, y) = \sqrt{.1 x^2 \cdot .1 y^2} + \arg(x + iy) $
	return {polar_r(x, y), polar_t(x, y)};
}

template <typename T>
static constexpr inline auto log_r(const T x, const T y) -> T {
	return .5 * log(x * x + y * y);
}

#if 0
// Some constants
constexpr float
    pi = kblib::pi<float>(),
    // Spiral term (shared)
    spiral_strength = 80, // How strong the spiral motif is
    // Spiral term (archimedian)
    spiral_tightness = 0.05f, // turning pitch of spiral
    spiral_sharpness = 1.40f, // how much to exaggerate peaks by (exponential)
    spiral_negative_bias = 1.2f, // how much to subtract from spiral term
    spiral_flatness = 1 / 1.25f, // how much to dampen peaks by (linear)
    // Plateau term
    plt_radius = 100, // Radius of the central plateau in ADU
    // plt_height = 24, //How high the central plateau is (logarithmic)
    plt_height_a = 48,   // How high the central plateau is (archimedian)
    plt_cliffiness = 20, // How smooth the slope to the central plateau is
    plt_cliffiness_x = 1 / plt_cliffiness,
    // plt_prescale = 1.15, //Scales the plateau term before normalization
    plt_flatness = 2.5f, // How much to flatten the central plateau by
    // plt_l10n_adj = plt_height_a/35-pi_down, //Reduces the large-scale
    // effect of the plateau term
    plateau_hinge = 120,
    // plt_l10n_adj2 = plt_l10n_adj/24,
    // Noise term
    noise_hinge = 155,   // Hinge which noise will be scaled away from
    noise_scale = 1.20f, // How much to scale noise by
    // Adjustment
    // sea_level = 135,
    sea_adjustment =
        noise_hinge - 8; // Overall factor to adjust the map height by
const float
    // Spiral term (log)
    spiral_curviness = 1.12f, // Base of the logarithm used for the spiral
    spiral_curviness_x = 1 / std::log(spiral_curviness);
#endif

auto getline_(std::istream& is, std::string& str) -> std::istream& {
	if (is.rdstate() == std::ios::eofbit) {
		is.setstate(std::ios::failbit);
	}
	if (is) {
		auto excflag = is.exceptions();
		is.exceptions(std::ios::goodbit);
		std::getline(is, str);
		is.clear(is.rdstate() & ~std::ios::failbit);
		is.exceptions(excflag);
	}
	return is;
}

auto operator>>(std::istream& is, params& p) -> std::istream& {
	std::string line;
	while (not is.eof() and getline_(is, line)) {
		std::istrstream ent(line.c_str(), line.length());
		// ent.exceptions(std::ios::failbit | std::ios::badbit);
		if (not (ent >> std::ws).eof() and ent.peek() != '#') {
			std::string name;
			try {
				ent >> name;

				const auto& r = params::entries.at(name);
				float val;
				if (ent >> kblib::expect('=') >> val) {
					p.*(r.var) = val;
				} else {
					std::cerr << "invalid parameter definition "
					          << kblib::quoted(line) << '\n';
					is.setstate(std::ios::failbit);
					return is;
				}
			} catch (std::out_of_range&) {
				if (is.exceptions() & std::ios::failbit) {
					is.setstate(std::ios::failbit);
					throw std::ios::failure(kblib::concat(
					    kblib::quoted(name), " is not a known parameter"));
				} else {
					std::cerr << "Unknown parameter: " << kblib::quoted(name)
					          << '\n';
					is.setstate(std::ios::failbit);
					return is;
				}
			} catch (std::ios::failure&) {
				if (is.exceptions() & std::ios::failbit) {
					is.setstate(std::ios::failbit);
					throw;
				} else {
					is.setstate(std::ios::failbit);
					return is;
				}
			}
		} else {
			;
		}
	}
	return is;
}

auto operator<<(std::ostream& os, params& p) -> std::ostream& {
	for (const auto& [name, ent] : params::entries) {
		os << name << " = " << p.*(ent.var) << '\n';
	}
	return os;
}

auto params::write_verbose(std::ostream& os) -> std::ostream& {
	for (const auto& [name, ent] : entries) {
		if (not ent.desc.empty()) {
			os << "# " << ent.desc << '\n';
		}
		os << name << " = " << this->*(ent.var) << '\n';
	}
	return os;
}

auto genSpiral(const bool archimedian, const params& p,
               vector<vector<float>>& buf, const vector<float>& x_vals,
               const vector<float>& y_vals) -> vector<vector<float>>& {

	// vector<vector<std::pair<float,float>>> pc(x_vals.size(),
	// vector<std::pair<float,float>>(y_vals.size()));
	// ts.push_back({steady_clock::now(),"Process setup"});
	for (size_t x = 0; x < x_vals.size(); ++x) {
		for (size_t y = 0; y < y_vals.size(); ++y) {
			auto r = polar_r(x_vals[x], y_vals[y]),
			     t = polar_t(x_vals[x], y_vals[y]);
			// float r = sqrt(x_vals[x]*x_vals[x]*.1f+y_vals[y]*y_vals[y]*.1f);
			// float t = atan2(x_vals[x],y_vals[y]);

			// $ buf_{x,y} = \left( buf_{x,y} - N_h \right) \cdot N_s +
			// sea\_adjustment $
			buf[x][y]
			    = (buf[x][y] - p.noise_hinge) * p.noise_scale + p.sea_adjustment();
			if (archimedian) {
				// $ buf_{x,y} = buf_{x,y} + \left(\left(\sin(t + S_t * r)
				// + 1.f\right)^{S_s}\cdot S_f - S_{nb}\right) $
				buf[x][y]
				    += p.spiral_strength
				       * (std::pow((std::sin(t + p.spiral_tightness * r) + 1.f),
				                   p.spiral_sharpness)
				              * p.spiral_flatness
				          - p.spiral_negative_bias);
			} else {
				// $ buf_{x,y} = buf_{x,y} + S_s\cdot\left(\left(\sin(t + S_{cx} *
				// \log(r)) + 1.f\right)^{S_s}\cdot S_f - S_{nb}\right) $
				buf[x][y]
				    += p.spiral_strength
				       * float(
				           std::pow((std::sin(t + p.spiral_curviness_x() * log(r))
				                     + 1.f),
				                    p.spiral_sharpness)
				               * p.spiral_flatness
				           - p.spiral_negative_bias);
			}
			// $$ {\tan^{-1}(P_{cx} \cdot (r - P_r)) \over \pi} + 0.5 $$
			auto plt
			    = std::atan(p.plt_cliffiness_x() * (r - p.plt_radius)) / (p.pi)
			      + 0.5f;
			// $$ buf_{x,y} = {(buf_{x,y} - P_h) \over (1 + P_f^2) - (P_f * plt)^2}
			// + P_{ha} \cdot (1 - plt) + P_h $$
			buf[x][y] = (buf[x][y] - p.plateau_hinge)
			                / ((1.f + square(p.plt_flatness))
			                   - square(p.plt_flatness * plt))
			            + p.plt_height_a * (1.f - plt) + p.plateau_hinge;
		}
	}
	return buf;
}

auto genSpiral(bool archimedian, const params& p, const row_array& x_vals,
               const col_array& y_vals) -> buffer_type {
	buffer_type buf(x_vals.size(), y_vals.size());
	// for (auto [row, yv] : kblib::zip(rows(buf), traverse(y_vals))) {
	for (const auto rx : kblib::range(y_vals.size())) {
		auto row = buf.row(rx);
		auto yv = y_vals[rx];
		const auto row_r
		    = x_vals.unaryExpr([yv = yv](float xv) { return polar_r(xv, yv); })
		          .eval();
		const auto row_t
		    = x_vals.unaryExpr([yv = yv](float xv) { return polar_t(xv, yv); })
		          .eval();

		if (archimedian) {
			// $ row = row + \left(\left(\sin(t + S_t * r) + 1.f\right)^{S_s}\cdot
			// S_f - S_{nb}\right) $
			row += p.spiral_strength
			       * (pow((sin(row_t + p.spiral_tightness * row_r) + 1.f),
			              p.spiral_sharpness)
			              * p.spiral_flatness
			          - p.spiral_negative_bias);
		} else {
			// $ row = row + S_s\cdot\left(\left(\sin(t + S_{cx} * \log(r))
			// + 1.f\right)^{S_s}\cdot S_f - S_{nb}\right) $
			row += p.spiral_strength
			       * (pow((sin(row_t + p.spiral_curviness_x() * log(row_r)) + 1.f),
			              p.spiral_sharpness)
			              * p.spiral_flatness
			          - p.spiral_negative_bias);
		}
		// $$ {\tan^{-1}(P_{cx} \cdot (r - P_r)) \over \pi} + 0.5 $$
		auto plt
		    = (p.plt_cliffiness_x() * (row_r - p.plt_radius)).atan() / params::pi
		      + 0.5f;
		// $$ row = {(row - P_h) \over (1 + P_f^2) - (P_f * plt)^2} + P_{ha} \cdot
		// (1 - plt) + P_h $$
		row = (row - p.plateau_hinge)
		          / ((1.f + square(p.plt_flatness))
		             - square(p.plt_flatness * plt))
		      + p.plt_height_a * (1.f - plt) + p.plateau_hinge;
	}
	return buf;
}
static inline auto absLess(float a, float b) -> bool {
	return std::abs(a) < std::abs(b);
}

template <typename T>
auto deref(std::pair<T, T> p) {
	return std::make_tuple(*p.first, *p.second);
}

template <typename T, typename F>
auto apply(const std::valarray<T>& v, F f) {
	using R = std::invoke_result_t<F, T>;
	std::valarray<R> ret(v.size());
	std::transform(begin(v), end(v), begin(ret), std::ref(f));
	return ret;
}

auto genWeather(vector<vector<float>>& buf, const vector<float>& x_vals,
                const vector<float>& y_vals) -> vector<vector<float>>& {
	assert(x_vals.size() == y_vals.size() and not x_vals.empty());

	double minx, miny, maxx, maxy;

	std::tie(minx, maxx)
	    = deref(std::minmax_element(x_vals.begin(), x_vals.end(), absLess));
	std::tie(miny, maxy)
	    = deref(std::minmax_element(y_vals.begin(), y_vals.end(), absLess));

	double minr{sqrt(square(minx) + square(miny))},
	    maxr{sqrt(square(maxx) + square(maxy)) + 1};

	clog << minr << ", " << maxr << endl;

	auto tempBuf = OlsenNoise::genNoise(0, floor(minr), 1, ceil(maxr - minr))[0];
	auto rainBuf
	    = OlsenNoise::genNoise(1024, floor(minr), 1, ceil(maxr - minr))[0];

	for (size_t x = 0; x < x_vals.size(); ++x) {
		for (size_t y = 0; y < y_vals.size(); ++y) {
			double r = std::sqrt(x_vals[x] * x_vals[x] + y_vals[y] * y_vals[y]);

			buf[x][y]
			    = rainBuf[round(r - minr)] / 1.2f - 60 * atan(r / 1000 - .2) + 40;
		}
	}
	return buf;
}

auto point(float x, float y) -> point_t {
	return (point_t{} << x, y).finished();
}

auto genWeather(const row_array& x_vals, const col_array& y_vals)
    -> buffer_type {

	assert(x_vals.size() == y_vals.size() and x_vals.size() > 0);

	buffer_type r = (x_vals.replicate(y_vals.size(), 1).square()
	                 + y_vals.replicate(1, x_vals.size()).square())
	                    .sqrt();
	float minr = floorf(r.minCoeff());
	float maxr = ceilf(r.maxCoeff()) + 1;
	/*
	   const auto minr =
	       (Eigen::AlignedBox2f(point(x_vals[0], y_vals[0]),
	                            point(x_vals.tail<1>()[0], y_vals.tail<1>()[0]))
	            .contains(point_t::Zero()))
	           ? 0.f
	           : polar_r(x_vals.abs().minCoeff(), y_vals.abs().minCoeff());
	   auto maxr = polar_r(x_vals.abs().maxCoeff(), y_vals.abs().maxCoeff());
	*/
	row_array tempBuf
	    = OlsenNoise::genNoise(0, floorf(minr), 1, ceilf(maxr - minr))[0]
	          .cast<float>();
	row_array rainBuf
	    = OlsenNoise::genNoise(1024, floorf(minr), 1, ceilf(maxr - minr))[0]
	          .cast<float>();

	buffer_type c_rain;
	c_rain.resizeLike(r);

	assert(r.minCoeff() >= minr);
	assert(r.maxCoeff() <= maxr);

	for (Eigen::Index y = 0; y < y_vals.size(); ++y) {
		for (Eigen::Index x = 0; x < x_vals.size(); ++x) {
			auto radius = r(x, y);
			auto i = lrintf(radius - minr);
			c_rain(x, y) = rainBuf(i);
		}
	}
	return c_rain / 1.2f - 60.f * atan(r / 1000.f - .2f) + 40.f;

	/*
	for (size_t y = 0; y < y_vals.size(); ++y) {
	   for (size_t x = 0; x < x_vals.size(); ++x) {
	      double r = polar_r(x_vals[x], y_vals[y]);

	      buf[y][x] =
	          rainBuf[round(r - minr)] / 1.2f - 60 * atan(r / 1000.f - .2f) + 40;
	   }
	}//*/
	/*
	for (const auto [row, yv] : kblib::zip(buf, y_vals)) {
	   const auto row_r = kblib::build_dy<std::valarray<float>>(
	       x_vals, [&yv = yv](auto xv) { return polar_r(xv, yv); });

	   const auto row_rain = kblib::build_copy<std::valarray<float>>(
	       rainBuf[apply(row_r, [=](float f) {
	          return static_cast<std::size_t>(roundf(f - minr));
	       })]);

	   row = row_rain / 1.2f - 60 * atan(row_r / 1000.f - .2f) + 40;
	}
	return buf;//*/
}

auto color(std::uint8_t r, std::uint8_t g, std::uint8_t b) -> pnm::rgb_pixel {
	return {r, g, b};
}

[[deprecated]] auto genColorMap(int type) -> palette {
	palette ret(256);
	if (type < 0 or kblib::to_unsigned(type) >= current_cmap.size()) {
		return ret;
	}
	switch (type) {
	case 0: // Greyscale
		for (unsigned char x : kblib::range(256)) {
			ret[x] = color(x, x, x);
		}
		break;
	case 1: // Natural
		for (unsigned char x : kblib::range(256)) {
			if (x < 132) { // ocean
				ret[x] = color(
				    0,
				    x / 3
				        + 4 * static_cast<char>(static_cast<double>(x) * x / 768.),
				    x + 70);
			} else if (x < 136) { // beach
				ret[x] = color(240, 240, 64);
			} else if (x < 200) { // grass, forests
				ret[x] = color(48 - (x - 168), 180 - (x - 168), (x - 100) / 2);
			} else if (x < 248) { // Prevent underflow of red
				ret[x] = color(0, 180 - (x - 168), (x - 160));
			} else if (x < 252) { // snowy forest
				ret[x] = color(x - 64, 248, x - 24);
			} else {                                  // snow
				ret[x] = color(x - 24, x - 12, x - 8); // Make snow bluish
			}
		}
		break;
	case 2: // Arctic
		// Broken
		for (unsigned char x : kblib::range(256)) {
			if (x < 149) {
				ret[x] = color(0, 0, x + 77);
			} else if (x == 149) {
				ret[x] = color(0, 128, 255);
			} else if (x < 155) {
				ret[x] = color(240 - 3 * (x - 155), 240, 64 + 4 * (x - 155));
			} else if (x < 162) {
				ret[x] = color(32 - (x - 162), 180 - (x - 162), 0 + (x - 162));
			} else if (x < 210) {
				ret[x] = color(0, 128, 32);
			} else {
				ret[x] = color(x - 4, x - 2, x); // Make snow bluish
			}
		}
		break;
	case 3: // Alien
		// Broken
		for (unsigned char x : kblib::range(256)) {
			if (x < 149) { // 0+
				ret[x] = color(0, 0, x + 77);
			} else if (x == 149) { // 149+
				ret[x] = color(0, 128, 255);
			} else if (x < 162) { // 150+
				ret[x] = color(240 - 3 * (x - 155), 240, 64 + 4 * (x - 155));
			} else if (x < 182) { // 162+
				ret[x] = color(32 - (x - 162), 180 - (x - 162), (x - 162));
			} else if (x < 210) {
				ret[x] = color(0, 128, 32);
			} else {                            // 182+
				ret[x] = color(x - 4, x - 2, x); // Make snow bluish
			}
		}
		break;
	case 4: // Purple, red (Natural * RGB=>GRB)
		for (unsigned char x : kblib::range(256)) {
			ret[x] = color(current_cmap[1].second[x].green(),
			               current_cmap[1].second[x].red(),
			               current_cmap[1].second[x].blue());
		}
		break;
	case 5: // Water ratio map
		for (unsigned char x : kblib::range(256)) {
			if (x < 132) { // ocean
				ret[x] = color(0, 0, 100);
			} else { // land
				ret[x] = color(48, 180, 48);
			}
		}
		break;
	case 6: // Elevation
		for (unsigned char x : kblib::range(256)) {
			if (x >= 132 && x <= 135) { // shore
				ret[x] = color(64, 64, 64);
			} else if (x < 132 && ((x - 132) % 10 == 0)) { // ocean depth marks
				ret[x] = color(128, 128, 128);
			} else if (x > 135 && ((x - 135) % 10 == 0)) { // land height marks
				ret[x] = color(192, 192, 192);
			} else {
				ret[x] = color(255, 255, 255);
			}
		}
		break;
	case 7: // Elevation-color
		for (unsigned char x : kblib::range(256)) {
			if (x >= 132 && x <= 135) { // shore
				ret[x] = color(64, 64, 64);
			} else if (x < 132 && ((x - 132) % 10 == 0)) { // ocean depth marks
				ret[x] = color(128, 128, 128);
			} else if (x > 135 && ((x - 135) % 10 == 0)) { // land height marks
				ret[x] = color(192, 192, 192);
			} else {
				ret[x] = color(255, 255, 255);
			}
		}
		break;
	case 8: // Turbo
		return {std::begin(turbo_srgb_bytes), std::end(turbo_srgb_bytes)};
		break;
	case 15: // Hash
	{
		unsigned h = -1u;
		for (unsigned char x : kblib::range(256)) {
			h = kblib::FNV32a_s(reinterpret_cast<const char*>(&h), sizeof(h));
			ret[x] = color(h & 255, (h >> 8) & 255, (h >> 16) & 255);
		}
	} break;
	default: // Fill cmap.h yourself and rebuild to get new maps
		return current_cmap[type].second;
		// for (std::uint16_t x = 0; x < 256; ++x) {
		// ret[x] =
		// color(cmap[type][x][0],cmap[type][x][1],cmap[type][x][2]);
		// }
		break;
	}
	return ret;
}

palette sRGB_to_cmap(gsl::span<const sRGB_tuple> input) noexcept {
	palette out;
	std::transform(input.begin(), input.end(), std::back_inserter(out),
	               [](const sRGB_tuple& tup) {
		               std::cerr << tup.red() << ',' << tup.green() << ','
		                         << tup.blue() << '\n';
		               return pnm::rgb_pixel{
		                   static_cast<std::uint8_t>(std::round(tup[0] * 255)),
		                   static_cast<std::uint8_t>(std::round(tup[1] * 255)),
		                   static_cast<std::uint8_t>(std::round(tup[2] * 255))};
	               });
	return out;
}