#ifndef mdg_grid_H
#define mdg_grid_H

/* Abstract multi-dimensional dyadic grids */
/* Last edited on 2009-08-21 18:25:24 by stolfilocal */

#define _GNU_SOURCE
#include <stdlib.h>
#include <stdint.h>

#include "bool.h"
#include <assert.h>
#include "sign.h"
#include "interval.h"
#include "box.h"
#include "vec.h"
#include "affirm.h"
#include "bool.h"
#include "indexing.h"
#include <stdlib.h>

#include "udg_grid.h"

/* 
   SPACE DIMENSIONS AND COORDINATE AXES */

#define mdg_MAX_GRID_DIM (4)
#define mdg_MAX_GRID_AXIS (mdg_MAX_GRID_DIM-1)
  /* Just to be safe. High grid dimensions {d} are dangerous because
    too many things have size {c^d} for some {c>1}.  In particular, 
    a {box_face_index_t} ranges from {0..3^d-1}. */

typedef box_dim_t mdg_dim_t; 
  /* Dimension of a box, space, etc.. */

typedef box_axis_t mdg_axis_t; 
  /* Identifies a coordinate axis of the enclosing space, from 0. */

typedef box_axis_set_t mdg_axis_set_t; 
  /* A packed subset of the coordinate axes. */

typedef box_axis_index_t mdg_axis_index_t; 
  /* Identifies a coordinate axis amonga set: 0 is the lowest. */

#define mdg_NO_AXIS (BOX_NO_AXIS)
  /* A {mdg_axis_t} value that means "no axis". */

#define mdg_AXIS_NOT_FOUND (BOX_AXIS_NOT_FOUND)
  /* A {mdg_axis_t} value that means "no such axis" */
  
box_axis_set_t mdg_axis_set_complement(mdg_dim_t d, mdg_axis_set_t A);
  /* The complement of {A} relative to the set {0..d-1}. */

/*
  MULTI-DIMENSIONAL GRIDS
  
  A /{d}-dimensional grid/ can be seen as the Cartesian product of {d}
  partitions that are unidimensional grids (see {udg_grid.h}). It is a
  partition {G} of {R^d} into boxes of various dimensions, the /items/
  of the grid, such that 
  
    (1) the intersection of the closure of any two items of {G} 
    is the union of items of {G};
    
    (2) for any axis {i}, the projections of all items parallel
    to {i} are open intervals with the same length {delta[i]}
    -- the /step/ of the grid along that  axis; and 
    
    (3) one of the items consists of a single point, the
    origin of {R^d}.
  
  A /cell/ of a grid is an item {C} with maximum dimension, i.e.
  {dim(C) == d}. Grid items with dimensions 0, 1, 2, 3 are called
  /vertices/,  /edges/, /faces/, and /blocks/ of the grid, respectively.
  
  It follows from this definition that all items with the same
  oriantation (same set of normal axes) are congruent translates of each
  other.  Indeed, any item {I} can be identified by its orientation and 
  the unique vertex that is the inferior corner of {I} (or, equivalently, by the 
  unique cell which has {I} as a inferior face).
  
  It folows also that every item with dimension {m < d} lies on the
  boundary of {2^{d-m}} items with dimension {m+1}.
  
  GRID POSITION VECTORS
  
  A vertex {v} on a fixed grid can be identified by an integer
  /position vector/ {pos[0..d-1]}, where {pos[i]} is the displacement
  from the origin to {v} along the coordinate axis {i}, divided by the
  corresponding grid step {delta[i]}. In particular, the vertex at the
  origin of {R^d} has position vector {(0,..0)}.
  
  Similarly, a cell can be identified by the position of its inferior
  corner. In particular, the cell whose inferior corner is the origin
  (the /base cell/ of the grid) has position vector {(0,..0)}.
  
  More generally, an item can be identified by its set of normal or
  spanning vectors, and by the position vector of its inferior corner. */

typedef udg_grid_pos_t mdg_grid_pos_t;
  /* Displacement of a grid item along an axis, in multiples of 
    the corresponding grid step.  Also relative displacement
    between two items. */
  
