/* Last edited on 2023-06-27 11:40:56 by stolfi */

/* ----------------------------------------------------------------------*/

  " {LO68} and {HI68} are the same for the 68.27% interval (that is," \
  " the 15.865% and  and 84.135% percentiles)"

#define dgr_quantile_NUM 5
#define dgr_quantile_FRACTIONS 0.02500,  0.15865,  0.50000,  0.84135,  0.97500
#define dgr_quantile_HEADERS   "lo 95%", "lo 68%", "median", "hi 68%", "hi 95%"
  /* Quantile fractions of runs that are to be written out to disk, and 
    their headers for {stderr}. */


        o->chProbs[s] = double_vec_new(10);
        int32_t nc = 0;
        while (argparser_next_is_number(pp))
          { double_vec_expand(&(o->chProbs[s]), nc);
            o->chProbs[s].e[nc] = argparser_get_next_double(pp, 0.0, 1.0);
            nc++;
          }
        double_vec_trim(&(o->chProbs[s]), nc);
    
void gdr_adjust_child_probs(double_vec_t *chProbs);
  /* Adjust the children probabilities {chProbs.e[0..c_max]} so that the average number of 
    children is 1 and the sum is 1; where {c_max} is {chProbs.ne-1}. 
    Will preserve the proportions of {chProbs.e[1..c_max]} but not that 
    of {chProbs.e[0]} to the rest. */

void gdr_adjust_child_probs(double_vec_t *chProbs)
  { 
    int32_t c_max = chProbs->ne - 1;  /* Child count varies in {0..c_max}. */
    
    /* Compute average num of children {nc_avg}: */
    double nc_avg = 0;
    for (int32_t c = 1; c <= c_max; c++) { nc_avg += c*chProbs->e[c]; }
    
    /* Scale {chProbs[1..c_max]} so that the average num of children is 1: */
    for (int32_t c = 1; c <= c_max; c++) { chProbs->e[c] /= nc_avg; }
    
    /* Adjust {chProbs[0]} so that sum of probs is 1: */
    double sum_p = 0;
    for (int32_t c = 1; c <= c_max; c++) { sum_p += chProbs->e[c]; }
    chProbs->e[0] = 1.0 - sum_p;
  }
  
/* ---------------------------------------------------------------------- */
    bool_t debug = FALSE;
    if (debug) { fprintf(stderr, "      > %s\n", __FUNCTION__); }

    demand(ysz_exp > 0, "invalid {ysz_exp}");
    demand(yCur >= 1, "should not be asked to adjust initial year");
    
    /* Compute histogram {chs_a} of child counts in previous year {y-1}: */
    int32_vec_t chs_a = gdr_sim_compute_child_count_histogram(y-1, fbi, par->e);
    
    int32_t i0 = fbi[y];  /* First individual of year {y}. */
    int32_t i1_raw = fbi[y+1]-1;       /* Current last individual in year {y}. */
    int32_t i1_exp = i0 + ysz_exp - 1; /* Desired last individual in year {y}. */
    
    int32_t ysz_raw = i1_raw - i0 + 1; /* Current size of year {yCur}. */
    assert(ysz_raw > 0);
    assert(ysz_exp > 0);
    if (ysz_raw < ysz_exp)
      { /* Randomly replicate individuals of year {yCur}. */
        if (debug) { fprintf(stderr, "        expanding year %d from %d to %d indivs\n", y, ysz_raw, ysz_exp); }
        int32_vec_expand(par, i1_exp);
        for (int32_t i = i1_raw+1; i <= i1_exp; i++)
          { /* Create individual {i} by twinning some other individual {ir} in year {y}: */
            int32_t ir = int32_abrandom(i0, i - 1); /* Choose individual to replicate. */
            int32_t ip = par->e[ir]; 
            par->e[i] = ip;
          }
      }
    else if (ysz_raw > ysz_exp)
      { /* Randomly delete individuals of year {y}: */
        if (debug) { fprintf(stderr, "        reducing year %d from %d to %d indivs\n", y, ysz_raw, ysz_exp); }
        int32_t i1 = i1_raw;
        while (i1 >= i1_exp)
          { /* Delete some random individual {ir} still in year {y}: */
            int32_t ir = int32_abrandom(i0, i1); /* Individual to delete. */
            if (ir < i1) { par->e[ir] = par->e[i1]; }
            i1--;
          }
      }
    else
      { /* No correction needed: */
        if (debug) { fprintf(stderr, "        year %d has desired size %d indivs\n", y, ysz_raw); }
      }
    /* Update the {fbi} table of year {y}: */
    fbi[y+1] = i1_exp + 1;

    /* Recompute histogram {chs_b} of child counts in previous year {y-1}: */
    int32_vec_t chs_b = gdr_sim_compute_child_count_histogram(y-1, fbi, par->e);

    if (debug) 
      { gdr_sim_compare_child_count_histograms("before and after adjustment", &chs_a, &chs_b); }
    
    /* Accumulate children counts in gloabal {chCounts}: */
    for (int32_t c = 0; c < chs_b.ne; c++)
      { int32_t co = chCounts->ne;
        int64_vec_expand(chCounts, c);
        while (co < chCounts->ne) { chCounts->e[co] = 0; co++; }
        chCounts->e[c] += chs_b.e[c];
      }
    free(chs_a.e);
    free(chs_b.e);
