/*

    ......... 2016 Ivan Mahonin

    This program is free software: you can redistribute it and/or modify

    it under the terms of the GNU General Public License as published by

    the Free Software Foundation, either version 3 of the License, or

    (at your option) any later version.

    This program is distributed in the hope that it will be useful,

    but WITHOUT ANY WARRANTY; without even the implied warranty of

    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License

    along with this program.  If not, see <http: licenses="" www.gnu.org="">.</http:>

*/

#include <cstdarg></cstdarg>

#include <cstdio></cstdio>

#include <cstring></cstring>

#include <cmath></cmath>

#include <iostream></iostream>

#include <iomanip></iomanip>

#include <fstream></fstream>

#include <sstream></sstream>

#include <vector></vector>

#include <set></set>

#include <map></map>

#include "surface.h"

using namespace std;

// options

const double initial_radius = 16.0;

const double min_radius = 1.0;

const double max_radius = 2048;

const double step = 0.1;

const double large_params_step = 0.05;

const double params_step = 0.025;

const double threshold = 0.012;

const char logfile[] = "results/coef.log";

const char final_coefs_file[] = "results/coef.hpp";

const char final_coefs_file2[] = "results/coef2.hpp";

#define COMPARE_HELPER(a) a < other.a ? true : other.a < a ? false :

#define COMPARE(a) a < other.a;

#define COMPARE2(a, b) COMPARE_HELPER(a) COMPARE(b)

#define COMPARE3(a, b, c) COMPARE_HELPER(a) COMPARE2(b, c)

#define COMPARE4(a, b, c, d) COMPARE_HELPER(a) COMPARE3(b, c, d)

const double pi = 3.14159265358979323846;

string strprintf(const char *format, ...)

{

	va_list args;

	va_start(args,format);

	static char buffer[1024*1024];

	vsprintf(buffer, format, args);

	return buffer;

}

template<typename t=""></typename>

inline bool less(const T &a, const T &b, bool def)

	{ return a < b ? true : b < a ? false : def; }

double gauss(double x, double radius) {

	static const double k = 1.0/sqrt(2.0*pi);

	return k/radius*exp(-0.5*x*x/(radius*radius));

}

class Params {

public:

	union {

		struct { double k0, k1, k2; };

		double k[3];

	};

	Params(): k0(), k1(), k2() { }

	Params(double k0, double k1, double k2): k0(k0), k1(k1), k2(k2) { }

	bool operator< (const Params &other) const

		{ return COMPARE3(k0, k1, k2); }

};

class Checker {

private:

	vector<double> blank;</double>

	vector<double> tmp;</double>

public:

	int offset;

	bool valid;

	double radius;

	Checker(double radius): offset(3), valid(), radius(radius) {

		blank.resize((int)ceil(radius*10));

		tmp.resize(blank.size() + 2*offset);

		for(int i = 0; i < (int)blank.size(); ++i)

			blank[i] = gauss(i, radius);

	}

	Params convert_params(const Params ¶ms) {

		double a = 1.0/params.k0;

		double b = a*cos(pi*params.k1);

		double c = 1.0/params.k2;

		Params p;

		p.k0 = (a*a + 2.0*c*b)/(c*a*a);

		p.k1 = -(c + 2.0*b)/(c*a*a);

		p.k2 = 1.0/(c*a*a);

		return p;

	}

	double check2(const Params ¶ms) {

		if ( abs(params.k0) > 0.99999999

		  //|| abs(params.k1) > 0.99999999

		  || abs(params.k2) > 0.99999999

		  || params.k0 < 1e-6

		  || abs(params.k2) < 1e-6 ) return INFINITY;

		return check(convert_params(params));

		//return check(params);

	}

