/* See dg_plot.h */
/* Last edited on 2009-08-21 15:31:44 by stolfilocal */

#include "dg_plot.h"
#include "mdg_grid.h"
#include "dg_locus.h"
#include "dg_tree.h"

#include "bz_basic.h"
#include "affirm.h"
#include "set32.h"
#include "jsmath.h"
#include "pswr.h"
#include "pswr_iso.h"
#include "pswr_color.h"
#include "vec.h"

#include <math.h>
#include <assert.h>
#include <values.h>

/* INTERNAL PROTOTYPES */

/* IMPLEMENTATIONS */

void dg_enum_locus_leaves
  ( mdg_dim_t d,            /* Dimension of grid. */
    dg_locus_t E,          /* Root locus. */
    mdg_rank_t r,           /* Rank of {E}. */
    bz_patch_t *b,         /* Shape of {E}, or NULL. */
    mdg_grid_size_t psz[],  /* Size vector of the tree pack {NE}. */
    dg_tree_vec_t *NE,     /* Tree star of {E}. */
    mdg_rank_t rMax,        /* Truncate the recursion just below this rank. */
    dg_locus_proc_t visit  /* Client face-processing procedure. */  
  )
  { 
    int nCells = NE->ne; /* Number of cells in star. */
    mdg_dim_t s = set32_count(E.norm);  /* Co-dimension of {E}. */
    mdg_dim_t m = d - s;                  /* Dimension of {E}. */
    bool_t leaf;
    
    /* Paranoia: */
    dg_tree_pack_check(NE, d, psz, mdg_cell_rank(E.cell));
    
    if (b != NULL) 
      { affirm(b->m == m, "inconsistent shape dimension"); }

    if (r <= rMax)
      { 
        fprintf(stderr, "+ dg_enum_locus_leaves(E = <%d:%lld> r = %d)\n", E.norm, E.cell, r);

        /* Check whether splitting-and-recursion is needed. */ 
        /* The face must be split if it is not flat, or if any adjacent cell
          has a split along an axis parallel to the face: */
        if (r >= rMax)
          { /* Pretend it is a leaf locus. */
            leaf = TRUE;
          }
        else
          { /* Check whether it is a leaf locus: */
            int i;
            /* dg_tree_pack_check(NE, d, psz, mdg_cell_rank(E.cell)); */
            leaf = TRUE; 
            /* Could be improved --- should skip levels with parallel splits. */
            for (i = 0; leaf && (i < nCells); i++)
              { dg_tree_node_t *NEi = NE->e[i];
                if (NEi != NULL)
                  { leaf &= (LOCH(*NEi) == NULL) || (HICH(*NEi) == NULL); }
              }
          }
        if (leaf) 
          { /* Leaf item, process it: */
            visit(d, E, r, b, NE);
          }
        else 
          { /* Recurse on children: */
            mdg_axis_t splax = mdg_split_axis(d, r);  /* Split axis. */
            if (set32_belongs(splax,E.norm))
              { /* The split axis is perpendicular to the item. */
                dg_locus_t E1 = dg_locus(E.norm,2*E.cell);  /* Same item in next level. */
                /* Shrink tree star around element: */
                dg_tree_vec_t NE1;
                assert(psz[splax] == 2);
                dg_tree_pack_child(d, psz, NE, splax, 1, &NE1);
                /* Recurse on same item in the next level, with same shape. */
                dg_enum_locus_leaves(d, E1, r+1, b, psz, &NE1, rMax, visit);
                free(NE1.e);
              }
            else
              { /* Axis {splax} is parallel to item {E}, split it: */
                /* Get the spanning axes of {E}: */
                mdg_axis_set_t SpnE = mdg_axis_set_complement(d, E.norm);
                /* Find relative position {splj} of {splax} among the {Spn(E)}: */
                mdg_axis_index_t splj = set32_index_from_elem(splax, SpnE);
                fprintf(stderr, "split on axis Spn(E)[%d] = %d\n", splj, splax);
                affirm((splj >= 0) &&(splj < d), "inconsistent axis sets");
                assert(psz[splax] == 1);

                /* Compute normal axes of splitting item: */
                mdg_axis_set_t NrmMD = set32_include(splax, E.norm); 
              
                /* Pieces of this locus: */
                dg_locus_t ELO = dg_locus(E.norm, 2*E.cell);
                dg_locus_t EHI = dg_locus(E.norm, 2*E.cell + 1);
                dg_locus_t EMD = dg_locus(NrmMD, 2*E.cell + 1);
                
                /* Shapes of the pieces: */
                bz_patch_t bMD = bz_patch_facet_new(m, b->n, b->g, splj);
                bz_patch_t bLO = bz_patch_new(m, b->n, b->g);
                bz_patch_t bHI = bz_patch_new(m, b->n, b->g);

                 /* Tree stars of the pieces: */
                dg_tree_vec_t NLO, NHI, NMD;
                
                /* Split the tree star: */
                assert(psz[splax] == 1);
                dg_tree_pack_child(d, psz, NE, splax, 0, &NLO);
                dg_tree_pack_child(d, psz, NE, splax, 1, &NHI);

                /* Split the item's shape: */
                bz_patch_split(b, splj, 0.5, &bLO, &bHI);
                bz_patch_get_facet(&bHI, splj, 0, &bMD);
                
                /* Recurse for each half of item {F}: */
                dg_enum_locus_leaves(d, ELO, r+1, &bLO, psz, &NLO, rMax, visit);
                dg_enum_locus_leaves(d, EHI, r+1, &bHI, psz, &NHI, rMax, visit);

                /* Recurse for the separating item {S} of {F}. */
                psz[splax] = 2;
                dg_tree_pack_split(d, psz, NE, splax, &NMD);
                dg_enum_locus_leaves(d, EMD, r+1, &bMD, psz, &NMD, rMax, visit);
                psz[splax] = 1;

                /* Reclaim temporary storage. */
                bz_patch_free(bMD); free(NMD.e);
                bz_patch_free(bLO); free(NLO.e);
                bz_patch_free(bHI); free(NHI.e);
              }
          }

        fprintf(stderr, "- dg_enum_locus_leaves(E = <%d:%lld> r = %d)\n", E.norm, E.cell, r);
      }
  }

