MODULE TestCurve EXPORTS Main;

(* Tests the "PZCurve" module. *)
(* Last edited on 1999-09-04 17:41:01 by hcgl *)

IMPORT Math, PSPlot, Wr, Fmt, LR3, Thread, Random;
IMPORT PZTypes, PZCurve;
FROM Stdio IMPORT stderr;
FROM PZTypes IMPORT LONG, LONGS, NAT, BOOL;

TYPE 
  Generator = PROCEDURE(t: LONG): LR3.T;
  Times = LONGS;
  Labels = LONGS;
  Points = ARRAY OF LR3.T;

CONST TwoPi = 2.0d0*3.141592653589793d0;

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

PROCEDURE Main() =
  BEGIN
    FOR k := 0 TO 2 DO
      FOR tJitter := FALSE TO TRUE DO
        FOR pJitter := FALSE TO TRUE DO
          DoTests(
            Circle, "Circle",
            uniform := NOT tJitter,
            exact := TwoPi,
            perturb := pJitter
          );
        END
      END
    END
  END Main;
  
PROCEDURE DoTests(
    generator: Generator; curveName: TEXT;
    uniform: BOOL;
    exact: LONG;
    perturb: BOOL;
  ) =
  VAR n: NAT;
  BEGIN
    WITH
      uniformName = ARRAY BOOL OF TEXT{"jittered", "uniform"}[uniform],
      perturbName = ARRAY BOOL OF TEXT{"exact", "perturbed"}[perturb],
      testName = curveName & "-" & uniformName & "-" & perturbName,
      noise = ARRAY BOOL OF LONG{0.00d0, 0.01d0}[perturb]
    DO
      Wr.PutText(stderr, "generator = " & curveName & "\n");
      Wr.PutText(stderr, uniformName & "\n");
      Wr.PutText(stderr, perturbName & "\n");
      Wr.PutText(stderr, 
        "exact length = " & 
        Fmt.Pad(Fmt.LongReal(exact, Fmt.Style.Fix, 6), 12) & "\n"
      );
      Wr.PutText(stderr, 
        "position noise = " & 
        Fmt.Pad(Fmt.LongReal(noise, Fmt.Style.Fix, 6), 12) & "\n"
      );
      Wr.PutText(stderr, "\n");
      n := 128;
      WHILE n > 0 DO
        WITH
          t = PickTimes(n, 1.0d0, uniform)^,
          p = PickPoints(generator, t, noise)^,
          c = PZCurve.New(n)
        DO
          c.setSamples(p);
          c.setTimes(t, 1.0d0);
          c.setLabels(t, 1.0d0);
          PrintTest(c);
          IF n = 8 THEN 
            WITH runName = testName & "-" & Fmt.Pad(Fmt.Int(n), 3, '0') DO
              PlotTest(runName, c, 100)
            END
          END;
        END;
        n := n DIV 2
      END;
    END
  END DoTests;

PROCEDURE PickTimes(n: NAT; T: LONG; uniform: BOOL): REF Times =
  (*
    Returns "n" increasing times equally spaced in "[0..T]" cyclic. *)
  BEGIN
    WITH
      rt = NEW(REF Times, n), t = rt^,
      step = T/FLOAT(n, LONG),
      rnd = NEW(Random.Default).init(fixed := TRUE)
    DO
      FOR i := 0 TO n-1 DO 
        t[i] := step*(FLOAT(i, LONG) + 0.5d0);
        IF NOT uniform THEN 
          t[i] := t[i] + step * rnd.longreal(-0.499d0, +0.499d0)
        END
      END;
      RETURN rt
    END
  END PickTimes;

PROCEDURE PickPoints(
    generator: Generator;
    READONLY t: Times; 
    noise: LONG;
  ): REF Points =
  (*
    Samples the generator "generator(t)" at times "t[0..LAST]",
    perturbs the samples with the given amount of noise. *)
  BEGIN
    WITH
      n = NUMBER(t),
      rp = NEW(REF Points, n), p = rp^,
      rnd = NEW(Random.Default).init(fixed := TRUE)
    DO
      FOR i := 0 TO n-1 DO 
        p[i] := LR3.Add(generator(t[i]), Jitter(rnd, noise))
      END;
      RETURN rp
    END;
  END PickPoints;
  