	double check(const Params ¶ms) {

		valid = false;

		//if (mode3a && fabs(params.k0 + (params.k1 + params.k2)/(params.k0 + params.k1 + params.k2)) >= 0.999999999)

		//	return INFINITY;

		if (fabs(params.k0) < 1e-10 && fabs(params.k1) < 1e-10 && fabs(params.k2) < 1e-10)

			return INFINITY;

		double k = 1.0 - params.k0 - params.k1 - params.k2;

		tmp[offset] = k;

		for(int i = offset+1, end = (int)tmp.size()-offset; i < end; ++i)

			tmp[i] = params.k0*tmp[i-1]

				   + params.k1*tmp[i-2]

				   + params.k2*tmp[i-3];

		for(int i = (int)tmp.size()-offset-1, end = offset-1; i > end; --i)

			tmp[i] = k*tmp[i]

				   + params.k0*tmp[i+1]

				   + params.k1*tmp[i+2]

				   + params.k2*tmp[i+3];

		double delta = 0.0;

		//double minv = tmp[offset];

		//double kblank = 1.0/blank[0];

		for(int i = 0; i < (int)blank.size(); ++i) {

			double t = tmp[i + offset];

			double d = blank[i] - t;

			delta += fabs(d*d);

			//if (t < minv) minv = t;

			//if (minv < 0) minv = 0;

			//delta += 2*fabs(t - minv);

			//delta += 0.1*fabs(d/(kblank*blank[i] + 0.1));

		}

		if (isinf(delta) || isnan(delta)) return INFINITY;

		valid = true;

		return sqrt(delta/blank.size())*blank.size();

	}

	void graph(const string &filename) {

		int base_size = 768;

		int count = blank.size();

		int scale = base_size/count + 1;

		int div   = count/768 + 1;

		int pad = max(16, count*scale/div/8);

		int gw = count*scale/div;

		int gh = gw;

		Surface::Color axis = Surface::convert_color(0, 0, 1, 1);

		Surface::Color ca = Surface::convert_color(1, 0, 0, 0.25);

		Surface::Color cb = Surface::convert_color(0, 0, 0, 0.75);

		Surface s(gw + 2*pad, gh + 2*pad);

		s.clear(1.0, 1.0, 1.0, 1.0);

		s.set_color(axis);

		s.move_to(pad, pad);

		s.line_by(0, gh);

		s.line_by(gw, 0);

		s.set_color(ca);

		s.move_to(pad, pad);

		for(int i = 0; i < count; ++i)

			s.line_to(pad + i*scale/div, pad + (int)round((1.0 - blank[i]/blank[0])*gh));

		if (valid) {

			s.set_color(cb);

			s.move_to(pad, pad);

			for(int i = 0; i < count; ++i)

				s.line_to(pad + i*scale/div, pad + (int)round((1.0 - tmp[i + offset]/blank[0])*gh));

		}

		s.save(filename);

	}

};

class Finder {

public:

	class Entry {

	public:

		double value;

		Params min;

		Params current;

		Params max;

		Finder *finder;

		Entry *parent;

		int level;

		int index;

		int max_variants;

		set<entry> variants;</entry>

		Entry(): value(INFINITY), finder(), parent(), level(), index(), max_variants() { }

		bool operator< (const Entry &other) const

			{ return COMPARE4(value, min, current, max); }

		void add_variant(const Entry &entry) {

			if (isinf(entry.value) || isnan(entry.value) || max_variants == 0) return;

			if ((int)variants.size() < max_variants || entry.value < variants.rbegin()->value) {

				variants.insert(entry);

				if (max_variants < (int)variants.size()) {

					set<entry>::iterator i = variants.end();</entry>

					variants.erase(--i);

				}

			}

		}

		void report_index() {

			if (parent) {

				parent->report_index();

				cout << "["  << index << "]";

			}

		}

		void report(bool open) {

			if (level < finder->max_report_level) {

				cout << "level " << level << " ";

				report_index();

				for(int i = level; i < finder->max_report_level; ++i)

					cout << "   ";

				cout << (open ? " begin" : " end  ")

					 << fixed

					 << setprecision(8)

					 << " " << current.k0

					 << " " << current.k1

					 << " " << current.k2

					 << " " << value

					 << " (" << finder->best << ")"

					 << endl;

			}

		}

