#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <map>
#include <cstring>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <iostream>
#include <memory>
#include <random>
#include <stdexcept>
#include <string>
#include <tuple>
#include <unistd.h>
#include <utility>
#include <vector>
#include <ctime>
#include <string>
#include <vector>
#include <sstream>
#include <iomanip>
#include <ostream>
#include <cstring>


class stop_watch
{
private:
public:
	struct timespec ts1_r;
	struct timespec ts2_r;

public:
	stop_watch()
	{
		memset(&ts1_r,0,sizeof(ts1_r));
		memset(&ts2_r,0,sizeof(ts2_r));
	}

	void start()
	{
		clock_gettime(CLOCK_REALTIME,&ts1_r);
	}

	void stop()
	{
		clock_gettime(CLOCK_REALTIME,&ts2_r);
	}

	long double get_micro_r() const
	{
		long double ret = (ts2_r.tv_sec - ts1_r.tv_sec) * 1000000.0L + (ts2_r.tv_nsec - ts1_r.tv_nsec) / 1000.0L;
		return ret;
	}

	long double get_seconds_r() const { return get_micro_r() / 1000000.0L; }
};



template<class C>
const typename C::value_type& random_element( C& container )
{
	static std::random_device rnd;
	auto it = container.begin();
	size_t num = rnd() % container.size();
	std::advance( it, num );

	return *it;
}

class bmp_stream
{
private:
	std::vector<char> buffer;
	std::ofstream out;
	size_t xres;
	size_t yres;
	size_t rowpos;
	size_t pixels;

	template<class T>
	void write( const T& t )
	{
		size_t spos = buffer.size();
		buffer.resize(buffer.size()+sizeof(t));
		char* dst = &buffer[spos];
		memcpy(dst,&t,sizeof(t));
	}

	template<class T>
	void write( const T& t, size_t num )
	{
		size_t spos = buffer.size();
		buffer.resize(buffer.size()+num);
		char* dst = &buffer[spos];
		memcpy(dst,&t,num);
	}

	void create_header( size_t _x, size_t _y, int bpp );

public:
	bmp_stream( size_t x, size_t y, const std::string& fname );

	~bmp_stream( );

	void flush( )
	{
		out.write(&buffer[0],buffer.size());
		buffer.resize(0);
		out.flush();
	}

	void add_pixel( int r, int g, int b );
};

template<size_t N>
struct rgb_pixel
{
	uint32_t r : N;
	uint32_t g : N;
	uint32_t b : N;

	rgb_pixel( uint32_t r_, uint32_t g_, uint32_t b_ ) :
		r(r_),
		g(g_),
		b(b_)
	{
	}
};



template<size_t N, class T>
rgb_pixel<N>& operator/=( rgb_pixel<N>& l, const T& t )
{
	l.r /= t;
	l.g /= t;
	l.b /= t;
	return l;
}

template<size_t N, class T>
rgb_pixel<N>& operator*=( rgb_pixel<N>& l, const T& t )
{
	l.r *= t;
	l.g *= t;
	l.b *= t;
	return l;
}



template<size_t N>
rgb_pixel<N>& operator-=( rgb_pixel<N>& l, const rgb_pixel<N>& r )
{
	if( l.r > r.r )
	{
		l.r -= r.r;
	}
	else
	{
		l.r = 0;
	}

	if( l.g > r.g )
	{
		l.g -= r.g;
	}
	else
	{
		l.g = 0;
	}

	if( l.b > r.b )
	{
		l.b -= r.b;
	}
	else
	{
		l.b = 0;
	}

	return l;
}

template<size_t N>
rgb_pixel<N>& operator+=( rgb_pixel<N>& l, const rgb_pixel<N>& r )
{
	if( l.r + r.r >= (1<<N) )
	{
		l.r = (1<<N)-1;
	}
	else
	{
		l.r += r.r;
	}

	if( l.g + r.g >= (1<<N) )
	{
		l.g = (1<<N)-1;
	}
	else
	{
		l.g += r.g;
	}

	if( l.b + r.b >= (1<<N) )
	{
		l.b = (1<<N)-1;
	}
	else
	{
		l.b += r.b;
	}

	return l;
}

class bmp
{
public:
	typedef rgb_pixel<8> pixel_type;
private:
	std::vector<pixel_type> pixels;
	size_t xdim;
	size_t ydim;
protected:
public:
	bmp ( size_t xs, size_t ys, const pixel_type& proto = pixel_type(0,0,0) );

