MODULE PZMultiScale EXPORTS Main;

(* multi-scale matching *)
(* Last edited on 1999-08-07 07:28:20 by hcgl *)

IMPORT PSPlot, LR3, LR4x4; 
IMPORT ParseParams, Process, FileRd, Wr, OSError, Thread, Stdio, Fmt, Rd;
IMPORT PZLRChain, PZSymbolChain, PZMatch, PZSymbol, PZPlot, PZGeo;
IMPORT PZLR3Chain;  

FROM Stdio IMPORT stderr;
FROM PZTypes IMPORT LONG, LONGS, INT, NAT, NATS, BOOL, BOOLS;
 
<* FATAL Wr.Failure, Thread.Alerted, OSError.E *>

TYPE  
  Interval = RECORD (* Identifies a segment of a filtered curve *)
      k : NAT;  (* Curve identifier *)
      i : LONG;  (* Fractional sample index of interval starting point *)
      j : LONG;  (* Fractional sample index of interval ending point *)
    END;
      (*
        The values "i" and "j" lie in the range "[0 _ m)" where
        "m" is the number of samples in the filtered curve.
        Note that "j" may be less than "i" in case o f warp-around.
        They are never equal. *)
   
  Candidate = RECORD 
      int:  ARRAY[0..1] OF Interval;
      dist:  LONG;
      count: NAT;   (* Number of matched steps *)
    END;              
      (*         
        A candidate matching between two curves. *)
   
  CandidateList = ARRAY OF Candidate; 
    
CONST 
  ExpectedMatchesPerCurve = 6;
  MaxMatchesPerCurve = 1000;
  Delta = 1.02d0;     (* Candidate widening factor *)
  Alpha = 1.2d0;      (* Corner broadening factor *)
  Lambda = 0.85d0;    (* MIN percentage of matching steps in match *)
  MaxPlotCands = 160; (* Max number of candidates to plot per band *)

TYPE
  MismatchOptions = PZMismatch.Options;
  Options = RECORD
      inName: TEXT;          (* Input CurvChain basis file *)
      nChains: NAT;     (* Input how many chains that will be read *)
      outName: TEXT;         (* Output float/curv chain file name (without ".???") *)
      bandMax: NAT;     (* maximum "band" (filtering wavelength parameter) *)
      bandMin: NAT;     (* minimum "band" (filtering wavelength parameter) *) 
      initDist: LONG;    (* Max average distance for initial candidate generation *)
      initLength: LONG;  (* Min length of initial candidates *)
      cutDist: REF CutDists; (* Max average matching distance for each band. *)
      cutLength: LONG;   (* Minimum length of matching segments *)
    END;
 
  CutDists = LONGS;

PROCEDURE Main() =
  VAR
    band : NAT;                          (* Filtering "band" parameter *)
    iBand: NAT;                          (* Band counter *)
    cand : REF CandidateList;                 (* Candidate matches *)
    chain : REF ARRAY OF REF PZSymbolChain.T; (* Coded curvature chains *)
    tPrev, t: REF ARRAY OF REF PZLRChain.T;   (* Time vector of each chain *)
    mPrev, m: REF NATS;          (* Number of samples in each chain *)
    step: LONG;                           (* Sampling stepsize for current chains *)
    minSize: NAT;                        (* Minimum number of letters in match *)
    extra: LONG;                          (* Retters to add on each side *)
  BEGIN
    WITH 
      o = GetOptions(),
      nChains = o.nChains
    DO
      band := o.bandMax; iBand := 0;
      chain := PZSymbolChain.ReadAll(o.inName, nChains, band, (*OUT*) step);
      m := GetNumSamples(chain^);
      
      minSize := FLOOR(o.initLength / step);
      cand := ComputeInitialCandidates(chain^, o.initDist, minSize);
      Wr.PutText(stderr, "step: " &  FLR(step,6,4) & "\n"); 
      PrintCandidates(chain^, cand^, m^, band, step);
      
      minSize := FLOOR(Lambda * (o.cutLength - Alpha * FLOAT(band, LONG)) / step);
      extra := 1.0d0;
      cand := RefineCandidates(cand^,chain^,o.cutDist[iBand],minSize,extra);
      PrintCandidates(chain^, cand^, m^, band, step);
      DrawCandidates(o.inName, nChains, band, cand^, m^);
      
      chain := NIL;
      t := ReadAllTMCFiles(o.inName, nChains, band);
      band :=  band DIV 2; INC(iBand);
      WHILE band >= o.bandMin AND NUMBER(cand^) > 0 DO 
        tPrev := t;
        mPrev := m;
        (* Read finer band chains *)
        chain := ReadAllCVCFiles(o.inName, nChains, band, (*OUT*) step);
        Wr.PutText(stderr, "step: " &  FLR(step,6,4) & "\n"); 
        m := GetNumSamples(chain^);
        t := ReadAllTMCFiles(o.inName, nChains, band);
        MapAllCandidates(cand^, tPrev^, mPrev^, t^, m^); 
        tPrev := NIL;
        (* General Match  *)
        Wr.PutText(stderr,"band :" & Fmt.Int(band) & "\n");
        (* PrintCandidates(chain^, cand^, m^, band, step); *)
        
        minSize := FLOOR(Lambda * (o.cutLength - Alpha * FLOAT(band, LONG)) / step);
        extra := 0.5d0 * Delta * Alpha * FLOAT(band, LONG) / step + 1.0d0;
        cand := RefineCandidates(cand^, chain^, o.cutDist[iBand], minSize, extra);
        Wr.PutText(stderr, Fmt.Int(NUMBER(cand^)) & " matches found\n");
        Wr.PutText(stderr," Final result: \n");
        PrintCandidates(chain^, cand^, m^, band, step);
        DrawCandidates(o.inName, nChains, band, cand^, m^);
        chain := NIL;
        band :=  band DIV 2; INC(iBand);
      END;
    END;
  END Main;          

      