		double find() {

			if (max_variants <= 0) return value;

			if (value < finder->best) finder->best = value;

			report(true);

			int i[3];

			for(i[0] = 0; i[0] < finder->sub_division; ++i[0])

			for(i[1] = 0; i[1] < finder->sub_division; ++i[1])

			for(i[2] = 0; i[2] < finder->sub_division; ++i[2]) {

				Params p;

				for(int j = 0; j < 3; ++j)

					p.k[j] = (max.k[j] - min.k[j])/finder->sub_division*(i[j] + 0.5) + min.k[j];

				finder->value(i[0], i[1], i[2]) = finder->checker.check2(p);

			}

			finder->clear_zones();

			for(i[0] = 1; i[0] < finder->sub_division-1; ++i[0])

			for(i[1] = 1; i[1] < finder->sub_division-1; ++i[1])

			for(i[2] = 1; i[2] < finder->sub_division-1; ++i[2])

				if (finder->is_local_minimum(i[0], i[1], i[2])) {

					Zone &z = finder->zone_for_value(i[0], i[1], i[2]);

					double &v = finder->value(i[0], i[1], i[2]);

					if (v < z.value) z.set(v, i);

				}

			for(i[0] = 0; i[0] < finder->zone_count; ++i[0])

			for(i[1] = 0; i[1] < finder->zone_count; ++i[1])

			for(i[2] = 0; i[2] < finder->zone_count; ++i[2]) {

				Zone &z = finder->zone(i[0], i[1], i[2]);

				if (!isinf(z.value)) {

					Entry sub;

					sub.finder = finder;

					sub.parent = this;

					sub.level = level + 1;

					sub.max_variants = max_variants/2;

					if (sub.max_variants < 1) sub.max_variants = 1;

					for(int j = 0; j < 3; ++j) {

						sub.min.k[j]     = (max.k[j] - min.k[j])/finder->sub_division*(z.index[j] + 0.5 - finder->zone_size) + min.k[j];

						sub.current.k[j] = (max.k[j] - min.k[j])/finder->sub_division*(z.index[j] + 0.5) + min.k[j];

						sub.max.k[j]     = (max.k[j] - min.k[j])/finder->sub_division*(z.index[j] + 0.5 + finder->zone_size) + min.k[j];

					}

					sub.value = z.value;

					add_variant(sub);

				}

			}

			if (level < finder->max_level) {

				vector<entry> v(variants.begin(), variants.end());</entry>

				variants.clear();

				for(vector<entry>::iterator i = v.begin(); i != v.end(); ++i) {</entry>

					i->index = i - v.begin();

					i->find();

					variants.insert(*i);

				}

			}

			if (!variants.empty() && variants.begin()->value < value) {

				value = variants.begin()->value;

				current = variants.begin()->current;

			}

			if (value < finder->best) finder->best = value;

			report(false);

			return value;

		}

	};

	struct Zone {

		double value;

		int index[3];

		Zone() { clear(); }

		void clear() { value = INFINITY; memset(index, 0, sizeof(index)); }

		void set(double value, int *i) { this->value = value; memcpy(index, i, sizeof(index)); }

	};

	Checker checker;

	int sub_division;

	int zone_count;

	int zone_size;

	int s0;

	int s1;

	int s2;

	int max_level;

	int max_report_level;

	vector<double> values;</double>

	vector<zone> zones;</zone>

	double best;

	Entry root;

	Finder(

		double radius,

		Params min,

		Params max,

		int zone_count,

		int zone_size,

		int max_variants,

		int max_level,

		int max_report_level

	):

		checker(radius),

		sub_division(zone_count*zone_size),

		zone_count(zone_count),

		zone_size(zone_size),

		s0(1),

		s1(sub_division),

		s2(sub_division*sub_division),

		max_level(max_level),

		max_report_level(max_report_level),

		values(sub_division*sub_division*sub_division),

		zones(zone_count*zone_count*zone_count),

		best(INFINITY)

