/* See dg_locus.h */
/* Last edited on 2009-08-21 15:28:04 by stolfilocal */

#include "dg_locus.h"
#include "mdg_grid.h"
#include "dg_tree.h"
#include "affirm.h"
#include "set32.h"
#include <math.h>

/* EXPORTED PROCS */

/* LOCUS IN A MULTIGRID */

mdg_dim_t dg_locus_dimension(mdg_dim_t d, dg_locus_t E)
  { return d - set32_count(E.norm); }

mdg_axis_t dg_locus_normal_axis(mdg_dim_t d, dg_locus_t E, mdg_axis_index_t j)
{ return set32_elem_from_index(j,E.norm); }

mdg_axis_t dg_locus_spanning_axis(mdg_dim_t d, dg_locus_t E, mdg_axis_index_t j)
  { mdg_axis_set_t ES = mdg_axis_set_complement(d, E.norm);
    return set32_elem_from_index(j,ES);
  }

dg_locus_t dg_locus_minimize_rank(mdg_dim_t d, dg_locus_t E)
  { mdg_rank_t r = mdg_cell_rank(E.cell);
    mdg_axis_set_t NrmE = E.norm;
    while (set32_belongs((r+d-1) % d, NrmE) && ((E.cell & 1LL) == 0))
      { E.cell >>= 1; r--; }
    return E;
  }

/* FACES AND SUB-FACES OF LOCI */

dg_locus_t dg_locus_facet(mdg_dim_t d, dg_locus_t E, mdg_axis_t ax, box_signed_dir_t dir)
  {
    mdg_axis_set_t da = 1 << ax;
    affirm(ax < d, "bad axis");
    if ((E.norm & da) != 0)
      { /* Locus {E} is already normal to axis {ax}: */
        return E;
      }
    else if (dir == box_NEG)
      { /* Just add {ax} to the normal set: */
        return dg_locus(E.norm | da, E.cell); 
      } 
    else if (dir == box_POS)
      { /* Must increment the cell index along axis {ax}: */
        mdg_cell_index_t mci = mdg_cell_axis_shift(d, E.cell, ax, 1);
        return dg_locus(E.norm | da, mci); 
      }
    else
      { /* Asked {dir == box_ZER}, means {E} itself: */
        return E; 
      } 
  }

dg_locus_t dg_locus_face(mdg_dim_t d, dg_locus_t E, box_face_index_t fi)
  { int ax = 0;
    mdg_axis_set_t ES = mdg_axis_set_complement(d, E.norm);
    while (ES != 0)
      { if (ES & 1)
          { box_signed_dir_t dir = ((fi + 1) % 3) - 1; 
            if (dir != box_ZER) { E = dg_locus_facet(d, E, ax, dir); }
            fi /= 3;
          }
        ax++; ES >>= 1;
      }
    return E;
  }

dg_locus_t dg_locus_extend(mdg_dim_t d, dg_locus_t E, mdg_axis_t ax, box_signed_dir_t dir)
  {
    mdg_axis_set_t da = 1 << ax;
    affirm(ax < d, "bad axis");
    if ((E.norm & da) == 0)
      { /* Locus {E} is already parallel to axis {ax}: */
        return E;
      }
    else if (dir == box_POS)
      { /* Just remove {ax} from the normal set: */
        return dg_locus(E.norm & (~da), E.cell); 
      } 
    else if (dir == box_NEG)
      { /* Must decrement the cell index along axis {ax}: */
        return dg_locus(E.norm & (~da), mdg_cell_axis_shift(d, E.cell, ax, -1)); 
      }
    else
      { /* Asked {dir == box_ZER}, means {E} itself: */
        return E; 
      } 
  }

dg_locus_t dg_locus_face_from_signature(mdg_dim_t d, dg_locus_t E, box_signed_dir_t dir[])
  { mdg_axis_index_t j = 0; mdg_axis_t ax = 0;
    mdg_axis_set_t ES = mdg_axis_set_complement(d, E.norm);
    while (ES != 0)
      { if (ES & 1) 
          { box_signed_dir_t dirj = dir[j];
            if (dirj != box_ZER) { E = dg_locus_facet(d, E, ax, dirj); }
            j++;
          }
        ax++; ES >>= 1;
      }
    return E;
  }
  
void dg_locus_box_root_relative(mdg_dim_t d, dg_locus_t E, interval_t B[])
  {
    int i;
    mdg_cell_box_root_relative(d, E.cell, B);
    for (i = 0; i < d; i++)
      { if (set32_belongs(i, E.norm)) 
          { HI(B[i]) = LO(B[i]); }
      }
  }

void dg_locus_star_size(mdg_dim_t d, dg_locus_t E, mdg_grid_size_t psz[])
  { int i;
    for (i = 0; i < d; i++)
      { psz[i] = (set32_belongs(i, E.norm) ? 2 : 1); }
  }

void dg_locus_star(mdg_dim_t d, dg_locus_t E, mdg_cell_index_t k[])
  {
    mdg_grid_pos_t pos[mdg_MAX_GRID_DIM];
    mdg_cell_position(d, E.cell, pos);
    int i, ix;
    k[0] = E.cell; ix = 1;
    for (i = 0; i < d; i++)
      { if (set32_belongs(i, E.norm))
          { /* Copy all cells already found, shift the rest down: */
            int ix0;
            for (ix0 = 0; ix0 < ix; ix0++)
              { k[ix + ix0] = k[ix0];
                k[ix0] = mdg_cell_axis_shift(d, k[ix0], i, -1);
              }
            ix = 2*ix;
          }
      }
  }

/* PRINTOUT */

void dg_locus_print(FILE *wr, mdg_dim_t d, dg_locus_t E)
  { 
    int i;
    mdg_grid_pos_t Ep[mdg_MAX_GRID_DIM];
    mdg_cell_position(d, E.cell, Ep);
    mdg_cell_print(wr, d, E.cell);
    fprintf(wr, ":");
    for (i = 0; i < d; i++) 
      { fprintf(wr, "%d", (set32_belongs(i, E.norm) ? 1 : 0)); }
  }

/* MANIPULATION OF LOCUS VECTORS */

vec_typeimpl(dg_locus_vec_t,dg_locus_vec,dg_locus_t);