PROCEDURE GetOptions(): Options =
  VAR o: Options;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO
      TRY
        pp.getKeyword("-inName");
        o.inName := pp.getNext();
        
        pp.getKeyword("-nChains");
        o.nChains := pp.getNextInt();
        
        pp.getKeyword("-outName");
        o.outName := pp.getNext();
        
        pp.getKeyword("-bandMax");
        o.bandMax := pp.getNextInt();
        
        pp.getKeyword("-bandMin");
        o.bandMin := pp.getNextInt();
        
        pp.getKeyword("-initLength");
        o.initLength := pp.getNextLongReal(0.0d0);
        
        pp.getKeyword("-initDist");
        o.initDist := pp.getNextLongReal(0.0d0);
        
        pp.getKeyword("-cutLength");
        o.cutLength := pp.getNextLongReal();
        
        pp.getKeyword("-cutDist");
        VAR nBands: NAT := 0; band: NAT := o.bandMax;
        BEGIN
          WHILE band >= o.bandMin DO band := band DIV 2; INC(nBands) END;
          o.cutDist := NEW(REF CutDists, nBands);
          FOR b := 0 TO nBands-1 DO
            o.cutDist[b] := pp.getNextLongReal(0.0d0);
          END
        END;
                 
        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: PZMultiScale \\\n");
          Wr.PutText(stderr, "  -inName   NAME \\\n");
          Wr.PutText(stderr, "  -nChains  NUMBER \\\n");

          Wr.PutText(stderr, "  -outName  NAME \\\n");
          Wr.PutText(stderr, "  -bandMax  NUMBER \\\n");
          Wr.PutText(stderr, "  -bandMin  NUMBER \\\n");
          Wr.PutText(stderr, "  -initLength NUMBER \\\n");
          Wr.PutText(stderr, "  -initDist   NUMBER \\\n");
          Wr.PutText(stderr, "  -cutLength  NUMBER \\\n");
          Wr.PutText(stderr, "  -cutDist    NUMBER ... NUMBER \n");
        Process.Exit(1);
      END;
    END;
    RETURN o
  END GetOptions;


