MODULE BITExtractSamples EXPORTS Main;
(* Last edited on 1999-12-22 14:27:01 by hcgl *)

(* Extracts matching segments of given chordlength from a set of candidates. *)

(*
  Refines a list of candidates (pairs of outline segments), using DP
  to find the optimum matching, and NL optimization to find the best
  relative rotation and displacement.

  From each refined candidate pair, the program extracts zero or more
  matching segments with specified number of steps.
*)

IMPORT ParseParams, FileRd, FileWr, Text, TextWr, FS, Wr, Fmt;
IMPORT Process, OSError, Thread, Stdio;
IMPORT
  PZGeo, PZShape, PZLR3Chain, PZLRChain,
  PZCandidate, PZMismatch, PZMatch, PZOptMatrix;

IMPORT LR3, LR4x4;

FROM Math IMPORT cos, sin, sqrt;
FROM LR3 IMPORT Dist;
FROM LR4x4 IMPORT Mul, Identity;
FROM PZShape IMPORT ComputeSpectrum, ComputePowerSpectrum;
FROM Minimizer IMPORT Value;
FROM PZTypes IMPORT LONG, LONGS, INT, NAT, BOOL;

FROM Stdio IMPORT stderr;

<* FATAL Wr.Failure, Thread.Alerted, OSError.E *>