/*
  TOROIDAL TOPOLOGY
  
  A (/canonical/) /{d}-torus/ is the Cartesian product of {d} flat
  circles (see {cp_pulse.h}). It is the Cartesian product of {d}
  intervals of the form {[0_wd[i])}, which is a {d}-dimensional box
  with inferior corner {(0,0,...)} plus all the inferior faces of that box.
  The topology is such that every point on any superior face of that box
  is identified with the corresponding point on the opposite inferior face.
  
  Each number {wd[i]} is called the /period/ or /extent/ of the torus
  along axis {i}.

  It can be understood as the quotient of {R^d} by an equivalence
  relation that equates any two points of {R^d} which, along every
  coordinate axis {i}, differ by an integer multiple of {wd[i]}.
 
  TOROIDAL GRIDS
  
  A /{d}-dimensional toroidal grid/ can be seen as the Cartesian
  product (in the partition sense) of {d} unidimensional circular
  grids (see {cp_pulse.h}). It consists of a finite subset {G} of
  items of an (infinite) {d}-dimensional grid, whose union is a
  {d}-dimensional torus {D} (the the /domain/ of the grid).
  
  The period {wd[i]} of the torus along axis {i} must be some integer
  multiple {gsz[i]*delta[i]} of the grid's step {delta[i]} along that
  axis. The integer {gsz[i]} is called the grid's /size/ along axis
  {i}. When {d==1}, we often say /circular/ instead of toroidal.
  
  Note that in a toroidal grid, as in an infinite grid, all items of a
  given dimension and orientation are topologically equivalent. */
  
typedef udg_grid_size_t mdg_grid_size_t;
  /* Size of a toroidal grid along an axis, in multiples of 
    the corresponding grid step. */
  
/*
  THE DYADIC MULTIGRID

  An (/infinite dyadic/) /{d}-dimensional multigrid/ {G*} is the union
  of an infinite collection of {d}-dimensional toroidal grids -- the
  /levels/ or /layers/ of {G*} -- each identified by a natural number,
  its /rank/. 

  By definition, level 0 has a single cell, the /root cell/,
  the {d}-dimensional cube {U^d = (0_1)^d}. Every cell in level {r} contains
  exactly two cells of rank {r+1}, its /children/; and, except for the
  root, every cell is contained in exactly one cell of rank {r-1}, its
  /parent/.  The two children of a cell are /sisters/ to each other,
  and the root cell is its own sister by definition.
  
  The children of a cell are obtained by splitting {C} in half,
  perpendicularly to axis {r \bmod d}. The half-cell that lies towards
  {-oo} along that axis is the /low child/ or /0-/child, and the other
  half is the /high child/ or /1-child/.
  
  Note that all levels of a multigrid {G*} cover the same domain {D =
  \U G* = [0_1)^d}, namely the union of the root cell {U^d} plus all
  its inferior faces, with toroidal topology. */

typedef udg_rank_t mdg_rank_t;
  /* Identifies a layer in a dyadic multigrid, or the depth of a node
    in a binary tree, etc. It may also be a relative rank (i.e.
    difference of two ranks). */

mdg_axis_t mdg_split_axis(mdg_dim_t d, mdg_rank_t r);
  /* The axis perpendicular to the cut that divides a {d}-dimensional
    cell of level {r} into its children cells; namely, axis {r \bmod d}. */

void mdg_grid_size(mdg_dim_t d, mdg_rank_t r, mdg_grid_size_t gsz[]);
  /* Stores in {gsz[0..d-1]} the number of grid cells along axis {i} in
    level {r} of the {d}-dimensional multigrid. */
    
mdg_grid_pos_t mdg_max_grid_pos(mdg_dim_t d, mdg_rank_t r, mdg_axis_t i);
  /* Maximum position along axis {i} of level {r} in a
    {d}-dimensional multigrid. */
    