/* ---------------------------------------------------------------------- */
  
#define gdr_probs_PREC 3 
  /* Number of decimal fraction digits to round the probabilities to. */
  
#define gdr_probs_freqs_PREC (gdr_probs_PREC + 2)
  /* Number of decimal fraction digits to round actual frequencies to. */

/* ---------------------------------------------------------------------- */
  
void gdr_draw_assign_rows_compact
  ( gdr_sim_state_t *st, 
    int32_t yMin, 
    int32_t yMax, 
    int32_t *nrows_P, 
    int32_t rdr[]
  );
  /* Chooses a row {rdr[i]} for each individual {i} in {st}. Also
    determines the number of rows {nrows} needed for the plot, and
    returns it in {*nrows_P} See {gdr_draw_plot_individuals} for the
    meaning of the other parameters.
    
    This procedure tries to produce the most compact plot possible. */
  
void gdr_draw_assign_rows_by_lineage
  ( gdr_sim_state_t *st, 
    int32_t yMin, 
    int32_t yMax, 
    int32_t *nrows_P, 
    int32_t rdr[]
  );
  /* Chooses a row {rdr[i]} for each individual {i} in {st}. Also
    determines the number of rows {nrows} needed for the plot, and
    returns it in {*nrows_P} See {gdr_draw_plot_individuals} for the
    meaning of the other parameters.
    
    This procedure tries to produce a planar plot where contemporanous
    lineages are separated. */


void gdr_draw_plot_individuals
  ( epswr_figure_t *eps, 
    gdr_sim_state_t *st, 
    int32_t ncols, 
    int32_t nrows, 
    double Xstep, 
    double Ystep, 
    int32_t yMin, 
    int32_t rdr[], 
    int32_t yRef
  );
  /* Plots onto {eps} the individuals in state {st}.
  
    The main plot area is conceptually divided into an array of cells
    with {ncols} columns and {nrows} rows, where each column corresponds
    to a year and has width {Xstep}, and each row has height {Ystep}.
    Column 0 corresponds to year {yMin}, probably negative.

    Rows are numbered {0..nrows-1} from bottom up, and
    columns are numbered {0..ncols-1} from right to left. 

    Each individual {i} with birth year {ybr[i]} and death year
    {ydt[i]} occupies columns {ybr[i]-yMin..ydt[i]-yMin} inter
    {0..ncols-1} of row {rdr[i]} of this matrix.

    (The individual is technically dead on year {ydt[i]}, but may have
    a child born on that year, so it may have stuff drawn on that
    column too.)

    If {yRef} is in {0..st.ny-1}, the procedure highlights the surviving
    lineages that count for that year. */

