/* Dyadic grids in arbitrary dimension */
/* Last edited on 2002-04-27 15:47:16 by stolfi */

#include <stdlib.h>
#include <math.h>
#include <js.h>
#include <dygrid.h>

nat8 dyg_depth(dyg_node_index k)
{
  int r = 1;
  while (k > 1) { k = (k/2); r++; }
  return r;
}

dyg_node *dyg_new_node(dyg_node *p)
{
  dyg_node *r = (dyg_node*)notnull(malloc(sizeof(dyg_node)), "no mem for dyg_node");
  r->c0 = NULL;
  r->c1 = NULL;
  r->data = NULL;
  r->index = 0;
  r->p = p;
  r->virtual = false;
  return r;
}

void dyg_free_tree_rec(dyg_node* node)
{
  if( node->c0 || node->c1 )
    {
      dyg_free_tree_rec(node->c0);
      dyg_free_tree_rec(node->c1);
    }
  node->c0 = NULL;
  node->c1 = NULL;
  dyg_free_node(&node);
}


void dyg_free_tree(dyg_header* header)
{
  if( header->root->c0 || header->root->c1 )
    {
      dyg_free_tree_rec(header->root->c0);
      dyg_free_tree_rec(header->root->c1);
    }
  header->root->c0 = NULL;
  header->root->c1 = NULL;
  dyg_free_node(&header->root);
  free(header->inicial_dimension);
}

void dyg_free_node(dyg_node **pp)
{
  dyg_node *p = (*pp);
  affirm(p->c0 == NULL, "node has c0");
  affirm(p->c1 == NULL, "node has c1");
  free(p->data);
  free(p);
  (*pp) = NULL;
}

dyg_header *dyg_new_tree(nat8 d, dyg_dimensions* dim, int inicial_axis)
{
  dyg_header *h = (dyg_header*)notnull(malloc(sizeof(dyg_header)), "no mem for dyg_header");
  dyg_node *r = dyg_new_node(NULL);
  h->d = d;
  h->root = r;

  h->inicial_dimension = (dyg_dimensions*) malloc(sizeof(dyg_dimensions));

  h->inicial_dimension->max[X] = dim->max[X];
  h->inicial_dimension->min[X] = dim->min[X];
  h->inicial_dimension->max[Y] = dim->max[Y];
  h->inicial_dimension->min[Y] = dim->min[Y];

  h->inicial_axis = inicial_axis;
  h->root->index = 1;
  return h;
}
  
void dyg_split_node(dyg_node *p)
{
  dyg_node *c0 = dyg_new_node(p);
  dyg_node *c1 = dyg_new_node(p);
  affirm(p->c0 == NULL, "node has c0");
  affirm(p->c1 == NULL, "node has c1");
  
  c0->index = 2 * p->index;
  c1->index = 2 * p->index + 1;
  p->c0 = c0;
  p->c1 = c1;
}
  
void dyg_unsplit_node(dyg_node *p)
{
  affirm(p->c0 != NULL, "node has no c0");
  affirm(p->c1 != NULL, "node has no c1");
  dyg_free_node(&(p->c0));
  dyg_free_node(&(p->c1));
}

void dyg_shatter_node(
  dyg_node *p, int maxlevel, 
  split_criterion split,
  data_adjust data,
  int* param,
  dyg_dimensions* dim,
  axis long_axis,
  bool virtual
)
/*
  Given a node with no children, replaces it 
  by a random subtree of depth at most maxlevel. */
{
  dyg_dimensions d;

  d.min[X] = dim->min[X];
  d.max[X] = dim->max[X];
  d.min[Y] = dim->min[Y];
  d.max[Y] = dim->max[Y];

  if ((maxlevel > 1) && (split(dim)))
    { 
      dyg_split_node(p);

      p->c0->virtual = virtual;
      p->c1->virtual = virtual;

      *param = data(p);

      if (long_axis == x_axis)
        { 
          double xmid = (dim->max[X] + dim->min[X])/2.0;
	  d.max[X] = xmid;
	  dyg_shatter_node(p->c0, maxlevel-1, split, data, param, &d, y_axis, virtual);
	  d.max[X] = dim->max[X];
	  d.min[X] = xmid;
	  dyg_shatter_node(p->c1, maxlevel-1, split, data, param, &d, y_axis, virtual);
        }
      else
        { 
          double ymid = (dim->max[Y] + dim->min[Y])/2.0;
      	  d.max[Y] = ymid;
          dyg_shatter_node(p->c0, maxlevel-1, split, data, param, &d, x_axis, virtual);
	  d.max[Y] = dim->max[Y];
	  d.min[Y] = ymid;
          dyg_shatter_node(p->c1, maxlevel-1, split, data, param, &d, x_axis, virtual);
        }
     }
  else
    *param = data(p);
} 
  
void dyg_insert_node(dyg_node_list** header, dyg_node* node)
{
  dyg_node_list* aux;
  dyg_node_list* tmp;

  affirm(node != NULL, "node is NULL");

  tmp = (dyg_node_list*)malloc(sizeof(dyg_node_list));
  
  tmp->next = NULL;
  tmp->node = node;

  if(*header == NULL )
    {
      *header = tmp;
      return;
    }
  
  for(aux = *header; aux->next != NULL; aux = aux->next);
  aux->next = tmp;
}