TYPE
  MismatchOptions = PZMismatch.Options;
  Options = RECORD
      (* Data identification parameters: *)
      input: TEXT;              (* Input candidate file name. *)
      output: TEXT;             (* Output file name prefix (usually same as input) *)
      chainDir: TEXT;           (* Directory where to find the chains *)
      chainPrefix: TEXT;        (* Chain file name prefix. *)
      band: NAT;                (* Nominal band width (lambda) for file names. *)
      extension: TEXT;          (* Extension of chain files, usually ".flc" *)

      (* Matching parameters: *)
      maxShift: LONG;           (* Max shift to try in alignment. *)
      mp: MismatchOptions;      (* Mismatch formula parameters. *)
      minLength: LONG;          (* Minimum length of matching segments. *)
      blurFactor: LONG;         (* Corner broadening factor. *)
      extraLength: LONG;        (* Extra length to add at each candidate end. *)

      (* Placement/optimization parameters: *)
      dontGuess : BOOL;         (* TRUE to use identity as the initial matrix  *)
      perturbMatrix: BOOL;      (* TRUE to perturb placement matrix *)
      maxEvals: NAT;            (* Max trials in evaluation (0 = dont' optimize). *)

      (* Segment extraction parameters: *)
      maxCands: NAT;            (* Number of candidates to generate *)
      nSteps: NAT;              (* Desired number of output sampling steps. *)
    END;

TYPE
  ChainData = PZLR3Chain.ReadData;

CONST
  ChainUnit = 0.1d0;  (* Coordinate quantization unit for chain files. *)
  ShapeUnit = 0.1d0;  (* Quantization unit for shape files. *)

PROCEDURE Main() =
  BEGIN
    WITH
      o = GetOptions(),
      candData = PZCandidate.Read(FileRd.Open(o.input & ".can")),
      cand = candData.c^,
      comment = candData.cmt,
      lambda = candData.lambda,
      needed = PZCandidate.GetCurves(cand)^,
      chData = PZLR3Chain.ReadAll(
        o.chainPrefix, o.band, o.extension, sel := needed, 
        dir := o.chainDir, headerOnly := FALSE
      ),
      ch = chData.chain^
    DO
      ProcessCandidates(cand, ch, o, lambda, comment)
    END;
  END Main;

PROCEDURE ProcessCandidates(
    READONLY cand: PZCandidate.List;
    READONLY ch: ARRAY OF ChainData;
    READONLY o: Options;
    lambda: LONG;
    comment: TEXT; 
  ) =
  VAR
    minAvgDist: REF LONGS := NIL;
    rCost: REF PZMatch.CostMatrix := NIL;
    nOutCands: NAT := 0;
    iCand: NAT := 0;
  BEGIN
    TRY FS.CreateDirectory(o.output) EXCEPT OSError.E => (*OK*) END;
    WITH
      nCands = NUMBER(cand),
      nCurves = NUMBER(ch),
      length = NEW(REF LONGS, nCurves)^
    DO
      FOR k := 0 TO nCurves-1 DO
        IF ch[k].c # NIL THEN length[k] := ChainLength(ch[k].c^) END
      END;
      WHILE iCand < nCands AND nOutCands < o.maxCands DO
        ProcessOneCandidate(
          iCand, cand[iCand], ch, length, o, lambda, comment,
          (*IO*)
          nOutCands := nOutCands,
          (*WORK*)
          minAvgDist := minAvgDist,
          rCost := rCost
        );
        INC(iCand)
      END;
    END
  END ProcessCandidates;

PROCEDURE ChainLength(READONLY c: PZLR3Chain.T): LONG =
  VAR s: LONG := 0.0d0;
  VAR j := LAST(c);
  BEGIN
    FOR i := 0 TO LAST(c) DO
      s := s + LR3.Dist(c[i], c[j]);
      j := i
    END;
    RETURN s
  END ChainLength;

PROCEDURE ProcessOneCandidate(
    iCand: NAT;
    READONLY cand: PZCandidate.T;
    READONLY ch: ARRAY OF ChainData;
    READONLY chLength: LONGS;
    READONLY o: Options;
    lambda: LONG;
    comment: TEXT;
    (*IO*)
    VAR nOutCands: NAT;                    (* Num candidates written so far *)
    (*WORK*)
    VAR minAvgDist: REF LONGS;
    VAR rCost: REF PZMatch.CostMatrix;
  ) =
  (*
    Let "a,b" be the two curves that constitute candidate "cand".
    This procedure rotates and translates curve "a" (through
    a 4×4 matrix "m") so as to optimize the match with curve "b".
    Then it extracts matching segments of prescribed length
    from the two candidates.
  *)
  VAR
    mt: REF PZMatch.T;
    mismatch: LONG;
    candLength: LONG;
    matchedLength: LONG;
    m: LR4x4.T;
    cam:  REF PZLR3Chain.T;
  BEGIN
    WITH
      aSeg = cand.seg[0],
      cha = ch[aSeg.cvx],
      ma = NUMBER(cha.c^),
      rna = aSeg.ns,
      aStep = chLength[aSeg.cvx] / FLOAT(ma, LONG),
      aExtra = CEILING(o.extraLength / aStep),
      na = MIN(rna + 2*aExtra, ma),
      aStretch = (na - rna) DIV 2,
      ia = (aSeg.ini - aStretch) MOD ma,
      ca = PZLR3Chain.Trim(cha.c^, ia, na)^,

      bSeg = cand.seg[1],
      chb = ch[bSeg.cvx],
      mb = NUMBER(chb.c^),
      rnb = bSeg.ns,
      bStep = chLength[bSeg.cvx] / FLOAT(mb, LONG),
      bExtra = CEILING(o.extraLength / bStep),
      nb = MIN(rnb + 2*bExtra, mb),
      bStretch = (nb - rnb) DIV 2,
      ib = (bSeg.ini - bStretch) MOD mb,
      cb = PZLR3Chain.Trim(chb.c^, ib, nb)^,

      step = 0.5d0 * (aStep + bStep)
    DO
      <* ASSERT na > 0 *>
      <* ASSERT nb > 0 *>
      IF aSeg.rev THEN PZLR3Chain.Reverse(ca) END;
      IF bSeg.rev THEN PZLR3Chain.Reverse(cb) END;
      ComputeMatrix(
        ca, cb, o,
        step := step,
        lambda := lambda,
        plotName := "",
        (*OUT*)
        m := m,
        mt := mt,
        mismatch := mismatch,
        length := candLength,
        matchedLength := matchedLength,
        (*WORK*)
        minAvgDist := minAvgDist,
        rCost := rCost
      );
      Wr.PutText(stderr, "# mismatch = " &  Fmt.LongReal(mismatch) & "\n");
      mt := PZMatch.RemoveDanglingEnds(mt^);
      cam := PZLR3Chain.Map(ca, m);
      WITH
        candName = o.output & "/" & Fmt.Pad(Fmt.Int(iCand), 6, '0'),
        candComment = comment & "\n" & CandComment(iCand, cand, ch)
      DO
        ExtractSampleSegments(
          cam^, cb, mt^, candName, 
          o.nSteps, 
          candComment,
          (*IO*) nOutCands
        )
      END;
      Wr.PutText(stderr, "\n");
    END
  END ProcessOneCandidate;

PROCEDURE ExtractSampleSegments(
    READONLY a, b: PZLR3Chain.T;  (* The matched curves. *)
    READONLY m: PZMatch.T;        (* The matching between "a" and "b" samples. *)
    output: TEXT;                 (* Output file name prefix, incl. cand num. *)
    nSteps: NAT;                  (* Requires number of sampling steps in output. *)
    candComment: TEXT;            (* Comments for output files *)
    (*IO*)
    VAR nOutCands: NAT;              (* Number of candidates written so far *)
  ) =
  (*
    Extracts from the given match one or more matching pairs of
    segments with the specified number of sampling steps, then writes each
    segment to disk.  The file names are "XXX-NN/a.flc" and
    "XXX-NN/b.flc" where "XXX" is "output" and "NN" is a 2-digit
    serial count from "00".
  *)
  VAR nChops: NAT := 0;
  BEGIN
    <* ASSERT nSteps MOD 2 = 0 *>
    WITH
      hSteps = nSteps DIV 2
    DO

      PROCEDURE Chop(i, j: NAT): BOOL =
        (*
          Tries to extract one or more segments from the sub-candidate
          defined by the pairs "m[i..j]". The segments will be numbered
          starting from "nChops", which is incremented accordingly.
          Returns TRUE if at least one segment was generated.
        *)
        BEGIN
          WITH
            k = Midpoint(m, i, j)
          DO
            IF m[k][0] - m[i][0] < hSteps 
            OR m[k][1] - m[i][1] < hSteps 
            OR m[j][0] - m[k][0] < hSteps 
            OR m[j][1] - m[k][1] < hSteps THEN
              RETURN FALSE
            ELSIF Chop(i,k) OR Chop(k,j) THEN
              RETURN TRUE
            ELSE
              WITH 
                nTag = Fmt.Pad(Fmt.Int(nChops), 2, '0'),
                name = output & "-" & nTag,
                sa = ExtractSegmentShape(a, m[k,0], nSteps, name, "a", candComment)^,
                sb = ExtractSegmentShape(b, m[k,1], nSteps, name, "b", candComment)^,
                sm = ComputeMeanSignal(sa, sb)^,
                sd = ComputeDiffSignal(sa, sb)^
              DO 
                ProcessSignal(sd, name, "d", candComment);
                ProcessSignal(sm, name, "m", candComment);
              END;
              INC(nChops);
              INC(nOutCands);
              RETURN TRUE
            END
          END
        END Chop;

      BEGIN
        IF NOT Chop(0,LAST(m)) THEN
          Wr.PutText(stderr, "# ExtractSampleSegments: curve too short\n")
        END
      END
    END
  END ExtractSampleSegments;

PROCEDURE Midpoint(
    READONLY m: PZMatch.T;
    i, j: NAT;
  ): NAT =
  (*
    Finds an index "k" such that "s[k] ~ (s[i]+s[j])/2"
    where "s[r] = (m[r][0] + m[r][1]).
  *)
  VAR kMin, kMax: NAT;
  BEGIN
    WITH
      sm = ((m[i][0] + m[i][1]) + (m[j][0] + m[j][1])) DIV 2
    DO
      kMin := i; kMax := j;
      WHILE kMin < kMax DO 
        WITH
          sMin = m[kMin][0] + m[kMin][1],
          sMax = m[kMax][0] + m[kMax][1],
          r = FLOAT(sm - sMin, LONG)/FLOAT(sMax - sMin, LONG),
          k = MIN(kMin + FLOOR(r*FLOAT(kMax - kMin, LONG)), kMax-1),
          sk0 = m[k][0] + m[k][1],
          sk1 = m[k+1][0] + m[k+1][1]
        DO
          <* ASSERT sMin <= sm *>
          <* ASSERT sMax >= sm *>
          IF sm < sk0 THEN
            kMax := k
          ELSIF sm > sk1 THEN
            kMin := k+1
          ELSIF sm - sk0 < sk1 - sm THEN
            RETURN k
          ELSE
            RETURN k+1
          END
        END
      END;
      <* ASSERT kMin = kMax *>
      RETURN kMin
    END 
  END Midpoint;
  
PROCEDURE ExtractSegmentShape(
    READONLY c: PZLR3Chain.T;  (* Mother chain *)
    ctr: NAT;                  (* Index of center sample *)
    nSteps: NAT;               (* Desired number of sampling steps *)
    output: TEXT;             (* Output file name prefix *)
    curveName: TEXT;           (* Curve name ("a" or "b"). *)
    candComment: TEXT;         (* Candidate's comments *)
  ): REF PZLRChain.T =
  (* 
    Extracts from curve "c" a segment of given "length", 
    centered on sample "c[ctr]".  Returns its shape function.
    Writes both to disk.
  *)
  BEGIN
    TRY FS.CreateDirectory(output) EXCEPT OSError.E => (*OK*) END;
    WITH
      hSteps = nSteps DIV 2,
      i = ctr - hSteps,
      j = ctr + hSteps,
      m = j-i+1,
      a = PZLR3Chain.Trim(c, i, m)^,
      rs = PZShape.Signal(a), s = rs^,
      comment = candComment & "\n" & SampleComment(curveName, ctr, i, j)
    DO
      WITH fileName = output & "/" & curveName & ".flc", wr = FileWr.Open(fileName) DO
        Wr.PutText(stderr, "# writing " & fileName & "\n");
        PZLR3Chain.Write(wr, comment, a, unit := ChainUnit);
        Wr.Close(wr)
      END;
      WITH fileName = output & "/" & curveName & ".fsh", wr = FileWr.Open(fileName) DO
        Wr.PutText(stderr, "# writing " & fileName & "\n");
        PZLRChain.Write(wr, comment, s, stride := 0.0d0, unit := ShapeUnit);
        Wr.Close(wr)
      END;
      RETURN rs
    END
  END ExtractSegmentShape;
  
PROCEDURE ProcessSignal(
    READONLY s: PZLRChain.T;
    output: TEXT;             (* Output file name prefix *)
    curveName: TEXT;           (* Curve name ("a" or "b"). *)
    candComment: TEXT;         (* Candidate's comments *)
  ) =
  (*
    Writes to disk the signal "s", and also its Fourier transform
    and power spectra.
  *)
  BEGIN
    TRY FS.CreateDirectory(output) EXCEPT OSError.E => (*OK*) END;
    WITH 
      fft = ComputeSpectrum(s)^,
      N = NUMBER(fft),
      FourierUnit = ShapeUnit / sqrt(FLOAT(N, LONG)),
      PowerUnit = FourierUnit*FourierUnit,
      pwr = ComputePowerSpectrum(fft)^, 
      comment = candComment & "\n" & SignalComment(curveName)
    DO
      WITH fileName = output & "/" & curveName & ".fsh", wr = FileWr.Open(fileName) DO
        Wr.PutText(stderr, "# writing " & fileName & "\n");
        PZLRChain.Write(wr, comment & "\nshape function", s, 
          stride := 0.0d0, unit := ShapeUnit/2.0d0
        );
        Wr.Close(wr)
      END;
      WITH fileName = output & "/" & curveName & ".fft", wr = FileWr.Open(fileName) DO
        Wr.PutText(stderr, "# writing " & fileName & "\n");
        PZLRChain.Write(wr, comment & "\nfourier transform (smart basis)",
          fft, stride := 0.0d0, unit := FourierUnit
        );
        Wr.Close(wr)
      END;
      WITH fileName = output & "/" & curveName & ".fpw", wr = FileWr.Open(fileName) DO
        Wr.PutText(stderr, "# writing " & fileName & "\n");
        PZLRChain.Write(wr, comment & "\nfourier transform (smart basis)",
          pwr, stride := 0.0d0, unit := PowerUnit
        );
        Wr.Close(wr)
      END;
    END;
  END ProcessSignal;

PROCEDURE ComputeMeanSignal(READONLY a, b: PZLRChain.T): REF PZLRChain.T =
  BEGIN
    WITH
      N = NUMBER(a),
      rm = NEW(REF LONGS, N), m = rm^
    DO   
      <* ASSERT NUMBER(b) = N *>
      FOR i := 0 TO N-1 DO
        m[i] := (a[i] + b[i])/2.0d0
      END;
      RETURN rm
    END
  END ComputeMeanSignal;

PROCEDURE ComputeDiffSignal(READONLY a, b: PZLRChain.T): REF PZLRChain.T =
  BEGIN
    WITH
      N = NUMBER(a),
      rd = NEW(REF LONGS, N), d = rd^
    DO   
      <* ASSERT NUMBER(b) = N *>
      FOR i := 0 TO N-1 DO
        d[i] := b[i] - a[i]
      END;
      RETURN rd
    END;
  END ComputeDiffSignal;

PROCEDURE ComputeMatrix(
    READONLY a, b: PZLR3Chain.T;   (* Trimmed segments for the two candidates *)
    READONLY o: Options;
    step: LONG;
    lambda: LONG;
    plotName: TEXT;
    (*OUT*)
    VAR m: LR4x4.T;
    VAR mt: REF PZMatch.T;
    VAR mismatch: LONG;
    VAR length: LONG;
    VAR matchedLength: LONG;
    (*WORK*)
    VAR minAvgDist: REF LONGS;
    rCost: REF PZMatch.CostMatrix;
  ) =

  VAR f, f1: Value;
      nEvals: NAT;
      aMap: REF PZLR3Chain.T;
  BEGIN
    WITH
      ma = NUMBER(a),
      mb = NUMBER(b),
      Ra = Dist(a[0], a[ma-1])/2.0d0,
      Rb = Dist(b[0], b[mb-1])/2.0d0,
      scale = sqrt(2.0d0/3.0d0) * (Ra + Rb)/2.0d0,
      minLength = o.minLength - o.blurFactor * 2.0d0 * lambda,
      minChainEdges = MAX(0, FLOOR((2.0d0 * minLength)/step))
    DO

      PROCEDURE GoalFunction(READONLY m: LR4x4.T): LONG =
        (*
          Computes the image "aMap" of chain "a" by "m",
          the best matching "mt" of "aMap" with "b",
          and returns the mismatch.
        *)
        BEGIN
          aMap := PZLR3Chain.Map(a, m);
          ComputeMatch(aMap^, b, o,
            minChainEdges := minChainEdges,
            step := step,
            mt := mt,
            mismatch := mismatch,
            length := length,
            matchedLength := matchedLength,
            minAvgDist := minAvgDist,
            rCost := rCost
          );
          RETURN mismatch
        END GoalFunction;

      BEGIN
        (* Compute an initial guess for the placement matrix: *)
        m := ComputeInitialMatrix(a, b, o);

        (* Compute initial mapped chain, matching, and mismatch: *)
        f := GoalFunction(m);
        <* ASSERT f = mismatch *>

        (* Optimize matrix *)
        PZOptMatrix.Optimize(
          func := GoalFunction,
          scaleX := scale, scaleY := scale,
          maxEvals := o.maxEvals,
          nEvals := nEvals,
          m := m,
          f := f,
          plotName := plotName
        );

        (* Compute final mapped chain, matching, and mismatch: *)
        f1 := GoalFunction(m);
        <* ASSERT f1 = f *>
        <* ASSERT f = mismatch *>
      END
    END
  END ComputeMatrix;

PROCEDURE ComputeInitialMatrix(
    READONLY a, b: PZLR3Chain.T;
    READONLY o: Options;
  ): LR4x4.T =
  VAR m: LR4x4.T;
  BEGIN
    WITH na = NUMBER(a), nb = NUMBER(b) DO
      IF o.dontGuess THEN
        m := Identity();
      ELSE
        m := ComputeTransformationMatrix(a[0], a[na-1], b[0], b[nb-1]);
      END;
      (* Perturbe initial matrix *)
      IF o.perturbMatrix THEN
        m := PerturbMatrix (m, dTheta := 0.1d0, dX := 1.0d0, dY := 1.0d0)
      END;
    END;
    RETURN m
  END ComputeInitialMatrix;

PROCEDURE ComputeTransformationMatrix(READONLY a1, b1, a2, b2 : LR3.T): LR4x4.T =
  BEGIN
    WITH
      c1 = PZGeo.SegMid(a1, b1),
      c2 = PZGeo.SegMid(a2, b2),
      u  = PZGeo.SegDir(a1, b1),
      v  = PZGeo.SegDir(a2, b2),
      T1 = PZGeo.Translation(LR3.Neg(c1)),
      T2 = PZGeo.Translation(c2),
      R =  PZGeo.Rotation(u, v)
    DO
      RETURN Mul(Mul(T1, R), T2);
    END;
  END ComputeTransformationMatrix;

PROCEDURE PerturbMatrix(VAR m: LR4x4.T; dTheta, dX, dY: LONG): LR4x4.T =
  VAR n: LR4x4.T := LR4x4.Identity();
  BEGIN
    WITH
      cs = cos(dTheta),
      sn = sin(dTheta)
    DO
      <* ASSERT cs # 0.0d0 OR sn # 0.0d0 *>
      n[1,1] := cs;
      n[1,2] := sn;
      n[2,1] := -sn;
      n[2,2] := cs;
      n[0,1] := dX;
      n[0,2] := dY;
    END;
    RETURN LR4x4.Mul(m, n)
  END PerturbMatrix;

PROCEDURE ComputeMatch(
    READONLY a, b: PZLR3Chain.T;
    READONLY o: Options;
    minChainEdges: NAT;
    step: LONG;
    (*OUT*)
    VAR mt: REF PZMatch.T;
    VAR mismatch: LONG;
    VAR length: LONG;
    VAR matchedLength: LONG;
    (*WORK*)
    VAR minAvgDist: REF LONGS;
    VAR rCost: REF PZMatch.CostMatrix;
  ) =
  BEGIN
    WITH
      NA = NUMBER(a),
      NB = NUMBER(b),
      ctr = PickCenter(a, b),
      mtOld = PZMatch.T{ctr}
    DO
      IF minChainEdges <= (NA-1)+(NB-1) THEN
        IF minAvgDist = NIL OR NUMBER(minAvgDist^) < NA + NB - 1 THEN
          minAvgDist := NEW(REF LONGS, NA + NB - 1)
        END;
        PZLR3Chain.RefineMatch(
          a, b,
          step := step,
          mtOld := mtOld,
          maxAdjust := CEILING(o.maxShift / step),
          critDistSqr := o.mp.critDistSqr,
          maxDistSqr := o.mp.maxDistSqr,
          (*OUT*)
          mt := mt,
          mismatch := mismatch,
          length := length,
          matchedLength := matchedLength,
          (*WORK*)
          minAvgDist := SUBARRAY(minAvgDist^, 0, NA+NB-1),
          rCost := rCost
        );
        (* If the matching has zero length, provide a useful "mismatch": *)
        IF mt = NIL OR NUMBER(mt^) < 2 OR length <= 0.0d0 THEN 
          <* ASSERT mismatch = 0.0d0 *>
          WITH
            ac = (NA-1) DIV 2,
            bc = IndexNearestPoint(b, a[ac])
          DO
            mismatch := LR3.DistSqr(a[ac], b[bc]);
          END
        END
      ELSE
        mismatch := LAST(LONG);
        length := 0.0d0;
        matchedLength := 0.0d0;
        mt := NEW(REF PZMatch.T, 1);
        mt[0] := ctr
      END
    END;
  END ComputeMatch;

PROCEDURE PickCenter(READONLY a, b: PZLR3Chain.T): PZMatch.Pair =
  BEGIN
    WITH
      ai = IndexNearestPoint(a, b[0]),
      af = IndexNearestPoint(a, b[LAST(b)]),
      bi = IndexNearestPoint(b, a[0]),
      bf = IndexNearestPoint(b, a[LAST(a)])
    DO
      RETURN PZMatch.Pair{(ai + af) DIV 2, (bi + bf) DIV 2}
    END
  END PickCenter;

PROCEDURE CandComment(
    iCand: NAT;
    READONLY cand: PZCandidate.T;
    READONLY ch: ARRAY OF ChainData;
  ): TEXT =
  BEGIN
    WITH
      wr = NEW(TextWr.T).init()
    DO
      Wr.PutText(wr, "Samples from candidate " & Fmt.Pad(Fmt.Int(iCand), 6, '0') & "\n");
      FOR j := 0 TO 1 DO 
        WITH 
          seg = cand.seg[j],
          chj = ch[seg.cvx],
          fin = (seg.ini + seg.ns - 1) MOD seg.tot,
          m = NUMBER(chj.c^),
          fm = FLOAT(m, LONG),
          pci = FLOAT(seg.ini, LONG) / fm,
          pcf = FLOAT(fin, LONG) / fm,
          pc = FLOAT(seg.ns, LONG) / fm
        DO
          Wr.PutText(wr, "  Segment [" & FI(j,1) & "] from chain " & FZI(seg.cvx,4));
          IF seg.rev THEN Wr.PutText(wr, " (reversed)") END;
          Wr.PutText(wr, "\n");
          Wr.PutText(wr, "    " & FI(seg.ns,1));
          Wr.PutText(wr, " samples [" & FI(seg.ini,4) & ".." & FI(fin,4) & "]\n"); 
          Wr.PutText(wr, "    " & FLR(pc,5,3));
          Wr.PutText(wr, " of the curve [" & FLR(pci,5,3) & "__" & FLR(pcf,5,3) & "]\n"); 
        END
      END;
      RETURN TextWr.ToText(wr)
    END
  END CandComment;

PROCEDURE SampleComment(
    curve: TEXT;    (* "a" or "b" *)
    ctr: NAT;       (* Relative index of center point *)
    i, j: NAT;      (* Range of relative indices *)
  ): TEXT = 
  BEGIN
    WITH
      wr = NEW(TextWr.T).init()
    DO
      Wr.PutText(wr, "  subsegment " & curve & "[" & FI(i,1) & ".." & FI(j,1) & "]");
      Wr.PutText(wr, " center = " & curve & "[" & FI(ctr,1) & "]");
      RETURN TextWr.ToText(wr)
    END
  END SampleComment;

PROCEDURE SignalComment(curveName: TEXT): TEXT = 
  BEGIN
    IF Text.Equal(curveName, "m") THEN
      RETURN "average of the two shape functions"
    ELSIF Text.Equal(curveName, "d") THEN
      RETURN "difference of the two shape functions"
    ELSE
      <* ASSERT FALSE *>
    END
  END SignalComment;
  
PROCEDURE IndexNearestPoint(READONLY b: PZLR3Chain.T; p: LR3.T ): NAT =
  VAR index : NAT;
      mindist :  LONG;
  BEGIN
     WITH
       m = NUMBER(b),
       mHalf = (m-1) DIV 2
     DO
       index := 0; mindist := Dist(b[0],p);
       FOR i:= mHalf TO 1 BY -1  DO
         WITH d = Dist(b[i],p) DO
           IF d <  mindist THEN index := i; mindist := d END;
         END
       END;
       FOR i:= mHalf+1 TO m-1 DO
         WITH d = Dist(b[i],p) DO
           IF d <  mindist THEN index := i; mindist := d END;
         END
       END;
     END;
     RETURN index
  END IndexNearestPoint;

PROCEDURE GetOptions(): Options =
  VAR o: Options;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO
      TRY
        (* Data identification parameters: *)

        pp.getKeyword("-input");
        o.input := pp.getNext();

        pp.getKeyword("-chainPrefix");
        o.chainPrefix := pp.getNext();

        IF pp.keywordPresent("-chainDir") THEN 
          o.chainDir := pp.getNext()
        ELSE
          o.chainDir := "."
        END;
        
        pp.getKeyword("-band");
        o.band := pp.getNextInt();

        IF pp.keywordPresent("-extension") THEN
          o.extension := pp.getNext()
        ELSE
          o.extension := ".flc"
        END;

        pp.getKeyword("-output");
        o.output := pp.getNext();

        (* Matching parameters: *)

        pp.getKeyword("-minLength");
        o.minLength := pp.getNextLongReal();

        pp.getKeyword("-blurFactor");
        o.blurFactor := pp.getNextLongReal(0.0d0, 5.0d0);

        pp.getKeyword("-extraLength");
        o.extraLength := pp.getNextLongReal(0.0d0);

        pp.getKeyword("-maxShift");
        o.maxShift := pp.getNextLongReal();

        o.mp := PZMismatch.ParseOptions(pp);
        
        (* Placement/optimization parameters: *)

        IF pp.keywordPresent("-dontGuess") THEN
          o.dontGuess := TRUE;
        ELSE
          o.dontGuess := FALSE;
        END;

        o.perturbMatrix := pp.keywordPresent("-perturbMatrix");

        IF pp.keywordPresent("-dontOptimize") THEN
          o.maxEvals := 0
        ELSE
          pp.getKeyword("-maxEvals");
          o.maxEvals := pp.getNextInt();
        END;

        (* Segment extraction parameters: *)

        pp.getKeyword("-nSteps");
        o.nSteps := pp.getNextInt(2,4096);
        
        IF pp.keywordPresent("-maxCands") THEN
          o.maxCands := pp.getNextInt(1)
        ELSE
          o.maxCands := LAST(NAT)
        END;

        pp.finish();
      EXCEPT
      | ParseParams.Error =>
          Wr.PutText(stderr, "Usage: BITExtractSamples \\\n");
          Wr.PutText(stderr, "  -input NAME \\\n");
          Wr.PutText(stderr, "  -chainPrefix NAME [ -chainDir DIR ] \\\n");
          Wr.PutText(stderr, "  -band BAND [ -extension EXT ] \\\n");
          Wr.PutText(stderr, "  -output  NAME \\\n");
          Wr.PutText(stderr, "  -maxShift NUM \\\n");
          Wr.PutText(stderr, "  " & PZMismatch.OptionsHelp & " \\\n");
          Wr.PutText(stderr, "  -minLength NUM \\\n");
          Wr.PutText(stderr, "  -blurFactor NUM -extraLength NUM \\\n");
          Wr.PutText(stderr, "  [ -dontGuess ] [ -perturbMatrix ] \\\n");
          Wr.PutText(stderr, "  { -dontOptimize | -maxEvals  NUM } \\\n");
          Wr.PutText(stderr, "  -nSteps NUMBER [ -maxCands NUMBER ]\n");
         Process.Exit(1);
      END;
    END;
    RETURN o
  END  GetOptions;

PROCEDURE FI(x: INT; w: NAT): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.Int(x), w)
  END FI;

PROCEDURE FZI(x: INT; w: NAT): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.Int(x), w, '0')
  END FZI;

PROCEDURE FLR(x: LONG; w, d: NAT): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(x, prec := d, style := Fmt.Style.Fix), w)
  END FLR;

BEGIN
  Main()
END BITExtractSamples.