#define PROG_NAME "msdemo"
#define PROG_DESC "plots demographic profiles for faculty career planning"
#define PROG_VERS "1.0"

#define msdemo_C_COPY \
  "Copyright © 2007 by the State University of Campinas (UNICAMP)."

/* Last edited on 2024-11-20 06:06:49 by stolfi */
/* See the authorship, rights and warranty notices in the PROG_INFO below. */

#define PROG_HELP \
  "  " PROG_NAME " \\\n" \
  "    -blabber {AMOUNT} \\\n" \
  "    " argparser_help_info_HELP " \\\n" \
  "    < {INFILE} \\\n" \
  "    > {OUTFILE}"

#define PROG_INFO \
  "NAME\n" \
  "  " PROG_NAME " - " PROG_DESC "\n" \
  "\n" \
  "SYNOPSIS\n" \
  PROG_HELP "\n" \
  "\n" \
  "DESCRIPTION\n" \
  "  The program bla bla bla" \
  " bla bla bla bla {X+Y} bla bla" \
  " bla {INFILE} bla \"foobar.ppm\" bla bla bla\n" \
  "\n" \
  "  Beware that bla bla bla BLEBBLE BLOB bla" \
  " bla bla bla bla.\n" \
  "\n" \
  "OPTIONS\n" \
  "  -blabber {AMOUNT}\n" \
  "    Blabbers for that {AMOUNT}. May also bla" \
  " bla bla bla bla bla bla bla bla bla bla bla bla" \
  " bla bla bla bla bla bla bla.\n" \
  "\n" \
  "DOCUMENTATION OPTIONS\n" \
  argparser_help_info_HELP_INFO "\n" \
  "\n" \
  "SEE ALSO\n" \
  "  ioplot(1).\n" \
  "\n" \
  "AUTHOR\n" \
  "  Created 2007-06-02 by Jorge Stolfi, IC-UNICAMP.\n" \
  "\n" \
  "MODIFICATION HISTORY\n" \
  "  Option bla bla added by J. Stolfi, IC-UNICAMP.\n" \
  "\n" \
  "WARRANTY\n" \
  argparser_help_info_NO_WARRANTY "\n" \
  "\n" \
  "RIGHTS\n" \
  "  " msdemo_C_COPY "\n" \
  "\n" \
  argparser_help_info_STANDARD_RIGHTS

#define _GNU_SOURCE
#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <math.h>

#include <bool.h>
#include <pswr.h>
#include <frgb.h>
#include <argparser.h>

/* DATA TYPES AND CONSTANTS */

typedef enum
  { st_3H  =  0,
    st_3E  =  1,
    st_3R  =  2,

    st_5H  =  3,
    st_5E  =  4,
    st_5R  =  5,

    st_6a  =  6,
    st_6b  =  7,
    st_6c  =  8,
    st_6d  =  9,
    
    st_6H  = 10,
    
    st_6p  = 11,
    st_6q  = 12,
    st_6r  = 13,
    st_6s  = 14,
    
    st_6E  = 15,
    
    st_6R  = 16,

    st_XX  = 17 
  } state_t;
  /* Codes of career states for a faculty member. */

#define msd_NS 18
  /* Number of career states. */

typedef struct stdesc_t
  { char *tag;     /* Two-letter code of state ("3H", "XX", etc.) */
    char *def;     /* verbose definition, for plot legend. */
    frgb_t color;  /* Color for plot. */
  } stdesc_t; 
  
/* COMMAND-LINE OPTIONS */

typedef struct options_t
  { bool_t op1;     /* A boolean argument. */
    char* op2;      /* A string-valued argument. */
  } options_t;

/* INTERNAL PROTOTYPES */

PSStream *msd_init_plot
  ( int aMax, 
    double pMax, 
    double hSize, 
    double vSize, 
    double *hScale, 
    double *vScale
  );
  /* Creates an Encapsulated Postscript stream for a demographic plot
    for maximum years since hiring {aMax}, and maximum relative cohort
    size {pMax}. Also draws axes, tics, scales, etc.. The plot proper
    (excluding axis arrows, tics, etc.) will measure {hSize} by
    {vSize} mm. 
    
    Returns in {*hScale} and {*vScale} the scales to use in 
    the horizontal and vertical axes. */