void dyg_flush_list(dyg_node_list* header)
{
  dyg_node_list* next;
  dyg_node_list* aux;

  affirm(header != NULL, "header is NULL");

  for(aux = header; aux != NULL; aux = next)
    {
      next = aux->next;
      free(aux);
    }
}

void dyg_get_cell_dimension(dyg_node_index index, dyg_dimensions* dim,
			   int inicial_axis, dyg_dimensions* ret_dim)
{
  int log2 = (int)(log(index) / log(2));

  ret_dim->min[X] = dim->min[X];
  ret_dim->max[X] = dim->max[X];
  ret_dim->min[Y] = dim->min[Y];
  ret_dim->max[Y] = dim->max[Y];

  if( log2 )
    log2 += (int)floor(index / pow(2, log2 + 1));
  for(; log2 > 0; log2--)
    {
      if( (index >> (log2 - 1)) % 2 )
	{
	  if( inicial_axis == x_axis )
	    ret_dim->min[X] = (ret_dim->min[X] + ret_dim->max[X]) * 0.5;
	  else
	    ret_dim->min[Y] = (ret_dim->min[Y] + ret_dim->max[Y]) * 0.5;
	}
      else
	{
	  if( inicial_axis == x_axis )
	    ret_dim->max[X] = (ret_dim->min[X] + ret_dim->max[X]) * 0.5;
	  else
	    ret_dim->max[Y] = (ret_dim->min[Y] + ret_dim->max[Y]) * 0.5;	  
	}
      inicial_axis = next_axis(inicial_axis);
    }
} 

dyg_node_list* dyg_find_adjacents_cells(
  dyg_node *h, 
  dyg_dimensions* dim,
  int naxis,
  int deep, int max_deep, double p[],
  bool virtual_cells)
{
  dyg_node_list* ret = 0, *next;
  dyg_dimensions d;
  int i;
  bool flag;
  double delta[DIM];

  if( deep > max_deep )
    return ret;

  for(i = 0; i < DIM; i++ )
    delta[i] = (dim->max[i] - dim->min[i]) * 0.5;

  flag = true;
  
  for(i = 0; (i < DIM - 1) && flag; i++ )
    if( LT(p[i], (dim->min[i] - delta[i])) || GT(p[i], (dim->max[i] + delta[i])) )
      flag = false;

  if( flag )
    {
      if( (h->c0 == NULL && h->c1 == NULL) )
	{
	  dyg_insert_node(&ret, h);
	  ret->next = 0;
	}
      else if(!virtual_cells && h->c0->virtual && h->c1->virtual)// É folha
	{
	  dyg_insert_node(&ret, h);
	  ret->next = 0;
	}
      else
	{
          for(i = 0; i < DIM && flag; i++ )
	    d.max[i] = dim->max[i];
          for(i = 0; i < DIM && flag; i++ )
            d.min[i] = dim->min[i];

          if( h->c1 != NULL && (virtual_cells || !h->c1->virtual) )
            {
              d.min[naxis] = (dim->min[naxis] + dim->max[naxis]) * 0.5;
	      ret = dyg_find_adjacents_cells(h->c1, &d, next_axis(naxis), deep + 1, max_deep, p, virtual_cells);
              d.min[naxis] = dim->min[naxis];           
            }

	  d.max[naxis] = (dim->min[naxis] + dim->max[naxis]) * 0.5;
          if( h->c0 != NULL && (virtual_cells || !h->c0->virtual) )
            {
	      if( !ret )
		ret = dyg_find_adjacents_cells(h->c0, &d, next_axis(naxis), deep + 1, max_deep, p, virtual_cells);
	      else
		{
		  next = ret;
		  while (next->next != NULL)
		    next = next->next;
		  next->next = dyg_find_adjacents_cells(h->c0, &d, next_axis(naxis), deep + 1, max_deep, p, virtual_cells);
		}
            }
	}
    } 
  return ret;
}

dyg_node* dyg_find_cell(
  dyg_node *h, 
  dyg_dimensions* dim,
  int naxis,
  double p[],
  bool virtual_cells)
{
  dyg_node* ret = 0;
  dyg_dimensions d;
  int i;
  bool flag;

  flag = true;
  
  for(i = 0; (i < DIM - 1) && flag; i++ )
    if( LT(p[i], dim->min[i]) || GT(p[i], dim->max[i]) )
    	flag = false;

  if( flag )
    {
      if( (h->c0 == NULL && h->c1 == NULL) || (h->virtual == virtual_cells) )// É folha
          return h;
      else
	      {
          for(i = 0; i < DIM && flag; i++ )
	          d.max[i] = dim->max[i];
          for(i = 0; i < DIM && flag; i++ )
            d.min[i] = dim->min[i];

          if( h->c1 != NULL || (h->c1->virtual == virtual_cells) )
            {
              d.min[naxis] = (dim->min[naxis] + dim->max[naxis]) * 0.5;
	            ret = dyg_find_cell(h->c1, &d, next_axis(naxis), p, virtual_cells);
              if( ret )
                return ret;
              d.min[naxis] = dim->min[naxis];           
            }

	        d.max[naxis] = (dim->min[naxis] + dim->max[naxis]) * 0.5;
          if( h->c0 != NULL || (h->c0->virtual == virtual_cells) )
     		      return dyg_find_cell(h->c0, &d, next_axis(naxis), p, virtual_cells);
	      }
    } 
  return 0;
}