/* SOMakeGrid -- Create dyadic grids for testing. */
/* Last edited on 2003-05-07 18:20:45 by ra014852 */

#include <SOGrid.h>
#include <SOParams.h>
#include <SOBasic.h>
#include <js.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* Creates grids for testing purposes.  The user may select
  between a few fixed geometries (a regular grid, a two-region
  one, etc.) or a random tree of a given maximum depth. */

/* INTERNAL PROTOTYPES */

typedef struct Options     /* Parsed command line options */
  { char *outName;
    SOGrid_Dim dim;        /* Grid dimension. */
    /* These are mutually exclusive: */
    bool regular;          /* TRUE asks for a regular grid. */
    bool circle;           /* TRUE for a grid that is denser within a circle. */
    bool front;            /* TRUE for a grid that is denser near a curved line. */
    bool line;             /* TRUE for a grid that is denser near the y = x line. */
    int maxRank;           /* Maximum rank of grid. */
  } Options;

Options *GetOptions(int argn, char **argc);

SOGrid_Tree *SOMakeGrid_Make
  ( SOGrid_Dim dim,        /* Grid dimension. */
    bool regular,          /* TRUE asks for a regular grid. */
    bool circle,           /* TRUE for a grid that is denser within a circle. */
    bool front,            /* TRUE for a grid that is denser near a curved line. */
    bool line,             /* TRUE for a grid that is denser near the y = x line. */
    int maxRank            /* Maximum rank of grid. */
  );
  
typedef bool SplitCriterion(SOGrid_Dim d, SOGrid_Rank r, double *lo, double *hi);
/* A predicate that tells whether the cell with rank {r},
  low corner {lo[0..d-1]}, and high corner {hi[0..d-1]} 
  should be subdivided. */

void SOMakeGrid_shatter_node
  ( SOGrid_Dim d,
    SOGrid_Node *p,
    SOGrid_Rank r,
    double *lo, 
    double *hi,
    SplitCriterion split
  );

/* IMPLEMENTATIONS */

int main(int argn, char **argc)
  {
    Options *o = GetOptions(argn, argc);
    SOGrid_Tree *tree = 
      SOMakeGrid_Make(o->dim, o->regular, o->circle, o->front, o->line, o->maxRank);
    FILE *wr = open_write(txtcat(o->outName, ".tree"));
    SOGrid_Tree_write(wr, tree);
    fclose(wr);
    return 0;
  }
    
SOGrid_Tree *SOMakeGrid_Make
  ( SOGrid_Dim dim,        /* Grid dimension. */
    bool regular,          /* TRUE asks for a regular grid. */
    bool circle,           /* TRUE for a grid that is denser within a circle. */
    bool front,            /* TRUE for a grid that is denser near a curved line. */
    bool line,             /* TRUE for a grid that is denser near the y = x line. */
    int maxRank            /* Maximum rank of grid. */
  )
  { SOGrid_Tree *tree = SOGrid_Tree_new(dim);
    SplitCriterion *toobig;
    double minSize = pow(0.5, ((double)maxRank)/((double)dim)); // Min cell size.
    double lo[dim], hi[dim];
    int i;
    
    auto SplitCriterion regular_criterion, circle_criterion, front_criterion, line_criterion;

    bool regular_criterion(SOGrid_Dim d, SOGrid_Rank r, double *lo, double *hi)
      { return (r < maxRank); }
  
    bool circle_criterion(SOGrid_Dim d, SOGrid_Rank r, double *lo, double *hi)
      { double d2 = 0;  // Distance from root cell center to cell center, squared
        double sz = 0;  // Maximum extent of cell.
        int i;
	if (r >= maxRank) { return FALSE; }
	//	if (r < 5) { return TRUE; }
        // Compute {d2}, {sz}:
        for (i = 0; i < d; i++)
          { double ci = (lo[i] + hi[i])/2.0 - 0.5;
            double szi = hi[i] - lo[i];
            d2 += ci*ci;
            if (szi > sz) { sz = szi; }
          }
        // Decide if cell is small enough:
	//        if(sqrt(d2) < 0.45 && r < 8) { return TRUE; }
        return (sz > (double) d2 * 1.8);
      }
  
    bool front_criterion(SOGrid_Dim d, SOGrid_Rank r, double *lo, double *hi)
      { double d2 = 0;  // Distance from root cell center to cell center, squared
        double sz = 0;  // Maximum extent of cell. 
        int i;
        if (r >= maxRank) { return FALSE; }
        // Compute {d2}, {sz}:
        for (i = 0; i < d; i++)
          { double ci = (lo[i] + hi[i])/2.0 - 0.5;
            double szi = hi[i] - lo[i];
            d2 += ci*ci;
            if (szi > sz) { sz = szi; }
          }
        // Decide if cell is small enough:
        { double dc2 = (d2 - 0.25)*(d2 - 0.25);
          return (sz > minSize + dc2/0.0625*(1.0 - minSize));
        }
      }

    bool line_criterion(SOGrid_Dim d, SOGrid_Rank r, double *lo, double *hi)
      { double dl = 1;  
        double sz = 1;
        int i;

        if (r <= maxRank) { return TRUE; }
        // Compute {dl}, {sz}:
        for (i = 0; i < d; i++) sz *= hi[i] - lo[i];

        dl = fabs(lo[0] - lo[1]) + fabs(hi[0] - hi[1]);

        // Decide if cell is small enough:
        return (sz > (double) dl * 0.02 && sz > 0.001);
      }
  
    if (regular)
      { toobig = &regular_criterion; }
    else if (circle)
      { toobig = &circle_criterion; }
    else if (front)
      { toobig = &front_criterion; }
    else if (line)
      { toobig = &line_criterion; }
    else
      { assert(FALSE, "must specify some split criterion"); }
      
    for (i = 0; i < dim; i++) { lo[i] = 0.0; hi[i] = 1.0; }

    SOMakeGrid_shatter_node(dim, tree->root, 0, lo, hi, toobig);
    
    return tree;
  }