void msd_plot_demography(PSStream *ps, double pIni, int aMax, double hScale, double vScale);
  /* Computes and plots the demography profile to stream {ps}. 
    The initial relative cohort size will be {pIni}.
    The simulation will run until {aMax} years since hiring. */

void msd_plot_cohort
  ( PSStream *ps, 
    double a, 
    stdesc_t *S, 
    double *fS, 
    int ns, 
    double hScale, 
    double vScale
  );
/* Plots the color bars for the cohort with {a} years since hiring.
  Assumes that there are {ns} states, with relative population sizes 
  {fS[0..ns-1]}.  The colorbar for state {s} is filled with {S[s].color}. */

double msd_trans_prob
  ( int a, 
    state_t r, 
    state_t s,
    double pD,
    double pE,
    double pR,
    double pL,
    double pT
  );
  /* Probability that a faculty member who is in state {r}, {a} years
    since hiring, will be in state {s} after one more year.
    The arguments {pD,pE,pR,pL,pT} are the probabilities of 
    events {D,E,R,L,T} happening, assuming that the faculty 
    member qualifies for them. */
  
double msd_prob_D(double a); /* Prob of a live faculty dying. */
double msd_prob_E(double a); /* Prob of ineligible for retirement becoming eligible. */
double msd_prob_R(double a); /* Prob of eligible actually retiring. */
double msd_prob_L(double a); /* Prob of an MS-3 passing to MS-5. */
double msd_prob_T(double a); /* Prob of an MS-5 passing to MS-6. */
  /* Each function {msd_prob_{X}} gives the probability of event {X}
    occurring to a faculty member who has been hired {a} years ago,
    irrespective of other events which may occur in the same year.  
    
    Each function assumes that the necessary requirement for the event
    have been met. In particular, {msd_prob_D} assumes that the person
    is currently alive; and all other functions assume that the person
    is alive and does not die during the year. So, for example, the
    probability of a faculty member currently in state {st_3H}
    retiring as MS-6, all in the same year, is

      {(1-msd_prob_D(a))*msd_prob_L(a)*msd_prob_T(a)*msd_prob_E(a)*msd_prob_R(a)}

    */

double msd_sigmoid(double x, double xMin, double xAvg, double xDev, double xMax);
  /* A sigmoid function of {x} that is 0 at {xMin} and 1 at {xMax},  has inflection point
     at {xAvg} and jerk points at {xAvg ± xDev}. */
     
double msd_humpoid(double x, double xMin, double xAvg, double xDev, double xMax);
  /* A bell-shaped function of {x} which is 0 at xMin, 1 at {xMax},
    has approximate maximum at {xAvg} and inflection points at
    approximately {xAvg ± xDev}.
    
    If {xMax-xAvg} and {xAvg-xMin} are equal, or both are greated than
    {5*xDev}, then then the position of the maximum and inflection
    points are exact. */

double msd_erf(double x, double xAvg, double xDev);
   /* The integral from {-oo} to {x} of a normal distribution with mean 
     {xAvg} and deviation {xDev}.  It is 0 at {-oo} and 1 at {+oo}. */

double msd_gauss(double x, double xAvg, double xDev);
   /* The Gaussian curve with mean {xAvg} and deviation {xDev},
     normalized so that it is 1 at {xAvg}. */

options_t *parse_options(int argc, char **argv);
  /* Parses the command line arguments and packs them as an {options_t}. */

int main(int argc,char** argv);

/* IMPLEMENTATIONS */

int main(int argc, char** argv)
  {
    options_t *o = parse_options(argc, argv);
    
    /* The horizontal coord is the number of years since hiring. */
    /* The vertical coord is the relative cohort size. */

    /* Assumes that people are hired at about 28 years of age. */
    
    int aMax = 65;  /* Max years since hiring. */
    double pIni = 1.00;  /* Initial relative cohort size. */
    double pMax = pIni;  /* Max relative cohort size. */

    /* Main plot area dimensions excl. arrows, tics, etc. (mm): */
    double hSize = 200.00; 
    double vSize = 100.00;
    
    double hScale, vScale;
    PSStream *ps = msd_init_plot(aMax, pMax, hSize, vSize, &hScale, &vScale);
    
    msd_plot_demography(ps, pIni, aMax, hScale, vScale);
    
    pswr_close_stream(ps);
    
    return 0;
  }