/* POSTSCRIPT PLOTTING OF 2D MULTIGRID ITEMS */

void dg_plot_2D_vertex(PSStream *ps, double *c, bz_patch_rdim_t n)
  {
    pswr_set_pen(ps, 0.0, 0.0, 0.0,   0.15, 0.0, 0.0);
    pswr_set_fill_color(ps, 0.0,0.0,0.0);
    pswr_dot(ps, c[0], c[1], 0.25,  TRUE, FALSE);
  }

void dg_plot_2D_edge(PSStream *ps, double c0[], double c1[], bz_patch_rdim_t n)
  {
    pswr_set_pen(ps, 0.0, 0.0, 0.0,   0.15, 0.0, 0.0);
    pswr_segment(ps, c0[0], c0[1], c1[0], c1[1]);
  }

#define dg_bilinear(x0,x1,c00,c01,c10,c11) \
  (1-(x0))*((1-(x1))*(c00) + (x1)*(c01)) + (x0)*((1-(x1))*(c10) + (x1)*(c11))
  /* Bilinear interpolation between values {c00,c01,c10,c11} with 
     parameters {x0,x1}, in [0 _ 1]. */

void dg_plot_2D_face_shade
  ( PSStream *ps,
    double c00[], 
    double c01[],
    double c10[], 
    double c11[],
    bz_patch_rdim_t n,
    int plotDepth,
    double fMin,
    double fMax
  )
  {
    int i0, i1, j;
    int steps = ipow(2, plotDepth);
    double p[bz_patch_MAX_RDIM], q[bz_patch_MAX_RDIM], r[bz_patch_MAX_RDIM], s[bz_patch_MAX_RDIM];
    for (i0 = 1; i0 <= steps; i0++)
      { for (i1 = 0; i1 <= steps; i1++)
          { double u0 = ((double)i0)/((double)steps);
            double u1 = ((double)i1)/((double)steps);
            double v0 = ((double)i0-1)/((double)steps);
            double v1 = ((double)i1)/((double)steps);
            for (j = 0; j < n; j++)
              { p[j] = dg_bilinear(u0,u1,c00[j],c01[j],c10[j],c11[j]);
                r[j] = dg_bilinear(v0,v1,c00[j],c01[j],c10[j],c11[j]);

              }
            if (i1 > 0) { dg_shade_quadrilateral(ps, p, q, r, s, n, fMin, fMax); }
            for (j = 0; j < n; j++) { q[j] = p[j]; s[j] = r[j]; }
          }
      }
  }