PROCEDURE Jitter(rnd: Random.T; noise: LONG): LR3.T =
  VAR v: LR3.T;
  BEGIN
    v[2] := 0.0d0;
    REPEAT
      v[0] := rnd.longreal(-1.0d0, +1.0d0);
      v[1] := rnd.longreal(-1.0d0, +1.0d0);
    UNTIL LR3.NormSqr(v) <= 1.0d0;
    v[0] := v[0] * noise;
    v[1] := v[1] * noise;
    v[2] := v[2] * noise;
    RETURN v
  END Jitter;

PROCEDURE PrintTest(
    READONLY c: PZCurve.T;  (* Curve to test. *)
  ) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    c. uniformize();
    WITH
      m = c.m,
      p = c.p^,
      pLen = PolyLength(p),
      cLen = c.tPeriod,
      mLen = 0.5d0 * (pLen + cLen)
    DO
      Wr.PutText(stderr, "m = " & Fmt.Pad(Fmt.Int(m), 5));
      Wr.PutText(stderr, "  curve length = " & 
        Fmt.Pad(Fmt.LongReal(cLen, Fmt.Style.Fix, 6), 12)
      );
      Wr.PutText(stderr, "  poly length = " & 
        Fmt.Pad(Fmt.LongReal(pLen, Fmt.Style.Fix, 6), 12)
      );
      Wr.PutText(stderr, "  average = " & 
        Fmt.Pad(Fmt.LongReal(mLen, Fmt.Style.Fix, 6), 12)
      );
      Wr.PutText(stderr, "\n")
    END
  END PrintTest;
  
PROCEDURE PlotTest(
    testName: TEXT;
    c: PZCurve.T;           (* The curve. *)
    ns: NAT;
  ) =
  BEGIN
    c.uniformize();
    WITH
      q = NEW(REF Points, ns)^,
      dq = NEW(REF Points, ns)^,
      r = NEW(REF Labels, ns)^,
      t = PickTimes(ns, c.tPeriod, uniform := TRUE)^,
      f = NEW(PSPlot.EPSFile).open(testName & ".eps")
    DO
      c.sample(t, q, dq, r);
      f.setScale(PSPlot.Axis.X, -2.0d0, +2.0d0);
      f.setScale(PSPlot.Axis.Y, -2.0d0, +2.0d0);
      f.setFillColor(PSPlot.Invisible);
      FOR i := 0 TO ns-1 DO 
        WITH a = q[i], b = q[(i+1) MOD ns] DO 
          f.segment(a[0], a[1], b[0], b[1])
        END;
      END;
      FOR i := 0 TO c.m-1 DO 
        WITH a = c.p[i], d = c.d[i] DO 
          f.segment(a[0], a[1], a[0]+d[0], a[1]+d[1])
        END;
      END;
      f.setLineColor(PSPlot.Red);
      FOR i := 0 TO c.m-1 DO 
        WITH a = c.p[i] DO 
          f.dot(a[0], a[1], 2.0)
        END;
      END;
      f.close()
    END
  END PlotTest;
  
PROCEDURE PolyLength(
    READONLY p: Points;     (* Corresponding sample points *)
  ): LONG =
  (*
    The length of the polygonal "p[0..LAST]". *)
  VAR s: LONG := 0.0d0;
  BEGIN
    WITH ns = NUMBER(p) DO
      FOR i := 0 TO ns-1 DO 
        WITH 
          Pa = p[i], 
          Pb = p[(i+1) MOD ns], 
          ds = LR3.Dist(Pa, Pb)
        DO
          s := s + ds
        END;
      END;
    END;
    RETURN s
  END PolyLength;

PROCEDURE Circle(t: LONG): LR3.T =
  BEGIN
    WITH
      theta = TwoPi*t,
      c = Math.cos(theta),
      s = Math.sin(theta)
    DO
      RETURN LR3.T{c, s, 0.0d0}
    END
  END Circle;

BEGIN
  Main()
END TestCurve.