void msd_plot_demography(PSStream *ps, double pIni, int aMax, double hScale, double vScale)
  {
    /* 
      Main events:

        H = hiring (at MS-3 level).
        L = promotion to MS-5 (Livre Docente).
        T = promotion to MS-6 (Titular).
        E = eligible for retirement. 
        R = actual retirement.
        D = death.

      Secondary events:

        M = mature (five years after T).
    */

    /* Main states: */
    int ns = msd_NS;
    stdesc_t S[ns];
    S[st_3H] = (stdesc_t){ "3H", "MS-3, não aposentável",       (frgb_t){{ 0.000f,0.600f,0.000f }} };
    S[st_3E] = (stdesc_t){ "3E", "MS-3, aposentável",           (frgb_t){{ 0.500f,0.800f,0.500f }} };
    S[st_3R] = (stdesc_t){ "3R", "MS-3, aposentado",            (frgb_t){{ 0.000f,0.300f,0.000f }} };

    S[st_5H] = (stdesc_t){ "5H", "MS-5, não aposentável",       (frgb_t){{ 0.900f,0.800f,0.000f }} };
    S[st_5E] = (stdesc_t){ "5E", "MS-5, aposentável",           (frgb_t){{ 0.950f,0.900f,0.500f }} };
    S[st_5R] = (stdesc_t){ "5R", "MS-5, aposentado",            (frgb_t){{ 0.450f,0.400f,0.000f }} };

    S[st_6a] = (stdesc_t){ "6a", "MS-6 há 1 anos, não aposen.", (frgb_t){{ 1.000f,0.200f,0.000f }} };
    S[st_6b] = (stdesc_t){ "6b", "MS-6 há 2 anos, não aposen.", (frgb_t){{ 1.000f,0.200f,0.000f }} };
    S[st_6c] = (stdesc_t){ "6c", "MS-6 há 3 anos, não aposen.", (frgb_t){{ 1.000f,0.200f,0.000f }} };
    S[st_6d] = (stdesc_t){ "6d", "MS-6 há 4 anos, não aposen.", (frgb_t){{ 1.000f,0.200f,0.000f }} };

    S[st_6H] = (stdesc_t){ "6H", "MS-6 antigo, não aposen.",    (frgb_t){{ 1.000f,0.200f,0.000f }} };
 
    S[st_6p] = (stdesc_t){ "6p", "MS-6 há 1 anos, aposentável", (frgb_t){{ 1.000f,0.600f,0.500f }} };
    S[st_6q] = (stdesc_t){ "6q", "MS-6 há 2 anos, aposentável", (frgb_t){{ 1.000f,0.600f,0.500f }} };
    S[st_6r] = (stdesc_t){ "6r", "MS-6 há 3 anos, aposentável", (frgb_t){{ 1.000f,0.600f,0.500f }} };
    S[st_6s] = (stdesc_t){ "6s", "MS-6 há 4 anos, aposentável", (frgb_t){{ 1.000f,0.600f,0.500f }} };

    S[st_6E] = (stdesc_t){ "6E", "MS-6 antigo, aposentável",    (frgb_t){{ 1.000f,0.600f,1.000f }} };

    S[st_6R] = (stdesc_t){ "6R", "MS-6 aposentado",             (frgb_t){{ 0.500f,0.100f,0.000f }} };

    S[st_XX] = (stdesc_t){ "XX", "falecido",                    (frgb_t){{ 0.750f,0.750f,0.750f }} };

    /* Fractions in steady-state profile: */
    double fS[ns];     /* Rel. num. of faculty in current cohort, per state. */
    double gS[ns];     /* Rel. num. of faculty in next cohort, per state. */

    int a, r, s;

    /* Initialize fractions for first (recently-hired) cohort: */  
    for (s = 0; s < ns; s++) { fS[s] = (s == 0 ? pIni : 0.000); }

    for (a = 0; a < aMax; a++)
      { /* Plot cohort fractions: */
        msd_plot_cohort(ps, a, S, fS, ns, hScale, vScale);

        /* Event probabilities for the year (assumed to be independent of level): */
        double pD = msd_prob_D(a); /* Prob of dying. */
        double pE = msd_prob_E(a); /* Prob of becoming eligible for retirement. */
        double pR = msd_prob_R(a); /* Prob of eligible actually retiring. */
        double pL = msd_prob_L(a); /* Prob of an MS-3 passing to MS-5. */
        double pT = msd_prob_T(a); /* Prob of an MS-5 passing to MS-6. */

        fprintf 
          ( stderr, "%4d  D = %6.4f  E = %6.4f  R = %6.4f   L = %6.4f  T = %6.4f\n", 
            a, pD, pE, pR, pL, pT
          );

        /* Compute fractions {gs} of next cohort: */
        for (s = 0; s < ns; s++)
          { /* Compute fraction {gS[s]} of state {s} in next cohort: */
            gS[s] = 0.0; 
            for (r = 0; r < ns; r++)
              { /* Compute transition prob {ptrans} from state {r} to state {s}: */
                double ptrans = msd_trans_prob(a, r, s, pD,pE,pR,pL,pT);
                /* Apply transition: */
                gS[s] += fS[r]*ptrans;
              }
          }

        /* Copy next fractions to current fractions: */
        for (s = 0; s < ns; s++) { fS[s] = gS[s]; }
      }
  }
  