PROCEDURE ComputeInitialCandidates(
   READONLY chain : ARRAY OF REF PZSymbolChain.T; 
   cutDist : LONG;   (* MAX average distance *)
   minSize: NAT;    (* Min number of letters in match *)
  ): REF CandidateList =
  
  VAR 
    c : REF CandidateList;  
    n : NAT := 0; 
  BEGIN 
    WITH
      nr = NUMBER(chain) 
    DO 
      c := NEW(REF CandidateList, nr * ExpectedMatchesPerCurve);
      FOR ka := 0 TO nr-1 DO 
        WITH
          a = chain[ka]^,
          ma = NUMBER(a) 
        DO
          FOR kb := ka+1 TO nr-1 DO 
            WITH 
              b = chain[kb]^,   mb = NUMBER(b),
              mxDiv = Mdc(ma,mb), 
              lenDiag = ma * mb DIV  mxDiv
            DO
              Wr.PutText(stderr, "ka: " & Fmt.Int(ka) &
               "   kb: " & Fmt.Int(kb) );
              Wr.PutText(stderr, "  ma: " & Fmt.Int(ma) &
               "   mb: " & Fmt.Int(mb)  & "    mxDiv : " & Fmt.Int(mxDiv));
              Wr.PutText(stderr,"  lenDiag : " & Fmt.Int(lenDiag) & "\n");
              FOR t := 0 TO mxDiv-1 DO 
                WITH d2 = NEW(REF NATS, lenDiag)
                DO   
                  FOR r := 0 TO lenDiag-1 DO 
                    d2[r] :=  PZSymbol.DistSqr(a[r MOD ma], b[(t-r) MOD mb],TRUE); 
                  END;
                  BreakDiagonal(t,d^,c,n,ma,mb,ka,kb,cutDist,minSize);
                END; 
              END;   
            END;
          END;
        END;
      END;
      WITH ct = NEW(REF CandidateList,n) DO
        ct^ :=  SUBARRAY(c^, 0, n);
        RETURN ct;
      END;
    END;
  END ComputeInitialCandidates;
 
PROCEDURE BreakDiagonal(
    t : NAT ;                    (* Index of diagonal *)
    READONLY d : NATS;
    VAR c: REF  CandidateList;
    VAR  n:  NAT; 
    ma : NAT;
    mb : NAT;
    ka : NAT ;
    kb : NAT;
    cutDist : LONG;
    minSize : NAT;
  ) =
  VAR 
    k: NAT;
    count,  kStart : NAT;
    dist : NAT;
    inx: NAT;
  BEGIN 
    WITH 
       nd = NUMBER(d), 
       dm  = NEW(REF BOOLS, nd)^
    DO 
      <* ASSERT nd MOD ma = 0 *>
      <* ASSERT nd MOD mb = 0 *>
      FOR r := 0 TO nd-1 DO 
        IF FLOAT(d[r] + d[(r+1) MOD nd], LONG) <= 2.0d0 * cutDist 
          THEN dm[r] := TRUE; 
          ELSE dm[r] := FALSE; 
        END;
      END;
      (* Look for a TRUE element "dm[kStart]" right after a FALSE element: *)
      k := 0;
      WHILE dm[k MOD nd] = TRUE AND k < nd DO INC(k) END;
      IF k = nd THEN 
        (* All TRUE *)
        kStart := 0
      ELSE
        inx := 0;
        WHILE dm[k MOD nd] = FALSE AND inx < nd DO INC(k); INC(inx) END;
        IF inx = nd THEN (* All FALSE *) RETURN END;
        kStart := k MOD nd; 
      END;
      (* Collect runs of TRUEs not shorter than "minSize": *)
      k := kStart;
      REPEAT 
        (* Here begins a new run *)
        dist := d[k MOD nd];
        count := 1;
        WITH 
          ia = k MOD ma,
          fb = (t - k) MOD mb
        DO 
          WHILE dm[(k+1) MOD nd] = TRUE AND count < nd DO
            dist := dist + d[(k+1) MOD nd];
            INC(count);
            INC(k);  
          END;
          IF count >= minSize THEN 
            IF n >= NUMBER(c^) THEN 
              WITH ct  = NEW(REF CandidateList, n * 2) DO 
                SUBARRAY(ct^, 0, NUMBER(c^)) := c^;
                c := ct; 
              END;
            END;
            WITH 
              fa = k MOD ma, 
              ib = (t - k) MOD mb,   
              s = c[n], sa = s.int[0], sb = s.int[1]
            DO 
              sa.k := ka; sb.k := kb;
              sa.i := ReduceSampleIndex(FLOAT(ia,LONG), ma); 
              sa.j := ReduceSampleIndex(FLOAT(fa,LONG), ma);
              sb.i := ReduceSampleIndex(FLOAT(ib,LONG), mb); 
              sb.j := ReduceSampleIndex(FLOAT(fb,LONG), mb);
              s.dist := FLOAT(dist,LONG) / FLOAT(count,LONG);
              s.count := count;
              INC(n);
            END;  
          END; 
        END;
        INC(k);
        WHILE dm[k MOD nd] = FALSE  DO INC(k) END;
      UNTIL (k MOD nd) = kStart ;
    END;  
  END BreakDiagonal;

  