#define mdg_MAX_RANK (udg_MAX_RANK)
  /* Maximum rank of any level in a multigrid. The limit 
    is such that a {mdg_cell_index_t} (below) will fit
    in one 64-bit word.

    In a 4-dimensional domain (for example, 3D space + time), the
    finest level has at least {2^15 = 32768} cells along each axis.
    That is enough for 1m resolution in a field 30 km wide, together
    with 5 min resolution in a 100-day simulation period.

    In a 3-dimensional domain (3D space, or 2D space + time), the
    finest level has {2^21 = 2097152} cells along each axis, which
    means 5 cm resolution in a field 100 km wide, together with 1 min
    resolution in a 1400-day simulation period.
    
    In a 2-dimensional domain (2D space, or 1D space + time), the
    finest level has at least {2^31 = 2147483648} cells along each
    axis, which means 1 mm resolution in a field 2000 km wide,
    together with 1 sec resolution in a 68-year simulation period.
    
    In a 1-dimensional domain (space or time), the finest level (rank
    62) has {2^62 = 4.6�10^18} cells, which means 10 picometer
    (sub-atomic) resolution over the Earth's circumference, or 1
    millisecond resolution in a 150,000-year simulation period. That
    should be enough for most practical purposes. */

/*
  CELL INDICES
  
  We can uniquely identify a cell of a fixed multigrid {G*}, at any
  level, by a a single /cell index/ {k}. By definition, the root cell
  in layer 0 has index 1; and the children of cell {k} are the cells
  with indices {2*k} (low) and {2*k + 1} (high). The parent of cell
  {k} has therefore index {floor(k/2)}. Note that these formulas are
  independent of the grid's dimension {d}.
  
  These rules assign consecutive indices to all cells of each layer of
  {G*}; in layer {r} the indices range from {2^r} to {2^{r+1}-1},
  inclusive. Therefore, the position of the leading 1 bit in the index
  of a cell is its rank. The remaining bits are the bits of the cell's
  position vector, interleaved: that is, bit {2^j} of {pos[i]} is bit
  {2^{d*j + i}} of the cell index {k}, for each axis {i} and each bit
  position {j} such that {d*j + i < r}.  */
    
typedef udg_cell_index_t mdg_cell_index_t;
  /* Identifies a cell in some multigrid. */

#define mdg_NO_CELL (udg_NO_CELL)
  /* A cell index that means /no such cell/. */

#define mdg_ROOT_CELL (1LL)
  /* The index of the root cell of any multigrid. */

#define mdg_cell_min_index(r) (1LL << (r))
#define mdg_cell_max_index(r) ((1LL << (r+1))-1LL)
  /* Cells in level {r} have conscutive indices in the range
    {[mdg_cell_min_index(r) .. mdg_cell_max_index(r)]}. The minimum index is
    the base cell of layer {r}, the maximum index is the cell
    adjacent to the superior corner of the domain. */

mdg_rank_t mdg_cell_rank(mdg_cell_index_t k);
  /* The rank of the cell whose index is {k}. */
  
#define mdg_check_cell_rank(k,r)   (((k)>>(r)) == 1LL)
  /* TRUE iff cell {k} belongs to level {r}, */

void mdg_cell_position(mdg_dim_t d, mdg_cell_index_t k, mdg_grid_pos_t pos[]);
  /* For each axis {i} in {0..d-1}, stores in {pos[i]} the 
    position of cell {k} within its layer. */
    
mdg_cell_index_t mdg_cell_from_position(mdg_dim_t d, mdg_rank_t r, mdg_grid_pos_t pos[]);
  /* Returns the index of the cell of level {r} whose position along axis {i}
    is {pos[i]}, for {i} in {0..d-1}. */
    
mdg_cell_index_t mdg_cell_axis_shift(mdg_dim_t d, mdg_cell_index_t k, mdg_axis_t ax, mdg_grid_pos_t dp);
  /* The index of the cell that is displaced from cell {k} by {dp}
    grid steps along axis {ax}. Takes into account the toroidal
    topology of each layer. */
    
