/* Filter a curve with "a-posteriori" parametrization. */
/* Last edited on 2023-10-15 03:33:24 by stolfi */

#define PROG_HELP \
  "    " PROG_NAME " \\\n" \
  "  -input {IN_NAME} \\\n" \
  "  -output {OUT_NAME} \\\n" \
  "  -lambdaMin {L_MIN} \\\n" \
  "  [ -nScales {N_SCALES} ] \\\n" \
  "  [ -nRepStages {NREPS} ] \\\n" \
  "  [ -lambdaRatio {L_RATIO}] \\\n" \
  "  { -parametric | \\\n" \
  "    -tolLength {TOL_LEN} -tolStep {TOL_STEP} [ -maxIter {MAX_IT} ] \\\n" \
  "  } \\\n" \
  "  [ -maxStep {MAX_STEP} ] \\\n" \
  "  [ -writeRaw ] \\\n" \
  "  [ -writeStep {WR_STEP} ] [ -outUnit {UNIT} ] \\\n" \
  "  [ -verbose ] \n"

#define PROG_DESC \
  "  This program takes an arbitrary curve {c} and outputs several " \
  "versions {d[0],d[1],...} of it, geometrically filtered with Gaussian " \
  "filters of scale {lambda}, {2*lambda}, {4*lambda}, etc.\n" \
  "\n" \
  "In signal-processing contexts, a Gaussian filter with " \
  "scale {sigma} is a filter whose kernel is a Gaussian distribution " \
  "having mean zero and standard deviation {sigma}, i.e.\n" \
  "\n" \
  "  G_sigma(t) ==  sigma/sqrt(2Pi) exp(-t^2/(2*sigma^2))\n" \
  "\n" \
  "`Geometric filtering' a curve {c} with a filter {F} means applying " \
  "{F} to {c}, after reparametrizing {c} in such a way that the filtered " \
  "curve has unit velocity.\n" \
 "\n" \ \
  "Actually, each output curve {d[i]} after {d[0]} is obtained by " \
  "filtering the preceding curve {d[i-1]} with a Gaussian filter whose " \
  "scale is computed so that the result will be roughly equivalent to " \
  "filtering {c} to scale {2^i*lambda}.  The purpose of this " \
  "optimization is to reduce the amount of work from {O(n log n)} to " \
  "{O(n)}, where {n} is the number of samples in curve {c[0]}."
  
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <values.h>
#include <assert.h>

#include <argparser.h>
#include <jsstring.h>
#include <jsmath.h>

#include <pz_types.h>
#include <vsignal.h>
#include <pz_r3_chain.h>
#include <pz_double_chain.h>

typedef
  struct options_t {
    char *input;         /* Input curve file (without ".flc"). */
    char *output;        /* Output curve file name (without ".flc"). */
    double lambdaMin;    /* Finest filtering scale length. */
    double lambdaRatio;  /* Ratio between successive scales. */
    unsigned nScales;    /* Number of scales (default +oo). */
    double tolLength;    /* Max change in length (rel. lambda) per iteration. */
    double tolStep;      /* Max change in interval size (rel lambda) per iter. */
    unsigned maxIter;    /* Maximum number of iterations (0 for parametric). */
    unsigned nRepStages; /* Number of reparametrization stages. */
    bool_t writeRaw;       /* TRUE to write raw approximation to input contour */
    unsigned writeStep;  /* Write the curve every this many filtering scales. */
    double maxStep;      /* Maximum sampling step for output (default +oo). */
    double outUnit;      /* Output scale factor. */
    bool_t verbose;        /* TRUE to make noises while working */
  } options_t;

pz_curve_t  pz_curve_geom_filter
  ( pz_curve_t *c,
    double lambdaOld,
    double lambdaNew,
    unsigned nRepStages,
    double tolLength,
    double tolStep,
    unsigned maxIter,
    unsigned *totIter, /* Total iterations performed for all stages */
    bool_t verbose
  );
 /* Takes a curve {c} that has already been geometrically filtered to
    scale {lambdaOld} and returns a copy {d} of {c} processed with a
    geometric Gaussian filter of scale {lambdaNew > lambdaOld},
    sampled at {n} uniformly spaced points, where {n} is chosen
    according to the Nyquist criterion.

    Upon entry, assumes {c} is parametrized by arclength of the
    {lambdaOld}-filtered version.  Returns the {lambdaNew}-filtered
    curve {d} with the natural parametrization, and reparametrizes {c}
    with the same parametrization.
  */