void gdr_draw_plot_year_line
  ( epswr_figure_t *eps, 
    int32_t ncols, 
    int32_t nrows, 
    double Xstep, 
    double Ystep, 
    int32_t yMin, 
    double R, 
    double G,
    double B,
    int32_t y
  );
  /* Plots a vertical line with the given color across the drawing, at the given year {y}.
    See {gdr_draw_plot_individuals} for the meaning of the other parameters. */


void gdr_draw_assign_rows_compact
  ( gdr_sim_state_t *st, 
    int32_t yMin, 
    int32_t yMax, 
    int32_t *nrows_P, 
    int32_t rdr[]
  )
  { bool_t debug = FALSE;
    if (debug) { fprintf(stderr, "    > %s\n", __FUNCTION__); }
    
    int32_t ni = st->ni;

    int32_t ncols = yMax - yMin + 1;
    int32_t nrows = ni; /* Tentative -- to be reduced later. */
    int32_t ncells = nrows*ncols; /* Number of cells in main plot area. */
    /* Use table {occ[row*ncols + col} to tell whether a cell is occupied:  */ 
    bool_t *occ = (bool_t *)notnull(malloc(ncells*sizeof(bool_t)), "no mem"); /* Cell occupancy */
    /* Table of displacements used by {find_free_row}: */
    int32_t nd = 3*nrows; /* Num of row displacements in table. */
    int32_t drow[nd];
    for (int32_t kc = 0; kc < nrows; kc++) 
     { assert(3*kc+2 < nd);
       drow[3*kc+0] = -2*kc-1; 
       drow[3*kc+1] = -2*kc-2; 
       drow[3*kc+2] = kc+1;
     } 
     assert(drow[3*nrows-1] >= nrows-1);

    /* Used by {assign_row}, see below. */
    double rChSum[ni];  

    auto int32_t assign_row(int32_t i);
      /* Assigns a row for the life of an individual {i}, so that the
        relevant cells are free with one cell to spare on each side.
        
        Must be called for {i} only after it has been called for all
        children of {i}.
        
        The procedure assumes that, for all {j} in {0..i}, {rdr[j]} is
        the row assigned to {j} in a previous iteration, or some
        arbitary row index for the first iteration; and {rChSum[j]} is
        the sum of {rdr[i]} for all children of {j} that have already
        been assigned.
        
        If {i} has any children, the ideal row for {i} is taken to be
        {(rChSum[i] + rdr[i])/(st->nch[i] + 1)}, rounded. If {i} has no
        children, but has a parent {p}, the ideal row is set to be
        {rdr[p]}. If {i} has no chidlren and no parent, its ideal row is
        the current {rdr[i]}.
        
        In any case, the procedure will to assign {i} to the row that
        has the required cells still free and is closest to that ideal
        row. */
    
    auto int32_t find_free_row(int32_t row_best, int32_t col0, int32_t col1);
      /* Finds the row that has all the cells {col0..col1} free
        and is closest to row {row_best}. */

    /* Assign a row {rdr[i]} to every individual {i}: */
    int32_t rowMax; /* Max row assigned. */
    int32_t niters = 4; /* Iterations of the placement loop. */
    /* Pretend the initial placement was row 0: */
    for (int32_t i = 0; i < ni; i++) { rdr[i] = 0; }
    for (int32_t it = 0; it < niters; it++)
      { if (debug) { fprintf(stderr, "      --- iteration %d ---\n", it); }
        rowMax = -1;
        /* Clear {occ,rChSum}: */
        for (int32_t k = 0; k < ncells; k++) { occ[k] = FALSE; }
        for (int32_t i = 0; i < ni; i++) { rChSum[i] = 0.0; }
        /* Assign rows in reverse chrono order: */
        for (int32_t i = ni-1; i >= 0; i--) 
          { int32_t row = assign_row(i);
            assert((row >= 0) && (row < nrows));
            if (row > rowMax) { rowMax = row; }
            rdr[i] = row;
            if (debug) { fprintf(stderr, "        i = %d row = %d\n", i, row); }
            /* Update the {rChSum} entry of the parent: */
            int32_t p = st->par.e[i];
            if (p != -1) 
              { assert((p >= 0) && (p < ni));
                rChSum[p] += (double)row;
              }
          }
      }
    /* Reduce plot grid: */
    nrows = rowMax+1;
    fprintf(stderr, "      new cell grid size = %d x %d\n", ncols, nrows);
    free(occ);

    (*nrows_P) = nrows;
    if (debug) { fprintf(stderr, "    < %s\n", __FUNCTION__); }
    return;

    int32_t assign_row(int32_t i)
      { bool_t debug = FALSE;
        if (debug) { fprintf(stderr, "      > %s\n", __FUNCTION__); }
    
        /* Get main parameters of individual {i}: */
        int32_t ybri = st->ybr.e[i];
        int32_t ydti = st->ydt.e[i];
        int32_t nchi = st->nch.e[i];
        if (nchi == -1) { assert(ydti == -1); nchi = 0; ydti = ybri; }
        
        /* Determine the column range of {i}'s life, including the death year: */
        int32_t col0 = ybri - yMin;
        int32_t col1 = ydti - yMin;
        
        /* Compute the ideal row {row_best}: */
        int32_t p = st->par.e[i];
        if (p == -1) { p = i; } else { assert(p < i); }
        if (debug) { fprintf(stderr, "        rChSum[i] = %.1f nch[i] = %d p = %d old rdr[p] = %d", rChSum[i], nchi, p, rdr[p]); }
        assert((p >= 0) && (p < ni));
        assert((rdr[p] >= 0) && (rdr[p] < nrows));
        int32_t row_best = (int32_t)floor(((double)rdr[p] + rChSum[i])/(nchi + 1.0) + 0.5);
        assert((row_best >= 0) && (row_best < nrows));
        if (debug) { fprintf(stderr, " row_best = %d", row_best); }
        
        /* Choose the row for {i}: */
        int32_t row = find_free_row(row_best, col0 - 1, col1 + 1);
        if (debug) { fprintf(stderr, " row = %d cols = {%d .. %d}\n", row, col0, col1); }
        assert((0 <= row) && (row < nrows));
        
        /* Mark cells {col0..col1} on that row as occupied: */
        bool_t *occr = &(occ[row*ncols]);
        for (int32_t col = col0; col <= col1; col++) { occr[col] = TRUE; }

        if (debug) { fprintf(stderr, "      < %s\n", __FUNCTION__); }
        return row;
      }

    int32_t find_free_row(int32_t row_best, int32_t col0, int32_t col1)
      { if (col0 < 0) { col0 = 0; }
        if (col1 >= ncols) { col1 = ncols-1; }
        int32_t row = row_best;
        int32_t kd = 0; /* Index into next row increment in {drow}. */
        for (int32_t kr = 0; kr < nrows; kr++)
          { /* Try on {row}: */
            bool_t *occr = &(occ[row*ncols]);
            bool_t ok = TRUE;
            for (int32_t col = col0; (col <= col1) && ok; col++)
              { if (occr[col]) { ok = FALSE; } }
            if (ok) { return row; }
            
            /* Get to the next closest valid row: */
            row = -1;
            while ((kd < nd) && ((row < 0) || (row >= nrows))) 
              { row = row_best + drow[kd]; kd++; }
            
            /* Must succeed eventually: */
            assert((row >= 0) && (row < nrows));
          }
        demand(FALSE, "failed to find a free row");
      }
  }    
    
