#include <stdio.h>
#include <math.h>

#include "tectonics.h"
#include "geometry.h"
#include "util.h"
#include "logging.h"

static void rebuild_tree(struct world_t *world)
{
   // remove any existing tree
   if (world->tree.nw != NULL) {
      quadtree_free(world->tree.nw);
      quadtree_free(world->tree.ne);
      quadtree_free(world->tree.sw);
      quadtree_free(world->tree.se);
      world->tree.nw = NULL;
      world->tree.ne = NULL;
      world->tree.sw = NULL;
      world->tree.se = NULL;
   }
   world->tree.id = -1;

   // build a new tree
   for (int i=0; i<world->n_points; i++) {
      quadtree_insert(&(world->tree), world->points, i);
   }
}


/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */

struct centroid_t {
   double x, y;
   int updates;
};

// monte-carlo lloyd's algorithm
static void relax_points(struct world_t *world, int iterations)
{
   // centroids to compute
   struct centroid_t *centroid =
      malloc(sizeof(struct centroid_t) * world->n_points);

   // initialize centroids with current locations
   for (int i=0; i<world->n_points; i++) {
      struct point_t *pt = world->points + i;
      struct centroid_t *c = centroid + i;
      c->x = pt->x;
      c->y = pt->y;
      c->updates = 1;
   }

   // approximate centroids
   for (int i=0; i<iterations; i++) {
      // generate a random point and find the closest
      // terrain point
      double x = rand01();
      double y = rand01();
      int closest = quadtree_get_closest(&(world->tree),
					 world->points,
					 (struct point_t){x, y});

      if (closest == -1) {
	 log_msg(WARN, "bad closest point for (%f, %f)", x, y);
      }
      else {
	 // average the centroid towards that random point
	 struct centroid_t *c = centroid + closest;
	 int u = c->updates;
	 c->x = (u*c->x + x)/(u+1);
	 c->y = (u*c->y + y)/(u+1);
	 c->updates += 1;
      }
   }

   // update positions from centroids
   for (int i=0; i<world->n_points; i++) {
      struct point_t *pt = world->points + i;
      struct centroid_t *c = centroid + i;
      pt->x = c->x;
      pt->y = c->y;
   }

   rebuild_tree(world);

   free(centroid);
}


/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */

void create_world(struct world_t *world, int n_points)
{
   world->n_points = n_points;
   world->points = malloc(sizeof(struct point_t) * n_points);
   if (world->points == NULL) {
      fprintf(stderr,
	      "ERROR: failed to allocate %d points\n",
	      n_points);
      return;
   }

   struct point_t center = { 0.5, 0.5 };
   world->tree = quadtree_new_node(center, 0.5);

   for (int i=0; i<n_points; i++) {
      struct point_t *pt = world->points + i;
      pt->x = rand01();
      pt->y = rand01();
   }

   rebuild_tree(world);

   for (int i=0; i<1; i++) {
      log_msg(INFO, "relaxing points - iteration %d", i+1);
      relax_points(world, 50*world->n_points);
   }
}


/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */

void free_world(struct world_t *world)
{
   if (world->tree.nw != NULL) {
      quadtree_free(world->tree.nw);
      quadtree_free(world->tree.ne);
      quadtree_free(world->tree.sw);
      quadtree_free(world->tree.se);
   }
   free(world->points);
}