mdg_cell_index_t mdg_cell_shift(mdg_dim_t d, mdg_cell_index_t k, mdg_grid_pos_t dp[]);
  /* Returns the index of the cell in a {d}-dimensional multigrid
    whose position vector is {dp[0..d-1]} relative to the inferior
    corner of cell {k}, within the layer {r} that contains {k}. Takes
    into account the toroidal topology of each layer. */
    
mdg_cell_index_t mdg_cell_index_add(mdg_dim_t d, mdg_cell_index_t ka, mdg_cell_index_t kb);
  /* The index of the cell whose position vector is the sum of the
    position vectors corresponding to indices {ka} and {kb}. The 
    two cells should have the same rank. */

mdg_cell_index_t mdg_cell_neighbor(mdg_dim_t d, mdg_cell_index_t k, box_face_index_t fi);
  /* The cell {k'} of rank as {k}, such that the highest-dimensional
    face shared by {k} and {k'} is face {fi} of {k}. In particular, if
    {fi == 0}, returns {k} itself. */

/* 
  ROOT-RELATIVE COORDINATES 
  
  The coordinates of a point {x} of {D} are said to be
  /root-relative/, because in them the inferior and superior corners of the
  root cell {U^d} are {(0,..0)} and {(1,..1}). Root-relative
  coordinates that differ by an integer are implicitly equivalent, and
  therefore they are often reduced modulo 1 to real fractions in the
  range {[0_1)}. */

void mdg_cell_box_root_relative(mdg_dim_t d, mdg_cell_index_t k, interval_t B[]);
  /* Stores in {B[0..d-1]} the inferior and superior coordinates of the cell
    {k}, relative to the root cell, in a {d}-dimensional grid. The
    inferior bounds will be in the range {[0_1)}; the superior bounds will
    be strictly greater than the inferior bounds, in the range
    {(0_1]}. */

/*
  THE CANONICAL GEOMETRY

  In the /canonical geometry/ of the dyadic hierarchy, 
  the grid of level 0 includes the /canonical root cell/ {C0},
  whose projection along axis {i} is the interval {C0_i = [0 _ s^i]}
  where {s = (1/2)^(1/d)}. It follows that every cell of any rank {r}
  is a copy of {C0}, uniformly scaled by a factor {s^r}, and with its
  axes cyclically permuted {r} times.
  
  The main purpose of the canonical geometry is to justify and
  visualize the topological features of the grid, in particular the
  cyclic permutation of the splitting axis as one descends in the
  multigrid. It is not necessarily the geometry used in an actual
  application, where the grid may be arbitrary distorted by some
  `shape function'. */

const double *mdg_cell_root_canonical_sizes(mdg_dim_t d);
  /* Returns a pointer to a (statically allocated) vector {sz[0..d-1]}
    such that {sz[i]} is the extent along axis {i} of the canonical
    root cell of the {d}-dimensional grid. Clients should treat that
    vector as read-only. */

void mdg_cell_box_canonical(mdg_dim_t d, mdg_cell_index_t k, interval_t B[]);
  /* Stores in {B[0..d-1]} the inferior and superior coordinates
    of the cell {k} in a {d}-dimensional grid. */

/*
  CELL PACKS
  
  A /cell pack/ is an array {\X} of cells in a given layer {r} of
  {G*}, defined by a position vector {pos[0..d-1]} in that layer and a /pack size
  vector/ {psz[0..d-1]} of positive integers. In the procedures below, a
  pack is usually specified by its size vector {psz} and the cell index
  {b} of the cell whose position is {pos}.
  
  Let {gsz[0..d-1]} be the size of layer {r} along each axis. A cell {C'}
  of layer {r} belongs to the pack {\X} if and only if its position
  vector {pos'} satisfies {pos'[i] = pos[i] + e[i] (mod gsz[i])} for
  every {i}, where each {e[i]} is some integer in the range {0..psz[i]-1}.
  The vector {e} is the /relative position/ of {C'} within the pack.

  DUPLICATE ENTRIES
  
  A cell pack is defined as an array, so it has exactly {N =
  psz[0]*..psz[d-1]} entries. This is true even when {psz[i] > gsz[i]} for
  some {i}. In this case the pack wraps around and overlap itself,
  and includes some cells two or more times, with different relative
  position vectors {e}.
  
  PACK INDICES
  
  Each entry of the cell pack is identified by a /slot index/ {pkx}
  in {0..N-1}, which is related to its relative position vector
  {e[0..d-1]} and the pack sizes {psz[0..d-1]} by the function
  {ix_pos_by_col}.
  
  The lower cell of the pack has relative position {e = (0,..0)} and
  slot index {0}. The cell with relative position {e[i] = psz[i]-1} for
  all {i} is the /upper cell/ of the pack, and has slot index {N-1}. */
  