void msd_plot_cohort
  ( PSStream *ps, 
    double a, 
    stdesc_t *S, 
    double *fS, 
    int ns, 
    double hScale, 
    double vScale
  )
  {
    /* Horizontal range covered by colored bars: */
    double xlo = (a + 0.10)*hScale;
    double xhi = xlo + 0.80*hScale;
    
    /* Scan states: */
    int s;
    double ylo = 0;  /* Bottom of current color bar: */
    for (s = 0; s < ns; s++)
      { 
        /* Top of color bar for state {s}: */
        double yhi = ylo + fS[s]*vScale;
        /* Select colors and pen width: */
        frgb_t *rgb = &(S[s].color);
        double R = rgb->c[0], G = rgb->c[1], B = rgb->c[2];
        pswr_set_fill_color(ps, R, G, B);
        pswr_set_pen(ps, 0,0,0, 0.10, 0,0);
        /* Draw the color bar: */
        pswr_rectangle(ps, xlo, xhi, ylo, yhi, TRUE, TRUE);
        /* Prepare for next state: */
        ylo = yhi;
      }
  }

double msd_trans_prob
  ( int a, 
    state_t r, 
    state_t s,
    double pD,
    double pE,
    double pR,
    double pL,
    double pT
  )
  {
            
    switch(r)
      {
        case st_3H:
          /* Possible transitions:
            D -> XX
            !D & L & E & R  -> 5R
            !D & L & E & !R -> 5E
            !D & L & !E -> 5H
            !D & !L & E & R -> 3R
            !D & !L & E & !R -> 3E
            !D & !L & !E -> 3H
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_5R)
              { return (1-pD)*pL*pE*pR; }
            else if (s == st_5E)
              { return (1-pD)*pL*pE*(1-pR); }
            else if (s == st_5H)
              { return (1-pD)*pL*(1-pE); }
            else if (s == st_3R)
              { return (1-pD)*(1-pL)*pE*pR; }
            else if (s == st_3E)
              { return (1-pD)*(1-pL)*pE*(1-pR); }
            else if (s == st_3H)
              { return (1-pD)*(1-pL)*(1-pE); }
            else
              { return 0; }
          }
        case st_3E:
         /* Possible transitions:
            D -> XX
            !D & L & R  -> 5R
            !D & L & !R -> 5E
            !D & !L & R -> 3R
            !D & !L & !R -> 3E
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_5R)
              { return (1-pD)*pL*pR; }
            else if (s == st_5E)
              { return (1-pD)*pL*(1-pR); }
            else if (s == st_3R)
              { return (1-pD)*(1-pL)*pR; }
            else if (s == st_3E)
              { return (1-pD)*(1-pL)*(1-pR); }
            else
              { return 0; }
          }
        case st_3R:
          /* Possible transitions:
            D -> XX
            !D -> 3R
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_3R)
              { return (1-pD); }
            else
              { return 0; }
          }
        case st_5H:
          /* Possible transitions:
            D -> XX
            !D & T & E & R  -> 6R
            !D & T & E & !R -> 6p
            !D & T & !E -> 6a
            !D & !T & E & R -> 5R
            !D & !T & E & !R -> 5E
            !D & !T & !E -> 5H
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD)*pT*pE*pR; }
            else if (s == st_6p)
              { return (1-pD)*pT*pE*(1-pR); }
            else if (s == st_6a)
              { return (1-pD)*pT*(1-pE); }
            else if (s == st_5R)
              { return (1-pD)*(1-pT)*pE*pR; }
            else if (s == st_5E)
              { return (1-pD)*(1-pT)*pE*(1-pR); }
            else if (s == st_5H)
              { return (1-pD)*(1-pT)*(1-pE); }
            else
              { return 0; }
          }
        case st_5E:
          /* Possible transitions:
            D -> XX
            !D & T & R  -> 6R
            !D & T & !R -> 6p
            !D & !T & R -> 5R
            !D & !T & !R -> 5E
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD)*pT*pR; }
            else if (s == st_6p)
              { return (1-pD)*pT*(1-pR); }
            else if (s == st_5R)
              { return (1-pD)*(1-pT)*pR; }
            else if (s == st_5E)
              { return (1-pD)*(1-pT)*(1-pR); }
            else
              { return 0; }
          }
        case st_5R:
          /* Possible transitions:
            D -> XX
            !D -> 5R
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_5R)
              { return (1-pD); }
            else
              { return 0; }
          }
        case st_6a:
          /* Possible transitions:
            D -> XX
            !D & E & R -> 6R
            !D & E & !R -> 6q
            !D & !E -> 6b
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD)*pE*pR; }
            else if (s == st_6q)
              { return (1-pD)*pE*(1-pR); }
            else if (s == st_6b)
              { return (1-pD)*(1-pE); }
            else
              { return 0; }
          }
        case st_6b:
          /* Possible transitions:
            D -> XX
            !D & E & R -> 6R
            !D & E & !R -> 6r
            !D & !E -> 6c
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD)*pE*pR; }
            else if (s == st_6r)
              { return (1-pD)*pE*(1-pR); }
            else if (s == st_6c)
              { return (1-pD)*(1-pE); }
            else
              { return 0; }
          }
        case st_6c:
          /* Possible transitions:
            D -> XX
            !D & E & R -> 6R
            !D & E & !R -> 6s
            !D & !E -> 6d
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD)*pE*pR; }
            else if (s == st_6s)
              { return (1-pD)*pE*(1-pR); }
            else if (s == st_6d)
              { return (1-pD)*(1-pE); }
            else
              { return 0; }
          }
        case st_6d:
          /* Possible transitions:
            D -> XX
            !D & E & R -> 6R
            !D & E & !R -> 6E
            !D & !E -> 6H
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD)*pE*pR; }
            else if (s == st_6E)
              { return (1-pD)*pE*(1-pR); }
            else if (s == st_6H)
              { return (1-pD)*(1-pE); }
            else
              { return 0; }
          }
        case st_6H:
          /* Possible transitions:
            D -> XX
            !D & E & R -> 6R
            !D & E & !R -> 6E
            !D & !E -> 6H
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD)*pE*pR; }
            else if (s == st_6E)
              { return (1-pD)*pE*(1-pR); }
            else if (s == st_6H)
              { return (1-pD)*(1-pE); }
            else
              { return 0; }
          }
        case st_6p:
          /* Possible transitions:
            D -> XX
            !D & R -> 6R
            !D & !R -> 6q
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD)*pR; }
            else if (s == st_6q)
              { return (1-pD)*(1-pR); }
            else
              { return 0; }
          }
        case st_6q:
          /* Possible transitions:
            D -> XX
            !D & R -> 6R
            !D & !R -> 6r
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD)*pR; }
            else if (s == st_6r)
              { return (1-pD)*(1-pR); }
            else
              { return 0; }
          }
        case st_6r:
          /* Possible transitions:
            D -> XX
            !D & R -> 6R
            !D & !R -> 6s
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD)*pR; }
            else if (s == st_6s)
              { return (1-pD)*(1-pR); }
            else
              { return 0; }
          }
        case st_6s:
          /* Possible transitions:
            D -> XX
            !D & R -> 6R
            !D & !R -> 6E
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD)*pR; }
            else if (s == st_6E)
              { return (1-pD)*(1-pR); }
            else
              { return 0; }
          }
        case st_6E:
          /* Possible transitions:
            D -> XX
            !D & R -> 6R
            !D & !R -> 6E
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD)*pR; }
            else if (s == st_6E)
              { return (1-pD)*(1-pR); }
            else
              { return 0; }
          }
        case st_6R:
          /* Possible transitions:
            D -> XX
            !D -> 6R
          */
          { if (s == st_XX) 
              { return pD; }
            else if (s == st_6R)
              { return (1-pD); }
            else
              { return 0; }
          }
        case st_XX:
          /* Possible transitions:
             -> XX
          */
          { if (s == st_XX) 
              { return 1; }
            else
              { return 0; }
          }
        default:
          demand(FALSE, "invalid current state");
      }
    assert(FALSE);
  }