double pz_filter_step
  ( pz_curve_t *cOrg,     /* Original curve, possibly reparametrized. */
    double lambda,        /* Standard deviation of convolution kernel. */
    pz_curve_t *cFil,     /* Filtered curve (preallocated). */
    unsigned minSamples,  /* Number of samples by Nyquist's limit */
    bool_t verbose,         /* TRUE to mumble while working. */
    /* WORK*/
    double_vec_t *ra,
    double_vec_t *ta,
    double_vec_t *tb      /* Working vector */
  );
  /* Convolves the curve {cOrg} with a Gaussian of unit area and
    characteristic duration {lambda}, i.e. {exp(-(t/lambda)^2/2)/A}
    where {A} is the unit-area normalization term.  Samples the
    filtered curve at about {minSamples} equally spaced times, and
    stores the result in {cFil}.

    Then reparametrizes {cOrg} so that the time lapse between two samples
    of {cOrg} is the length of {cFil} between the points corresponding
    to those samples. (This step is omitted if the filtered curve
    {cFil} has length less than lambda.)

    Reallocates {cFil} if it hasn't been initialized, or has too few or
    too many points.  Returns the largest (in absolute value) signed
    change in the spacing between consecutive sampling times of {cOrg}. */
  

void pz_debug(char *msg, double_vec_t *x);
  /* Prints {msg} and {x} to {stderr}. */

double pz_reparametrize
  ( pz_curve_t *cOrg,
    pz_curve_t *cFil,
    double minLength,
    bool_t verbose,
     /* (WORK) */
    double_vec_t *ra,
    double_vec_t *ta,
    double_vec_t *tb
  );
  /* Given a curve {cOrg}, and a parametrically filtered and resampled
    version {cFil} of {cOrg}, this procedure reparametrizes both by
    the arc length of {cFil}.

    More precisely, it sets the time {cFil.t[j]} to be the length of
    the curve {cFil(t)} from sample {cFil.p[0]} to sample {cFil.p[j]}.

    Then it sets the time {cOrg.t[i]} associated with each sample
    {cOrg.p[i]} to the value {t} such that {cFil(t)} is the point of
    {cFil} that corresponds to {cOrg.p[i]}.

    The routine then computes the time intervals between the
    consecutive samples of {cOrg}, and compares those intervals with
    the corresponding intervals of old times.  The returned value is
    the maximum absolute change in those intervals.

    As a special case, if the length of {cFil} is less than
    {minLength}, the curve {cOrg} is not reparametrized, and the
    procedure returns {minLength}. */

typedef struct pz_curve_data_t
  { char *cmt;
    pz_curve_t curve;
  } pz_curve_data_t;

pz_curve_data_t pz_curve_read(char *fileName, bool_t verbose);
  /* Reads a {pz_r3_chain_t} from a file named "{fileName}",
    using {pz_r3_chain_read}, and returns it packaged as
    a {pz_curve_t}.  Returns also the comments read from the 
    file. */

void pz_curve_write_raw
  ( options_t *o,
    char *cmt,
    pz_curve_t *c,
    double unit
  );
  /* Writes curve {c} to files "{name}000.lbl", "{name}000.flc",
    "{name}000.lambda" (as per {pz_curve_write}). Here {name} is
    {o->output}. */

void pz_curve_write_sampled
  ( options_t *o,
    char *cmt,
    pz_curve_t *c,
    double lambda,
    double lambdaMin,
    unsigned totIter,
    double unit
  );
  /* Resamples curve {c} and writes the result to files
    "{name}{NNN}.lbl", "{name}{NNN}.flc", "{name}{NNN}.lambda" (as per
    {pz_curve_write}).  Here {name} is {o->output}, and
    {NNN} is {lambda/lambdaMin} rounded to integer and 
    zero-padded to 3 digits. */

pz_curve_t pz_curve_sample(pz_curve_t c, unsigned n);
  /* Resamples a curve at {n} points equally spaced in time. */

void pz_curve_write
  ( char *name,
    char *cmt,
    pz_curve_t *c,
    double lambda,
    double unit
  );
  /* Writes the curve {c} to files called "{name}.lbl" (labels)
    and "{label}.flc" (points).  Also writes a file called 
    "{name}.lambda" with the value of {lambda}. */

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