void gdr_draw_assign_rows_by_lineage
  ( gdr_sim_state_t *st, 
    int32_t yMin, 
    int32_t yMax, 
    int32_t *nrows_P, 
    int32_t rdr[]
  )
  { bool_t debug = FALSE;
    if (debug) { fprintf(stderr, "    > %s\n", __FUNCTION__); }
    
    int32_t ni = st->ni;

    int32_t ncols = yMax - yMin + 1;
    int32_t nrows = ni; /* Tentative -- to be reduced later. */
    
    /* The general strategy is to recursively assign rows to the lineage
      of each individual. First some children are
      considered in increasing order of birth year, and their lineages
      get assigned recursively. Then the individual is assigned to the
      lowest possible row. Then the remaining children are considered
      in reverse chrono order, and their lineages of the recursively
      assigned. Then the individual is reassigned, if possible, to a
      better position, still above the frist set of chidlren and below the
      second set ones, if any.
      
      Each lineage is deposited on top of all other individuals already
      assigned. For that we keep a /horizon/ table {hir[0..ncols-1]}
      with the higest occupied row on each column. */
      
    /* The horizon table: */
    int32_t hir[ncols];
    for (int32_t col = 0; col < ncols; col++) { hir[col] = -1; }
    
    /* Initialize {rdr[i]} to {-1} for all {i}.  That will mean "not assigned yet". */
    for (int32_t i = 0; i < ni; i++) { rdr[i] = -1; }
    
    auto int32_t assign_lineage_rows(int32_t i);
      /* Sets {rdr[i]}, and {rdr[j]} for all descendants of {i}.
        Updates the horizon {hir[0..ncols-1]} as needed. 
        Returns the number of individuals assigned.*/
    
    /* Scan individuals in chrono order, assign lineages of unassigned ones: */
    int32_t rowMax = -1;
    for (int32_t i = 0; i < ni; i++)
      { if (rdr[i] == -1)
          { /* Still unassigned: */
            assert(st->par.e[i] == -1); /* Otherwise should have been assigned. */
            int32_t nass = assign_lineage_rows(i);
            if (debug) { fprintf(stderr, "        root indiv %d assigned to row %d, %d descendants\n", i, rdr[i], nass); }
            assert(rdr[i] >= 0);
          }
        if (rdr[i] > rowMax) { rowMax = rdr[i]; }
      }
    
    /* Reduce plot grid: */
    nrows = rowMax+1;
    fprintf(stderr, "      new cell grid size = %d x %d\n", ncols, nrows);

    (*nrows_P) = nrows;
    if (debug) { fprintf(stderr, "    < %s\n", __FUNCTION__); }
    return;
    
    auto int32_t assign_to_lowest_free_row(int32_t col0, int32_t col1);
      /* Finds the lowest row where it is possible to place a 
        a lifeline that spans grid columns {col0...col1} so that it is 
        above all previous lifelines that enter those columns
        and has at least one free cell on each side.
        
        Also updates {hir[col0..col1]} to that row. */
        
    auto bool_t child_goes_below(int32_t r, int32_t nchi);
      /* Decides whether child {r} out of {nchi} children 
        should go below or above the parent: */
    
    int32_t assign_lineage_rows(int32_t i)
      { bool_t debug = FALSE;
        if (debug) { fprintf(stderr, "      > %s\n", __FUNCTION__); }
    
        /* Get main parameters of individual {i}: */
        int32_t ybri = st->ybr.e[i];
        int32_t ydti = st->ydt.e[i];
        int32_t nchi = st->nch.e[i];
        if (nchi == -1) { assert(ydti == -1); nchi = 0; ydti = ybri; }
        
        /* Determine the column range of {i}'s life, including the death year: */
        int32_t col0 = ybri - yMin;
        int32_t col1 = ydti - yMin;
        
        /* Assign its lineage: */
        int32_t nass = 1; /* Number of individuals in lineage. */
        int32_t rowi = -1;
        if (nchi == 0)
          { /* Just place {i} as low as possible. */
            rowi = assign_to_lowest_free_row(col0, col1);
          }
        else
          { /* Find all the chidlren. We need not scan too far since life is short. */
            int32_t chi[nchi];
            int32_t kchi = 0; /* Counts children found. */
            for (int32_t j = 0; (j < ni) && (kchi < nchi); j++)
              { int32_t p = st->par.e[j];
                if (p == i) { chi[kchi] = j; kchi++; }
              }
            assert(kchi == nchi);

            /* Sort children in order of birth year: */
            for (int32_t r = 1; r < nchi; r++)
              { int32_t jr = chi[r];
                int32_t ybrr = st->ybr.e[jr];
                int32_t s = r; 
                while ((s > 0) && (st->ybr.e[chi[s-1]] > ybrr))
                  { chi[s] = chi[s-1]; s--; }
                chi[s] = jr;
              }
              
            /* Assign some children in order of birth: */
            for (int32_t r = 0; r < nchi; r++)
              { if (child_goes_below(r, nchi))
                  { int32_t nassj = assign_lineage_rows(chi[r]);
                    nass += nassj;
                  }
              }

            /* Tentatively assign {i} as low as possible: */
            rowi = assign_to_lowest_free_row(col0, col1);
            
            /* Assign the remaining children in reverse order: */
            for (int32_t r = nchi-1; r >= 0; r--)
              { if (! child_goes_below(r, nchi))
                  { int32_t nassj = assign_lineage_rows(chi[r]);
                    nass += nassj;
                  }
              }
              
            /* Raise {i} to midway of its children: */
            fprintf(stderr, "!! individual %d -- raising not implemented\n", i);
          }
        rdr[i] = rowi;
        if (debug) { fprintf(stderr, "          indiv %d assigned to row %d, %d children, %d descendants \n", i, rowi, nchi, nass); }
        if (debug) { fprintf(stderr, "      < %s\n", __FUNCTION__); }
        return nass;
      }
      
    bool_t child_goes_below(int32_t r, int32_t nchi)
      { /* Let's try to put 1, even below: */
        return (r == 1) || ((r % 2) == 0);
      }
    
    int32_t assign_to_lowest_free_row(int32_t col0, int32_t col1)
      { assert((0 <= col0) && (col0 <= col1) && (col1 < ncols));
        /* Find maximum of {hir[col0-1..col1+1} to leave space on each side: */
        int32_t hirMax = -1;
        for (int32_t col = col0-1; col <= col1+1; col++)
          { if ((col >= 0) && (col < ncols))
              { if (hir[col] > hirMax) { hirMax = hir[col]; } }
          }

        /* Assign to row above: */
        int32_t row = hirMax + 1;
        assert((row >= 0) && (row < nrows));

        /* Update horizon on cols {col0..col1}: */
        for (int32_t col = col0; col <= col1; col++) { hir[col] = row; }
        
        return row;
      }
  }

