/* See dglocus.h */
/* Last edited on 2004-08-24 01:01:29 by stolfi */

#include <dglocus.h>

#include <dgbasic.h>
#include <dgtree.h>
#include <js.h>
#include <math.h>

/* EXPORTED PROCS */

dg_Locus_vec_t dg_find_critical_loci(dg_Dim d, dg_Node *root, bool vtxOnly, bool maximal)
  { dg_Locus_vec_t crit = dg_Locus_vec_new(100);
    int nCrit = 0;

    auto void visit
      ( dg_Locus E, 
        dg_Rank r, 
        dg_NodeStar *NE, 
        dg_LocusData *dtE, 
        bool maxE
      );

    void visit
      ( dg_Locus E, 
        dg_Rank r, 
        dg_NodeStar *NE, 
        dg_LocusData *dtE, 
        bool maxE
      )
      { if ((maximal == maxE) && ((! vtxOnly) || (dg_locus_dimension(d, E) == 0)))
          { if (nCrit >= crit.nel) { dg_Locus_vec_expand(&crit, nCrit); }
            crit.el[nCrit] = E;
            nCrit++;
          }
      }
      
    dg_enum_critical_loci(d, root, visit, NULL, NULL, NULL, NULL);
    dg_Locus_vec_trim(&crit, nCrit);
    return crit;
  }