double msd_prob_D(double a)
  {
    return msd_sigmoid(a, 0, 62, 10, 72);
  }
  
double msd_prob_E(double a)
  {
    return msd_sigmoid(a, 22, 37, 3, 42);
  }
  
double msd_prob_R(double a)
  {
    return msd_sigmoid(a, 27, 42, 5, 52);
  }
  
double msd_prob_L(double a)
  {
    double aref = 37.0; /* Academic age of reference. */
    double rref = 4.00; /* Determines the (MS5+MS6)/MS3 ratio at age {aref}. */
    /* A curve {s} with a hump at 10 years and a plateau afterwards: */
    double s1 = msd_humpoid(a, 3, 10, 5, 72);
    double s2 = msd_sigmoid(a, 3, 10, 5, 72);
    double s = 0.67*s1 + 0.33*s2;
    return (rref/aref)*s;
  }
  
double msd_prob_T(double a)
  {
    double aref = 37.0; /* Academic age of reference. */
    double rref = 2.50; /* Determines the MS6/MS5 ratio at age {aref}. */
    /* A curve {s} with a hump at 15 years and a plateau afterwards: */
    double s1 = msd_humpoid(a, 6, 15, 5, 72);
    double s2 = msd_sigmoid(a, 6, 15, 5, 72);
    double s = 0.50*s1 + 0.50*s2;
    return (rref/aref)*s;
  }
  