int main(int argc, char **argv)
  { pz_curve_t cOld, cNew;
    double lambdaOld, lambdaNew;
    unsigned totIter;
    unsigned  iScale = 0;
    double MaxRelativeLambda = 0.45; /* Stop when lambda/L is grater than this. */

    options_t *o = pz_get_options(argc, argv);
    char *inputFileName = txtcat(o->input , ".flc");
    pz_curve_data_t rp = pz_curve_read(inputFileName, o->verbose);
    pz_curve_t c = rp.curve;
    char *cCmt = txtcat4(rp.cmt, "pz_filter:\n  input file = ", inputFileName, "\n");
    
    if (o->writeRaw) { pz_curve_write_raw(o, cCmt, &c, o->outUnit); }

    cOld = c;
    lambdaOld = 0.0;
    lambdaNew = o->lambdaMin;
    iScale = 0;
    /* Filter for successive scales: */
    while (1)
      { if ((lambdaNew > MaxRelativeLambda * cOld.tPeriod) || (iScale >= o->nScales)) { break; }
        cNew = pz_curve_geom_filter
          ( &cOld, lambdaOld, lambdaNew,
            o->nRepStages, o->tolLength, o->tolStep, o->maxIter, &totIter, o->verbose
          );
        if (cNew.tPeriod <= lambdaNew) { break; }
        if ((iScale == o->nScales) || (iScale % o->writeStep == 0))
          { pz_curve_write_sampled(o, cCmt, &cNew, lambdaNew, o->lambdaMin, totIter, o->outUnit); }
        cOld = cNew;
        lambdaOld = lambdaNew;
        lambdaNew = lambdaNew * o->lambdaRatio;
        iScale++;
      }
    fprintf(stderr, "maximum lambda = %8.5f\n", lambdaOld);
  }

pz_curve_t  pz_curve_geom_filter
  ( pz_curve_t *c,
    double lambdaOld,
    double lambdaNew,
    unsigned nRepStages,
    double tolLength,
    double tolStep,
    unsigned maxIter,
    unsigned *totIter, /* Total iterations performed for all stages */
    bool_t verbose
  )
  { int MaxIntermStageIter = 3; /* Max iterations in intermediate stages. */
    pzcurve_t d; 
    double change;
    unsigned stageIter;       /* Iteration count per stage. */
    unsigned maxStageIter;    /* Max iterations in current stage. */
    double maxLength;         /* Max length of {lambdaNew}-filtered curve. */
    double_vec_t ra = double_vec_new(0);  /* Work vectors for FilterStep. */
    double_vec_t ta = double_vec_new(0);  /* Work vectors for FilterStep. */
    double_vec_t tb = double_vec_new(0);  /* Work vectors for FilterStep. */

    maxLength = c->tPeriod;
    totIter = 0;
    int k;
    for (k = 1; k <= nRepStages ; k++)
      {
        stageIter = 0;
        if (k < nRepStages)
          { maxStageIter = MaxIntermStageIter; }
        else
          { maxStageIter = UINT_MAX; }

        double exp = ((double)nRepStages - k)/imax(1, nRepStages - 1);
        double lambdaStg = lambdaNew * pow(2.0, exp);
        int minSamples = (int)ceil(4.0 * maxLength/lambdaStg);
        /* double exp = ((double)k)/((double)nRepStages); */
        /* double lambdaStg = lambdaOld * pow(lambdaNew/lambdaOld, exp); */
        while (1)
          { (*totIter)++;
            stageIter++;
            if (verbose) { fprintf(stderr, "\n"); }

            { /* Compute the parameter {lambdaInc} of a convolution
                that applied to {c} has the same effect as applying
                convolution {lambdaStg} to the original curve {p}: */
              double lambdaInc = sqrt(lambdaStg*lambdaStg - lambdaOld*lambdaOld);
              double TOld = c->tPeriod + 0.0;
              if (verbose)
                { fprintf(stderr, "iteration %d", (*totIter));
                  fprintf(stderr, "lambdaStg = %8.5f", lambdaStg);
                  fprintf(stderr, "lambdaInc = %8.5f\n", lambdaInc);
                }
              change = pz_filter_step(c, lambdaInc, d, minSamples, verbose, ra,ta,tb);
              affirm((d.tPeriod < lambdaInc) || (c->tPeriod == d.tPeriod), "period bug");
              if ((d.tPeriod <= lambdaInc) || (d.m < 2))
                { fprintf(stderr, "! curve collapsed, stage abandoned !\n"); break; }
              if ((fabs(change) <= tolStep*lambdaStg) && (fabs(TOld - d.tPeriod) <= tolLength*lambdaStg))
                { fprintf(stderr, "! stage converged !\n"); break; }
              if ((*totIter) >= maxIter)
                { if (maxIter > 1)
                    { fprintf(stderr, "! too many iterations, stopped !\n"); }
                  break;
                }
              if (stageIter >= maxStageIter)
                { fprintf(stderr, "! too many iterations, interm stage stopped !\n"); break; }
            }
        }
        
      }
    return d;
  }