	/**
	 * Writes the bmp image to a file
	 * @param fname is the name of the file to write
	 */
	void write( const std::string& fname ) const;
	
	void fill( const pixel_type& col )
	{
		for (size_t x = 0; x < xdim; ++x)
		{
			for (size_t y = 0; y < ydim; ++y)
			{
				this->raw(x,y) = col;
			}
		}
	}

	pixel_type& raw( size_t x, size_t y )
	{
		assert(x < xdim);
		assert(y < ydim);
		return pixels[x + xdim*y];
	}

	const pixel_type& raw( size_t x, size_t y ) const
	{
		return pixels[x + xdim*y];
	}

	void raw_line_add( int32_t x0, int32_t y0, int32_t x1, int32_t y1, const pixel_type& col )
	{
		const int32_t dx = std::abs(x1-x0);
		const int32_t dy = std::abs(y1-y0);

		const int32_t sx = (x0<x1) ? 1 : -1;
		const int32_t sy = (y0<y1) ? 1 : -1;

		int32_t err = dx-dy;

		while(true)
		{
			raw(x0,y0) += col;
			if( x0 == x1 && y0 == y1 ) break;
			int32_t e2 = 2*err;
			
			if( e2 > -dy )
			{
				err -= dy;
				x0 += sx;
			}
			if( e2 < dx )
			{
				err += dx;
				y0 += sy;
			}
		}
	}

	void raw_line( int32_t x0, int32_t y0, int32_t x1, int32_t y1, const pixel_type& col )
	{
		const int32_t dx = std::abs(x1-x0);
		const int32_t dy = std::abs(y1-y0);

		const int32_t sx = (x0<x1) ? 1 : -1;
		const int32_t sy = (y0<y1) ? 1 : -1;

		int32_t err = dx-dy;

		while(true)
		{
			raw(x0,y0) = col;
			if( x0 == x1 && y0 == y1 ) break;
			int32_t e2 = 2*err;

			if( e2 > -dy )
			{
				err -= dy;
				x0 += sx;
			}
			if( e2 < dx )
			{
				err += dx;
				y0 += sy;
			}
		}
	}

	void max_contrast_all( )
	{
		uint32_t min = std::numeric_limits<uint32_t>::max();
		uint32_t max = 0;

		for( auto& p : pixels )
		{
			min = std::min(min,uint32_t(p.r));
			min = std::min(min,uint32_t(p.g));
			min = std::min(min,uint32_t(p.b));

			max = std::max(max,uint32_t(p.r));
			max = std::max(max,uint32_t(p.g));
			max = std::max(max,uint32_t(p.b));
		}

		uint32_t range = max - min;
		double factor = 255.0 / range;
		for( auto& p : pixels )
		{
			p.r -= min;
			p.g -= min;
			p.b -= min;

			p *= factor;
		}
	}

	std::vector<std::pair<int32_t,int32_t>> raw_line_pixels( int32_t x0, int32_t y0, int32_t x1, int32_t y1 ) const
	{
		std::vector<std::pair<int32_t,int32_t>> ret;
		const int32_t dx = std::abs(x1-x0);
		const int32_t dy = std::abs(y1-y0);

		const int32_t sx = (x0<x1) ? 1 : -1;
		const int32_t sy = (y0<y1) ? 1 : -1;

		int32_t err = dx-dy;

		while(true)
		{
			ret.push_back({x0,y0});
			if( x0 == x1 && y0 == y1 ) break;
			int32_t e2 = 2*err;
			
			if( e2 > -dy )
			{
				err -= dy;
				x0 += sx;
			}
			if( e2 < dx )
			{
				err += dx;
				y0 += sy;
			}
		}
		return ret;
	}

	std::vector<pixel_type>& get_pixels( ) { return pixels; }
	const std::vector<pixel_type>& get_pixels( ) const { return pixels; }
};