double msd_sigmoid(double x, double xMin, double xAvg, double xDev, double xMax)
  {
    if (x < xMin) { return 0; }
    if (x > xMax) { return 1; }
    double s = msd_erf(x, xAvg, xDev);
    double sMin = msd_erf(xMin, xAvg, xDev);
    double sMax = msd_erf(xMax, xAvg, xDev);
    return (s - sMin)/(sMax - sMin);
  }

double msd_erf(double x, double xAvg, double xDev)
  {
    return 0.5*(1 + erf((x - xAvg)/(M_SQRT2*xDev)));
  }

double msd_humpoid(double x, double xMin, double xAvg, double xDev, double xMax)
  {
    if (x < xMin) { return 0; }
    if (x > xMax) { return 0; }
    /* Get a gaussian {g} and its values {gMin,gMax} at {xMin,xMax}: */
    double g = msd_gauss(x, xAvg, xDev);
    double gMin = msd_gauss(xMin, xAvg, xDev);
    double gMax = msd_gauss(xMax, xAvg, xDev);
    /* Get a erfc {e} and its values {eMin,eMax} at {xMin,xMax}: */
    double s = msd_erf(x, xAvg, xDev);
    double sMin = msd_erf(xMin, xAvg, xDev);
    double sMax = msd_erf(xMax, xAvg, xDev);
    /* Subtract from {s} the sigmoid {e} shifted so as to cancel at {xMin,xMax}: */
    double d = gMin + (gMax - gMin)*(s - sMin)/(sMax - sMin);
    return g - d;
  }