double pz_filter_step
  ( pz_curve_t *cOrg,     /* Original curve, possibly reparametrized. */
    double lambda,        /* Standard deviation of convolution kernel. */
    pz_curve_t *cFil,     /* Filtered curve (preallocated). */
    unsigned minSamples,  /* Number of samples by Nyquist's limit */
    bool_t verbose,         /* TRUE to mumble while working. */
    /* WORK*/
    double_vec_t *ra,
    double_vec_t *ta,
    double_vec_t *tb      /* Working vector */
  )
  {
    double change;
    if ((cFil == NULL) || (cFil->m < minSamples) || (cFil->m > (12*minSamples) / 10))
      { cFil = pz_curve_new(minSamples);
        if (verbose)
          { fprintf(stderr, "resizing to %d samples\n", cFil->m); }
      }
    if (verbose)
      { fprintf(stderr, "period = %8.5f", cOrg->tPeriod );
        fprintf(stderr, "filtering...\n");
      }
    pz_curve_setToFiltered(cFil, cOrg, lambda);
    assert(cFil->tPeriod + 0.0 == cFil->tPeriod );  /* To catch NaN */
    change = pz_reparametrize(cOrg, cFil, lambda, verbose, ra, ta, tb);
    affirm((cFil->tPeriod < lambda) || (cOrg->tPeriod == cFil->tPeriod), "period bug");
    return change;
  }

void pz_debug(char *msg, double_vec_t *x)
  {
    fprintf(stderr, msg);
    int i;
    for (i = 0; i < x->nel; i++)
      { if ((i > 0) && (i % 5 == 0))
          { fprintf(stderr, "\n"); }
        else
          { fprintf(stderr, " "); }
        fprintf(stderr, "%12.8f", x->el[i]);
      }
    fprintf(stderr, "\n");
  }

double pz_reparametrize
  ( pz_curve_t *cOrg,
    pz_curve_t *cFil,
    double minLength,
    bool_t verbose,
     /* (WORK) */
    double_vec_t *ra,
    double_vec_t *ta,
    double_vec_t *tb
  )
  {
    double maxChange, minStep, maxStep;
    double tao, tan;
    unsigned maxChangeIndex;
    assert(cFil->tPeriod == cOrg->tPeriod );
    assert(cFil->tPeriod + 0.0 == cFil->tPeriod ); /* To catch NaNs*/
    assert(cOrg->tPeriod + 0.0 == cOrg->tPeriod ); /* To catch NaNs*/
    double_vec_expand(ra, cOrg->m-1);
    double_vec_expand(ta, cFil->m-1);
    double_vec_expand(tb, cOrg->m-1);

    double_vec_t tOrg = &(cOrg->t);
    double_vec_t tFil = &(cFil->t);
    double TOld = cOrg->tPeriod;
    double_vec_t tFilOld = SUBARRAY(ta, 0, cFil->m); /* Original times of filtered samples. */
    double_vec_t tOrgNew = SUBARRAY(tb, 0, cOrg->m); /* New times of original samples. */
    double_vec_t arg = SUBARRAY(ra, 0, cOrg->m); /* Work array for reparametrization. */
    /* Save the original times of filtered samples: */
    tFilOld = tFil;

    /* Set the times of filtered samples to be length since {cFil.p[0]}: */
    pz_curve_uniformize(cFil);

    double TNew = cFil->tPeriod;

    if (verbose)
      { fprintf(stderr, "filtered curve length = %8.5f", TNew);
        fprintf(stderr, "change = %8.5f", TNew-TOld);
      }
    /* Apply the reparametrization {tFilOld -> tFil} to {tOrg} too: */
    if ((cOrg->m < 2) || (TNew < minLength))
      { if (verbose)
          { fprintf(stderr, "filtered curve collapsed, ");
            fprintf(stderr, "length = %8.5f", TNew);
            fprintf(stderr, "- original times not changed \n");
          }
        return minLength;
      }
    else
      { /* Compute {arg[i]} == fractional index of {tOrgOld[i]} in {tFil}: */
        pz_double_chain_args_from_values(&tFilOld, TOld, tOrg, &arg);
        assert(arg[0] + 0.0 == arg[0]);
        pz_double_chain_values_from_args(tFil, TNew, arg, tOrgNew);
      }
    cOrg.tPeriod = TNew;

    /* Compute variation in time intervals from {tOrg} to {tOrgNew}: */
    maxChange = 0.0; maxChangeIndex = 0;
    minStep = +INF;
    maxStep = -INF;
    tao = tOrg[cOrg.m-1] - TOld;
    tan = tOrgNew[cOrg.m-1] - TNew;
    for (i = 0; i < cOrg.m; i++)
      { double tbo = tOrg[i];
        double tbn = tOrgNew[i];
        double stepSize = tbn - tan;
        double stepChange = (tbn-tan)-(tbo-tao);
        assert(tan <= tbn );
        assert(tao <= tbo );
        if (stepSize > maxStep){ maxStep = stepSize; }
        if (stepSize < minStep){ minStep = stepSize; }
        if (fabs(stepChange) > fabs(maxChange))
          { maxChange = stepChange;  maxChangeIndex = i; }
        tao = tbo; tan = tbn;
      }
    if (verbose)
      { fprintf(stderr, "new time steps: min = %8.5f", minStep);
        fprintf(stderr, "max = %8.5f", maxStep);
        fprintf(stderr, "max change = %8.5f", maxChange);
        fprintf(stderr, "at sample %d\n", maxChangeIndex);
      }

    /* Now set the times in {cOrg} appropriately: */
    cOrg.setTimes(tOrgNew, TNew);
    assert(cOrg.tPeriod + 0.0 == cOrg.tPeriod );  /* To catch NaN */;
    
    return maxChange;
  }