std::pair<size_t, size_t> load_bmp( const std::string& fname, bmp& bmap )
{
	std::cout << "Loading " << fname << "..." << std::flush;
	std::ifstream in(fname);
	std::string p6;
	uint32_t xs;
	uint32_t ys;
	std::string cpp;

	in >> p6;
	in >> ys;
	in >> xs;
	in >> cpp;

	bmap = bmp(ys,xs);

	char dummy;
	in.read(&dummy,1);

	for (size_t i = 0; i < xs; ++i)
	{
		for (size_t j = 0; j < ys; ++j)
		{
			unsigned char r,g,b;
			in.read(reinterpret_cast<char*>(&r),1);
			in.read(reinterpret_cast<char*>(&g),1);
			in.read(reinterpret_cast<char*>(&b),1);
			bmap.raw(j,i) = {r,g,b};
		}
	}

	std::cout << "done" << std::endl;
	return { ys, xs };
}

size_t distance( ssize_t x0, ssize_t x1, ssize_t y0, ssize_t y1 )
{
	ssize_t dx = x0-x1;
	ssize_t dy = y0-y1;

	return dx*dx+dy*dy;
}

ssize_t rgb_distance( ssize_t r0, ssize_t r1, ssize_t g0, ssize_t g1, ssize_t b0, ssize_t b1 )
{
	ssize_t dr = r0-r1;
	ssize_t dg = g0-g1;
	ssize_t db = b0-b1;

	return dr*dr+dg*dg+db*db;
}


struct environment
{
	ssize_t x = 0;
	ssize_t y = 0;

	size_t lines = 0;

	// minimal line length (pixels squared)
	ssize_t m = 0;
	// maximum line length (pixels squared)
	ssize_t M = 1'000'000;
};

struct line
{
	int16_t x0 = 0;
	int16_t y0 = 0;
	int16_t x1 = 0;
	int16_t y1 = 0;

	uint8_t r;
	uint8_t g;
	uint8_t b;
	uint8_t skip_next = 0;
	bool active = true;
	bool auto_deactivate = false;

	static uint64_t count;

	line( int x0_, int y0_, int x1_, int y1_, int r_, int g_, int b_ ) :
		x0(x0_),
		y0(y0_),
		x1(x1_),
		y1(y1_),
		r(r_),
		g(g_),
		b(b_)
	{
		count++;
	}

	line( ) 
	{
		++count;
	}

	line( const line& other ) :
		x0(other.x0),
		y0(other.y0),
		x1(other.x1),
		y1(other.y1),
		r(other.r),
		g(other.g),
		b(other.b),
		skip_next(other.skip_next),
		active(other.active),
		auto_deactivate(other.auto_deactivate)
	{
		++count;
	}

	line( line&& other ) :
		x0(other.x0),
		y0(other.y0),
		x1(other.x1),
		y1(other.y1),
		r(other.r),
		g(other.g),
		b(other.b),
		skip_next(other.skip_next),
		active(other.active),
		auto_deactivate(other.auto_deactivate)
	{
		++count;
	}


	line& operator=( const line& other )
	{
		x0 = other.x0;
		y0 = other.y0;
		x1 = other.x1;
		y1 = other.y1;
		r = other.r;
		g = other.g;
		b = other.b;
		skip_next = other.skip_next;
		active = other.active;
		auto_deactivate = other.auto_deactivate;
		return *this;
	}

	~line( )
	{
		--count;
	}

	void randomize( const environment& env, const bmp& orig  )
	{
		static std::random_device rnd;
		static std::uniform_int_distribution<ssize_t> dist(0,255);

		ssize_t d = 0;

//		r = g = b = 128;
//		r = dist(rnd);
//		g = dist(rnd);
//		b = dist(rnd);

		do
		{

			x0 = rnd() % env.x;
			x1 = rnd() % env.x;

			y0 = rnd() % env.y;
			y1 = rnd() % env.y;

			d = distance( x0,x1,y0,y1);
		}
		while( d < env.m || d > env.M );

		auto pixels = orig.raw_line_pixels( x0,y0,x1,y1 );
		size_t pr = 0;
		size_t pg = 0;
		size_t pb = 0;

		for (size_t i = 0; i < pixels.size(); ++i)
		{
			pr += orig.raw(pixels[0].first,pixels[0].second).r;
			pg += orig.raw(pixels[0].first,pixels[0].second).g;
			pb += orig.raw(pixels[0].first,pixels[0].second).b;
		}

		r = pr / pixels.size();
		g = pg / pixels.size();
		b = pb / pixels.size();
	}

