/* See {Heuristics.h} */
#include <Heuristics.h>

#include <Octf.h>
#include <Triangulation.h>
#include <r4.h>
#include <r4x4.h>
#include <Random.h>
#include <math.h>
#include <Debug.h>
#include <rn.h>
#define _GNU_SOURCE
#include <stdio.h>
 
CONST debug == FALSE;

TYPE
  Node == Triangulation.Node;
  Wall == Triangulation.Wall;
  @{Edge->?} == Triangulation.edge;
  @{EDGE->?}S == Edge_vec_t;
  WALLS == Wall_vec_t;
  Walls == REF WALLS;
  DRF == REF uint_vec_t;

  /*
r4_t Dir(v: r4_t)
  /* Same as r4_Dir(v), but returns random unit vector if "v" is zero. */
  { ??? s = r4_Norm(v);
    {
      if (s == 0.0)
      {
        return (r4_t){1.0, 0.0, 0.0, 0.0};
      }
      else
      {
        return r4_Scale(1.0/s, v);
      }
    }    
  } /* END Dir */

Node ChooseReferenceNode(
     Place_t @p; 
     bool_vec_t *variable; 
     rand_t *coins
 )
  /*
    If any of the edges out of "OrgV(a)" has a fixed tip, picks 
    one such @{edge->?}, with probability proportional to the
    number of "ONext"s until the next fixed @{edge->?}.
    Otherwise picks a random arc out of "OrgV(a)".  */
      v : REF Node_vec_t;
  { ??? star = Triangulation.NeighborNode(a,top);
    ??? nv = Triangulation->numberNeighborNode(star);
    ??? k = coins.integer(0, nv-1);
    {
      Node_vec_t v = Node_vec_new(nv);
      for (j = 0; j < nv; j++)
      {
        v[j] = star[j];
      }
      for (i = 0; i < k; i++)
      {
        if (! variable[v[i]->num]){ return a; }
        a = OPrev(a) 
      }
      return a;
    }
  } /* END ChooseReferenceArc */
*/


@{EDGE->?}S *Collect@{Edge->?}s(*u: Node; ElemTableRec_t *top)==
  uint NT = 0;
  { ??? NE = top->NE;
    ??? t = Edge_vec_new(NE)^;
    {
      for (i = 0; i < NE; i++)
      {
        ??? e = top->edge[i];
          ??? eun = e.node[0]->num;
        ??? evn = e.node[1]->num;
        {
          if (((eun == u->num) || (evn == u->num)))
          {
            t[NT] = e;
            NT++;
          } 
        }
      }

      ??? r = Edge_vec_new(NT);
    {
        r^ = SUBARRAY(t,0,NT);
        return r;
      }
    }
  } /* END Collect@{Edge->?}s */

WALLS *CollectWalls(*e: @{Edge->?}; ElemTableRec_t *top)==
  uint NT = 0;
    uint ct;
  { int NF = top->wall.nel;
    ??? t = Wall_vec_new(NF)^;
    {
      for (i = 0; i < NF; i++)
      {
        ct = 0;
        ??? f = top->wall[i];
          ??? fun = f.node[0]->num;
        ??? fvn = f.node[1]->num;
        ??? fwn = f.node[2]->num;
        ??? eun = e.node[0]->num;
        ??? evn = e.node[1]->num;
        {
          if (((fun == eun) || (fun == evn))){ ct++; }
          if (((fvn == eun) || (fvn == evn))){ ct++; }
          if (((fwn == eun) || (fwn == evn))){ ct++; }
           if (ct == 2)
           {
             t[NT] = f;
             NT++;
           }
         }
       }
       ??? r = Wall_vec_new(NT);
    {
         r^ = SUBARRAY(t,0,NT);
         return r;
       }
     }
   } CollectWalls;

/* UNUSED */
DRF ComputeWallRingDegree(*walls: REF Wall_vec_ts, ElemTableRec_t *top)
  { ??? drf = NEW(DRF, top->NE);
    {
      for (l = 0; l < top->NE; l++)
      {
        ??? fie = walls[l];
        {
          drf^[l] =  ((fie^).nel);
        }
      }
      return drf;
    }
  } /* END ComputeWallRingDegree */

/* UNUSED */
Wall_vec_ts *CollectWallsRing(ElemTableRec_t *top)==
  { ??? walls = NEW(REF Wall_vec_ts, top->NE);
    {
      for (l = 0; l < top->NE; l++)
      {
        ??? e = top->edge[l];
        ??? fie = CollectWalls(e,top);
        {
          walls[l] =  fie;
        }
      }
      return walls;
    }
  } /* END CollectWallsRing */

/* UNUSED */
double ComputeWallRingDegreeTotal(*drf: DRF)
      n : uint = 0;
  { for(l = 0; l < drf.nel; l++){
      ??? d = drf^[l];
      {
         n = n + d;
       }
    }
    return ((double)n);
  } /* END ComputeWallRingDegreeTotal */

PROCEDURE EquaAngle(Place_t @p, Coords_t *c, bool_vec_t *variable, 
    ElemDataRec_t *top;
 )
  { Node_t u = OrgV(a);
    ??? @{edge->?}s = Collect@{Edge->?}s(u,top);
    {
      for (j = 0 TO LAST(@{edge->?}s^))
      {
        ??? e =  @{edge->?}s^[j];
        ??? fe = CollectWalls(e,top);
        { 
          EqualizeDiedralAngle(e,fe,c,u);
        }
      }   
    }
  } /* END EquaAngle */