pz_curve_data_t pz_curve_read(char *fileName, bool_t verbose)
  {
    fprintf(stderr, "reading file %s ...\n", fileName );

    FILE *rd = open_read(fileName, TRUE);
    pz_r3_chain_t chd = pz_r3_chain_read(rd, /*header_only:*/ FALSE, /*recenter:*/ pz_ctr_SMPS);
    r3_vec_t *p = &(chd.c);
    int m = p->nel;
    pz_curve_t c = pz_curve_new(m);
    double_vec_t s = double_vec_new(m);
    pz_curve_setSamples(&c, p);
    /* Initial times are lengths along the polygonal: */
    double L = pz_r3_chain_linear_lengths(p, s);
    if (verbose)
      { fprintf(stderr, "polygonal length = %8.5f\n", L); }
    pz_curve_setTimes(&c, s, L);
    /* Recompute lengths along interpolated curve: */
    pz_curve_uniformize(&c);
    if (verbose) { fprintf(stderr, "adjusted curve length = %8.5f\n", c->tPeriod); }
    /* Use the natural parameter values of the original curve as labels: */
    pz_curve_setLabels(c, c->t, c->tPeriod);
    return (pz_curve_data_t){cmt = chd.cmt, curve = c};
  }

void pz_curve_write_raw
  ( options_t *o,
    char *cmt,
    pz_curve_t *c,
    double unit
  )
  {
    pz_curve_write(o->output & "000", cmt, c, 0.0, unit);
  }

void pz_curve_write_sampled
  ( options_t *o,
    char *cmt,
    pz_curve_t *c,
    double lambda,
    double lambdaMin,
    unsigned totIter,
    double unit
  )
  {

    char *FmtMaxStep(double maxStep)
      { if (maxStep != +INF)
          { return " maxStep = %8.5f", maxStep }
        else
          { return ""; }
      }

    char *name = o->output & Fmt.Pad(Fmt.Int(ROUND(lambda/lambdaMin)), 3, '0');
    double L = c->tPeriod;
    double step = MIN(o->maxStep, lambda/4.0);
    int n = CEILING(L/step);
    pz_curve_write(
      name,
      cmt
      & "  filtered to scale lambda = %8.5f", lambda & " (pixels)\n"
      & "  nRepStages = %d", o->nRepStages
      & "  tolLength = %8.5f", o->tolLength
      & "  tolStep = %8.5f", o->tolStep
      & "  maxIter = %d", o->maxIter & "\n"
      & "  sampled at %d", n & " equally spaced points"
      & FmtMaxStep(o->maxStep) & "\n"
      & "  totIter = %d", totIter & ")\n",
      pz_curve_sample(c, n),
      lambda,
      unit
    );

  }