PROCEDURE Mdc(a : INT; b : INT) : INT =
  VAR m : NAT := 1;
      i : NAT := 2;
  BEGIN 
    WHILE  i <= MIN(a,b) DO 
      IF a MOD i = 0 AND b MOD i = 0 THEN 
        a := a DIV i; b := b DIV i;
        m := m * i;
        ELSE  INC(i);  
      END;
    END;
    RETURN m;
  END Mdc;
 

PROCEDURE MapAllCandidates(
    VAR cand : CandidateList;
    READONLY tPrev : ARRAY OF REF PZLRChain.T;
    READONLY mPrev : NATS;
    READONLY t : ARRAY OF REF PZLRChain.T;
    READONLY m: NATS;
  ) =
  BEGIN 
    WITH nCand = NUMBER(cand)  DO 
      FOR i := 0 TO nCand-1  DO 
        FOR j := 0 TO 1 DO 
          WITH 
            s = cand[i].int[j],
            tca = tPrev[s.k]^,
            ma = mPrev[s.k],
            tcb = t[s.k]^,
            mb = m[s.k]
          DO
            s.i := MapSampleIndex(tca, ma, tcb, mb, s.i);
            s.j := MapSampleIndex(tca, ma, tcb, mb, s.j);
          END
        END;
      END;
    END;
  END MapAllCandidates;
 

PROCEDURE MapSampleIndex(
    READONLY tca : PZLRChain.T;
    ma : NAT;
    READONLY tcb : PZLRChain.T;
    mb: NAT;
    x: LONG;
  ): LONG =
  (*
    Returns the index "y" of a point in a filtered curve "fb", with "mb"
    equally-spaced samples, given its index "x" in some other filtered curve "fa", 
    with "ma" equally-spaced samples, and the corresponding time
    vectors "tca" and "tcb". *)
  BEGIN
    WITH
      La = tca[LAST(tca)] - tca[0],
      ta = tca[0] + x * La / FLOAT(ma, LONG),
      tb = PZLRChain.MapTime(tca, tcb, ta),
      Lb = tcb[LAST(tcb)] - tcb[0],
      y = (tb - tcb[0])/Lb * FLOAT(mb, LONG)
    DO
      <* ASSERT y >= 0.0d0 *>
      <* ASSERT y < FLOAT(mb, LONG) *>
      RETURN y
    END
  END MapSampleIndex;

 
PROCEDURE RefineCandidates(
    READONLY cand: CandidateList;
    READONLY chain: ARRAY OF REF PZSymbolChain.T; 
    cutDist: LONG;
    minSize: NAT;
    extra: LONG; 
  ) : REF CandidateList =
  VAR 
    match : REF PZMatch.T;
    dist  : LONG;
    n, nMatched : NAT := 0;
    c : REF CandidateList;
  BEGIN 
    WITH
      nCand = NUMBER(cand)
    DO
      (* Wr.PutText(stderr,"RefineCandidates: nCand = " & Fmt.Int(nCand) & "\n"); *)
      c := NEW(REF CandidateList, nCand);
      FOR i := 0 TO nCand-1 DO 
        WITH 
          s = cand[i],
          sa = s.int[0],
          ma = NUMBER(chain[sa.k]^),
          ia = FLOOR(sa.i - extra) MOD ma,
          na = (CEILING(sa.j + extra) - ia) MOD ma + 1,
          ca = PZSymbolChain.Trim(chain[sa.k]^, ia, na)^,
          sb = s.int[1],
          mb = NUMBER(chain[sb.k]^),
          ib = FLOOR(sb.i - extra) MOD mb,
          nb = (CEILING(sb.j + extra) - ib) MOD mb + 1,
          cb = PZSymbolChain.Trim(chain[sb.k]^, ib, nb)^
        DO 
          IF na > 0 AND nb > 0  THEN 
            <* ASSERT NUMBER(ca)>0 *>
            <* ASSERT NUMBER(cb)>0 *>
            PZSymbolChain.ReverseAndComplement(cb);
            ComputeCurvChainAvgDist(ca, cb, TRUE, cutDist, dist, nMatched, match);
            PZSymbolChain.ReverseAndComplement(cb);
            IF dist <= cutDist AND nMatched > minSize THEN 
              (* Wr.PutText(stderr, "dist : " &  Fmt.LongReal(dist) & "\n"); *)
              WITH
                r = c[n],
                ra = r.int[0],
                rb = r.int[1],
                nm = NUMBER(match^),
                mia = match[0,0], 
                mfa = match[nm-1,0],
                mib = nb - 1 - match[nm-1,1], 
                mfb = nb - 1 - match[0,1]  
              DO
                r.dist := dist;
                r.count := nMatched;
                ra.i := ReduceSampleIndex(FLOAT(ia+mia, LONG), ma);
                ra.j := ReduceSampleIndex(FLOAT(ia+mfa, LONG), ma);
                ra.k := sa.k;
                rb.i := ReduceSampleIndex(FLOAT(ib+mib, LONG), mb);
                rb.j := ReduceSampleIndex(FLOAT(ib+mfb, LONG), mb);
                rb.k := sb.k;
              END;
              INC(n);
            END;
          END;
        END;
      END;
      SortAndPruneCandidates(c^, n, NUMBER(chain), MaxMatchesPerCurve);
      WITH ct = NEW(REF CandidateList, n) DO
        ct^ := SUBARRAY(c^,0,n);
        RETURN ct
      END;
    END
  END RefineCandidates;
 
 