void SOMakeGrid_shatter_node
  ( SOGrid_Dim d,
    SOGrid_Node *p,
    SOGrid_Rank r,
    double *lo, 
    double *hi,
    SplitCriterion *split
  )
  /*
    Given a node with no children, replaces it
    by a random subtree of depth at most maxRank. */
  {
    if (split(d, r, lo, hi))
      { double lo_child[d], hi_child[d];
        int i;
        SOGrid_Node_split(p);
        // Compute low and high coordinates of low child:
        for (i = 0; i < d; i++)
          { lo_child[i] = lo[i]; 
            hi_child[i] = (i == (r%d) ? (lo[i]+hi[i])/2 : hi[i]);
          }
        SOMakeGrid_shatter_node(d, p->c[LO], r+1, lo_child, hi_child, split);
        
        // Compute low and high coordinates of high child:
        for (i = 0; i < d; i++)
          { lo_child[i] = (i == (r%d) ? (lo[i]+hi[i])/2 : lo[i]); 
            hi_child[i] = hi[i];
          }
        SOMakeGrid_shatter_node(d, p->c[HI], r+1, lo_child, hi_child, split);
      }
  }

#define PPUSAGE SOParams_SetUsage

Options *GetOptions(int argn, char **argc)
  {
    Options *o = (Options *)notnull(malloc(sizeof(Options)), "no mem");
    SOParams_T *pp = SOParams_NewT(stderr, argn, argc);

    PPUSAGE(pp, "SOMakeGrid \\\n");
    PPUSAGE(pp, "  -outName NAME \\\n");
    PPUSAGE(pp, "  -dim D \\\n");
    PPUSAGE(pp, "  -maxRank NUM \\\n");
    PPUSAGE(pp, "  { -regular | -circle | -front | -line }");

    SOParams_GetKeyword(pp, "-outName");                               
    o->outName = SOParams_GetNext(pp);  

    SOParams_GetKeyword(pp, "-dim");                               
    o->dim = SOParams_GetNextInt(pp, 1, 4);  

    SOParams_GetKeyword(pp, "-maxRank");                               
    o->maxRank = SOParams_GetNextInt(pp, 0, SOGRID_MAX_RANK);  

    o->regular = FALSE;
    o->circle = FALSE;
    o->front = FALSE;
    o->line = FALSE;
    if (SOParams_KeywordPresent(pp, "-regular"))
      { o->regular = TRUE; }
    else if (SOParams_KeywordPresent(pp, "-circle"))
      { o->circle = TRUE; }
    else if (SOParams_KeywordPresent(pp, "-front"))
      { o->front = TRUE; }
    else if (SOParams_KeywordPresent(pp, "-line"))
      { o->line = TRUE; }
    else
      { SOParams_Error(pp, "must specify type of grid"); }
    
    SOParams_Finish(pp);
    return o;
  }