pz_curve_t pz_curve_sample(pz_curve_t c, unsigned n)
  {
    double tau = pz_double_chain_new(n)^;
    double p = pz_r3_chain_new(n)^;
    double v = pz_r3_chain_new(n)^;
    double r = pz_double_chain_new(n)^;
    double T = c->tPeriod;
    double R = c->rPeriod;
    pz_curve_t d = pz_curve_new(n);
    for (i = 0; i <= n-1 ; i++)
      { tau[i] = ((double)i)/((double)n) * T; }
    c->sample(tau, p, v, r);
    d.setSamples(p);
    d.setTimes(tau, T);
    d.setLabels(r, R);
    return d;
  }

void pz_curve_write
  ( char *name,
    char *cmt,
    pz_curve_t *c,
    double lambda,
    double unit
  )
  {
    char *fileName;
    
    fileName = txtcat(name , ".lbl");
    fprintf(stderr, "writing file %s ...\n", fileName );

    double adjUnit = pz_double_chain_adjust_unit(unit, c->r^);
    pz_double_chain_write(open_write(fileName, TRUE), cmt, c->r^, c->rPeriod, unit = adjUnit);

    fileName = txtcat(name , ".flc");
    fprintf(stderr, "writing file %s ...\n", fileName );

    double adjUnit = pz_r3_chain_adjust_unit(unit, c->p^);
    pz_r3_chain_write(open_write(fileName, TRUE), cmt, c->p^, unit = adjUnit);

    fileName = txtcat(name , ".lambda");
    fprintf(stderr, "writing file %s ...\n", fileName );

    FILE *wr = open_write(fileName, TRUE);
    fprintf(wr, "%12.6f", lambda & "\n");
    fclose(wr);
  }

options_t *pz_get_options(int argc, char **argv)
  {
    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);
    
    options_t *o = (options_t *)notnull(malloc(sizeof(options_t)), "out of mem");

    argparser_get_keyword(pp, "-input");
    o->input = argparser_get_next(pp);

    argparser_get_keyword(pp, "-output");
    o->output = argparser_get_next(pp);

    argparser_get_keyword(pp, "-lambdaMin");
    o->lambdaMin = argparser_get_next_double(pp);

    if (argparser_keyword_present(pp, "-nScales"))
      { o->nScales = argparser_get_next_int(pp, 0, (unsigned.ne - 1)); }
    else
      { o->nScales = (unsigned.ne - 1); }

    if (argparser_keyword_present(pp, "-parametric"))
      { o->maxIter = 0;
        o->tolLength = 0.0;
        o->tolStep = 0.0
      }
    else
      { argparser_get_keyword(pp, "-tolLength");
        o->tolLength = argparser_get_next_double(pp, 0.0);
        argparser_get_keyword(pp, "-tolStep");
        o->tolStep = argparser_get_next_double(pp, 0.0);

        if (argparser_keyword_present(pp, "-maxIter"))
          { o->maxIter = argparser_get_next_int(pp, 0, (unsigned.ne - 1)); }
        else
          { o->maxIter = (unsigned.ne - 1); }
      }

    if (( argparser_keyword_present(pp, "-nRepStages") )){
      o->nRepStages = argparser_get_next_int(pp, 1, (unsigned.ne - 1))
    }else{
      o->nRepStages = 1;
    }

    if (argparser_keyword_present(pp, "-maxStep"))
      { o->maxStep = argparser_get_next_double(pp, 1.0e-6, (double.ne - 1)); }
    else
      { o->maxStep = (double.ne - 1); }

    if (argparser_keyword_present(pp, "-lambdaRatio"))
      { o->lambdaRatio =  argparser_get_next_double(pp); }
    else
      { o->lambdaRatio = 2.0; }

    o->writeRaw = argparser_keyword_present(pp, "-writeRaw");

    if (( argparser_keyword_present(pp, "-writeStep") )){
      o->writeStep = argparser_get_next_int(pp, 1, (unsigned.ne - 1))
    }else{
      o->writeStep = 1;
    }

    if (argparser_keyword_present(pp, "-outUnit"))
      { o->outUnit = argparser_get_next_double(pp, 0.1e-10, 0.1e10); }
    else
      { o->outUnit = 0.001; }

    o->verbose = argparser_keyword_present(pp, "-verbose");

    argparser_finish(pp);

    return o;
  }