	{

		root.finder = this;

		root.min = min;

		root.max = max;

		root.max_variants = max_variants;

	}

	double find() {

		for(int j = 0; j < 3; ++j)

			root.current.k[j] = 0.5*(root.max.k[j] - root.min.k[j]) + root.min.k[j];

		root.value = checker.check2(root.current);

		best = root.value;

		return root.find();

	}

	double& value(int i0, int i1, int i2)

		{ return values[i0*s0 + i1*s1 + i2*s2]; }

	bool is_local_minimum(int i0, int i1, int i2) {

		double *vp = &value(i0, i1, i2);

		double v = *vp;

		return v < *(vp-s0) && v <= *(vp+s0)

			&& v < *(vp-s1) && v <= *(vp+s1)

			&& v < *(vp-s2) && v <= *(vp+s2);

	}

	Zone& zone_for_value(int i0, int i1, int i2)

		{ return zones[(i2/zone_size*zone_count + i1/zone_size)*zone_count + i0/zone_size]; }

	Zone& zone(int i0, int i1, int i2)

		{ return zones[(i2*zone_count + i1)*zone_count + i0]; }

	void clear_zones() {

		for(vector<zone>::iterator i = zones.begin(); i != zones.end(); ++i)</zone>

			i->clear();

	}

};

double find_near(Checker &checker, double step, Params ¶ms) {

	double best_value = checker.check2(params);

	Params best_params = params;

	while(true) {

		bool found = false;

		for(int i = 0; i < 3; ++i) {

			for(int j = -1; j <= 1; j += 2) {

				Params p = params;

				p.k[i] += j*step;

				double v = checker.check2(p);

				if (v < best_value) { best_value = v; best_params = p; }

			}

		}

		if (found) params = best_params; else break;

	}

	return best_value;

}

double find_near2(Checker &checker, double step, Params ¶ms) {

	Finder finder(

		checker.radius,

		Params(params.k0-step, params.k1-step, params.k2-step),

		Params(params.k0+step, params.k1+step, params.k2+step),

		4,

		4,

		1,

		16,

		1

	);

	finder.find();

	params = finder.root.current;

	return finder.root.value;

}

double find_near3(Checker &checker, double step, Params ¶ms) {

	Finder finder(

		checker.radius,

		Params(params.k0-step, params.k1-step, params.k2-step),

		Params(params.k0+step, params.k1+step, params.k2+step),

		8,

		8,

		1,

		10,

		4

	);

	finder.find();

	params = finder.root.current;

	return finder.root.value;

}

double initial_find(Checker &checker, Params ¶ms) {

	Finder finder(

		checker.radius,

		Params(0, 0, -1),

		Params(1, 1, 1),

		16,

		16,

		1,

		4,

		8

	);

	finder.find();

	params = finder.root.current;

	return finder.root.value;

}

map< double, pair<params, double=""> > coefs;</params,>

double fix_radius(double r)

	{ return round(r/step)*step; }

void log_begin() {

	ofstream f(logfile, ios_base::app);

	f << endl << "double coefs[][5] = {" << endl;

}

void log_params(double radius, const Params ¶ms, double value) {

	ofstream f(logfile, ios_base::app);

	f << "    { " << radius

	  << setprecision(20)

	  << ", "     << params.k0

	  << ", "     << fabs(params.k1)

	  << ", "     << params.k2

	  << ", "     << value

	  << " },"    << endl;

	cout << endl;

	coefs[fix_radius(radius)].first = params;

	coefs[fix_radius(radius)].second = value;

}

void log_params(const Checker &checker, const Params ¶ms, double value) {

	log_params(checker.radius, params, value);

}

void log_params(const Finder &finder) {

	log_params(finder.checker, finder.root.current, finder.root.value);

}

void log_end() {

	ofstream f(logfile, ios_base::app);

	f << "};" << endl << endl;

}

