MODULE PZSymbolChain;

IMPORT Rd, Wr, Thread, FileFmt, NPut, NGet, FPut, FGet, Fmt, OSError, FileRd;
IMPORT PZMatch, PZProc, PZSymbol;
FROM Math IMPORT sqrt;
FROM Stdio IMPORT stderr;

  
PROCEDURE Trim(READONLY c: T; start, length: CARDINAL): REF T =
  BEGIN 
    WITH 
      N = NUMBER(c),
      ini = start MOD N,
      fin = ini + length - 1
    DO 
      WITH s = NEW(REF T, length) DO
        IF fin < N  THEN 
          s^ := SUBARRAY(c, ini, length)
        ELSE 
          SUBARRAY(s^, 0, N-ini) := SUBARRAY(c, ini, N-ini);
          SUBARRAY(s^, N-ini, (fin+1) MOD N):= SUBARRAY(c, 0, (fin+1) MOD N);
        END;
        RETURN s
      END;
    END;
  END Trim;


PROCEDURE ReverseAndComplement(VAR c: T) =
  VAR aux: PZSymbol.T;
  BEGIN 
    WITH N = NUMBER(c), L = N-1 DO 
      FOR i := 0 TO (N DIV 2)-1  DO 
        aux := c[i]; c[i] := c[L-i]; c[L-i] := aux
      END;
      FOR i := 0 TO N-1 DO c[i] := PZSymbol.Complement(c[i]) END
    END;  
  END ReverseAndComplement; 
 
CONST FileVersion = "97-02-03";

PROCEDURE Write(
    wr :Wr.T; 
    cmt: TEXT;
    READONLY c: T;
    length: LONGREAL;
    epsilon, delta: LONGREAL;
  ) =
   <* FATAL Thread.Alerted, Wr.Failure *>
  BEGIN
    WITH  
      n = LAST(c)
    DO 
      FileFmt.WriteHeader(wr, "PZSymbolChain.T", FileVersion);
      FileFmt.WriteComment(wr, cmt, '|');
      NPut.Int(wr, "samples", NUMBER(c));        FPut.EOL(wr);
      NPut.LongReal(wr, "length", length);       FPut.EOL(wr);
      NPut.LongReal(wr, "epsilon", epsilon);     FPut.EOL(wr);
      NPut.LongReal(wr, "delta", delta);         FPut.EOL(wr);
      FOR i := 0 TO n  DO 
        FPut.Char(wr, c[i]); FPut.EOL(wr)
      END;
      FPut.EOL(wr);
      FileFmt.WriteFooter(wr, "PZSymbolChain.T");
      Wr.Flush(wr);
     END;
  END Write;  


PROCEDURE Read(rd :Rd.T; headerOnly: BOOLEAN := FALSE): ReadData  =
  VAR d: ReadData;
  BEGIN
    FileFmt.ReadHeader(rd,"PZSymbolChain.T", FileVersion);
    d.cmt := FileFmt.ReadComment(rd, '|');
    d.samples := NGet.Int(rd, "samples"); FGet.EOL(rd);
    WITH 
      n = d.samples 
     DO 
      d.length := NGet.LongReal(rd, "length"); FGet.EOL(rd);
      d.epsilon := NGet.LongReal(rd, "epsilon"); FGet.EOL(rd);
      d.delta := NGet.LongReal(rd, "delta"); FGet.EOL(rd);
      IF headerOnly THEN 
        d.c := NIL
      ELSE
        d.c := NEW(REF T, n);
        FOR i := 0 TO n-1 DO 
          FGet.Skip(rd, FGet.Blanks);
          d.c[i] := FGet.Char(rd);
        END;
        FileFmt.ReadFooter(rd, "PZSymbolChain.T");
      END;
      RETURN d;
    END;
  END Read;  


PROCEDURE ReadAll(
    puzzle: TEXT; 
    band: CARDINAL;
    extension: TEXT;
    READONLY sel: ARRAY OF BOOLEAN;
    headerOnly: BOOLEAN := FALSE;
  ): REF ARRAY OF ReadData =
  CONST 
    NoData = ReadData{
      cmt := "", samples := 0, length := 0.0d0,
      epsilon := 0.0d0, delta := 0.0d0, c := NIL
    };
  <* FATAL Wr.Failure, Thread.Alerted, Rd.Failure *>
  BEGIN
    Wr.PutText(stderr, "reading symbolic chains:\n");
    WITH 
      rr = NEW(REF ARRAY OF ReadData, NUMBER(sel)), r = rr^
    DO
      FOR k := 0 TO LAST(sel) DO
        IF sel[k] THEN 
          WITH 
            chainTag  = Fmt.Pad(Fmt.Int(k), 4, '0'),
            bandTag = Fmt.Pad(Fmt.Int(band), 3, '0'),
            fileName =
              chainTag & "/" & puzzle & "-" & chainTag & "-f" & bandTag & extension,
            rd = OpenRead(fileName)
          DO
            IF rd # NIL THEN
              Wr.PutText(stderr, fileName & "\n");
              r[k] := Read(rd, headerOnly);
              Rd.Close(rd)
            ELSE
              r[k] := NoData
            END
          END  
        ELSE
          r[k] := NoData
        END
      END;
      RETURN rr
    END
  END ReadAll;

PROCEDURE OpenRead(name: TEXT): Rd.T =
  BEGIN
    TRY
      RETURN FileRd.Open(name)
    EXCEPT
      OSError.E => RETURN NIL
    END
  END OpenRead;