	void mutate( const environment& env )
	{
		static std::random_device rnd;
		static std::uniform_real_distribution<double> prob_dist(0,1);
		static std::uniform_int_distribution<ssize_t> dist(-10,10);

		x0 += dist(rnd);
		x1 += dist(rnd);
		y0 += dist(rnd);
		y1 += dist(rnd);

		x0 = std::max(int16_t(0),std::min(x0,int16_t(env.x-1)));
		x1 = std::max(int16_t(0),std::min(x1,int16_t(env.x-1)));
		y0 = std::max(int16_t(0),std::min(y0,int16_t(env.y-1)));
		y1 = std::max(int16_t(0),std::min(y1,int16_t(env.y-1)));
		
		if( prob_dist(rnd) < 0.1 )
		{
			active = !active;
		}
		if( prob_dist(rnd) < 0.1 )
		{
			skip_next += dist(rnd);
		}

		// let them do some wraparound
		r += dist(rnd);
		g += dist(rnd);
		b += dist(rnd);

		ssize_t d = distance( x0, x1, y0, y1 );

		auto_deactivate = ( d < env.m || d > env.M );
	}

	bool is_active( ) const { return !auto_deactivate && active; }
};

uint64_t line::count = 0;

struct chromo
{
	std::vector<line> lines;
	// per line (and one for the probability itself)
	double mutation_probability = 0.00005;

	// = 0 means "never calculated"
	size_t score = 0;

	std::unique_ptr<chromo> mate( const chromo& other, const environment& env, bool clone )
	{
		static std::random_device rnd;
		static std::uniform_int_distribution<size_t> block_dist(1,5);
		static std::uniform_real_distribution<double> mut_dist(0,1);
		static std::uniform_real_distribution<double> mut_mut_dist(0.5,2);

		std::unique_ptr<chromo> ret(new chromo);
		ret->mutation_probability = (other.mutation_probability + mutation_probability)/2.0;
		ret->lines.reserve(other.lines.size());

		const std::vector<line>* src = &other.lines;
		for( size_t i = 0; i < src->size(); )
		{
			if( mut_dist(rnd) > 0.5 && !clone )
			{
				src = &other.lines;
			}
			else
			{
				src = &lines;
			}
			ssize_t block = block_dist(rnd);
			block = std::min( block, ssize_t(src->size()) - ssize_t(i) );

			size_t copies = 1;
			double dr = mut_dist(rnd);
			if( dr < mutation_probability )
			{
				copies = 0;
			}
			else if( dr > (1-mutation_probability) )
			{
				copies = block_dist(rnd);
			}

			for (size_t k = 0; k < copies; ++k)
			{
				for (ssize_t j = 0; j < block; ++j)
				{
					ret->lines.push_back( (*src)[i+j] );
					if( mut_dist(rnd) < mutation_probability )
					{
						ret->lines.back().mutate(env);
					}
				}
			}
			i += block;
		}
		ret->lines.resize(env.lines * ( 1 + (.2 * mut_dist(rnd))));
		if( mut_dist(rnd) < mutation_probability )
		{
			ret->mutation_probability *= mut_mut_dist(rnd);
		}
		return ret;
	}
};


void paint_remaining( const std::vector<line>& lines, bmp& bmap, size_t start )
{
	for (size_t i = start; i < lines.size(); ++i)
	{
		if( lines[i].is_active() )
		{
			if( lines[i].skip_next > 0 )
			{
				i+= lines[i].skip_next;
				continue;
			}
			const auto& l = lines[i];
			bmap.raw_line(l.x0,l.y0,l.x1,l.y1,{l.r,l.g,l.b});
		}
	}
}

void paint( const std::vector<line>& lines, bmp& bmap )
{
	bmap.fill({255,255,255});

	for (size_t i = 0; i < lines.size(); ++i)
	{
		if( lines[i].is_active() )
		{
			if( lines[i].skip_next > 0 )
			{
				i+= lines[i].skip_next;
				continue;
			}
			const auto& l = lines[i];
			bmap.raw_line(l.x0,l.y0,l.x1,l.y1,{l.r,l.g,l.b});
		}
	}
}


void paint( const chromo& chr, bmp& bmap )
{
	paint(chr.lines,bmap);
}