typedef uint64_t mdg_pack_slot_index_t;
  /* Identifies a cell slot in the pack of some locus {E}. */

#define mdg_NOT_SLOT ((mdg_pack_slot_index_t)-1)
  /* The {mdg_pack_slot_index_t} of a cell that is not in the pack. */

uint64_t mdg_pack_slot_count(mdg_dim_t d, mdg_grid_size_t psz[]);
  /* Returns the number of slots in a cell pack with sinze
    {psz[i]} along each axis {i}, namely  {PROD { psz[i] : i \in 0..d-1 }}. */

void mdg_pack_cells
  ( mdg_dim_t d, 
    mdg_cell_index_t b,
    mdg_grid_size_t psz[], 
    mdg_cell_index_t k[]
  );
/* Returns in {k[0..N-1]} the indices of all cells in the pack whose
  lower cell is {b} and whose size along axis {i} is {psz[i]}, for {i}
  in {0..d-1}; where {N = mdg_pack_slot_count(d,psz)}. */

mdg_pack_slot_index_t mdg_pack_find_cell
  ( mdg_dim_t d, 
    mdg_cell_index_t k, 
    mdg_cell_index_t b,
    mdg_grid_size_t psz[]
  );
  /* Returns the slot index {pkx} of cell {k} in the cell pack with
    lower cell {b} and size vector {psz[0..d-1]}. If {k} occurs
    multiple times in the pack, returns the first occurrence. If {k}
    is not part of the pack, returns {mdg_NOT_SLOT} */
   
bool_t mdg_pack_intersects_cell
  ( mdg_dim_t d, 
    mdg_cell_index_t k, 
    mdg_cell_index_t b,
    mdg_grid_size_t psz[]
  );
  /* Returns TRUE if cell {k} intersects some cell of 
    the pack with lower cell {b} and size {psz[0..d-1]}. */
   
void mdg_pack_box_root_relative
  ( mdg_dim_t d, 
    mdg_cell_index_t b,
    mdg_grid_size_t psz[], 
    interval_t B[]
  );
  /* Stores in {B[0..d-1]} the lower and upper coordinates of the box
    of {R^d} that corresponds to the cell pack with lower cell {b}
    and size vector {psz[0..d-1]}.
    
    The lower endpoints will be in the range {[0_1)}; the upper
    endpoint will be strictly greater. Due to the toroidal geometry,
    the upper endpoint {HI(B[i])} will be 1 or more whenever the
    position of {b} along axis {i} is {gsz[i]-psz[i]} or more,
    where {gsz[i]} is the grid's size along that axis. */

/* 
   MISCELLANEOUS */

sign_t mdg_rank_cmp(mdg_rank_t *x, mdg_rank_t *y);
  /* Returns {-1} if {*x < *y}, {0} if {*x == *y}, {+1} if {*x > *y}.
    Useful for sorting things by rank. */ 

/* 
  PRINTOUT */

void mdg_cell_print(FILE *wr, mdg_dim_t d, mdg_cell_index_t k);
  /* Prints the cell {k} in the format 
    "{RANK}:({POS[0]},..{POS[d-1]})" */

vec_typedef(mdg_cell_index_vec_t,mdg_cell_index_vec,mdg_cell_index_t);
     /* A vector of {mdg_cell_index_t}. */

#endif