void gdr_draw_plot_individuals
  ( epswr_figure_t *eps, 
    gdr_sim_state_t *st, 
    int32_t ncols, 
    int32_t nrows, 
    double Xstep, 
    double Ystep, 
    int32_t yMin, 
    int32_t rdr[], 
    int32_t yRef
  )
  { bool_t debug = FALSE;
    if (debug) { fprintf(stderr, "    > %s\n", __FUNCTION__); }
    
    int32_t ni = st->ni; /* Number of individuals created. */
    
    double lifw = 0.50*Ystep; /* Life line width. */

    double dbrr = 0.40*lifw;  /* Radius of birth dot. */

    double arlw = 0.33*lifw;  /* Arrow line width. */
    double alen = 4.00*arlw;  /* Arrowhead length. */
    double awid = 1.50*arlw;  /* Arrowhead width. */
    double halw = 0.25*arlw;  /* Width of arrows white halo. */
    double datr = 0.50*arlw;  /* Radius of arrow tail dot. */
 
    auto void draw_indiv_pass_0(int32_t i);
    auto void draw_indiv_pass_1(int32_t i);
      /* Draws parts of an individual {i} that belong to layer 0 and layer 1,
        respectively.  */
    
    auto void draw_life(int32_t col0, int32_t col1, int32_t row);
      /* Draws the life of an individual on the given {row},  
        spanning columns {col0..col1} (corresponding to its 
        years of birth and death). */
        
    auto void draw_arrow(int32_t col, int32_t row0, int32_t row1);
      /* Draws a vertical arrow from row {row0}
        to row {row1} on column {col}.  Meant to be 
        from some point on an individual life to the 
        start of the life of one of its children. */
      
    auto void draw_tail_dot(int32_t col, int32_t row);
      /* Draws a dot at the given row and column, meant to be
        the point on the life of some individial when it had a child. */
      
    auto void draw_birth_dot(int32_t col, int32_t row, bool_t exp);
      /* Draws a dot at the given column and row, meant to indicate 
        the birth of an individual. */
    
    /* Plot the individuals: */
    for (int32_t pass = 0; pass < 2; pass++) 
      { /* Pass 0 draws lives, pass 1 draws arrows. */
        for (int32_t i = 0; i < ni; i++)
          { if (pass == 0)
              { draw_indiv_pass_0(i); }
            else if (pass == 1)
              { draw_indiv_pass_1(i); }
          }
      }

    if (debug) { fprintf(stderr, "    < %s\n", __FUNCTION__); }
    return;

    void draw_indiv_pass_0(int32_t i)
      { int32_t ybri = st->ybr.e[i];
        int32_t ydti = st->ydt.e[i];
        int32_t nchi = st->nch.e[i];
        if (nchi != -1)
          { /* Expanded. Draw fat lifetime segment: */
            assert(ydti >= ybri);
            int32_t col0 = ybri - yMin;
            int32_t col1 = ydti - yMin;
            int32_t row = rdr[i];
            if (debug) 
              { fprintf(stderr, "      life of i = %d years {%d .. %d}", i, ybri, ydti);
                fprintf(stderr, " cols = {%d .. %d} row = %d\n", col0, col1, row);
              }
            draw_life(col0, col1, row);
          }
      }

    void draw_indiv_pass_1(int32_t i)
      { int32_t ybri = st->ybr.e[i]; /* Year of birth. */
        int32_t ydti = st->ydt.e[i]; /* Year of death. */
        int32_t nchi = st->nch.e[i]; /* Number of children. */
        if (nchi == -1) { assert(ydti == -1); ydti = ybri; }
      
        int32_t p = st->par.e[i]; /* Parent index. */
        if (p != -1)
          { /* Draw parenting arrow: */
            assert((p >= 0) && (p < ni));
            int32_t col = ybri - yMin;
            int32_t row0 = rdr[p];
            int32_t row1 = rdr[i];
            draw_arrow(col, row0, row1);
            draw_tail_dot(col, row0);
          }

        { /* Draw birth dot: */
          bool_t exp = (nchi != -1);
          int32_t col = ybri - yMin;
          int32_t row = rdr[i];
          draw_birth_dot(col, row, exp);
        }
      }

    void draw_life(int32_t col0, int32_t col1, int32_t row)
      { if (debug) { fprintf(stderr, "      %s: columns {%d..%d}\n", __FUNCTION__, col0, col1); }
        demand((0 <= col0) && (col0 <= col1) && (col1 < ncols), "bad column range");
        demand((0 <= row) && (row < nrows), "bad row");
        double X0 = Xstep*(ncols - 1 - col0 + 0.5);
        double Y0 = Ystep*(row + 0.5);

        double X1 = Xstep*(ncols - 1 - col1 + 0.5);
        double Y1 = Y0;

        epswr_set_pen(eps, 0.800,0.500,0.200, lifw, 0,0);
        epswr_segment(eps, X0, Y0, X1, Y1);
      }

    void draw_arrow(int32_t col, int32_t row0, int32_t row1)
      { demand((0 <= col) && (col < ncols), "bad column");
        demand((0 <= row0) && (row0 < nrows), "bad row {row0}");
        demand((0 <= row1) && (row1 < nrows), "bad row {row1}");
        demand(row0 != row1, "row collision");
        double X0 = Xstep*(ncols - 1 - col + 0.5);
        double Y0 = Ystep*(row0 + 0.5);
        double X1 = X0;
        double Y1 = Ystep*(row1 + 0.5);
        Y1 += (Y1 > Y0 ? -dbrr : + dbrr);
        
        /* Draw the white shadow: */
        epswr_set_pen(eps, 1.000,1.000,1.000, arlw + 2*halw, 0,0);
        epswr_segment(eps, X0, Y0, X1, Y1);
        epswr_arrowhead(eps, X0, Y0, X1, Y1, awid, alen, 1.0, TRUE, TRUE);

        /* Draw the arrow proper: */
        epswr_set_pen(eps, 0.000,0.000,0.000, arlw, 0,0);
        epswr_segment(eps, X0, Y0, X1, Y1);
        epswr_arrowhead(eps, X0, Y0, X1, Y1, awid, alen, 1.0, TRUE, TRUE);
      }

    void draw_birth_dot(int32_t col, int32_t row, bool_t exp)
      { double X = Xstep*(ncols - 1 - col + 0.5);
        double Y = Ystep*(row + 0.5);

        epswr_set_pen(eps, 0.000,0.000,0.000, arlw, 0,0);
        epswr_dot(eps, X, Y, dbrr, exp, TRUE);
      }

    void draw_tail_dot(int32_t col, int32_t row)
      { double X = Xstep*(ncols - 1 - col + 0.5);
        double Y = Ystep*(row + 0.5);

        epswr_set_pen(eps, 0.000,0.000,0.000, arlw, 0,0);
        epswr_dot(eps, X, Y, datr, TRUE, TRUE);
      }
  }
    