PROCEDURE SortAndPruneCandidates(
    VAR c: CandidateList;
    VAR nCands: NAT;
    nChains : NAT;
    maxCandsPerCurve: NAT;
  ) =
  
  PROCEDURE RemoveDuplicates() =
    VAR m: NAT;
    BEGIN 
      m := 1;
      FOR i := 1 TO nCands-1 DO
        IF (c[i].int[0] # c[m-1].int[0]) 
        OR (c[i].int[1] # c[m-1].int[1])
        THEN
          IF m # i THEN c[m] := c[i] END;
          INC(m)
        END
      END;
      nCands := m;
    END RemoveDuplicates;

  PROCEDURE RemoveExcessCandidates() =
    VAR m: NAT := 0;
    BEGIN 
      WITH candsPerCurve = NEW(REF NATS, nChains)^ DO
        Wr.PutText(stderr,"RemoveExcessCandidates ...nCands:" & Fmt.Int(nCands) & "\n"); 
        FOR i := 0 TO nChains-1 DO candsPerCurve[i] := 0 END;
        FOR i := 0 TO nCands-1 DO 
          WITH
            k0 = c[i].int[0].k,
            ct0 = candsPerCurve[k0],
            k1 = c[i].int[1].k,
            ct1 = candsPerCurve[k1]
          DO
            IF ct0 < maxCandsPerCurve OR ct1 < maxCandsPerCurve THEN
              INC(ct0); INC(ct1);
              IF i > m THEN c[m] := c[i] END;
              INC(m)
            END
          END
        END; 
        nCands := m;
        Wr.PutText(stderr,"final nCands:" & Fmt.Int(nCands) & "\n");
      END;
    END RemoveExcessCandidates; 
  
  PROCEDURE Better(READONLY a, b: Candidate): BOOL =
    BEGIN
      IF a.dist # b.dist THEN RETURN a.dist < b.dist END;
      IF a.count # b.count THEN RETURN a.count > b.count END;
      IF a.int[0].k # b.int[0].k THEN RETURN a.int[0].k < b.int[0].k END;
      IF a.int[1].k # b.int[1].k THEN RETURN a.int[1].k < b.int[1].k END;
      IF a.int[0].i # b.int[0].i THEN RETURN a.int[0].i < b.int[0].i END;
      IF a.int[1].i # b.int[1].i THEN RETURN a.int[1].i < b.int[1].i END;
      IF a.int[0].j # b.int[0].j THEN RETURN a.int[0].j < b.int[0].j END;
      IF a.int[1].j # b.int[1].j THEN RETURN a.int[1].j < b.int[1].j END;
      RETURN FALSE
    END Better;
  
  VAR j: NAT;
      s: Candidate;
  BEGIN
    IF nCands = 0 THEN RETURN END;
    FOR i := 1 TO nCands-1 DO
      IF Better(c[i], c[i-1]) THEN
        s := c[i];
        j := i;
        WHILE j > 0 AND Better(s, c[j-1]) DO 
          c[j] := c[j-1]; DEC(j)
        END;
        c[j] := s
      END
    END;
    RemoveDuplicates();
    RemoveExcessCandidates(); 
    (* Juntar candidatos sobrepostos? *)
  END SortAndPruneCandidates;


PROCEDURE ReduceSampleIndex(x: LONG; m: NAT): LONG =
  (*
    Reduces the sample index "x" modulo the sample count "m". *)
  BEGIN
    WITH
      fm = FLOAT(m, LONG)
    DO
      WHILE x < 0.0d0 DO x := x + fm END;
      WHILE x >= fm DO x := x - fm END;
      RETURN x
    END
  END ReduceSampleIndex;


PROCEDURE GetNumSamples(READONLY chain: ARRAY OF REF PZSymbolChain.T): REF NATS =
  BEGIN
    WITH 
      nChains = NUMBER(chain),
      m = NEW(REF NATS, nChains)
    DO 
      FOR i := 0 TO nChains-1 DO 
        m[i] := NUMBER(chain[i]^)
      END;
      RETURN m
    END
  END GetNumSamples;
  

PROCEDURE PrintCandidates (
   READONLY chain: ARRAY OF REF  PZSymbolChain.T;
   READONLY cand : CandidateList; 
   READONLY  m: NATS; 
   band: NAT;
   step: LONG;
  )=
  BEGIN 
    Wr.PutText(stderr, "=== band = " & Fmt.Int(band));
    Wr.PutText(stderr, " step = " & Fmt.LongReal(step) & " === \n");
    FOR i := 0 TO LAST(cand)  DO
      Wr.PutText(stderr, "curves = (");
      Wr.PutText(stderr, Fmt.Int(cand[i].int[0].k) & "," & Fmt.Int(cand[i].int[1].k)); 
      Wr.PutText(stderr, ") ");
      Wr.PutText(stderr,"dist = " & FLR(cand[i].dist, 6, 2)); 
      Wr.PutText(stderr, " ");
      Wr.PutText(stderr,"nMatched = " & Fmt.Int(cand[i].count) & "\n"); 
      FOR j := 0 TO 1 DO 
        WITH 
          c = cand[i].int[j],
          mc = m[c.k],
          ns  = FLOAT(mc,LONG),
          ic = FLOOR(c.i) MOD mc,
          nc = (CEILING(c.j) - ic) MOD mc + 1,
          pci = FLOAT(ic, LONG) / ns,
          pcj = FLOAT(ic + nc, LONG) / ns
        DO 
          Wr.PutText(stderr,
            "cand[" & Fmt.Int(i) & "," & Fmt.Int(j) & "]: k = " & Fmt.Int(c.k)); 
          Wr.PutText(stderr, " i = " & FLR(c.i, 6, 2));
          Wr.PutText(stderr, "(" & FLR(pci, 3, 2));  
          Wr.PutText(stderr, ") j = " & FLR(c.j, 6, 2)); 
          Wr.PutText(stderr, "(" & FLR(pcj, 3, 2) & ") "  );
          IF j = 0 THEN 
            FOR r := 0 TO nc-1  DO 
              WITH ch = chain[c.k][(ic + r) MOD mc] DO
                Wr.PutText(stderr, Fmt.Char(ch))
              END; 
            END;
          ELSE
            FOR r := nc-1 TO 0 BY -1 DO 
              WITH ch = chain[c.k][(ic + r) MOD mc] DO
                Wr.PutText(stderr, Fmt.Char(PZSymbol.Complement(ch)))
              END; 
            END;
          END;  
          Wr.PutText(stderr, "\n");
        END;
      END; 
      FOR i := 1 TO 80 DO Wr.PutChar(stderr, '-') END;
      Wr.PutText(stderr, "\n");
    END;  
  END PrintCandidates; 
 

PROCEDURE FLR(x: LONG; w, p: NAT): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Fix, prec := p), w)
  END FLR;
  
  