void graph(Checker &checker, const Params ¶ms) {

	int ri = roundf(checker.radius);

	if (fabs(checker.radius - ri) > 1e-5) return;

	bool is_power_of_two = false;

	for(int i = 1; i <= 4096; i *= 2)

		if (ri == i) is_power_of_two = true;

	if (!is_power_of_two && ri%16 != 0)

		return;

	checker.check2(params);

	checker.graph(strprintf("iir_%04d.tga", ri));

}

void test() {

	Checker checker(16);

	Params p(0.9304775714874, -0.02728819847107, 0.9257968425751);

	cout << fixed << setprecision(8) << checker.check2(p) << endl;

	checker.graph("test.tga");

}

void walk(

	double current_radius,

	double target_radius,

	double step,

	double large_params_step,

	double params_step,

	Params params

) {

	cout << "walk from " << current_radius << " to "<< target_radius << endl;

	step = fabs(step);

	if (target_radius < current_radius) step = -step;

	Params p = params;

	{ Checker checker(current_radius); find_near3(checker, large_params_step, p); }

	for(double r = current_radius; (target_radius - r)/step > -0.5; r += step) {

		Checker checker(r);

		cout << "   " << r << endl;

		double value = find_near2(checker, params_step, p);

		if (value > threshold) {

			value = find_near3(checker, large_params_step, p);

			value = find_near2(checker, params_step, p);

		}

		if (value > threshold) {

			value = initial_find(checker, p);

			value = find_near3(checker, large_params_step, p);

			value = find_near2(checker, params_step, p);

		}

		value = checker.check2(p);

		graph(checker, p);

		log_params(checker, p, value);

		cout << "  " << r

			 << " " << p.k0

			 << " " << p.k1

			 << " " << p.k2

			 << " " << value

			 << endl;

	}

}

void find_initial_params() {

	cout << "find initial params" << endl;

	Finder finder(

		initial_radius,

		Params(0, 0, -1),

		Params(1, 1,  1),

		16,

		16,

		1,

		4,

		8

	);

	finder.find();

	cout << finder.checker.check2(finder.root.current) << endl;

	cout << find_near3(finder.checker, large_params_step, finder.root.current) << endl;

	graph(finder.checker, finder.root.current);

	log_params(finder);

}

void save_coefs(const string &filename, const map< double, pair<params, double=""> > &coefs) {</params,>

	ofstream f(filename.c_str(), ios_base::out | ios_base::trunc);

	f << endl << "double coefs[][5] = {" << endl;

	for(map< double, pair<params, double=""> >::const_iterator i = coefs.begin(); i != coefs.end(); ++i)</params,>

		f << "    { " << i->first

		  << setprecision(20)

		  << ", "     << i->second.first.k0

		  << ", "     << fabs(i->second.first.k1)

		  << ", "     << i->second.first.k2

		  << ", "     << i->second.second

		  << " },"    << endl;

	f << "};" << endl << endl;

}

void save_coefs2(const string &filename, const map< double, pair<params, double=""> > &coefs) {</params,>

	ofstream f(filename.c_str(), ios_base::out | ios_base::trunc);

	f << endl;

	f << "double min_radius = " << min_radius << endl;

	f << "double max_radius = " << max_radius << endl;

	f << "double step = " << step << endl;

	f << "double coefs[][3] = {" << endl;

	for(double r = min_radius; r < max_radius + 0.5*step; r += step) {

		Params p = coefs.count(fix_radius(r)) ? coefs.find(fix_radius(r))->second.first : Params(0.5, 0.0, 0.5);

		f << "    " << setprecision(20)

		  << "{ "   << p.k0

		  << ", "   << fabs(p.k1)

		  << ", "   << p.k2

		  << " },"  << endl;

	}

	f << "};" << endl << endl;

}