void dg_plot_2D_face_bands
  ( PSStream *ps,
    double c00[], 
    double c01[],
    double c10[], 
    double c11[],
    bz_patch_rdim_t n,
    int plotDepth,
    double fStart,  /* Synchronize isolines with this level. */
    double fStep,   /* isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax,       /* Maximum isoline index. */
    double *R, double *G, double *B
  )
  {
    demand(n == 3, "can't plot bands with more than one function value"); 
    int i0, i1, j;
    int steps = ipow(2, plotDepth);
    double p[bz_patch_MAX_RDIM], q[bz_patch_MAX_RDIM], r[bz_patch_MAX_RDIM], s[bz_patch_MAX_RDIM];
    for (i0 = 1; i0 <= steps; i0++)
      { for (i1 = 0; i1 <= steps; i1++)
          { double u0 = ((double)i0)/((double)steps);
            double u1 = ((double)i1)/((double)steps);
            double v0 = ((double)i0-1)/((double)steps);
            double v1 = ((double)i1)/((double)steps);
            for (j = 0; j < n; j++)
              { p[j] = dg_bilinear(u0,u1,c00[j],c01[j],c10[j],c11[j]);
                r[j] = dg_bilinear(v0,v1,c00[j],c01[j],c10[j],c11[j]);

              }
            if (i1 > 0) 
              { dg_bands_in_quadrilateral
                  (ps, p, q, r, s, n, fStart, fStep, kMin, kMax, R, G, B);
              }
            for (j = 0; j < n; j++) { q[j] = p[j]; s[j] = r[j]; }
          }
      }
  }

void dg_plot_2D_face_isolines
  ( PSStream *ps,
    double c00[], 
    double c01[],
    double c10[], 
    double c11[],
    bz_patch_rdim_t n,
    int plotDepth,
    double fStart,  /* Synchronize isolines with this level. */
    double fStep,   /* Isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax        /* Maximum isoline index. */
  )
  {
    demand(n == 3, "can't plot isolines with more than one function value"); 
    int i0, i1, j;
    int steps = ipow(2, plotDepth);
    double p[bz_patch_MAX_RDIM], q[bz_patch_MAX_RDIM], r[bz_patch_MAX_RDIM], s[bz_patch_MAX_RDIM];
    for (i0 = 1; i0 <= steps; i0++)
      { for (i1 = 0; i1 <= steps; i1++)
          { double u0 = ((double)i0)/((double)steps);
            double u1 = ((double)i1)/((double)steps);
            double v0 = ((double)i0-1)/((double)steps);
            double v1 = ((double)i1)/((double)steps);
            for (j = 0; j < n; j++)
              { p[j] = dg_bilinear(u0,u1,c00[j],c01[j],c10[j],c11[j]);
                r[j] = dg_bilinear(v0,v1,c00[j],c01[j],c10[j],c11[j]);

              }
            if (i1 > 0) 
              { dg_isolines_in_quadrilateral(ps, p, q, r, s, n, fStart, fStep, kMin, kMax); }
            for (j = 0; j < n; j++) { q[j] = p[j]; s[j] = r[j]; }
          }
      }
  }

void dg_shade_quadrilateral
  ( PSStream *ps,
    double p00[],   /* Coords and values at corner [0,0]. */
    double p01[],   /* Coords and values at corner [0,1]. */
    double p10[],   /* Coords and values at corner [1,0]. */
    double p11[],   /* Coords and values at corner [1,1]. */
    bz_patch_rdim_t n,      /* Dimension of range space. */
    double fMin,    /* Minimum function value. */
    double fMax     /* Maximum function value. */
  )
  {
    double ctr[bz_patch_MAX_RDIM];
    int j;
    for (j = 0; j < n; j++) 
      { ctr[j] = (p00[j]+p01[j]+p10[j]+p11[j])/4.0; }
    /* Paint triangles: */
    dg_shade_triangle(ps, p00, p01, ctr, n, fMin, fMax);
    dg_shade_triangle(ps, p01, p11, ctr, n, fMin, fMax);
    dg_shade_triangle(ps, p11, p10, ctr, n, fMin, fMax);
    dg_shade_triangle(ps, p10, p00, ctr, n, fMin, fMax);
  }

void dg_bands_in_quadrilateral
  ( PSStream *ps,
    double p00[],   /* Coords and values at corner [0,0]. */
    double p01[],   /* Coords and values at corner [0,1]. */
    double p10[],   /* Coords and values at corner [1,0]. */
    double p11[],   /* Coords and values at corner [1,1]. */
    bz_patch_rdim_t n,      /* Dim of point vectors (incl. function values). */
    double fStart,  /* Synchronize isolines with this level. */
    double fStep,   /* Isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax,       /* Maximum isoline index. */
    double *R, double *G, double *B
  )
  {
    demand(n == 3, "can't plot bands with more than one function value");
    pswr_bands_in_quadrilateral
      ( ps,
        p00[0], p00[1], p00[2],
        p01[0], p01[1], p01[2],
        p10[0], p10[1], p10[2],
        p11[0], p11[1], p11[2],
        fStart, fStep, kMin, kMax,
        R, G, B
      );
  }