void gdr_draw_plot_year_line
  ( epswr_figure_t *eps, 
    int32_t ncols, 
    int32_t nrows, 
    double Xstep, 
    double Ystep, 
    int32_t yMin, 
    double R, 
    double G,
    double B,
    int32_t y
  )
  {
    int32_t col = y - yMin;
    double X = Xstep*(ncols - 1 - col + 0.5);
    double Y0 = Ystep*(-2 + 0.5);
    double Y1 = Ystep*(nrows + 1 + 0.5);

    double linw = 0.25*Ystep;
    
    epswr_set_pen(eps, R,G,B, linw, 0,0);
    epswr_segment(eps, X, Y0, X, Y1);
  }

/* ---------------------------------------------------------------------- */

epswr_figure_t *gdr_draw_create_eps_figure(char *prefix, double Xsize, double Ysize);
  /* Creates an EPS figure object that writes to the file "{prefix}-evol.eps"
    
    The figure will measure {Xsize} by {Ysize} millimeters, with a 
    small margin all around. */

epswr_figure_t *gdr_draw_create_eps_figure(char *prefix, double Xsize, double Ysize)
  {
    double pt_per_mm = 72.0/25.4;
    
    /* Device plot window size and margins (pt): */
    double hPlotSize = Xsize*pt_per_mm;
    double vPlotSize = Ysize*pt_per_mm;
    double figMrg = 2.0*pt_per_mm;  /* Margin outside device plot window. */
    
    bool_t eps_verbose = FALSE;

    char *name = jsprintf("%s-evol", prefix);
    epswr_figure_t *eps = epswr_new_named_figure
      ( NULL, NULL, name, -1, NULL, 
        hPlotSize, vPlotSize, figMrg, figMrg, figMrg, figMrg, 
        eps_verbose
      );
    free(name);
    
    /* Set the client plot window: */
    double cliMrg = 2.0; /* Margin between device window and client window (mm): */
    double xMin = -cliMrg, xMax = Xsize + 2*cliMrg;
    double yMin = -cliMrg, yMax = Ysize + 2*cliMrg;
    epswr_set_client_window(eps, xMin, xMax, yMin, yMax);
    
    return eps;
  }
  
  {
    /* Find range of {y} that actually occur in {st}: */
    int32_t yMin = INT32_MAX;
    int32_t yMax = INT32_MIN;
    for (int32_t i = 0; i < ni; i++)
      { int32_t ybri = st->ybr.e[i]; /* Year of birth. */
        int32_t ydti = st->ydt.e[i]; /* Year of death. */
        int32_t nchi = st->nch.e[i];
        if (nchi == -1) { assert(ydti == -1); ydti = ybri; }
        assert(ydti >= ybri);
        if (ybri < yMin) { yMin = ybri; }
        if (ydti > yMax) { yMax = ydti; }
      }
    assert(yMin <= yMax);
    fprintf(stderr, "    full year range {%d .. %d}\n", yMin, yMax);
  }