void EqualizeDiedralAngle(e : @{Edge->?}; *f : REF WALLS;
                               Coords_t *c; *u : Node)
  { ??? eu = e.node[0];
    ??? ev = e.node[1];
    ??? m == ((f^).nel); 
    ??? angle == 2.0 * ((double)Math.Pi)/((double)m);
    {
      for (i = 0; i < m; i++)
      {
        ???;
      }
    }
  } /* END EqualizeDiedralAngle */


/* UNUSED */
PROCEDURE CollectWingPoint(*ff: REF WALLS;
                           @{Edge->?} *e;
                          ) : REF uint_vec_t ==
  /* Retorna para cada aresta o conjunto de nodes dobradicas */
  CONST INIT_STACK_SIZE == 10000;
      stack = NEW(REF uint_vec_t ,INIT_STACK_SIZE);
      uint nstack = 0;

  bool_t Present(n : uint)
    uint nstack1 = nstack;
    VAR {
      while (nstack1 > 0)
      {
        nstack1 = nstack1 - 1;
        if (stack.el[nstack1] == n){ return TRUE; }
      }
      return FALSE;
    } Present;

  { ??? fi = ff^[0];
    ??? wingi = WingPoint(e, fi);
    {
      stack.el[nstack] = wingi; nstack++;
    }

    for(i = 1; i < ff.nel; i++)
    {
      ??? f = ff^[i];
      ??? wing = WingPoint(e, f);
      {
        if ((! Present(wing))){ stack.el[nstack] = wing; nstack++; }
      }
    }
    ??? r = NEW(REF uint_vec_t, nstack);
    {
      r^ = SUBARRAY(stack^, 0, nstack); 
      return r;
    }
  } /* END CollectWingPoint */


uint WingPoint(*e: @{Edge->?}; *f: Wall)
  uint *wp;
  { ??? eun = e.node[0]->num;
    ??? evn = e.node[1]->num;
    ??? fn = f.node;
    {
      for (i = 0; i < 3; i++)
      {
        if (((fn[i]->num!=eun) && (fn[i]->num!=evn)
        )){ wp = fn[i]->num;
        }
      }
    }
    return wp; 
  } /* END WingPoint */

  /*
void FlattenNode(
    Arc a; 
    Coords_t *c;
    bool_vec_t *variable; 
    <*UNUSED); rand_t *coins;
 ) ==

      somaki: double = 0.0;
      somakiti: double = 0.0;
      uint n = 0;

  void CalcTi(*bar: r3_t; e: Arc)
    /* 
      Computes the pseudo-force "ti" and its weight "ki"
      and accumulates them in "somatiki", "somaki"
    */
    {
      /* Check if hinge exists and has a spring: */
      ??? b = OPrev(Sym(e));
      ??? be = NARROW(b.edge, Triang.edge);
      ??? bl = LeftF(b);
      ??? br = LeftF(Sym(b));
      {
        if ((! (bl->exists) && (br->exists) && (be.spring))){ return;}
      }
      /* Do the computations: */
      ??? f = ONext(e);
        Node_t p = OrgV(Sym(e))->num;
      Node_t q = OrgV(Sym(f))->num;
      Node_t w = OrgV(Sym(ONext(ONext(Sym(f)))))->num;

      ??? di = r3_sub(c[q], c[p]);
      ??? ni = LR3Extras.Cross(di, r3_sub(bar, c[q]));
      ??? mi = LR3Extras.Cross(r3_sub(c[w], c[p]), di);

      ??? dimod = LR3.Norm(di);
      ??? nimod = LR3.Norm(ni);
      ??? mimod = LR3.Norm(mi);

      ??? sinTheta = LR3x3.Det(LR3x3.T{ni, mi, di}) / (nimod*mimod*dimod);
      ??? cosTheta = LR3.Cos(ni, mi);
      ??? theta = Math.atan2(sinTheta, cosTheta);
      ??? feta = ThetaFunc(theta);

      ??? ki = dimod;
      ??? ti = feta * nimod / dimod;
      {
        if (debug)
        {
          Debug.LongReal("CalcTi: ni == ", ni);
          Debug.LongReal("        mi == ", mi);
          Debug.LongReal("        di == ", di);
          Debug.LongReal("        sin, cos == ", LRN.T{sinTheta, cosTheta});
          Debug.LongReal("        theta, feta == ", LRN.T{theta, feta});
          Debug.LongReal("        ki, ti == ", LRN.T{ki, ti});
        }
        somaki = somaki + ki;
        somakiti = somakiti + ki * ti;
        N = N+1;
      }
    } CalcTi;

  double ThetaFunc(Theta: double)
    CONST PPi == ((double)Math.Pi);
          PiCube == PPi * PPi * PPi;
          B == 0.5;
          A == (1.0 - B*PPi) / PiCube;        
    {
      ??? ThetaCube = Theta * Theta * Theta;
      {
        return A * ThetaCube +  B * Theta;
      }
    } ThetaFunc;

  Arc *an;
  { Node_t vx = OrgV(a);
    ??? v = vx->num;
    ??? bar = Triang.NeighborBarycenter(a,c);
    ??? n = Triang.NodeNormal(a, c);
    {
      if ((! variable[v]) && (vx->exists)){ return; }
      if (debug){ Debug.LongReal("   bar ==  ", bar); }
      an = a;
      do 
      {
        CalcTi(bar, an); 
        an = ONext(an) 
      } while (an != a);
      ??? t = somakiti / somaki;
      Node_t ao = OrgV(a)->num;
      ??? u = r3_add(bar, r3_scale(t, n));
      {
        if (debug)
        {
          Debug.LongReal("   t ==    ", LRN.T{t}); 
          Debug.LongReal("   n ==    ", n);
          Debug.LongReal("   u ==    ", u);
        }
        c[ao] = u;
      }
    }
  } /* END FlattenNode */
*/

{
} Heuristics.