/* Dyadic grids in arbitrary dimension */
/* Last edited on 2002-04-27 13:11:20 by stolfi */

#ifndef dygrid_H
#define dygrid_H

#include <stdlib.h>


#define DIM 3

typedef unsigned char nat8;
typedef unsigned long nat64;

typedef enum{X, Y, Z, T} dimensions;

typedef enum{false, true} bool;

#define NUM_AXIS 2
typedef enum{x_axis, y_axis} axis;
#define next_axis(P) ((P+1)%NUM_AXIS)

#define PRECISION 1e-10
#define EQUAL(A,B) (fabs(A - B) < PRECISION)
#define GT(A,B) ((A - B) > PRECISION)
#define LT(A,B) ((B - A) > PRECISION)

/* NODE INDICES */

typedef nat64 dyg_node_index;
/* 
  Every node of the infinite tree has an integer index.
  The root node has index 1.  If $k$ is the index of a node, 
  then $2*k$ and $2*k+1$ are the indices of its children. */
  
nat8 dyg_depth(dyg_node_index k);
/*
  The depth of a node with index $k$.  The root's depth is 1. */
  
#define DYG_MAX_INDEX (1152921504606846975L)
#define DYG_MAX_DEPTH(d) (60/d)

typedef struct dyg_dimensions {
  double min[DIM];
  double max[DIM];
}dyg_dimensions;

/*
  The maximum node index is $2^60-1$, meaning that in any concrete
  tree the maximum depth of an existing node is 60 (counting the depth
  of the root as 1).  
  
  Thus the finest achievable 4-dimensional (3D+time) grid has
  $2^15 = 32768$ cells in each direction. That is enough for 1m
  resolution in a 30 km wide 3-d field, and 10 min resolution in a
  200-day simulation period.
  
  The finest achievable 3-dimensional grid has $2^20 = 1048576$ cells
  in each direction, which means 10 cm resolution in a 100 km wide
  field, and 1 min resolution in a 700-day simulation period. */

/* THE TREE STRUCTURE */

/* It is not obvious whether the tree structure must be stored at all. */

typedef struct dyg_node {
  struct dyg_node  *p;   /* Parent node */
  struct dyg_node  *c0;  /* "Low" children node */
  struct dyg_node  *c1;  /* "High" children node */
  dyg_node_index index;
  bool virtual;
  void*  data; /* "Extra" data */
} dyg_node;

typedef struct dyg_header {
  nat8  d;         /* Dimension */
  dyg_node *root;  /* Root node of the tree */
  dyg_dimensions* inicial_dimension;
  int    inicial_axis;
} dyg_header;

typedef struct dyg_node_list {
  dyg_node *node;
  dyg_dimensions dim;
  struct dyg_node_list *next;
}dyg_node_list;

typedef bool split_criterion(dyg_dimensions*);
typedef int data_adjust(dyg_node*);
/*
  A procedure that returns 1 (true) if the given box needs splitting. */

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
);

dyg_header *dyg_new_tree(nat8 d, dyg_dimensions* dim, int inicial_axis);
/*
  Creates a new $d$-dimensional tree with a single node (the root). */

void dyg_free_tree(dyg_header* header);
void dyg_free_tree_rec(dyg_node* node);


dyg_node *dyg_new_node(dyg_node *p);
/*
  Creates a new node with no children and parent $p$. */
  
void dyg_free_node(dyg_node **pp);
/*
  Reclaims the storage used by node $*pp$, which must not have
  any children, and sets $(*pp) = NULL$. */
  
void dyg_split_node(dyg_node *p);
/*
  The node $p$ must have no children.  Creates two new nodes and
  appends them as children of $p$. */
  
void dyg_unsplit_node(dyg_node *p);
/*
  The node $p$ must have two children but no grandchildren.
  Deletes and reclaims the two children of $p$. */

void dyg_insert_node(dyg_node_list** header, dyg_node* node);

void dyg_flush_list(dyg_node_list* header);


void dyg_get_cell_dimension(
  dyg_node_index index, dyg_dimensions* dim,
  int inicial_axis, dyg_dimensions* ret_dim);

dyg_node_list* dyg_find_adjacents_cells(
  dyg_node *h, 
  dyg_dimensions* dim,
  int inicial_axis,
  int deep, int max_deep, double p[], 
  bool virtual_cells );

dyg_node* dyg_find_cell(
  dyg_node *h, 
  dyg_dimensions* dim,
  int naxis,
  double p[],
  bool virtual_cells );

#endif