double msd_gauss(double x, double xAvg, double xDev)
  { 
    double z = (x - xAvg)/xDev;
    return exp(-z*z/2);
  }

PSStream *msd_init_plot
  ( int aMax, 
    double pMax, 
    double hSize, 
    double vSize, 
    double *hScale, 
    double *vScale
  )
  {
    double mm = epswr_pt_per_mm;  /* One mm in Postscript points. */
    
    /* Lower bounds for meaningful client coordinates: */
    int aMin = 0;       /* years since hiring. */
    double pMin = 0.0;  /* Relative population size. */
    
    /* Client-to-plot scale factors: */
    *hScale = hSize/(aMax - aMin);  /* Hor scale (mm per year). */
    *vScale = vSize/(pMax - pMin);  /* Ver scale (mm per unit of population). */

    /* Extra margin for axis arrows, tics, etc.: */
    double hMrg =  5.0;  /* Hor margin (mm). */
    double vTop =  5.0;  /* Top margin (mm). */
    double vBot = 10.0;  /* Bot margin (mm). */
    
    /* Client plot ranges (mm): */
    double xMin = aMin*(*hScale) - hMrg, xMax = aMax*(*hScale) + hMrg;
    double yMin = pMin*(*vScale) - vBot, yMax = pMax*(*vScale) + vTop;
    
    /* Overall figure dimensions (pt): */
    double fMrg = 4.0; /* Extra margin around figure. */
    double hPageSize = (xMax - xMin)*mm + 2*fMrg;
    double vPageSize = (yMax - yMin)*mm + 2*fMrg;
    
    /* Plot window dimensions (pt): */
    double hMin = fMrg, hMax = hPageSize - fMrg;
    double vMin = fMrg, vMax = vPageSize - fMrg;
    
    /* Start plot stream: */
    PSStream *ps = pswr_new_stream("out/msdemo", NULL, TRUE, "doc", "letter", FALSE, hPageSize, vPageSize);
    pswr_new_canvas(ps, "1");
    pswr_set_window
      ( ps,
        xMin, xMax, yMin, yMax,
        hMin, hMax, vMin, vMax
      );
    pswr_set_pen(ps, 0,0,0, 0.25, 0,0);
    pswr_set_fill_color(ps, 0,0,0);
    pswr_set_label_font(ps, "TimesRoman", 10.0);
    
    pswr_axis(ps, HOR, 0,  xMin, xMax);
    pswr_arrowhead(ps, xMin,0, xMax,0, 1.5,3.0, 1.0, TRUE, TRUE);
    int a;
    for (a = 0; a < aMax; a += 5)
      { double x = a*(*hScale);
        pswr_tic(ps, HOR, x,0, 1.5, 0.5);
        char *lab = jsprintf("%d", a);
        pswr_label(ps, lab, x,-2.0, 0.0, 0.5,1.0);
        free(lab);
      }

    pswr_axis(ps, VER, 0,  yMin, yMax); 
    pswr_arrowhead(ps, 0,yMin, 0,yMax, 1.5,3.0, 1.0, TRUE, TRUE);
    
    return ps;
  }

options_t *parse_options(int argc, char **argv)
  {
    /* Initialize argument parser: */
    argparser_t *pp = argparser_new(stderr, argc, argv);
    argparser_set_help(pp, PROG_NAME " version " PROG_VERS ", usage:\n" PROG_HELP);
    argparser_set_info(pp, PROG_INFO);
    argparser_process_help_info_options(pp);
    
    /* Allocate the command line argument record: */
    options_t *o = (options_t *)malloc(sizeof(options_t)); 
    
    /* Parse keyword parameters: */
    
    /* Parse the boolean option {op1}: */
    o->op1 = argparser_keyword_present(pp, "-op1");
    
    /* Parse the string option {op2}: */
    if (argparser_keyword_present(pp, "-op2"))
      { o->op2 = argparser_get_next(pp); }
    else
      { o->op2 = "NONE"; }

    /* Parse positional arguments: */
    argparser_skip_parsed(pp);
    
    /* Check for spurious arguments: */
    argparser_finish(pp);
    
    return o;
  }