void dg_enum_critical_loci
  ( dg_Dim d, 
    dg_Node *root,
    dg_LocusVisitProc *visit, 
    dg_RootLocusDataProc *root_data,
    dg_ShrinkLocusDataProc *shrink_data,
    dg_SplitLocusDataProc *split_data,
    dg_FreeLocusDataProc *free_data
  )
  { 
    auto void dg_enum_critical_sub_loci
      ( dg_Locus E, 
        dg_Rank r, 
        dg_NodeStar *NE, 
        dg_LocusData *dtE
      );

    void dg_enum_critical_sub_loci
      ( dg_Locus E, 
        dg_Rank r, 
        dg_NodeStar *NE, 
        dg_LocusData *dtE
      )
      {
        /* fprintf(stderr, "%*s [enter rank = %d locus = ", 2*r, "", r); */
        /* dg_locus_print(stderr, d, E); */
        affirm(dg_node_star_is_valid(r, NE), "invalid locus"); 
        if (dg_node_star_is_leaf(NE))
          { /* Leaf locus. */
            /* fprintf(stderr, " leaf"); */
            visit(E, r, NE, dtE, FALSE);
            /* fprintf(stderr, "]\n"); */
          }
        else
          { /* Locus is not a leaf locus; split cells of {K(E)} and recurse. */
            /* fprintf(stderr, " nonleaf: "); */
            dg_Axis splax = dg_split_axis(d, r); /* Splitting axis. */
            /* Generate subtrees: */
            if (dg_axis_belongs(splax, E.norm))
              { 
                /* The split axis is perpendicular to the item. */
                dg_AxisIndex splj = dg_find_axis(E.norm, splax);
                affirm(splj < d, "inconsistent axis sets");
                /* fprintf(stderr, " shrink on axis %d = Nrm(E)[%d]\n", splax, splj); */                /* fprintf(stderr, "]\n"); */

                /* Get locus {E1} that is the same item in the next level: */
                dg_Locus E1 = dg_locus(E.norm, 2*E.cell);  

                /* Shrink node star around element: */
                dg_NodeRef_vec_t NE1;
                dg_node_star_shrink(NE, splj, &NE1);
                dg_LocusData *dtE1 = 
                  (shrink_data == NULL ? NULL : shrink_data(E, r, splax, NE, dtE));
                /* Recurse on same item in the next level, with same shape. */
                dg_enum_critical_sub_loci(E1, r+1, &NE1, dtE1);
                if (free_data != NULL) { free_data(E1, r+1, &NE1, dtE1); }
                free(NE1.el); 
              }
            else
              { /* The split axis is parallel to item {E}, split it: */
                /* fprintf(stderr, " bisect locus:"); */

                /* Compute normal axes of splitting item: */
                dg_Axes NrmMD = dg_axis_include(splax, E.norm); 

                /* Find rel position {splj} of {splax} among the {Nrm(EMD)}: */
                dg_AxisIndex splj = dg_find_axis(NrmMD, splax);
                affirm(splj < d, "inconsistent axis sets");
                /* fprintf(stderr, " split on axis %d = Nrm(EMD)[%d]\n", splax, splj); */
                /* fprintf(stderr, "]\n"); */

                /* Pieces of this locus: */
                dg_Locus ELO = dg_locus(E.norm, 2*E.cell);
                dg_Locus EHI = dg_locus(E.norm, 2*E.cell + 1);
                dg_Locus EMD = dg_locus(NrmMD, 2*E.cell + 1);

                /* Node stars of the pieces: */
                dg_NodeRef_vec_t NLO, NHI, NMD;
                dg_LocusData *dtLO, *dtHI, *dtMD;

                dg_node_star_split(NE, splj, &NLO, &NHI, &NMD);
                if (split_data != NULL)
                  { split_data(E, r, splax, NE, dtE, &dtLO, &dtHI, &dtMD); }
                else
                  { dtLO = dtHI = dtMD = NULL; }

                /* Recurse for each half of item {F}: */
                visit(ELO, r+1, &NLO, dtLO, TRUE);
                dg_enum_critical_sub_loci(ELO, r+1, &NLO, dtLO);

                visit(EHI, r+1, &NHI, dtHI, TRUE);
                dg_enum_critical_sub_loci(EHI, r+1, &NHI, dtHI);

                /* Recurse for the separating item {S} of {F}. */
                visit(EMD, r+1, &NMD, dtMD, TRUE);
                dg_enum_critical_sub_loci(EMD, r+1, &NMD, dtMD);

                /* Reclaim temporary storage. */
                if (free_data != NULL) 
                  { free_data(EMD, r+1, &NMD, dtMD); 
                    free_data(ELO, r+1, &NLO, dtLO);  
                    free_data(EHI, r+1, &NHI, dtHI);
                  }
                free(NMD.el);
                free(NLO.el);
                free(NHI.el);
              }
          }
        /* fprintf(stderr, "%*s [exit rank = %d locus = ", 2*r, "", r); */
        /* dg_locus_print(stderr, d, E); */
        /* fprintf(stderr, "]\n"); */
      }

    affirm(root != NULL, "tree must be non-null"); 
    affirm(root->index = 1, "must start at root");
    /* Start recursion with root cell and its faces: */
    { /* Allocate a root {dg_NodeStar} of sufficient size: */
      int max_ncells = (1 << d);
      dg_WingIndex ix;
      dg_NodeRef_vec_t NE = dg_NodeRef_vec_new(max_ncells);
      for (ix = 0; ix < max_ncells; ix++) { NE.el[ix] = root; }
      
      /* Enumerate the faces of the root cell: */
      dg_Axes nrm;
      for (nrm = 0; nrm < (1<<d); nrm++)
        { dg_Locus E = dg_locus(nrm,1);
          int m = dg_axes_count(nrm); /* Co-dimension of {E}. */
          NE.nel = (1 << m);
          dg_LocusData *dtE = (root_data == NULL ? NULL : root_data(E, root));
          visit(E, 0, &NE, dtE, TRUE);
          dg_enum_critical_sub_loci(E, 0, &NE, dtE);
          if (free_data != NULL) { free_data(E, 0, &NE, dtE); }
        }
      free(NE.el);
    }
  }

bool dg_node_star_is_valid(dg_Rank r, dg_NodeStar *N)
  { int i;
    for (i = 0; i < N->nel; i++)
      { dg_Node *Ni = N->el[i];
        if ((Ni == NULL) || (dg_cell_rank(Ni->index) != r))
          { return FALSE; }
      }
    return TRUE;
  }

bool dg_node_star_is_leaf(dg_NodeStar *N)
  { int i;
    for (i = 0; i < N->nel; i++)
      { dg_Node *Ni = N->el[i];
        if (Ni->nodes == 1) { return TRUE; }
      }
    return FALSE;
  }