/*
 * Copyright (c) 2014 Daniel Hartmeier <daniel@benzedrine.ch>
 * All rights reserved.
 *
 * COPYING, REDISTRIBUTION AND USE IN ANY FORM ARE PROHIBITED WITHOUT AN
 * EXPLICIT WRITTEN PERMISSION.
 */

#include <ode/ode.h>

class World {
public:
	World();
	World(const World &o);
	~World();
	void read();
	void randomize();
	void run();
	void shift();
	void print() const;
	void bounds(dVector3 &v, int t, bool print) const;
private:
	void init();
	void mass(float t);
	float collision() const;
	bool tighten(dReal *b);
	void loosen(float d);
	static void near(void *data, dGeomID o1, dGeomID o2);
	dWorldID world;
	dSpaceID space;
	dJointGroupID jg;
	dGeomID wall[6], item[101];
};

static int id[100];
static const float gravity = 1.5;
static World *best = NULL;
static float min_surface = 999999.9;
static bool refine = false;

void World::init()
{
	world = dWorldCreate();
	space = dSweepAndPruneSpaceCreate(0, dSAP_AXES_XYZ);
	jg = dJointGroupCreate(0);
	dWorldSetGravity(world, -gravity, -gravity, -gravity);
	dWorldSetQuickStepNumIterations(world, 10);
	wall[0] = wall[1] = wall[2] = 0;
	wall[3] = dCreatePlane(space, 0.0, 0.0, 1.0, 0.0);
	wall[4] = dCreatePlane(space, 0.0, 1.0, 0.0, 0.0);
	wall[5] = dCreatePlane(space, 1.0, 0.0, 0.0, 0.0);
	for (int i = 0; i < 101; ++i)
		item[i] = 0;
}

World::World()
{
	init();
}

World::World(const World &o)
{
	init();
	for (int i = 0; o.item[i]; ++i) {
		const dReal *p = dGeomGetPosition(o.item[i]);
		const dReal *R = dGeomGetRotation(o.item[i]);
		dVector3 d;
		dGeomBoxGetLengths(o.item[i], d);
		dBodyID body = dBodyCreate(world);
		item[i] = dCreateBox(space, d[0], d[1], d[2]);
		dGeomSetBody(item[i], body);
		dGeomSetPosition(item[i], p[0], p[1], p[2]);
		dGeomSetRotation(item[i], R);
	}
	for (int i = 0; i < 3; ++i)
		if (o.wall[i]) {
			dVector4 p;
			dGeomPlaneGetParams(o.wall[i], p);
			wall[i] = dCreatePlane(space, p[0], p[1], p[2], p[3]);
		}
}

World::~World()
{
	dJointGroupDestroy(jg);
	for (int i = 0; item[i]; ++i) {
		dBodyDestroy(dGeomGetBody(item[i]));
		dGeomDestroy(item[i]);
	}
	for (int i = 0; i < 6; ++i)
		if (wall[i])
			dGeomDestroy(wall[i]);
	dSpaceDestroy(space);
	dWorldDestroy(world);
}

void World::read()
{
	char s[256];
	float l, w, h;

	for (int n = 0; fgets(s, sizeof(s), stdin); ++n) {
		if (n == 100 ||
		    sscanf(s, "%d %f %f %f", &id[n], &l, &w, &h) != 4) {
			printf("Error\n");
			exit(0);
		}
		dBodyID body = dBodyCreate(world);
		item[n] = dCreateBox(space, l, w, h);
		dGeomSetBody(item[n], body);
	}
}

void World::randomize()
{
	const int d = 800000;
	for (int i = 0; item[i]; ++i) {
		dGeomSetPosition(item[i],
		    (float)arc4random_uniform(d) / 1000.0 + 25.0,
		    (float)arc4random_uniform(d) / 1000.0 + 25.0,
		    (float)arc4random_uniform(d) / 1000.0 + 25.0);
		dMatrix3 R;
		for (int j = 0; j < 12; ++j)
			R[j] = (float)arc4random_uniform(2001) / 1000.0 - 1.0;
		dGeomSetRotation(item[i], R);
	}
}

void World::mass(float t)
{
	for (int i = 0; item[i]; ++i) {
		dVector3 v;
		dGeomBoxGetLengths(item[i], v);
		dMass mass;
		dMassSetBoxTotal(&mass, t, v[0], v[1], v[2]);
		dBodySetMass(dGeomGetBody(item[i]), &mass);
	}
}