void paint_layers( const std::vector<line>& lines, bmp& bmap )
{
	bmap.fill({0,0,0});

	for (size_t i = 0; i < lines.size(); ++i)
	{
		if( lines[i].is_active() )
		{
			if( lines[i].skip_next > 0 )
			{
				i+= lines[i].skip_next;
				continue;
			}
			const auto& l = lines[i];
			auto pixels = bmap.raw_line_pixels(l.x0,l.y0,l.x1,l.y1);
			for( auto& p : pixels )
			{
				bmap.raw(p.first,p.second) += bmp::pixel_type(2,2,2);
				
			}
		}
	}
}

void paint_layers( const chromo& chr, bmp& bmap )
{
	paint_layers( chr.lines, bmap );
}

template<class PX>
ssize_t generate_score( const std::vector<PX>& a, const std::vector<PX>& b )
{
	ssize_t score = 0;
	for (size_t i = 0; i < a.size(); ++i)
	{
		score += rgb_distance( a[i].r, b[i].r, a[i].g, b[i].g, a[i].b, b[i].b );
	}
	return score;
}

ssize_t evaluate( const std::vector<line>& lines, const bmp& orig )
{
	ssize_t score = 0;

	bmp bmap(orig);

	const auto& op = orig.get_pixels();
	const auto& np = bmap.get_pixels();

	paint(lines,bmap);

	for (size_t i = 0; i < op.size(); ++i)
	{
		score += rgb_distance( op[i].r, np[i].r, op[i].g, np[i].g, op[i].b, np[i].b );
	}
	return score;
}

ssize_t evaluate( const chromo& chr, const bmp& orig )
{
	if( chr.score != 0 ) return chr.score;
	ssize_t score = evaluate(chr.lines,orig);
	return score;
}

void seed( const bmp& /*orig*/, const std::vector<line>& proto_lines, std::vector<std::unique_ptr<chromo>>& population, size_t individuals, size_t /*linelength*/, const environment& env )
{
	std::cout << "Start seeding out of " << proto_lines.size() << " proto lines ..." << std::flush;
	population.reserve(individuals);

	std::unique_ptr<chromo> ni( new chromo);
	ni->lines.reserve(proto_lines.size());
	for (size_t i = 0; i < proto_lines.size(); ++i)
	{
		ni->lines.push_back(proto_lines[i]);
	}

	population.emplace_back( new chromo(*ni) );
	for (size_t i = 1; i < individuals; ++i)
	{
		// some bad clones
		population.emplace_back(  ni->mate(*population[i-1],env,true) );
	}
	std::cout << "done" << std::endl;
}


void evolve( const std::string& fname, const std::vector<line>& proto_lines )
{
	bmp orig(1,1);
	size_t x;
	size_t y;
	std::tie(x,y) = load_bmp( fname, orig );

	environment env;
	env.x = x;
	env.y = y;

//	env.lines = 0.25*x*y;
	env.lines = 0.3*x*y;

	env.m = 10*10;
	env.M = 50*50;

	std::vector<std::unique_ptr<chromo>> population;
	const size_t individuals = 400;

	seed(orig, proto_lines, population,10,env.lines,env);

	const size_t generations = 100000;
	const size_t reproducers = 10;
	const size_t kidnum = 10;
	const size_t clones = 10;

	std::cout << "Target line amount is " << env.lines << "\n";
	for (size_t j = 0; j < generations; ++j)
	{
		std::cout << "Evaluating up to " << population.size() << " individuals..." << std::flush;
		double score = 0.0;
		for (size_t i = 0; i < population.size(); ++i)
		{
			population[i]->score = evaluate(*population[i],orig);
			score += population[i]->score;
		}
		score /= population.size();
		std::cout << "done\n";
		std::cout << "Average score "<< score << "\n",
		std::cout << "Sorting them..." << std::flush;
		std::sort(population.begin(), population.end(), []( auto& l, auto& r ) { return l->score < r->score; } );

		population.resize(std::min(population.size(),individuals));
		std::cout << "done" << std::endl;

		static chromo* last_painted = nullptr;

		if( last_painted != population[0].get() )
		{
			last_painted = population[0].get();
			std::ostringstream on;
			on << "generation_" << std::setw(5) << std::setfill('0') << j << ".bmp";

			bmp b(x,y);
			paint(*population[0],b);
			b.write(on.str());
			paint_layers(*population[0],b);
			b.max_contrast_all();
			b.write("layers_" + on.str());
		}
		std::cout << "At generation " << j << " we have " << line::count << " lines (" << double(line::count)/population.size() << " per chromo)\n";

		std::vector<std::unique_ptr<chromo>> kids;
		std::cout << "Running " << reproducers << " reproduction cycles..." << std::flush;

		for (size_t i = 0; i < reproducers; ++i )
		{
			for (size_t k = 0; k < kidnum; ++k)
			{
				kids.push_back( population[i]->mate(*random_element(population),env,false ));
			}
			for (size_t k = 0; k < clones; ++k)
			{
				kids.push_back( population[i]->mate(*population[0],env,true ));
			}
		}
		std::cout << "done\n";
		std::cout << "Merging " << kids.size() << " offsprings..." << std::flush;
		population.reserve( population.size() + kids.size() );
		for (size_t i = 0; i < kids.size(); ++i)
		{
			population.push_back( std::move(kids[i]));
		}
		std::cout << "done" << std::endl;
	}
}


