/* See SOGrid.h */
/* Last edited on 2003-07-08 00:22:35 by stolfi */

#include <SOGrid.h>
#include <js.h>
#include <filefmt.h>
#include <fget.h>
#include <nget.h>
#include <stdlib.h>

/* INTERNAL PROTOTYPES */

void SOGrid_subtree_write(FILE *wr, SOGrid_Node *p);
  /* Writes the subtree rooted at {p} to {wr}. If {p} is 
    NULL, writes "*"; else writes "(A,B)" where "A" and "B"
    are the low and high subtrees of {p}. */

SOGrid_Node *SOGrid_subtree_read(FILE *rd, SOGrid_Node *p, SOGrid_Index k);
  /* Reads a subtree from {rd}. If the subtree is not NULL, its root
    will have parent {p} and index {k}. */

/* IMPLEMENTATIONS */

SOGrid_Rank SOGrid_rank(SOGrid_Index k)
  { int r = 1;
    while (k > 1) { k = (k/2); r++; }
    return r;
  }

SOGrid_Index SOGrid_min_index(SOGrid_Rank r)
  { return ((SOGrid_Index)1) << r; }

SOGrid_Index SOGrid_max_index(SOGrid_Rank r)
  { return (((SOGrid_Index)1) << (r + 1)) - 1; }

SOGrid_Axis SOGrid_longest_axis(SOGrid_Dim d, SOGrid_Rank r)
  { return r % d; }

void SOGrid_Node_split(SOGrid_Node *n)
  { assert(n->c[LO] == NULL, "node has c[LO]");
    assert(n->c[HI] == NULL, "node has c[HI]");
    n->c[LO] = SOGrid_Node_new(n, 2 * n->index);
    n->c[HI] = SOGrid_Node_new(n, 2 * n->index + 1);
    n->count = 3;
  }

void SOGrid_Node_unsplit(SOGrid_Node *n)
  { assert(n->c[LO] != NULL, "node has no c[LO]");
    assert(n->c[HI] != NULL, "node has no c[HI]");
    SOGrid_Node_free(n->c[LO]); n->c[LO] = NULL; 
    SOGrid_Node_free(n->c[HI]); n->c[HI] = NULL;
  }

/* ALLOCATION */

SOGrid_Tree *SOGrid_Tree_new(SOGrid_Dim d)
  { void *v = notnull(malloc(sizeof(SOGrid_Tree)), "no mem for SOGrid_Tree");
    SOGrid_Tree *h = (SOGrid_Tree*)v;
    h->d = d;
    h->root = SOGrid_Node_new(NULL, 1);
    return h;
  }
  
void SOGrid_subtree_free(SOGrid_Node* root)
  { 
    if (root != NULL)
      { 
        SOGrid_subtree_free(root->c[LO]); 
        root->c[LO] = NULL;
        SOGrid_subtree_free(root->c[HI]); 
        root->c[HI] = NULL;
        SOGrid_Node_free(root);
      }
  }

void SOGrid_Tree_free(SOGrid_Tree *t)
  { if (t != NULL)
    { SOGrid_subtree_free(t->root);
      free(t);
    }
  }

#define SOGrid_Tree_FileFormat "2003-05-07"
  
void SOGrid_Tree_write(FILE *wr, SOGrid_Tree *t)
  { filefmt_write_header(wr, "SOGrid_Tree", SOGrid_Tree_FileFormat);
    fprintf(wr, "dim = %d\n", t->d);
    SOGrid_subtree_write(wr, t->root);
    filefmt_write_footer(wr, "SOGrid_Tree");
    fflush(wr);
  }
  
void SOGrid_subtree_write(FILE *wr, SOGrid_Node *p)
  { if (p == NULL)
      { fputc('*', wr); }
    else
      { fputc('(', wr);
        SOGrid_subtree_write(wr, p->c[LO]); 
        fputc(',', wr);
        SOGrid_subtree_write(wr, p->c[HI]); 
        fputc(')', wr);
      }
  }    
  
SOGrid_Tree *SOGrid_Tree_read(FILE *rd)
  { int dim;
    SOGrid_Tree *tree;
    filefmt_read_header(rd, "SOGrid_Tree", SOGrid_Tree_FileFormat);
    dim = nget_int(rd, "dim"); fget_eol(rd);
    assert(dim >= 0, "bad dimension");
    tree = SOGrid_Tree_new(dim);
    fget_skip_formatting_chars(rd);
    tree->root = SOGrid_subtree_read(rd, NULL, 1);
    fget_skip_formatting_chars(rd);
    filefmt_read_footer(rd, "SOGrid_Tree");
    if (rd != stdin) { fclose(rd); }
    return tree;
  }
  
SOGrid_Node *SOGrid_subtree_read(FILE *rd, SOGrid_Node *p, SOGrid_Index k)
  { int c = fget_char(rd);
  if (c == '*')
    { return NULL; }
  else if (c == '(')
    { SOGrid_Node *this = SOGrid_Node_new(p, k);
      // Read low subtree:
      this->c[LO] = SOGrid_subtree_read(rd, this, 2*k);
      // Check for comma:
      fget_skip_formatting_chars(rd);
      c = fget_char(rd);
      assert(c == ',', "comma expected");
      fget_skip_formatting_chars(rd);
      // Read high subtree:
      this->c[HI] = SOGrid_subtree_read(rd, this, 2*k + 1);
      // Check for close parenthesis:
      fget_skip_formatting_chars(rd);
      c = fget_char(rd);
      assert(c == ')', " expected ')'");
      return this;
    }
  else
    { assert(FALSE, "expected '*' or '('");
      return NULL;
    }
  }

SOGrid_Node *SOGrid_Node_new(SOGrid_Node *p, SOGrid_Index k)
  { void *v = notnull(malloc(sizeof(SOGrid_Node)), "no mem for SOGrid_Node");
    SOGrid_Node *n = (SOGrid_Node*)v;
    n->p = p;
    n->c[LO] = NULL;
    n->c[HI] = NULL;
    n->index = k;
    n->count = 1;
    return n;
  }

void SOGrid_Node_free(SOGrid_Node *node)
  { 
    assert(node->c[LO] == NULL, "node has c[LO]");
    assert(node->c[HI] == NULL, "node has c[HI]");
    free(node);
  }

void SOGrid_cell_coords(
  SOGrid_Index k,
  SOGrid_Rank r,
  SOGrid_Dim d,
  double *min,
  double *max
)
{ int fulldivs = r/d;  /* Number of whole division cycles. */
  int j;
  for(j = 0; j < d; j++)
    { int jrank = fulldivs + ((r%d) > j ? 1 : 0); /* Rank of cell along axis {j} */

      double jwidth = 1.0/(double)(1 << jrank);   /* Width of cell (rel root). */

      /* Find index {jindex} of cell along axis {j}: */
      int jindex = 0, z = (k - (1 << r))>> j, p = 1;
      int ax = (r+d+j-1)%d; /* The {+d} is needed if {r == 0}. */
      while (z != 0) { if (z & 1) { jindex = jindex | p; } z >>= d; p <<= 1; }
      min[ax] = (double)jindex * jwidth; 
      max[ax] = min[ax] + jwidth;     
    }	  
}

/*
SOGrid_NodeRefVec SOGrid_NodeRefVec_new(nat nel)
  {  // This is not a macro only because gcc does not allow cast of struct: 
    vec_t v = vec_new(nel, sizeof(SOGrid_Node *));
    SOGrid_NodeRefVec r;
    r.nel = v.nel; r.el = (SOGrid_Node **)v.el;
    return r;
  }
*/