void World::run()
{
	dVector3 a;
	int state = 0;
	int idle = 0;
	int limit = 30;

	mass(refine ? 0.001 : 0.5);
	for (int pass = 0; pass < 100000; ++pass) {

		for (int i = 0; i < (refine ? 1 : 3); ++i) {
			dSpaceCollide(space, this, &near);
			dWorldQuickStep(world, refine ? 0.005 : 0.05);
			dJointGroupEmpty(jg);
		}

		dVector3 b;
		bounds(b, 2, false);
		if (b[0] > 20000.0 || b[1] > 20000.0 || b[2] > 20000.0)
			break;
		float surface = 2.0 * (b[0] * b[1] + b[0] * b[2] +
		    b[1] * b[2]);
		float w = collision();
		float d = dCalcPointsDistance3(a, b);
		a[0] = b[0]; a[1] = b[1]; a[2] = b[2];
		if (refine) {
			if (w >= 0.005)
				loosen(0.001);
			else if (d > -0.001 && d < 0.001)
				return;
			continue;
		}
		if (!best || surface < min_surface) {
			if (best)
				delete best;
			best = new World(*this);
			min_surface = surface;
		}
		if (tighten(b))
			limit++;
		if (d > -0.1 && d < 0.1)
			idle++;
		else
			idle = 0;
		if (idle < 3)
			continue;
		if (state >= limit)
			return;
		bounds(b, 1, false);
		if (state == 0) {
			wall[0] = dCreatePlane(space, -1.0, 0.0, 0.0, -b[0]);
			dWorldSetGravity(world, gravity, -gravity, -gravity);
		} else if (state == 1) {
			wall[1] = dCreatePlane(space, 0.0, -1.0, 0.0, -b[1]);
			dWorldSetGravity(world, gravity, gravity, -gravity);
		} else if (state == 2) {
			wall[2] = dCreatePlane(space, 0.0, 0.0, -1.0, -b[2]);
			dWorldSetGravity(world, gravity, gravity, gravity);
		} else {
			dVector3 v;
			dWorldGetGravity(world, v);
			switch (arc4random_uniform(3)) {
			case 0:
				dWorldSetGravity(world, -v[0], v[1], v[2]);
				break;
			case 1:
				dWorldSetGravity(world, v[0], -v[1], v[2]);
				break;
			case 2:
				dWorldSetGravity(world, v[0], v[1], -v[2]);
				break;
			}
		}
		state++;
		idle = 0;
	}
}

void World::shift()
{
	dVector3 b;

	bounds(b, 0, false);
	for (int i = 0; item[i]; ++i) {
		const dReal *p = dGeomGetPosition(item[i]);
		dGeomSetPosition(item[i], p[0] - b[0], p[1] - b[1],
		    p[2] - b[2]);
	}
}

void World::print() const
{
	dVector3 b;

	bounds(b, 1, false);
	printf("%.2f %.2f %.2f\n", b[0], b[1], b[2]);
	bounds(b, 1, true);
}

void World::bounds(dVector3 &v, int t, bool print) const
{
	const int sign[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 },
	    { -1, 1, -1 }, { -1, -1 , 1 }, { 1, -1, 1 }, { 1, 1, 1 },
	    { -1, 1, 1 } };
	dVector3 min, max;
	int n = 0;

	for (int j = 0; j < 3; ++j)
		min[j] = max[j] = 0.0;
	for (int i = 0; item[i]; ++i) {
		if (print)
			printf("%d", id[n++]);
		dVector3 d;
		dGeomBoxGetLengths(item[i], d);
		for (int k = 0; k < 8; ++k) {
			dVector3 v, p;
			for (int j = 0; j < 3; ++j)
				v[j] = sign[k][j] * d[j] / 2.0;
			dGeomGetRelPointPos(item[i], v[0], v[1], v[2], p);
			if (print)
				printf(" (%.2f, %.2f, %.2f)",
				    p[0], p[1], p[2]);
			for (int j = 0; j < 3; ++j)
				if (p[j] < min[j])
					min[j] = p[j];
				else if (p[j] > max[j])
					max[j] = p[j];
		}
		if (print)
			printf("\n");
	}
	for (int j = 0; j < 3; ++j)
		v[j] = (t == 0 ? min[j] : (t == 1 ? max[j] : max[j] - min[j]));
}

void World::near(void *data, dGeomID o1, dGeomID o2)
{
	World *world = (World *)data;
	const int N = 8;
	dContact contact[N];

	int n = dCollide(o1, o2, N, &contact[0].geom, sizeof(dContact));
	for (int i = 0; i < n; ++i) {
		contact[i].surface.mode = dContactSoftERP | dContactSoftCFM;
		contact[i].surface.mu = refine ? 0.5 : 0;
		contact[i].surface.soft_erp = refine ? 0.1 : 0.9;
		contact[i].surface.soft_cfm = refine ? 0.1 : 0.05;
		dJointID c = dJointCreateContact(world->world, world->jg,
		    contact + i);
		dJointAttach(c, dGeomGetBody(o1), dGeomGetBody(o2));
	}
}

float World::collision() const
{
	const int N = 8;
	dContactGeom contact[N];
	float max = 0.0;

	for (int i = 0; item[i]; ++i)
		for (int j = i + 1; item[j]; ++j) {
			int r = dCollide(item[i], item[j], N, &contact[0],
			    sizeof(dContactGeom));
			for (int k = 0; k < r; ++k)
				if (contact[k].depth > max)
					max = contact[k].depth;
		}
	return max;
}

bool World::tighten(dReal *b)
{
	bool r = false;
	for (int i = 0; i < 3 && wall[i]; ++i) {
		dVector4 p;
		dGeomPlaneGetParams(wall[i], p);
		if (-p[3] > b[i]) {
			dGeomPlaneSetParams(wall[i], p[0], p[1], p[2], -b[i]);
			r = true;
		}
	}
	return r;
}

void World::loosen(float d)
{
	for (int i = 0; i < 3 && wall[i]; ++i) {
		dVector4 p;
		dGeomPlaneGetParams(wall[i], p);
		dGeomPlaneSetParams(wall[i], p[0], p[1], p[2], p[3] - d);
	}
}

int main(int argc, char *argv[])
{
	dInitODE();
	World *first = new World;
	first->read();
	for (time_t end = time(NULL) + 120; time(NULL) < end; ) {
		World *world = new World(*first);
		world->randomize();
		world->run();
		delete world;
	}
	delete first;
	if (best) {
		refine = true;
		best->run();
		best->shift();
		best->print();
		delete best;
	}
	dCloseODE();
}