std::pair<bmp,std::vector<line>>  paint_useful( const bmp& orig, bmp& clone, std::vector<line>& retlines, bmp& layer, const std::string& outprefix, size_t x, size_t y )
{
	const size_t pixels = (x*y);
	const size_t lines = 0.3*pixels;
//	const size_t lines = 10000;

//	const size_t start_accurate_color = lines/4;

	std::random_device rnd;

	std::uniform_int_distribution<size_t> distx(0,x-1);
	std::uniform_int_distribution<size_t> disty(0,y-1);
	std::uniform_int_distribution<size_t> col(-15,15);
	std::uniform_int_distribution<size_t> acol(0,255);

	const ssize_t m = 1*1;
	const ssize_t M = 50*50;

	retlines.reserve( lines );

	for (size_t i = retlines.size(); i < lines; ++i)
	{
		size_t x0;
		size_t x1;

		size_t y0;
		size_t y1;

		size_t dist = 0;
		do
		{
			x0 = distx(rnd);
			x1 = distx(rnd);

			y0 = disty(rnd);
			y1 = disty(rnd);

			dist = distance(x0,x1,y0,y1);
		}
		while( dist > M || dist < m );

		std::vector<std::pair<int32_t,int32_t>> points = clone.raw_line_pixels(x0,y0,x1,y1);

		ssize_t r = 0;
		ssize_t g = 0;
		ssize_t b = 0;

		for (size_t i = 0; i < points.size(); ++i)
		{
			r += orig.raw(points[i].first,points[i].second).r;
			g += orig.raw(points[i].first,points[i].second).g;
			b += orig.raw(points[i].first,points[i].second).b;
		}

		r += col(rnd);
		g += col(rnd);
		b += col(rnd);

		r /= points.size();
		g /= points.size();
		b /= points.size();
		
		r %= 255;
		g %= 255;
		b %= 255;

		r = std::max(ssize_t(0),r);
		g = std::max(ssize_t(0),g);
		b = std::max(ssize_t(0),b);

//		r = acol(rnd);
//		g = acol(rnd);
//		b = acol(rnd);

//		if( i > start_accurate_color )
		{
			ssize_t dp = 0; // accumulated distance of new color to original
			ssize_t dn = 0; // accumulated distance of current reproduced to original
			for (size_t i = 0; i < points.size(); ++i)
			{
				dp += rgb_distance(
										orig.raw(points[i].first,points[i].second).r,r,
										orig.raw(points[i].first,points[i].second).g,g,
										orig.raw(points[i].first,points[i].second).b,b
									);

				dn += rgb_distance(
										clone.raw(points[i].first,points[i].second).r,orig.raw(points[i].first,points[i].second).r,
										clone.raw(points[i].first,points[i].second).g,orig.raw(points[i].first,points[i].second).g,
										clone.raw(points[i].first,points[i].second).b,orig.raw(points[i].first,points[i].second).b
									);

			}

			if( dp > dn ) // the distance to original is bigger, use the new one
			{
				--i;
				continue;
			}
			// also abandon if already too bad
//			if( dp > 100000 )
//			{
//				--i;
//				continue;
//			}
		}
		
		layer.raw_line_add(x0,y0,x1,y1,{1u,1u,1u});
		clone.raw_line(x0,y0,x1,y1,{(uint32_t)r,(uint32_t)g,(uint32_t)b});
		retlines.push_back({ (int)x0,(int)y0,(int)x1,(int)y1,(int)r,(int)g,(int)b});

		static time_t last = 0;
		time_t now = time(0);
		if( i % (lines/100) == 0 )
		{
			std::ostringstream fn;
			fn << outprefix + "perc_" << std::setw(3) << std::setfill('0') << (i/(lines/100)) << ".bmp"; 
			clone.write(fn.str());
			bmp lc(layer);
			lc.max_contrast_all();
			lc.write(outprefix + "layer_" + fn.str());
		}

		if( (now-last) > 10 )
		{
			last = now;
			static int st = 0;
			std::ostringstream fn;
			fn << outprefix + "inter_" << std::setw(8) << std::setfill('0') << i << ".bmp";
			clone.write(fn.str());

			++st;
		}
	}
	clone.write(outprefix + "clone.bmp");
	return { clone, retlines };
}