PROCEDURE Match(
    READONLY a: T;
    READONLY b: T;
    maxDist: LONGREAL;
    removeUnmatchedEnds: BOOLEAN;
    VAR (*OUT*) avgCost: LONGREAL;
    VAR (*OUT*) nMatched: CARDINAL; 
    VAR (*OUT*) match: REF PZMatch.T;
  ) =
  
  TYPE Pair = PZMatch.Pair;
  
  BEGIN

    PROCEDURE StepCostSqr(READONLY p, q: Pair): LONGREAL =
      BEGIN
        IF p[0] = q[0] OR p[1] = q[1] THEN 
          RETURN LAST(LONGREAL)
        ELSE
          RETURN 
            PZSymbol.IntgDistSqr(
              a1 := a[p[0]], a2 := a[q[0]],  
              b1 := b[p[1]], b2 := b[q[1]],  
              complement := FALSE
            )
        END
      END StepCostSqr;

    VAR totCostSqr: LONGREAL;
    BEGIN
      WITH
        NA = NUMBER(a),
        NB = NUMBER(b),
        maxDistSqr = maxDist*maxDist
      DO
        PZMatch.Match(NA, NB, 
          stepCost := StepCostSqr, 
          maxCost := maxDistSqr, 
          (*OUT*) 
          totCost := totCostSqr, 
          match := match
        );
        IF removeUnmatchedEnds THEN
          match := PZMatch.RemoveUnmatchedEnds(match^, StepCostSqr, maxDistSqr);
          totCostSqr := PZMatch.MatchCost(match^, StepCostSqr, maxDistSqr)
        END;
        WITH n = NUMBER(match^) DO
          <* ASSERT n > 0 *>
          IF n = 1 THEN
            avgCost := maxDist
          ELSE
            avgCost := sqrt(totCostSqr/FLOAT(n-1, LONGREAL))
          END
        END;
        nMatched := PZMatch.NMatched(match^, StepCostSqr, maxDistSqr)
      END;
    END
  END Match;


PROCEDURE OutMatch(
    READONLY a: T;
    READONLY b: T;
    loc: LONGREAL;
    align: LONGREAL;
    maxShift: LONGREAL;
    maxDist: LONGREAL;
    cutDist: LONGREAL; 
    step: LONGREAL;
    minChainSteps: CARDINAL;
    VAR (*OUT*) mismatch: LONGREAL;
    VAR (*OUT*) nMatched: CARDINAL; 
    VAR (*OUT*) match: REF PZMatch.T;
    VAR (*OUT*) xCenter, yCenter: CARDINAL;
    VAR (*OUT*) minAvgDist: REF ARRAY OF LONGREAL;
    VAR (*WORK*) XYCost: REF CostMatrix;
  ) =
  
  TYPE Pair = PZMatch.Pair;
  
  BEGIN

    WITH 
      maxDistSqr = maxDist*maxDist
    DO

      PROCEDURE StepCostSqr(READONLY p, q: Pair): LONGREAL =
        BEGIN
          IF p[0] = q[0] OR p[1] = q[1] THEN 
            RETURN LAST(LONGREAL)
          ELSE
            RETURN 
              PZSymbol.IntgDistSqr(
                a1 := a[p[0]], a2 := a[q[0]],  
                b1 := b[p[1]], b2 := b[q[1]],  
                complement := FALSE
              )
          END
        END StepCostSqr;
        
      PROCEDURE ShiftCostSqr(shift: INTEGER): LONGREAL =
        BEGIN
          IF FLOAT(ABS(shift), LONGREAL) > maxShift THEN
            RETURN LAST(LONGREAL)
          ELSE
            RETURN 0.0d0
          END;
        END ShiftCostSqr;

      PROCEDURE ComputeMismatch(
          totCostSqr: LONGREAL; 
          totChainSteps: CARDINAL;
        ): LONGREAL =
        BEGIN
          IF totChainSteps < minChainSteps THEN
            RETURN LAST(LONGREAL)
          ELSE
            WITH 
              length = step * FLOAT(totChainSteps, LONGREAL) / 2.0d0,
              intgCostSqr = step * totCostSqr
            DO
              RETURN PZProc.Mismatch(intgCostSqr, length, cutDist*cutDist)
            END
          END;
        END ComputeMismatch;
        
      BEGIN
        WITH
          NA = NUMBER(a),
          NB = NUMBER(b),
          locInt = ROUND(loc),
          alignInt = ROUND(align),
          maxShiftInt = CEILING(maxShift/2.0d0) * 2
        DO
          <* ASSERT NA > 0 *>
          <* ASSERT NB > 0 *>
          PZMatch.OutMatch(
            NA, NB, 
            stepCost := StepCostSqr, 
            maxCost := maxDistSqr,
            L := locInt,
            A := (alignInt DIV 2) * 2 + (locInt MOD 2),
            shiftMin := -maxShiftInt,
            shiftMax := +maxShiftInt,
            shiftCost := ShiftCostSqr,
            mismatch := ComputeMismatch,
            (*OUT*) mBest := mismatch,
            (*OUT*) match := match,
            (*OUT*) xCenter := xCenter,
            (*OUT*) yCenter := yCenter,
            (*OUT*) minTotCost := minAvgDist,
            (*WORK*) XYCost := XYCost
          );
          <* ASSERT NUMBER(match^) > 0 *>
          nMatched := PZMatch.NMatched(match^, StepCostSqr, maxDistSqr);
          <* ASSERT  minAvgDist[0] = 0.0d0 *>
          FOR k := 1 TO LAST(minAvgDist^) DO 
            minAvgDist[k] := sqrt(minAvgDist[k]/(FLOAT(k, LONGREAL)/2.0d0))
          END;
        END;
      END
    END
  END OutMatch;

BEGIN
END PZSymbolChain.