/* See figstree.h. */
/* Last edited on 2009-08-21 14:32:04 by stolfilocal */
/* Last edited on 2009-08-14 21:55:17 by douglas */

  
#include "figstree.h"

void fig_tree
  ( PlotOptions *o, 
    char *tag,
    dg_tree_node_t *root
  )
  {
    int nlevels = 0, nnodes = 0;

    /* Count nodes and levels */
    auto void count_nodes(dg_tree_node_t *p, mdg_rank_t r);
    void count_nodes(dg_tree_node_t *p, mdg_rank_t r)
      { if (p != NULL)
        { nnodes++; 
          if (r >= nlevels) { nlevels = r+1; }
          count_nodes(LOCH(*p), r+1);
          count_nodes(HICH(*p), r+1);
        }
      }

    fprintf(stderr, "counting nodes...\n");
    count_nodes(root, 0);
    fprintf(stderr, "%d nodes found in %d levels.\n", nnodes, nlevels);
    
    /* Nodes per level, their positions, children, and ancestors. */
    double xc[nnodes]; 
    dg_tree_node_t *nd[nnodes];
    mdg_cell_index_t ix[nnodes];
    int father[nnodes], child[2*nnodes];
    mdg_rank_t rk[nnodes];

    int beg[nlevels+1];

    /* Gather nodes per level: */
    int i, j, k;

    fprintf(stderr, "gathering node data...\n");
    nd[0] = root; ix[0] = 1; rk[0] = 0;
    father[0] = -1; child[0] = -1; child[1] = -1;
    beg[0] = 0;
    j = 1;
    for (k = 1; k <= nlevels; k++) { beg[k] = -1; }
    i = 0;
    while (i < j)
      { mdg_rank_t r = rk[i];
        dg_tree_node_t *p = nd[i];
        int ich;
        for (ich = 0; ich < 2; ich++)
          { dg_tree_node_t *ch = p->ch[ich];
            if (ch != NULL)
              { nd[j] = ch;
                ix[j] = ch->index;
                rk[j] = r+1;
                father[j] = i;
                child[2*j] = -1; child[2*j+1] = -1;
                child[2*i+ich] = j;
                if (beg[r+1] == -1) { beg[r+1] = j; }
                fprintf(stderr, " %d(%d)", (int)(ix[j]), rk[j]);
                j++;
              }
            else
              { child[2*i+ich] = -1; }
          }
        i++;
      }
    affirm(j == nnodes, "bug nnodes");
    beg[nlevels] = nnodes;
    fprintf(stderr, "finished gathering nodes.\n");
    
    /* Compute the X coordinates */
    fprintf(stderr, "computing X coordinates...\n");
    for (i = 0; i < nnodes; i++) { xc[i] = NAN; }
    
    auto void relax_level(mdg_rank_t k);
    void relax_level(mdg_rank_t k)
      { 
        fprintf(stderr, "  recomputing level %d ...\n    ", k);
        /* Recompute {xc[i]} from neighbors */
        double w = 5.0; /* Weight of child rel to parent. */
        for(i = beg[k]; i < beg[k+1]; i++)
          { double xtot = 0.0, xq, wtot = 0.0;
            xq = (father[i] != -1 ? xc[father[i]] : NAN);
            if (! isnan(xq)) { xtot += xq + (double)(ix[i] % 2) - 0.5; wtot += 1.0; }
            xq = (child[2*i] != -1 ? xc[child[2*i]] : NAN);
            if (! isnan(xq)) { xtot += w*xq; wtot += w; }
            xq = (child[2*i+1] != -1 ? xc[child[2*i+1]] : NAN);
            if (! isnan(xq)) { xtot += w*xq; wtot += w; }
            xc[i] = xtot / (wtot == 0.0 ?  1.0 : wtot);
            fprintf(stderr, " %d:%.2f", i, xc[i]);
          }
        fprintf(stderr, "\n");
        
        fprintf(stderr, "  fixing overlaps ...\n    ");
        /* Ensure minimum spacing of 1 between successive nodes: */
        for (i = beg[k]+1; i < beg[k+1]; i++)
          { double overlap = (xc[i-1] + 0.5) - (xc[i] - 0.5);
            if (overlap > 0.0)
              { /* Find start of solid block ending with node {i-1}: */
                int r = i-1;
                double shift = overlap/2.0, totgap = 0.0;
                while (r > beg[k])
                  { double gap = (xc[r] - 0.5) - (xc[r-1] + 0.5);
                    if ( gap >= shift) { break; }
                    r--; 
                    totgap += gap;
                    shift = (overlap - totgap)/((double)i-r+1);
                  }
                fprintf(stderr, " %d..%d:%+.2f", r,i,shift);
                xc[r] -= shift; 
                while (r < i) { r++; xc[r] = xc[r-1] + 1.0; }
              }
          }
        fprintf(stderr, "\n");
        fprintf(stderr, "  range = [%6.2f _ %6.2f].\n", xc[beg[k]], xc[beg[k+1]-1]);
     }

    int npasses = 60;
    for (j = 0; j < npasses; j++)
      { for (k = 0; k < nlevels; k++) { relax_level(k); }
        for (k = nlevels-1; k >= 0; k--) { relax_level(k); }
      }
    fprintf(stderr, "finished computing the coordinates.\n");
      
    fprintf(stderr, "computing scales...\n");
    /* Plot scale: */
    double dy = 2.0; /* Relative spacing between levels */
    
    /* Now draw tree: */
    interval_t tr[2];   /* Bounding box of tree, slightly widened */
    interval_t win[2];  /* Plot window (in client coords). */
    
    tr[0] = (interval_t){{ INF, -INF }};
    for (i = 0; i < nnodes; i++)
      { interval_t xbi = (interval_t){{ xc[i]-0.5, xc[i]+0.5 }};
        tr[0] = interval_join(&(tr[0]), &xbi);
      }
    tr[1] = (interval_t){{ -dy*(nlevels-1)-0.5, +0.5 }};
    fprintf(stderr, 
      "tree bounding box = [%6.2f _ %6.2f] � [%6.2f _ %6.2f]\n", 
      LO(tr[0]), HI(tr[0]), LO(tr[1]), HI(tr[1])
    );
    
    double xsz = HI(tr[0]) - LO(tr[0]);
    double xscale = 8.0*25.4/xsz;
    double ysz = HI(tr[1]) - LO(tr[1]);
    double yscale = 6.0*25.4/ysz;
    double scale = (xscale < yscale ? xscale : yscale);
    int ax;
    for (ax = 0; ax < 2; ax++)
      { /* Widen box slightly: */
        LO(tr[ax]) += -0.05;
        HI(tr[ax]) += +0.05;

        /* Plotting window (client units): */
        LO(win[ax]) = LO(tr[ax]) - 0.05;
        HI(win[ax]) = HI(tr[ax]) + 0.05;
      }
    fprintf(stderr, "done computing scales...\n");

    fprintf(stderr, "starting figure...\n");
    start_figure(o, tag, tag, win, scale);
    
    /* Plot the edges: */
    fprintf(stderr, "plotting edges...\n");
    pswr_set_pen(ps, 0.0,0.0,0.0, 0.10, 0.0, 0.0);
    
    for (i = 0; i < nnodes; i++)
      { int j = father[i];
        if (j != -1)
          { double xi = xc[i], yi = -dy*rk[i];
            double xj = xc[j], yj = -dy*rk[j];
            pswr_segment(ps, xi, yi, xj, yj); 
          }
      }
    
    /* Plot the nodes: */
    fprintf(stderr, "plotting nodes...\n");
    pswr_set_pen(ps, 0.0,0.0,0.0, 0.20, 0.0, 0.0);
    pswr_set_fill_color(ps, 1.0,1.0,1.0);
    for (i = 0; i < nnodes; i++)
      { double xi = xc[i], yi = -dy*rk[i];
        pswr_dot(ps, xi, yi, 1.0, TRUE, TRUE);
      }

    /* Close figure: */
    finish_figure(o);
    fprintf(stderr, "done.\n");
  }