/* See plot1d.h. */
/* Last edited on 2021-08-20 16:34:53 by stolfi */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <values.h>

#include <affirm.h>
#include <pswr.h>

#include <dg_grid.h>
#include <dg_locus.h>
#include <dg_pulse.h>
#include <dg_tent.h>

#include <plot1d.h>
  
#define DIM 1
#define MAXTAGSZ 200 
  
#define DRAW_COORD_LINES FALSE
#define DRAW_FRAME TRUE

void compute_1d_pulse_ranges
  ( dg_pulse_kind_t pkind, /* Kind of mother spline. */
    dg_cont_t c,           /* Continuity class of spline. */
    dg_degree_t g,         /* Degree of spline. */
    dg_grid_size_t gsz,     /* Period of spline (#intervals). */
    double wd,             /* Width of root cell. */
    interval_t fr[]        /* (OUT) Nominal graph range (integers). */
  )
  { 
    /* Initialize intervals to empty: */
    fr[0] = (interval_t){{ +INF, -INF }};
    fr[1] = (interval_t){{ +INF, -INF }};
    
    /* Get the maximum pulse index for this pulse family: */
    dg_pulse_family_t fam = dg_pulse_family(pkind, c, g);
    int maxpix = fam.nmp - 1;
    
    /* Unite ranges for all pulse indices: */
    dg_pulse_mother_index_t pix;
    for (pix = 0; pix <= maxpix; pix++)
      { /* Get the support count: */
        dg_grid_size_t m = dg_pulse_supp_count(&fam, pix, gsz);
        /* Compute support for zero shift: */
        double xMin = 0, xMax = m*wd/gsz;
        /* Unite with {fr[0]}: */
        if (xMin < LO(fr[0])) { LO(fr[0]) = xMin; }
        if (xMax > HI(fr[0])) { HI(fr[0]) = xMax; }
        /* Estimate range of spline values from Bézier coeffs: */
        fr[1] = get_pulse_range(pkind, c, g, pix, gsz);
      }
    /* Round ranges to integers: */
    if (LO(fr[1]) > 0.0)
      { LO(fr[1]) = 0.0; }
    else
      { LO(fr[1]) = floor(LO(fr[1])); }
    if (HI(fr[1]) < 0.0)
      { HI(fr[1]) = 0.0; }
    else
      { HI(fr[1]) = ceil(HI(fr[1])); }
 } 


void compute_1d_common_pulse_ranges
  ( dg_cont_t c,          /* Continuity class of spline. */
    dg_grid_size_t gsz,    /* Period of spline (#intervals). */
    double wd,            /* Width of root cell. */
    interval_t fr[]       /* (OUT) Nominal graph range. */
  )
  { fr[0] = (interval_t) {{ +INF, -INF }};
    fr[1] = (interval_t) {{ +INF, -INF }};
    dg_pulse_kind_t pkind;
    for (pkind = 0; pkind <= dg_pulse_kind_MAX; pkind++)
      { dg_degree_t gmax = 2*c + 1; /* Should be enough? */
        dg_degree_t g;
        for (g = 0; g <= gmax; g++)
          { interval_t frp[2];  /* Ranges for pulses {pkind,g,c}. */
            compute_1d_pulse_ranges(pkind, c, g, gsz, wd, /* OUT */ frp);
            box_join(2, fr, frp, fr);
          }
      }
  } 

void compute_1d_plot_scales
  ( PSStream *ps,
    PlotOptions *o,
    interval_t fr[],   /* (IN) Nominal function range. */
    double scale,       /* Plot scale (one client unit in mm). */
    double  mag[],      /* Graph scaling factors (function to client). */
    interval_t gr[],   /* (OUT) Range of function coords for plot. */
    interval_t win[]   /* (OUT) Plot window in client coords. */
  )
  {
    /* Extra function coords margin for graph draw/clip: */
    double fMrg[2]; 
    fMrg[0] = 0.05; 
    fMrg[1] = 0.05;
    
    /* Extra margin for plotting window (client units): */
    double wMrg[2];
    wMrg[0] = 2.0/scale/mag[0]; 
    wMrg[1] = 1.5/scale/mag[1];
    
    int i;
    for (i = 0; i < 2; i++)
      { /* Function coords range for evaluation and graph clipping: */
        LO(gr[i]) = LO(fr[i])-fMrg[i];
        HI(gr[i]) = HI(fr[i])+fMrg[i];

        /* Plotting window (client units): */
        LO(win[i]) = mag[i]*(LO(gr[i])-wMrg[i]);
        HI(win[i]) = mag[i]*(HI(gr[i])+wMrg[i]);
      }

    /* Ensure bottom space for tick labels: */
    double wBot = 5.0/scale/mag[1] - wMrg[1] + LO(fr[1]);
    if (wBot > 0.0) { LO(win[1]) -= mag[1]*wBot; }
  }