void erase_bad( std::vector<line>& lines, const bmp& orig )
{
	ssize_t current_score = evaluate(lines,orig);

	std::vector<line> newlines(lines);

	uint32_t deactivated = 0;
	std::cout << "current_score = " << current_score << "\n";
	for (size_t i = 0; i < newlines.size(); ++i)
	{
		newlines[i].active = false;
		ssize_t score = evaluate(newlines,orig);
		if( score > current_score )
		{
			newlines[i].active = true;
		}
		else
		{
			current_score = score;
			++deactivated;
		}
		if( i % 1000 == 0 )
		{
			std::ostringstream fn;
			fn << "erase_" << std::setw(6) << std::setfill('0') << i << ".bmp";
			bmp tmp(orig);
			paint(newlines,tmp);
			tmp.write(fn.str());
			paint_layers(newlines,tmp);
			tmp.max_contrast_all();
			tmp.write("layers_" + fn.str());
			std::cout << "\r i = " << i << std::flush;
		}
	}
	std::cout << "\n";
	std::cout << "current_score = " << current_score << "\n";
	std::cout << "deactivated = " << deactivated << "\n";
	
	
	bmp tmp(orig);
	
	paint(newlines,tmp);
	tmp.write("newlines.bmp");
	lines.clear();
	for (size_t i = 0; i < newlines.size(); ++i)
	{
		if( newlines[i].is_active() )
		{
			lines.push_back(newlines[i]);
		}
	}
}

std::pair<bmp,std::vector<line>>  try_ggate( const std::string& fname )
{
	bmp orig(100,100);
	size_t x,y;

	std::tie(x,y) = load_bmp(fname,orig);
	orig.write("orig.bmp");

	size_t ir = 0;
	size_t ig = 0;
	size_t ib = 0;
	for( auto& p : orig.get_pixels() )
	{
		ir += p.r;
		ig += p.g;
		ib += p.b;
	}
	ir /= orig.get_pixels().size();
	ig /= orig.get_pixels().size();
	ib /= orig.get_pixels().size();

	bmp clone(orig);
//	clone.fill({255,255,255});
	clone.fill({(uint32_t)ir,(uint32_t)ig,(uint32_t)ib});

	bmp layer(x,y);
	layer.fill({0,0,0});
	std::vector<line> retlines;
	stop_watch sw;

	sw.start();
	paint_useful( orig, clone, retlines, layer, "standard_", x, y );
	sw.stop();
	std::cout << "Initial distribution took " << sw.get_seconds_r() << "s to write " << retlines.size() << " lines \n";
	sw.start();
	erase_bad(retlines,orig);
	sw.stop();
	std::cout << "Removing unnecessary took " << sw.get_seconds_r() << " and left us with " << retlines.size() <<" lines \n";

	sw.start();
	paint_useful( orig, clone, retlines, layer, "second_", x, y );
	sw.stop();
	std::cout << "Second layer of useful lines took " << sw.get_seconds_r() << " and we now have " << retlines.size() <<" lines \n";

	return { clone, retlines };
}