PROCEDURE ComputeCurvChainAvgDist(
    READONLY a, b: PZSymbolChain.T;
    trimEnds : BOOL;
    cutDist: LONG;
    VAR (*OUT*) avgDist : LONG;
    VAR (*OUT*) nMatched: NAT;
    VAR (*OUT*) match: REF PZMatch.T;
  ) =
  BEGIN
    PZSymbolChain.Match(
      a, b, 
      skipCost := 2.0d0 * cutDist, (* !!! Deveria ser outra formula... *)
      removeUnmatchedEnds := trimEnds,
      avgCost := avgDist, nMatched := nMatched, match := match
    );
  END ComputeCurvChainAvgDist;


PROCEDURE DrawCandidates(
   inName : TEXT;
   nChains, band : NAT;
   READONLY cand: CandidateList;
   READONLY cvcSize: NATS;
  )=
  VAR w: LR3.T; 
  CONST 
    PlotSize = 40.0d0;
    NPlotsPerRow = 4; 
    NRows = 4;
    NPlotsPerPage = (NRows * NPlotsPerRow);
    MaxFragmentWidth = 300.0d0;
    PlotXMargin = 0.5d0 * (PSPlot.LetterXSize - FLOAT(NPlotsPerRow,LONG)*PlotSize);
    PlotYMargin = 0.5d0 * (PSPlot.LetterYSize - FLOAT(NRows,LONG)*PlotSize);
  BEGIN 
    w[0] := 0.0d0; w[2] := 0.0d0;  
    Wr.PutText(stderr,"DrawCandidates...\n");
    WITH 
      nCands = NUMBER(cand),
      name = inName & "-f" & Fmt.Pad(Fmt.Int(band), 3, '0') & "-cand.ps",
      flc = ReadAllFLCFiles(inName, nChains, band)^,
      ps = NEW(PSPlot.PSFile).open(name)
    DO
      Wr.PutText(stderr, name & "\n");  
      FOR i := 0 TO MIN(MaxPlotCands, nCands)-1  DO 
        WITH 
          plotNum = i MOD NPlotsPerPage,
          plotX = plotNum MOD NPlotsPerRow,
          plotY = plotNum DIV NPlotsPerRow,
          xSize = PlotSize,
          ySize = PlotSize,
          xCenter = PlotSize * (FLOAT(plotX, LONG) + 0.5d0) + PlotXMargin,
          yCenter = PlotSize * (FLOAT(plotY, LONG) + 0.5d0) + PlotYMargin
        DO 
          Wr.PutText(stderr, "plot(" & Fmt.Int(plotX) & "," & Fmt.Int(plotY) & ")\n");
          IF plotNum = 0 THEN ps.beginPage() END;
          ps.beginDrawing(xSize, ySize, xCenter, yCenter);
          ps.setScale(PSPlot.Axis.X, -MaxFragmentWidth, +MaxFragmentWidth);
          ps.setScale(PSPlot.Axis.Y, -MaxFragmentWidth, +MaxFragmentWidth);
          ps.setLineWidth(0.2);
          ps.setLineColor(PSPlot.Black);
          ps.frame();
          FOR j := 0 TO 1 DO 
            IF j = 0 THEN w[1]:= +1.0d0 ELSE w[1]:= -1.0d0 END;
            WITH 
              interval = cand[i].int[j],
              chainNum = interval.k,
              chain = flc[chainNum]^,
              n = NUMBER(chain),
              m = cvcSize[chainNum],
              scale = FLOAT(n, LONG)/FLOAT(m, LONG),
              ini = FLOOR(interval.i * scale), 
              fin = CEILING(interval.j * scale),
              p = chain[ini MOD n],
              q = chain[fin MOD n],
              o = PZGeo.SegMid(p, q), 
              u = PZGeo.SegDir(p, o),
              t = PZGeo.Translation(LR3.Neg(o)),
              r = PZGeo.Rotation(u, w),
              m = LR4x4.Mul(t,r),
              mappedChain = PZLR3Chain.Map(chain, m)^
            DO 
              DrawCandidate(ps, mappedChain, ini, fin);
            END;
          END;
          ps.endDrawing();
          IF plotNum = NPlotsPerPage-1  THEN ps.endPage() END;
        END  
      END;
      ps.close()  
    END   
  END DrawCandidates;
 
   