void plot_1d_graph_decoration
  ( PSStream *ps,
    PlotOptions *o,
    double scale,
    double mag[],
    interval_t fr[],  /* Function domain and range, for ticks &c. */
    double tstep[],    /* Space between major (labeled) ticks. */
    double sstep[],    /* Space between minor (unlabeled) ticks. */
    interval_t gr[],  /* Function domain and range for graph/clip. */
    interval_t win[]  /* Plotting area (client units). */
  )
  { double ht = (o->eps ? 0.5 : 1.0)/scale; /* Half-length of major ticks. */
    double hs = 0.5*ht; /* Half-length of minor ticks. */
    
    /* Draw coordlines and tick marks: */
    /* Minor ticks */
    plot_1d_ticks_and_lines
      ( ps, o, 0, hs, sstep[0], FALSE, FALSE, &(fr[0]), mag[0], &(gr[1]), mag[1], win );
    plot_1d_ticks_and_lines
      ( ps, o, 1, hs, sstep[1], FALSE, FALSE, &(fr[1]), mag[1], &(gr[0]), mag[0], win );
    /* Major ticks */
    plot_1d_ticks_and_lines
      ( ps, o, 0, ht, tstep[0], TRUE,  FALSE, &(fr[0]), mag[0], &(gr[1]), mag[1], win );
    plot_1d_ticks_and_lines
      ( ps, o, 1, ht, tstep[1], TRUE,  FALSE, &(fr[1]), mag[1], &(gr[0]), mag[0], win );
    
    /* Draw coordinate axes: */
    pswr_set_pen(ps, 0.0,0.0,0.0, 0.10, 0.0, 0.0);
    plot_1d_axes(ps, o, &(gr[0]), mag[0], &(gr[1]), mag[1]);
    
    if (DRAW_FRAME)
      { /* Draw frame around plot: */
        pswr_set_pen(ps, 0.0,0.0,0.0, 0.15, 0.0, 0.0);
        pswr_frame(ps);
      }
  }
  
void plot_1d_tree_leaves
  ( PSStream *ps,
    PlotOptions *o, 
    dg_tree_node_t *G,
    interval_t rootCell[],
    interval_t xr,
    double xmag,
    interval_t yr, /* Range along Y axis. */
    double ymag     /* Extra scale factor for Y axis. */
  )
  {
    double ylo = ymag*LO(yr);
    double yhi = ymag*HI(yr);
    
    auto void do_plot(dg_rank_t r, dg_tree_node_t *n, interval_t cell[]);
    
    void do_plot(dg_rank_t r, dg_tree_node_t *n, interval_t cell[])
      { if (n == NULL)
          { /* Do nothing */ }
        else if ((LOCH(*n) == NULL) || (HICH(*n) == NULL))
          { /* Leaf cell, plot: */
            affirm((LOCH(*n) == NULL) && (HICH(*n) == NULL), "ch");
            int i;
            for (i = 0; i < 2; i++)
              { double x = cell[0].end[i];
                if ((x >= LO(xr)) && (x <= HI(xr)))
                  { pswr_segment(ps, xmag*x, ylo, xmag*x, yhi); }
              }
          }
        else
          { /* Non-leaf cell, recurse: */
            interval_t save = cell[r % DIM];
            affirm((LOCH(*n) != NULL) && (HICH(*n) != NULL), "ch");
            cell[r % DIM] = interval_split(save, 0);
            do_plot(r+1, LOCH(*n), cell);
            cell[r % DIM] = interval_split(save, 1);
            do_plot(r+1, HICH(*n), cell);
            cell[r % DIM] = save;
          }
      }

    pswr_set_fill_color(ps, o->fill.c[0], o->fill.c[1], o->fill.c[2]);
    do_plot(0, G, rootCell);
  }