int main(int argc, const char *argv[])
{
	std::string f = "ggate.ppm";

	if( argc > 1 )
	{
		f = argv[1];
	}
	
	stop_watch sw;

	std::vector<std::tuple<double,double,double>> allpixels;

	const size_t num = 1;
	bmp b(1,1);
	std::vector<line> lines;
	for (size_t i = 0; i < num; ++i)
	{
		std::tie( b, lines ) = try_ggate(f);

		if( allpixels.empty() )
		{
			allpixels.resize(b.get_pixels().size(),std::make_tuple(0.0,0.0,0.0));
		}
		for (size_t j = 0; j < b.get_pixels().size(); ++j)
		{
			std::get<0>(allpixels[j]) += b.get_pixels()[j].r;
			std::get<1>(allpixels[j]) += b.get_pixels()[j].g;
			std::get<2>(allpixels[j]) += b.get_pixels()[j].b;
		}
		evolve(f,lines);
	}
	for (size_t i = 0; i < allpixels.size(); ++i)
	{
		b.get_pixels()[i].r = std::get<0>(allpixels[i]) / num;
		b.get_pixels()[i].g = std::get<1>(allpixels[i]) / num;
		b.get_pixels()[i].b = std::get<2>(allpixels[i]) / num;
	}
	b.write("merged.bmp");
}


bmp_stream::bmp_stream( size_t x, size_t y, const std::string& fname ) :
	xres(x),
	yres(y),
	rowpos(0),
	pixels(0)
{
	out.open(fname.c_str());
	if(!out)
	{
		throw std::runtime_error("Could not open " + fname);
	}
	create_header(x,y,24); // 8 bit per pixel per color
}

bmp_stream::~bmp_stream( )
{
	// pad the rest of the file with black pixels
	while( pixels < (xres*yres) )
	{
		add_pixel(0,0,0);
	}
	flush();
//	out.write(&buffer[0],buffer.size());
}

void bmp_stream::add_pixel( int r, int g, int b )
{
	++pixels;
	++rowpos;

	unsigned char col[3];
	col[0] = b;
	col[1] = g;
	col[2] = r;
	buffer.insert(buffer.end(),{ (char)col[0], (char)col[1], (char)col[2] });
//	out.write(reinterpret_cast<const char*>(col),3);

	if( rowpos == xres )
	{
		// the position in the row means we are finished with the row. Close it
		// by rounding its bytes up to 4-byte boundary and start with the next
		// row
		rowpos = 0;
		char dummy = 0xcc;
		switch( xres*3 % 4 )
		{
			case 1:
				buffer.push_back(dummy);
//				write(dummy);
			case 2:
				buffer.push_back(dummy);
//				write(dummy);
			case 3:
				buffer.push_back(dummy);
//				write(dummy);
			case 0:
			default:
				break;
		}
	}
}

void bmp_stream::create_header( size_t _x, size_t _y, int bpp )
{
	// BMP Header
	const char magic[] = { 'B', 'M' };
	write(magic,2);
	// XXX (#31) approximation, must be corrected by header and stuff
	uint32_t fsize = _x*_y*(bpp/8) + 0x100;
	write(fsize);
	// write 2 times 2 implementation defined bytes
	write(magic,2);
	write(magic,2);

	// beginning of the actual bitmap data...
	uint32_t start = 0x100;
	write(start);

	// DIB Header
	uint32_t hsize = 40; // as defined by the format
	write(hsize);

	int32_t x = _x;
	int32_t y = - _y;
	write(x);
	write(y);
	int16_t colplanes = 1; // fixed by the format
	write(colplanes);
	uint16_t bitdepth = bpp;
	write(bitdepth);
	uint32_t compmethod = 0; // no compression
	write(compmethod);
	uint32_t rawbmpsize = (x+1)*y*3;
	write(rawbmpsize);
	int32_t hres = 240;
	int32_t vres = 240;
	write(hres);
	write(vres);
	int32_t colpal = 0;
	write(colpal);
	write(colpal);
	char dummy[0x100] = { 0 };
	size_t padbytes = 0x100 - buffer.size();//out.tellp();
//	out.write(dummy,padbytes);
	write(dummy,padbytes);
	buffer.reserve(1.1*_x*_y+buffer.size());
}

bmp::bmp ( size_t xs, size_t ys, const pixel_type& proto ) :
	xdim(xs),
	ydim(ys)
{
	pixels.resize(xdim*ydim,proto);
}

void bmp::write( const std::string& fname ) const
{
	bmp_stream bs(xdim,ydim,fname);
	for( size_t i = 0; i < pixels.size(); i++)
	{
		const pixel_type& pix = pixels[i];
		bs.add_pixel( pix.r, pix.g, pix.b );
	}
}

// vim: tabstop=4 shiftwidth=4 noexpandtab