PROCEDURE DrawCandidate(
   READONLY ps: PSPlot.PSFile;
   READONLY c: PZLR3Chain.T;
   ini, fin: INT;
  )=
  BEGIN
    Wr.PutText(stderr, "range = [" & Fmt.Int(ini) & ".." & Fmt.Int(fin) & "]\n");
    WITH
      n = NUMBER(c),
      in = ini MOD n,
      fn = fin MOD n
    DO
      ps.setLineSolid();
      ps.setFillColor(PSPlot.Invisible);
      IF in < fn THEN 
        ps.setLineWidth(0.1);
        ps.setLineColor(PSPlot.Black);
        PZPlot.LR3Chain(ps, SUBARRAY(c, 0, in+1), FALSE);
        PZPlot.LR3Chain(ps, SUBARRAY(c, fn, n-fn), FALSE);
        ps.setLineWidth(0.25);
        ps.setLineColor(PSPlot.Red);
        PZPlot.LR3Chain(ps, SUBARRAY(c, in, fn-in+1), FALSE);
        ps.segment(c[0][0], c[0][1], c[n-1][0], c[n-1][1]);
      ELSE
        ps.setLineWidth(0.1);
        ps.setLineColor(PSPlot.Black);
        PZPlot.LR3Chain(ps, SUBARRAY(c, fn, in-fn+1), FALSE);
        ps.setLineWidth(0.25);
        ps.setLineColor(PSPlot.Red);
        PZPlot.LR3Chain(ps, SUBARRAY(c, 0, fn+1), FALSE);
        PZPlot.LR3Chain(ps, SUBARRAY(c, in, n-in), FALSE);
        ps.segment(c[0][0], c[0][1], c[n-1][0], c[n-1][1]);
      END;
    END
  END DrawCandidate;