void plot_1d_function
  ( PSStream *ps,
    PlotOptions *o,
    double func (double x),
    double xstep,
    interval_t xr,
    double xmag,
    interval_t yr,
    double ymag
  )
  { 
    double x1 = NAN; /* Previous sample point, or {NAN} if none */
    double y1;       /* Function value at {x1}. */
        
    auto void plot_sample(double x);
    
    void plot_sample(double x)
      { double x0 = x;
        double y0 = func(x0);
        if 
          ( (! isnan(x1)) &&
            ((y0 >= LO(yr)) || (y0 <= HI(yr))) &&
            ((y1 >= LO(yr)) || (y1 <= HI(yr)))
          )
          { pswr_segment(ps, xmag*x1, ymag*y1, xmag*x0, ymag*y0); }
        x1 = x0; y1 = y0;
      }

    int N = (int)(1.0/xstep);
    int Nmin = (int)ceil(LO(xr)*N);
    int Nmax = (int)floor(HI(xr)*N);
    double h = 0.00001*(HI(xr) - LO(xr));
    int i;
    fprintf(stderr, "xr = [%7.4f _ %7.4f] xstep = %7.4f i = {%d..%d}\n", 
            LO(xr), HI(xr), xstep, Nmin, Nmax);
    for (i = Nmin; i <= Nmax; i++)
      { double x = ((double) i)*xstep;
        if ((i % N) == 0)
          { plot_sample(x-h);
            plot_sample(x+h);
          }
        else
          { plot_sample(x); }
      }
    fflush(ps->file);
  }

void plot_1d_locus
  ( PSStream *ps,
    PlotOptions *o, 
    dg_locus_t E,
    dg_tree_node_t *G,
    interval_t rootCell[],
    interval_t xr,
    double xmag
  )
  {
    dg_dim_t dE = dg_locus_dimension(DIM, E);
    interval_t box[DIM]; /* Box of {E}. */
    dg_locus_box_root_relative(DIM, E, box);
    box_box_map(DIM, box, rootCell, box);
    switch (dE)
      { 
        case 0:
          { double x = LO(box[0]);
            if ((x >= LO(xr)) && (x <= HI(xr)))
              { pswr_set_fill_color(ps,  0.0, 0.0, 0.0);
                pswr_dot
                  ( ps, 
                    xmag*x, 0.0,
                    0.5,
                    TRUE, FALSE      
                  );
              }
          }
          break;
        case 1:
          { double x0 = LO(box[0]);
            double x1 = HI(box[0]);
            if ((x1 >= LO(xr)) && (x0 <= HI(xr)))
              { pswr_set_pen(ps,  0.0,0.0,0.0, 0.5,  0.0, 0.0);
                pswr_segment
                  ( ps, 
                    xmag*LO(box[0]), 0.0, 
                    xmag*HI(box[0]), 0.0
                  );
              }
          }
          break;
        default:
          affirm(FALSE, "bad locus dimension");
      }
  }

void plot_1d_axes
  ( PSStream *ps,
    PlotOptions *o,
    interval_t *xr,
    double xmag,
    interval_t *yr,
    double ymag
  )
  { 
    if ((LO(*xr) < 0.0) && (HI(*xr) > 0.0))
      { double tlo = ymag*LO(*yr);
        double thi = ymag*HI(*yr);
        pswr_segment(ps, 0.0, tlo, 0.0, thi);
      }
    if ((LO(*yr) < 0.0) && (HI(*yr) > 0.0))
      { double tlo = xmag*LO(*xr);
        double thi = xmag*HI(*xr);
        pswr_segment(ps, tlo, 0.0, thi, 0.0);
      }
  }

void plot_1d_coord_lines
  ( PSStream *ps,
    PlotOptions *o,
    dg_axis_t ax,      /* Axis PERPENDICULAR to the lines. */
    int nsteps,      /* Number of steps (num of lines minus 1) */
    interval_t *nr, /* Cordinate range along axis {ax}  */
    double nmag,     /* Extra scale factor along axis {ax}. */
    interval_t *tr, /* Range along the lines. */
    double tmag      /* Extra scale factor along the lines. */
  )
  { /* Plot coord lines: */
    int i;
    double tlo = tmag*LO(*tr);
    double thi = tmag*HI(*tr);
            
    for (i = 0; i <= nsteps; i++)
      { double t = (nsteps == 0 ? 0.5 : ((double)i)/nsteps);
        double s = (nsteps == 0 ? 0.5 : ((double)(nsteps-i))/nsteps);
        double c = s*LO(*nr) + t*HI(*nr);
        double cp = nmag*c;  /* Client coordinate for plotting */
        if (ax == 0)
          { pswr_segment(ps, cp, tlo, cp, thi); }
        else
          { pswr_segment(ps, tlo, cp, thi, cp); }
      }
  }