void dg_isolines_in_quadrilateral
  ( PSStream *ps,
    double p00[],   /* Coords and values at corner [0,0]. */
    double p01[],   /* Coords and values at corner [0,1]. */
    double p10[],   /* Coords and values at corner [1,0]. */
    double p11[],   /* Coords and values at corner [1,1]. */
    bz_patch_rdim_t n,      /* Dimension of range space. */
    double fStart,   /* Synchronize isolines with this level. */
    double fStep,   /* Isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax        /* Maximum isoline index. */
  )
  {
    demand(n == 3, "can't plot isolines with more than one function value");
    pswr_isolines_in_quadrilateral
      ( ps,
        p00[0], p00[1], p00[2],
        p01[0], p01[1], p01[2],
        p10[0], p10[1], p10[2],
        p11[0], p11[1], p11[2],
        fStart, fStep, kMin, kMax
      );
  }

void dg_shade_triangle
  ( PSStream *ps,
    double a[],     /* First vertex. */
    double b[],     /* Second vertex. */
    double c[],     /* Third vertex. */
    bz_patch_rdim_t n,      /* Dimension of range space. */
    double fMin,    /* Minimum function value. */
    double fMax     /* Maximum function value. */
  )
  { double R, G, B;
    double f[bz_patch_MAX_RDIM];
    int j;
    affirm (n >= 2, "invalid point coordinates");
    /* Compute mean function values normalized to [-1 _ +1]: */
    for (j = 0; j < n-2; j++)
      { /* Compute mean function value: */
        f[j] = (a[j+2]+b[j+2]+c[j+2])/3.0;
        /* Clip to {[fMin _ fMax]}: */
        if (f[j] < fMin) { f[j] = fMin; }
        if (f[j] > fMax) { f[j] = fMax; }
        /* Normalize to [-1 _ +1]: */
        if (f[j] > 0.0) 
          { f[j] /= fMax; }
        else if (f[j] < 0.0) 
          { f[j] /= -fMin; }
      }
    /* Compute color {R,G,B} for white-background plotting: */
    double Rs = 1.0, Gs = 0.342, Bs = 0.000;
    double Y0 = 0.95;
    double Y1 = 0.30;
    if (n == 2)
      { R = 1.0; G = 0.9; B = 0.8; }
    else if (n == 3)
      { pswr_color_scale_1(f[0], Rs, Gs, Bs, Y0, &R, &G, &B); }
    else if (n == 4)
      { pswr_color_scale_2(f[0], f[1], Rs, Gs, Bs, Y0, &R, &G, &B); }
    else if (n >= 5)
      { pswr_color_scale_3(f[0], f[1], f[2], Y0, Y1, &R, &G, &B); }
    pswr_set_fill_color(ps, R,G,B);
    pswr_triangle(ps, a[0], a[1], b[0], b[1], c[0], c[1], TRUE, FALSE);
  }

PSStream *dg_plot_init
  ( char *name,
    bool_t epsformat,
    char *paperSize,
    char *caption,
    interval_t bbox[]
  )
  { double hpt, vpt;
    if (epsformat)
      { hpt = 5.0*72.0; vpt = hpt * 4.0/3.0; }
    else
      { pswr_get_paper_dimensions(paperSize, &hpt, &vpt); }
    PSStream *ps = pswr_new_stream(name, NULL, epsformat, "doc", paperSize, hpt+6, vpt+8);
    pswr_set_canvas_layout(ps, 3.0, 4.0, TRUE, 0.0, 0.0, (caption == NULL ? 0 : 1), 1,1);
    double dx = 0.1*(HI(bbox[0]) - LO(bbox[0]))/2;
    double dy = 0.1*(HI(bbox[1]) - LO(bbox[1]))/2;
    pswr_new_picture
      ( ps, 
        LO(bbox[0])-dx, HI(bbox[0])+dx, 
        LO(bbox[1])-dy, HI(bbox[1])+dy
      );
    pswr_set_pen(ps, 0,0,0, 0.15, 0,0);
    if (caption != NULL)
      { pswr_add_caption(ps, caption, 0.0); }
    pswr_set_pen(ps, 1,0,0, 0.25, 0,0);
    return ps;
  }

void dg_plot_finish ( PSStream *ps )
  { pswr_close_stream(ps); }