void check_coefs() {

	if (coefs.size() < 3) return;

	cout << "check coefs:" << endl;

	map< double, pair<params, double=""> >::const_iterator i = coefs.begin();</params,>

	Params p2 = i->second.first; ++i; p2.k1 = fabs(p2.k1);

	Params p1 = i->second.first; ++i; p1.k1 = fabs(p1.k1);

	while(i != coefs.end()) {

		double r = i->first;

		Params p0 = i->second.first; ++i; p0.k1 = fabs(p0.k1);

		for(int i = 0; i < 3; ++i)

			if (fabs(p0.k[i] - p1.k[i]) > 2.0*fabs(p1.k[i] - p2.k[i]) && fabs(p0.k[i] - p1.k[i]) > 0.0000005)

				cout << "  big step for k" << i << " for radius " << r << ": " << fabs(p0.k[i] - p1.k[i]) << endl;

		p2 = p1;

		p1 = p0;

	}

	cout << "check coefs complete" << endl;

}

void load_coefs(const string &filename, map< double, pair<params, double=""> > &coefs) {</params,>

	cout << "loading..."<< endl;

	ifstream f(filename.c_str());

	while(f) {

		string line;

		getline(f, line);

		stringstream s(line);

		string w[6];

		double radius = 0, value = 0;

		Params p;

		s >> w[0] >> radius

		  >> w[1] >> p.k0

		  >> w[2] >> p.k1

		  >> w[3] >> p.k2

		  >> w[4] >> value

		  >> w[5];

		if ( w[0] == "{"

		  && w[1] == ","

		  && w[2] == ","

		  && w[3] == ","

		  && w[4] == ","

	  	  && w[5] == "},"

	  	  && radius + 0.5*step > min_radius

		  && radius - 0.5*step < max_radius

		  && fabs(radius - fix_radius(radius)) < 1e-5 )

		{

			radius = fix_radius(radius);

			Checker checker(radius);

			double v = checker.check2(p);

			if (fabs(value - v) > 1e-5)

				cout << "  loading: wrong value for radius " << radius << " (" << value << " vs " << v << ")"<< endl;

			if (!coefs.count(radius) || v < coefs[radius].second) {

				coefs[radius].first = p;

				coefs[radius].second = v;

			}

		}

	}

	cout << "done"<< endl;

}

bool is_valid_coefs(double radius) {

	return coefs.count(fix_radius(radius)) && coefs[fix_radius(radius)].second < threshold;

}

void build_diapasones_down(map<double, double=""> &diapasones) {</double,>

	bool valid = true;

	double begin = initial_radius;

	double end = begin;

	for(double r = initial_radius; r > min_radius - 0.5*step; r -= step) {

		bool v = is_valid_coefs(r);

		if (!valid && v)

			diapasones[begin] = end;

		if (v) begin = r; else end = r;

		valid = v;

	}

	if (!valid) diapasones[begin] = end;

}

void build_diapasones_up(map<double, double=""> &diapasones) {</double,>

	bool valid = true;

	double begin = initial_radius;

	double end = begin;

	for(double r = initial_radius; r < max_radius + 0.5*step; r += step) {

		bool v = is_valid_coefs(r);

		if (!valid && v)

			diapasones[begin] = end;

		if (v) begin = r; else end = r;

		valid = v;

	}

	if (!valid) diapasones[begin] = end;

}

void process() {

	log_begin();

	load_coefs(logfile, coefs);

	save_coefs(final_coefs_file, coefs);

	save_coefs2(final_coefs_file2, coefs);

	check_coefs();

	if (!is_valid_coefs(initial_radius))

		find_initial_params();

	cout << "diapasones to find:" << endl;

	map<double, double=""> diapasones;</double,>

	build_diapasones_down(diapasones);

	build_diapasones_up(diapasones);

	for(map<double, double="">::const_iterator i = diapasones.begin(); i != diapasones.end(); ++i)</double,>

		cout << "  " << i->first << " " << i->second << endl;

	cout << "process diapasones:" << endl;

	for(map<double, double="">::const_iterator i = diapasones.begin(); i != diapasones.end(); ++i)</double,>

		walk(i->first, i->second, step, large_params_step, params_step, coefs[fix_radius(i->first)].first);

	log_end();

	save_coefs(final_coefs_file, coefs);

	save_coefs2(final_coefs_file2, coefs);

	check_coefs();

}

int main() {

	test();

	process();

}