void plot_1d_ticks_and_lines
  ( PSStream *ps,
    PlotOptions *o,
    dg_axis_t ax,       /* Axis along which ticks are to be placed. */
    double hsize,     /* Half-length of tick marks (client units). */
    double step,      /* Spacing between ticks (func units). */
    bool_t labeled,     /* TRUE to print label at left/under each tick. */
    bool_t lines,       /* TRUE draws coord lines at tick positions. */
    interval_t *cr,  /* Only place ticks within this coordinate range. */
    double cmag,      /* Extra scale factor along axis {ax}. */
    interval_t *gr,  /* Trim coord lines to this coordinate range. */
    double gmag,      /* Extra scale factor across axis {ax}. */
    interval_t win[] /* Plot window in client units. */
  )
  { double wd = HI(*cr) - LO(*cr);
    int nt;     /* Number of tick intervals. */
    nt = (floor)(wd/step + 0.5);
    char *fmt;  /* Format for major tick labels in each axis. */
    if (labeled) 
      { /* Use sign if negative; use fraction if small step. */
        fmt = plot_1d_tick_format(wd/((double)nt), cr);
      }
    else
      { fmt = NULL; }

    if (lines)
      { /* Draw coord lines: */
        pswr_set_pen(ps, 0.0,0.0,0.0, 0.075, 0.50, 0.50);
        plot_1d_coord_lines(ps, o, ax, nt, cr, cmag, gr, gmag);
      }

    pswr_set_pen(ps, 0.0,0.0,0.0, 0.10, 0.0, 0.0);
    pswr_set_label_font(ps, "TimesRoman", (o->eps ? 8.0 : 12.0));
    plot_1d_ticks(ps, o, ax, hsize, fmt, nt, cr, cmag, win);
  }    

char *plot_1d_tick_format(double step, interval_t *fr)
  { /* If the range extends into negative numbers, print the sign too: */
    bool_t wsign = (LO(*fr) < -0.9999*step);
    /* Compute the number {nd} of decimals needed: */
    int nd = 0;
    while ((nd <= 2) && (fabs(floor(step + 0.5) - step) > 1.0e-5 * step))
      { step *= 10.0; nd++; }
    /* If the range extends into negative numbers, print the sign too: */
    if (wsign) 
      { return (nd == 0 ? "%+.0f" : (nd == 1 ? "%+.1f" : "%+.2f")); }
    else
      { return (nd == 0 ? "%.0f" : (nd == 1 ? "%.1f" : "%.2f")); }
  }

void plot_1d_ticks
  ( PSStream *ps,
    PlotOptions *o,
    dg_axis_t ax,       /* Axis along which ticks are to be placed. */
    double hsize,     /* Half-length of tick marks (client units). */
    char *fmt,        /* Format for tick labels, or NULL, */
    int nsteps,       /* Number of tick steps in {cr}. */
    interval_t *cr,  /* Only place ticks within this coordinate range. */
    double mag,       /* Scale fator from func units to client units. */
    interval_t win[] /* Plot window in client units. */
  )
  { int i;
    double tickmrg = 2.5*hsize;
    char buf[300];
    /* ?? Should check whether 0 is well inside range {win[1-ax]}. */
    /* ?? However that requires knowing the character's bbox. */
    for (i = 0; i <= nsteps; i++)
      { double t = (nsteps == 0 ? 0.5 : ((double)i)/nsteps);
        double s = (nsteps == 0 ? 0.5 : ((double)(nsteps-i))/nsteps);
        double c = s*LO(*cr) + t*HI(*cr);
        double cp = mag*c;  /* Client coordinate for plotting */
        if 
          ( (cp > LO(win[ax]) + tickmrg) && 
            (cp <= HI(win[ax]) - tickmrg) && 
            (mag*fabs(c) > tickmrg)
          )
          { if (fmt != NULL) { sprintf(buf, fmt, c); }
            if (ax == 0)
              { pswr_segment(ps, cp, -hsize, cp, +hsize);
                if (fmt != NULL) 
                  { pswr_label(ps, buf, cp, -hsize, 0.0, 0.5, 1.3); }
              }
            else if (ax == 1) 
              { pswr_segment(ps, -hsize, cp, +hsize, cp);
                if (fmt != NULL) 
                  { pswr_label(ps, buf, -hsize, cp, 0.0, 1.3, 0.5); }
              }
            else
              { affirm(FALSE, "bad axis"); }
          }
      }
    fflush(ps->file);
  }