PROCEDURE ReadAllCVCFiles(
    inName : TEXT;
    nChains, band: NAT; 
    VAR step: LONG;
  ): REF ARRAY OF REF PZSymbolChain.T =
  VAR epsilon, delta: LONG;
  <* FATAL Rd.Failure *>
  BEGIN 
   WITH 
     chain = NEW(REF ARRAY OF REF  PZSymbolChain.T, nChains )
   DO 
     FOR i := 0 TO nChains-1 DO  
        WITH 
          chainTag  = Fmt.Pad(Fmt.Int(i), 4, '0'),
          bandTag = Fmt.Pad(Fmt.Int(band), 3, '0'),
          fileName = inName & "-" & chainTag & "-f" & bandTag & "-c-u.cvc", 
          rd = FileRd.Open(fileName),
          data = PZSymbolChain.Read(rd)
        DO
          chain[i] := data.c;
          IF i = 0 THEN 
            step := data.length / FLOAT(data.samples, LONG);
            epsilon := data.epsilon; 
            delta := data.delta
          ELSE
            <* ASSERT epsilon = data.epsilon *>
            <* ASSERT delta = data.delta *>
          END;
          Rd.Close(rd);
        END;
      END;
     RETURN chain;
    END;  
  END ReadAllCVCFiles;
  
  
PROCEDURE ReadAllTMCFiles(
    inName : TEXT;
    nChains, band: NAT; 
  ): REF ARRAY OF REF PZLRChain.T = 
  <* FATAL Rd.Failure *>
  BEGIN  
    WITH chain = NEW(REF ARRAY OF REF  PZLRChain.T, nChains )
    DO 
      FOR i := 0 TO nChains-1 DO  
        WITH 
          chainTag  = Fmt.Pad(Fmt.Int(i), 4, '0'),
          bandTag = Fmt.Pad(Fmt.Int(band), 3, '0'),
          fileName = inName & "-" & chainTag & "-f" & bandTag & ".tmc",  
          rd =FileRd.Open(fileName),
          data = PZLRChain.Read(rd)
        DO
          chain[i] := data.c;
          Rd.Close(rd); 
        END;
      END;
      RETURN chain;
    END;
  END  ReadAllTMCFiles; 
  
  
PROCEDURE ReadAllFLCFiles(
    inName : TEXT;
    nChains, band: NAT; 
  ): REF ARRAY OF REF PZLR3Chain.T =
  <* FATAL Rd.Failure *>
  BEGIN 
   WITH 
     chain = NEW(REF ARRAY OF REF PZLR3Chain.T, nChains)
   DO 
     FOR i := 0 TO nChains-1 DO  
       WITH 
         chainTag  = Fmt.Pad(Fmt.Int(i), 4, '0'),
         bandTag = Fmt.Pad(Fmt.Int(band), 3, '0'),
         fileName = inName & "-" & chainTag & "-f" & bandTag & "-c-u.flc",
         rd = FileRd.Open(fileName),
         data = PZLR3Chain.Read(rd)
       DO
         Wr.PutText(stderr, fileName & "\n");
         chain[i] := data.c;
         Rd.Close(rd);
       END;
     END;
     RETURN chain;
    END;  
  END ReadAllFLCFiles;
 
 
BEGIN
  Main()
END PZMultiScale.