MODULE MREvolver2D EXPORTS Main;

IMPORT Texture2DGray AS Tx, FeatureMapper2D AS Mapper, PSPlot;
IMPORT CsiFeatures2D, CrossFeatures2D, PhiFeatures2D, 
       FunGradFeatures2D, FunSlopeFeatures2D;
IMPORT ParseParams, Process, Wr, Thread, Text, Fmt, FileWr, OSError;
IMPORT GaussRandom, LR2;
IMPORT Random, Math;
FROM Stdio IMPORT stderr;

CONST
  Debug = FALSE;
  SqrtH = 0.70710678118654752440d0;

TYPE
  LONG = LONGREAL;
  LONGS = ARRAY OF LONG;

  NAT = CARDINAL;
  INT = INTEGER;
  BOOL = BOOLEAN;

  Pixel = Tx.Pixel;
  LRPixel = Tx.LongRealPixel;
  
  ThickTx = ARRAY OF Tx.T;
  MaskedTx = ARRAY OF Tx.T;  (* Last layer is a boolean mask *)

TYPE
  Options = RECORD
      inName: TEXT;    (* Sample texture filename prefix. *)
      outName: TEXT;   (* Output filename prefix. *)
      random: LONG;    (* Amount of random noise to inject in each layer *)
      fold: BOOL;      (* TRUE constrains pixels to [-1 _ +1] *)
      norm: BOOL;      (* TRUE to normalize the gradient locally *)
      UX, VY: NAT;     (* Synthetic texture size *)
      VX: NAT;         (* Brick-row displacement of synthetic texture *)
      map: Mapper.T;   (* The local feature mapper *)
      cmp: Comparator; (* The feature comparator *)
      maxEvals: NAT;   (* Max number of evals at finer scale *)
      initStep: LONG;  (* Initial step length *)
      fuzz: LONG;      (* Fuzz constant for charge distribution *)
    END;

  Nhood = Tx.PixelNeighborhood;
  Nhoods = ARRAY OF Nhood;

  Features = ARRAY OF LONG;

  FeatureDistr = ARRAY OF Features;

  Comparator = OBJECT
      NF: NAT;           (* Number of local features *)
    METHODS
      print(wr: Wr.T);
        (*
          Prints a description of the comparator to "wr". *)
          
      eval(
          READONLY A: Features;
          qA: LONG;
          VAR dEdA: Features;
          READONLY B: Features;
          qB: LONG;
          VAR dEdB: Features;
          VAR E: LONGREAL;
        );
        (*
          Adds to "E" the energy due to the charge pair "(A,qA)" and "(B,qB)".
          Also adds to "dEdA", "dEdB" the corresponding derivatives. *)
    END;  
    
  TxEvaluator = OBJECT
      NL: NAT;        (* Number of layers expected by "eval" and "grad" *)
      flat: BOOL;     (* TRUE if evaluator is degenerate (constant) *)
    METHODS
      print(wr: Wr.T);
        (*
          Prints a description of the evaluator to "wr". *)
 
      eval(READONLY A: ThickTx; (*OUT*) VAR E: LONG; fName: TEXT := "");
        (*
          Evaluates texture "A", returning its energy in "E".
          If "fName" is not empty, writes debugging info
          to a file called "fName". *)

      grad(READONLY A: ThickTx; (*OUT*) READONLY dEdA: ThickTx);
        (*
          Returns in "dEdA" the gradient of the energy, computed 
          for the texture "A" --- which must have been the argument
          of the last "eval". *)
    END;

CONST
  LayerName = ARRAY OF TEXT{"A", "B", "C", "D", "E", "F", "G", "H"};
  ZeroPixel = LRPixel{0.0d0, ..};
  UnitPixel = LRPixel{1.0d0, ..};

PROCEDURE Main() =
  BEGIN
    WITH
      o = GetOptions(),
      S = ReadSample(o.inName, o.map.NL)^,
      A = MakeThickTx(o.map.NL, o.UX, o.VX, o.VY)^,
      rnd = NEW(Random.Default).init(fixed := TRUE)
    DO
      <* ASSERT Tx.NChannels = 1 *>
      (* Generate synthetic texture: *)
      EVAL GenTexture(
        (*OUT*) A, (* IN *) S, o, rnd, o.maxEvals, 
        logScale := 0, xOrg := -0.5d0, yOrg := -0.5d0
      );
      (* Write final state and derivatives: *)
      WriteComments(o.outName, o);
      WriteTextures(A, o.map.NL, o.outName & "-final", wd := o.UX, ht := o.VY);
    END
  END Main;

PROCEDURE ReadSample(
    name: TEXT; 
    NL: NAT;
  ): REF MaskedTx =
  (*
    Reads an "NL"-layered image to use as the sample. The image is
    stored as a texture only for technical reasons (a texture is
    easier to shrink by 1/2 than an image).
    
    Layers 0 thru "NL-1" are the image proper. Layer "NL" is a mask that
    is 0 ("undefined") along the top and left borders, and 1 
    ("defined") elsewhere.  Note that the top row and rightmost column
    of the image are lost.
  *)
  BEGIN
    WITH
      S = NEW(REF MaskedTx, NL+1)
    DO
      FOR i := 0 TO NL-1 DO
        WITH
          fName = name & "-" & LayerName[i] & ".pgm"
        DO
          S[i] := Tx.Read(fName, VX := 0, a := ZeroPixel, s := UnitPixel);
          IF i > 0 THEN
            <* ASSERT S[i].UX = S[0].UX *>
            <* ASSERT S[i].VY = S[0].VY *>
            <* ASSERT S[i].VX = S[0].VX *>
          END
        END
      END;
      S[NL] := SampleMask(S[0].UX, S[0].VX,  S[0].VY);
      RETURN S
    END;
  END ReadSample;

PROCEDURE MakeThickTx(NL, UX, VX, VY: NAT): REF ThickTx =
  BEGIN
    WITH
      A = NEW(REF ThickTx, NL)
    DO
      FOR i := 0 TO NL-1 DO
        A[i] := Tx.New(UX, VX, VY, ZeroPixel, ZeroPixel);
      END;
      RETURN A
    END;
  END MakeThickTx;

PROCEDURE SampleMask(UX, VX, VY: NAT): Tx.T =
  BEGIN
    WITH
      M = Tx.New(UX, VX, VY, ZeroPixel, ZeroPixel)
    DO
      FOR y := 0 TO VY-1 DO
        FOR x := 0 TO UX-1 DO
          IF x = UX-1 OR y = VY-1 THEN
            M.st[y,x] := 0.0
          ELSE
            M.st[y,x] := 1.0
          END
        END
      END;
      RETURN M
    END;
  END SampleMask;

PROCEDURE WriteComments(fName: TEXT; READONLY o: Options) =
  <* FATAL Thread.Alerted, Wr.Failure, OSError.E *>
  BEGIN
    WITH
      wr = FileWr.Open(fName & ".txt")
    DO
      Wr.PutText(wr, "--- feature mapper ----------------------------------\n");
      o.map.print(wr);
      Wr.PutText(wr, "--- distribution comparator -------------------------\n");
      o.cmp.print(wr);
      Wr.PutText(wr, "--- other parameters --------------------------------\n");
      Wr.PutText(wr, "  maxEvals = " & Fmt.Int(o.maxEvals) & "\n");
      Wr.PutText(wr, "  initStep = " & FLR(o.initStep, 6) & "\n");
      Wr.PutText(wr, "  random   = " & FLR(o.random, 6) & "\n");
      Wr.PutText(wr, "  fold     = " & FB(o.fold) & "\n");
      Wr.PutText(wr, "  norm     = " & FB(o.norm) & "\n");
      Wr.PutText(wr, "\n");
      FOR i := 0 TO o.map.NL-1 DO
        Wr.PutText(wr, "  layer " & LayerName[i] & " range = [-1 _ +1]\n");
      END;
      Wr.Close(wr)
    END
  END WriteComments;

PROCEDURE WriteTextures(
    READONLY A: ThickTx; 
    NL: NAT;
    name: TEXT;
    wd, ht: NAT;
    logScale: NAT := 0;
    xOrg, yOrg: LONG := -0.5d0;
  ) =
  (*
    Writes each layer of texture "A" as 
      (1) a simple grayscale image, in a format suitable for reading it back; 
      (2) as a greyscale image, magnified "2^logScale" times in area 
          and replicated twice each direction.
      (3) a Postscript file, where the texture is rotated by "logScale"
          times 45 degrees, magnified by "sqrt(2)^logScale",
          and clipped to a rectangle of size "wd" by "ht".
    
    For layers 0 thru "NL-1", -1.5 is mapped to black, +1.5 to white.
    For layer "NL", 0 is mapped to black, 1 to white.
  *)
  VAR M: Tx.T;
      a, s: LRPixel;
  CONST BasePixelSize = 1.0d0;  (* Unmagnified pixel size, mm. *)
  BEGIN
    WITH
      masked = (NL < NUMBER(A)),
      mag = Pow(2, logScale DIV 2),
      extraX = ORD(A[0].UX <= A[0].VY),
      replX = 1 + extraX * (logScale MOD 2),
      replY = 1 + (1 - extraX) * (logScale MOD 2),
      step = Dir45(logScale),
      pixelScale = Math.pow(2.0d0, FLOAT(logScale, LONG)/2.0d0),
      pixelSize = BasePixelSize * pixelScale,
      xOrg = xOrg, yOrg = yOrg
    DO
      IF masked THEN M := A[NL] ELSE M := Tx.Null END;
      FOR i := 0 TO LAST(A) DO
        IF i < NL THEN
          a := LRPixel{0.0d0, ..};
          s := LRPixel{1.0d0, ..};
        ELSE
          a := LRPixel{0.5d0, ..};
          s := LRPixel{0.5d0, ..};
        END;
        Tx.Write(
          name & "-" & LayerName[i] & "-r.pgm",
          A[i], a := a, s := s,
          NX := replX, NY := replY, mag := mag
        );
        Tx.Write(
          name & "-" & LayerName[i] & "-u.pgm",
          A[i], a := a, s := s,
          NX := 1, NY := 1, mag := 1
        );
        IF i # NL THEN 
          WritePSLayer(
            A[i], M,
            name := name & "-" & LayerName[i] & ".eps", 
            cos := step[0], sin := step[1],
            xOrg := xOrg - 0.25d0 * (step[0] + step[1]), 
            yOrg := yOrg - 0.25d0 * (step[0] - step[1]),
            wd := FLOAT(wd, LONG)/pixelScale + 0.5d0, 
            ht := FLOAT(ht, LONG)/pixelScale + 0.5d0,
            pixelSize := pixelSize,
            black := LRPixel{-1.0d0, ..}, 
            white := LRPixel{+1.0d0, ..}
          )
        END
      END;
    END;
  END WriteTextures;

PROCEDURE WritePSLayer(
    READONLY t: Tx.T;   (* Texture layer *)
    READONLY m: Tx.T;   (* Mask or "Tx.Null" *)
    name: TEXT;         (* File name prefix *)
    xOrg, yOrg: LONG;   (* Texture point at lower left corner of picture. *)
    wd, ht: LONG;       (* Picture dimensions in texture pixels. *)
    cos, sin: LONG;     (* Cosine and sine of texture rotation. *)
    pixelSize: LONG;    (* Pixel side in mm *)
    black: LRPixel;     (* Pixel value to map to black. *)
    white: LRPixel;     (* Pixel value to map to white. *)
  ) =
  (*
    Writes an Encapsulated Postscript description of a texture layer "t".
    The parameters "name", "xOrg", "yOrg", "wd", "ht", "cos",
    "sin", and "pixelSize" are simply passed on to "Tx.WritePS". 
    
    The pixel colors are 
    
    The texture "m" is interpreted as a binary mask, that tells where
    the texture 't" is defined ("m=1") or undefined ("m=0").
  *)
  BEGIN
    <* ASSERT Tx.NChannels = 1 *>
    WITH
      masked = (m.st # NIL),
      black = black[0],
      white = white[0],
      pixelScale = 1.0d0/(white - black)
    DO

      PROCEDURE Exists(x, y: INT): BOOL = 
        (*
          Given reduced pixel indices "x" and "y", returns TRUE iff
          the pixel exists. *)
        VAR c: Pixel;  
        BEGIN
          IF NOT masked THEN
            RETURN TRUE
          ELSE
            Tx.Reduce(m, x, y);
            <* ASSERT 0 <= x AND x < m.UX *>
            <* ASSERT 0 <= y AND y < m.VY *>
            Tx.GetPixel(m, x, y, c);
            RETURN c > 0.5
          END
        END Exists;

      PROCEDURE PixelToGray(READONLY p: Pixel): Tx.Color = 
        BEGIN
          WITH
            p = FLOAT(p, LONG),
            g = MAX(0.0d0, MIN(1.0d0, pixelScale * (p - black)))
          DO
            RETURN PSPlot.Gray(FLOAT(g))
          END;
        END PixelToGray;

      BEGIN 
        Tx.WritePS(
          name := name, 
          t := t,
          xOrg := xOrg, yOrg := yOrg,
          wd := wd, ht := ht,
          cos := cos, sin := sin,
          pixelSize := pixelSize,
          pixelColor := PixelToGray,
          outline := TRUE,
          defined := Exists
        )
      END
    END;
  END WritePSLayer;
  
PROCEDURE Dir45(rot: NAT): LR2.T =
  BEGIN
    CASE rot MOD 8 OF
    | 0 => RETURN LR2.T{+1.0d0, 00.0d0};
    | 1 => RETURN LR2.T{+SqrtH, +SqrtH};
    | 2 => RETURN LR2.T{00.0d0, +1.0d0};
    | 3 => RETURN LR2.T{-SqrtH, +SqrtH};
    | 4 => RETURN LR2.T{-1.0d0, 00.0d0};
    | 5 => RETURN LR2.T{-SqrtH, -SqrtH};
    | 6 => RETURN LR2.T{00.0d0, -1.0d0};
    | 7 => RETURN LR2.T{+SqrtH, -SqrtH};
    ELSE <* ASSERT FALSE *>
    END
  END Dir45;

PROCEDURE GenTexture(
    READONLY A: ThickTx;   (* OUT: Synthetic texture. *)
    READONLY S: MaskedTx;  (* IN: Sample and mask. *)
    READONLY o: Options;   (* IN: Command line parameters *)
    rnd: Random.T;         (* I/O: Randomness source *)
    topEvals: NAT;         (* IN: Max energy evaluations in top layer *)
    logScale: NAT;         (* IN: Log2(original_area/current_area) *)
    xOrg, yOrg: LONG;      (* IN: reference point for texture plotting *)
  ): BOOLEAN =
  (*
    
    Fills "A" with a synthetic texture resembling "S", by multiscale
    annealing.  Returns TRUE on success, FALSE on failure --- in
    particular, if "S" has no pixels.
    
    Texture "S" should have "map.NL+1" layers.  Layers 0 thru
    "NL-1" are the sample proper, and correspond to layers of "A";
    layer "NL" is interpreted as a mask that is 1 where "S" is
    defined, 0 where "S" is undefined.
    
    Two textures are assumed to be similar if they have similar
    distributions of local statistics, as computed by the 
    featre mapper "o.map" and compared by "o.cmp".
    
    For output purposes, the textures "A" and "S" are assumed to 
    have been shrunk by a factor of "2^logScale" in area.
  *)
  BEGIN
    WITH
      NL = NUMBER(A),
      AUX = A[0].UX, AVX = A[0].VX, AVY = A[0].VY,
      txName = o.outName & "-" & FI(logScale,2,'0'),
      V = NEW(StatEvaluator).init(
        S, AUX, AVX, AVY, 
        o.map, o.cmp, o.fold, o.norm, 
        fName := txName & "-sampl"
      )
    DO
      <* ASSERT NUMBER(S) = NL+1 *>
      (* Write sample texture: *)
      WriteTextures(
        S, NL, txName & "-sampl", 
        wd := o.UX, ht := o.VY, logScale := logScale,
        xOrg := xOrg, yOrg := yOrg
      );
      (* Compute size of reduced texture and sample: *)
      IF V.flat THEN
        RETURN FALSE
      ELSE
        (* Compute starting texture by recursion on scale: *)
        ComputeStartingGuess(
          (*OUT*) A, (*IN*) S, o, rnd, topEvals, 
          logScale := logScale, 
          xOrg := xOrg, yOrg := yOrg
        );
        (* Write starting texture: *)
        WriteTextures(
          A, NL, txName & "-start", 
          wd := o.UX, ht := o.VY, logScale := logScale,
          xOrg := xOrg, yOrg := yOrg
        );
        (* Optimize synthetic texture: *)
        OptimizeTexture(A, V, topEvals, o.initStep, rnd);
        (* Write annealed texture: *)
        WriteTextures(
          A, NL, txName & "-synth", 
          wd := o.UX, ht := o.VY, logScale := logScale,
          xOrg := xOrg, yOrg := yOrg
        );
        RETURN TRUE
      END;
    END
  END GenTexture;
  
PROCEDURE ComputeStartingGuess(
    READONLY A: ThickTx;   (* OUT: synthetic texture *)
    READONLY S: MaskedTx;  (* IN: Sample and mask *)
    READONLY o: Options;   (* IN: Command line options *)
    rnd: Random.T;         (* I/O: Randomness source *)
    topEvals: NAT;         (* IN: Max energy evaluations in current layer *)
    logScale: NAT;         (* IN: Log2(original_area/current_area) *) 
    xOrg, yOrg: LONG;      (* IN: reference point for texture plotting *)
  ) =
  VAR RUX, RVX, RVY, BUX, BVX, BVY: NAT;
  CONST QuantumNoiseVar = 1.0d-12;
  BEGIN
    WITH
      NL = NUMBER(A),
      AUX = A[0].UX, AVX = A[0].VX, AVY = A[0].VY, 
      SUX = S[0].UX, SVX = S[0].VX, SVY = S[0].VY
    DO
      BUX := AUX; BVX := AVX; BVY := AVY;
      Tx.HalfShrunkPeriods(BUX, BVX, BVY);
      RUX := SUX; RVX := SVX; RVY := SVY;
      Tx.HalfShrunkPeriods(RUX, RVX, RVY);
      
      IF 2*BUX*BVY = AUX*AVY THEN
        (* Can shrink; recurse on coarser scale: *)
        WITH
          R = NEW(REF MaskedTx, NL+1)^,
          B = NEW(REF ThickTx, NL)^
        DO
          (* Shrink sample texture and mask: *)
          R[NL] := Tx.New(RUX, RVX, RVY, ZeroPixel, ZeroPixel);
          Tx.HalfShrink(S[NL], R[NL], combine := Tx.MinCross);
          FOR k := 0 TO NL-1 DO
            R[k] := Tx.New(RUX, RVX, RVY, ZeroPixel, ZeroPixel);
            Tx.HalfShrink(S[k], R[k], combine := Tx.AvgCross);
            MaskOutTexture(R[k], R[NL], UnitPixel);
          END;
          (* Allocate the shrunken synthetic texture: *)
          FOR k := 0 TO NL-1 DO
            B[k] := Tx.New(BUX, BVX, BVY, ZeroPixel, ZeroPixel);
          END;
          (* Try to generate the shrunken synthetic texture: *)
          IF GenTexture(
              (*OUT*) B, (*IN*) R, o, rnd, 
              topEvals := 2*topEvals, logScale := logScale + 1,
              xOrg := 0.5d0 * (yOrg - xOrg),
              yOrg := 0.5d0 * (xOrg + yOrg)
            )
          THEN
            (* Expand it to original size: *)
            FOR k := 0 TO NL-1 DO
              Tx.HalfExpand(B[k], A[k], Tx.Avg4x4);
            END
          ELSE
            (* Copy pixels from "S" into "A"; *)
            SprinklePixels(S, A, rnd);
          END
        END
      ELSE
        (* Copy pixels from "S" into "A"; *)
        SprinklePixels(S, A, rnd);
      END;
      (* Equalize variance: *)
      FOR k := 0 TO NL-1 DO
        VAR varA, avgA: LRPixel;
            varS, avgS: LRPixel;
        BEGIN
          Tx.MeanAndVariance(A[k], avgA, varA);
          MaskedMeanAndVariance(S[k], S[NL], avgS, varS);
          WITH 
            avg = avgS[0] - avgA[0],
            dev = o.random * (Math.sqrt(QuantumNoiseVar + MAX(0.0d0, varS[0] - varA[0])))
          DO
            Message("Adding noise with avg = " & FLR(avg) & " dev = " & FLR(dev));
            Tx.NRandom(A[k], rnd, LRPixel{avg,..}, LRPixel{dev,..}, clear := FALSE)
          END
        END
      END;
    END 
  END ComputeStartingGuess;

PROCEDURE OptimizeTexture(
    READONLY A: ThickTx;  (* I/O: Texture to optimize *)
    V: TxEvaluator;       (* IN: Evaluation (``energy'') function. *)
    maxEvals: NAT;        (* IN: Max iterations *)
    initStep: LONG;       (* IN: Initial step *)
    rnd: Random.T;        (* I/O: Randomness source *)
  ) = 
  CONST HugeEnergy = 1.0d+20;
        TinyImprovement = 1.0d-6;
        MaxBadSteps = 5;
  VAR dt, E, EOld, dtOld, dtStart: LONG;
      nEvals: NAT := 0;
      badSteps: NAT := 0;
  BEGIN
    <* ASSERT maxEvals > 0 *>
    WITH
      NL = NUMBER(A),
      AUX = A[0].UX, AVX = A[0].VX, AVY = A[0].VY, 
      dEdA = NEW(REF ThickTx, NL)^
    DO
      TRY
        FOR k := 0 TO NL-1 DO
          dEdA[k] := Tx.New(AUX, AVX, AVY, ZeroPixel, ZeroPixel);
        END;
        StartOptimization(A, V, initStep, (*OUT*) E, dEdA, dt, (*IO*) nEvals);
        dtStart := dt;

        LOOP
          IF nEvals >= maxEvals THEN RETURN END;
          EOld := E;
          dtOld := dt;
          DownhillStep((*IN*) V, maxEvals, (*IO*) A, E, dEdA, dt, nEvals);
          Message("");
          IF E < -HugeEnergy THEN 
            RETURN
          ELSIF dt = 0.0d0 
          OR SQR(EOld-E)/(SQR(EOld) + SQR(E)) < SQR(TinyImprovement) THEN
            INC(badSteps);
            IF badSteps > MaxBadSteps THEN RETURN END;
            Message("(downhill step failed, retrying in random direction...)");
            PerturbGradient(dEdA, rnd, initStep, dt)
          ELSE
            badSteps := 0;        
            dt := 2.0d0 * dt
          END;
        END;
      FINALLY
        Message(" final energy = " & ELR(E, 6))
      END;
    END;
  END OptimizeTexture;
  
PROCEDURE StartOptimization(
    READONLY A: ThickTx;       (* IN: the texture to optimize *)
    V: TxEvaluator;            (* IN: the texture evaluator *)
    initStep: LONG;            (* IN: Initial step length *)
    VAR E: LONG;               (* OUT: the energy of "A" *)
    READONLY dEdA: ThickTx;    (* OUT: the energy gradient at "A" *)
    VAR dt: LONG;              (* OUT: time step *)
    VAR nEvals: NAT;           (* IN/OUT: counts energy evaluations *)
  ) =
  CONST TinyDT = 1.0d-10;
  VAR gradNorm: LONG;
  BEGIN
    WITH
      NL = V.NL
    DO
      (* Compute initial energy and gradient for "A":*)
      V.eval(A, E);
      V.grad(A, dEdA);
      INC(nEvals);
      (* Estimate initial step: *)
      gradNorm := 0.0d0;
      FOR k := 0 TO NL-1 DO
        gradNorm := MAX(gradNorm, Tx.InfNorm(dEdA[k]))
      END;
      dt := MAX(TinyDT, initStep/gradNorm);
      Message(
        "  starting probe: " &
        "  E = " & FLR(E, 6) & 
        "  dt = " & ELR(dt, 6) & 
        "  grad = " & FLR(gradNorm, 6)
      );
    END
  END StartOptimization;

PROCEDURE DownhillStep(
    V: TxEvaluator;            (* IN: The texture evaluator. *)
    maxEvals: NAT;             (* IN: Max energy evaluations *)
    READONLY A: ThickTx;       (* IN/OUT: the texture to optimize *)
    VAR E: LONG;               (* IN/OUT: the energy of "A" *)
    READONLY dEdA: ThickTx;    (* IN/OUT: the energy gradient at "A" *)
    VAR dt: LONG;              (* IN/OUT: time step *)
    VAR nEvals: NAT;           (* IN/OUT: counts energy evaluations *)
  ) =
  CONST TinyDT = 1.0d-10;
  VAR EOld, dsqr, h: LONG;
  BEGIN
    (* Search from "A" along the gradient "dEdA" for a lower point: *)
    WITH
      NL = V.NL
    DO
      EOld := E;
      h := -dt;
      REPEAT 
        dsqr := 0.0d0;
        FOR k := 0 TO NL-1 DO
          dsqr := dsqr + Tx.Modify(A[k], h, dEdA[k]);
        END;
        (* Invariant: the current "A" is the original "A" minus "dt" times "dEdA". *)
        V.eval(A, E);
        INC(nEvals);
        WITH
          step = Math.sqrt(dsqr)
        DO
          Message(
            "  iteration " & Fmt.Pad(Fmt.Int(nEvals), 6) &
            "  E = " & FLR(E, 6) & 
            "  dt = " & ELR(dt, 6) &
            "  step = " & ELR(step, 6)
          );
        END;
        IF E >= EOld THEN 
          IF nEvals >= maxEvals OR dt <= TinyDT THEN
            (* Go back to where we started: *)
            FOR k := 0 TO NL-1 DO
              EVAL Tx.Modify(A[k], dt, dEdA[k])
            END;
            dt := 0.0d0;
            E := EOld;
            RETURN
          END;
          h := 0.61803398875d0 * dt;
          dt := dt - h;
        END;
      UNTIL E < EOld;
      (* Downhill step succeeded, compute new gradient: *)
      V.grad(A, dEdA);
    END
  END DownhillStep;
  
PROCEDURE PerturbGradient(
    READONLY dEdA: ThickTx; 
    rnd: Random.T;
    initStep: LONG;            (* IN: Initial step length *)
    VAR dt: LONG;              (* OUT: time step *)
  ) =
  (*
    Randomizes gradient values, and adjusts time step so that 
    it jumps by "initStep". *)
  CONST TinyDT = 1.0d-10;
  VAR sumSqr: LONG := 1.0d-20;
      s: LRPixel;
  BEGIN
    WITH  
      NL = NUMBER(dEdA),
      UX = dEdA[0].UX, VY = dEdA[0].VY
    DO
      sumSqr := 0.0d0;
      FOR k := 0 TO NL-1 DO
        Tx.SumSquares(dEdA[k], s);
        sumSqr := sumSqr + s[0]
      END;
      WITH
        gradNorm = Math.sqrt(sumSqr),
        amp = gradNorm/FLOAT(NL*UX*VY, LONG)
      DO
        dt := MAX(TinyDT, initStep/gradNorm);
        FOR k := 0 TO NL-1 DO
          WITH t = dEdA[k] DO
            FOR y := 0 TO VY-1 DO
              WITH ty = t.st[y] DO
                FOR x := 0 TO UX-1 DO 
                  WITH c = FLOAT(ty[x], LONG) DO
                    ty[x] := FLOAT(c + amp*GaussRandom.LongReal(rnd))
                  END
                END
              END
            END
          END
        END
      END
    END
  END PerturbGradient;

PROCEDURE MaskedMeanAndVariance(
    READONLY t: Tx.T;
    READONLY m: Tx.T; 
    VAR a, v: LRPixel;
  ) =
  (*
    Like Tx.MeanAndVariance, but considers only pixels of "t" where "m" is 1.
  *)
  VAR c: LONG;
      N: NAT := 0;
  BEGIN
    <* ASSERT NUMBER(a) = 1 *>
    WITH
      UX = t.UX, VY = t.VY
    DO
      <* ASSERT m.UX = UX AND m.VY = VY *>
      (* Compute means: *)
      a := LRPixel{0.0d0, ..};
      FOR y := 0 TO VY-1 DO
        WITH ty = t.st[y], my = m.st[y] DO
          FOR x := 0 TO UX-1 DO 
            IF my[x] > 0.5 THEN 
              INC(N);
              c := FLOAT(ty[x], LONG);
              WITH ak = a[0] DO ak := ak + c END
            END
          END
        END
      END;
      WITH fN = FLOAT(N, LONGREAL) DO
        IF N < 1 THEN a[0] := 0.0d0 ELSE a[0] := a[0]/fN END
      END;
      
      (* Compute variances: *)
      v := LRPixel{0.0d0, ..};
      FOR y := 0 TO VY-1 DO
        WITH ty = t.st[y], my = m.st[y] DO
          FOR x := 0 TO UX-1 DO 
            IF my[x] > 0.5 THEN 
              c := FLOAT(ty[x], LONG);
              WITH vk = v[0], dk = c - a[0] DO vk := vk + dk*dk END
            END
          END
        END
      END;
      WITH fN = FLOAT(N - 1, LONGREAL) DO
        IF N < 2 THEN v[0] := 0.0d0 ELSE v[0] := v[0]/fN END
      END;
    END
  END MaskedMeanAndVariance;

PROCEDURE MaskOutTexture(
    READONLY t: Tx.T;       (* A texture *)
    READONLY m: Tx.T;       (* A boolean mask (0=undefined, 1=defined) *)
    READONLY fill: LRPixel; (* The fill value *)
  ) =
  (*
    Sets "t[y,x]" to "fill" whenever "m[y,x]" is "0".
  *)
  BEGIN
    <* ASSERT NUMBER(fill) = 1 *>
    WITH
      UX = t.UX, VY = t.VY,
      f = FLOAT(fill[0])
    DO
      <* ASSERT m.UX = UX AND m.VY = VY *>
      FOR y := 0 TO VY-1 DO
        WITH ty = t.st[y], my = m.st[y] DO
          FOR x := 0 TO UX-1 DO 
            IF my[x] < 0.5 THEN 
              ty[x] := f
            END
          END
        END
      END;
    END
  END MaskOutTexture;

PROCEDURE SprinklePixels(
    READONLY S: MaskedTx;
    READONLY A: ThickTx;
    rnd: Random.T
  ) =
  (*
    Fills "A" with randomly selected pixels from the sample texture
    "S[0..NL-1]" where "S[NL]" is 1.
  *)
  BEGIN
    WITH
      NL = NUMBER(A),
      AUX = A[0].UX, AVY = A[0].VY,
      SUX = S[0].UX, SVY = S[0].VY,
      DX = (SUX DIV 2) - ((SUX DIV 2) MOD 2) + 1,
      DY = (SVY DIV 2) - ((SVY DIV 2) MOD 2) + 1,
      M = S[NL]
    DO
      FOR yA := 0 TO AVY-1 DO
        FOR xA := 0 TO AUX-1 DO
          VAR xS, yS: NAT;
          BEGIN 
            xS := rnd.integer(0, SUX-1);
            yS := rnd.integer(0, SVY-1);
            REPEAT
              xS := (xS + DX) MOD SUX;
              IF xS < DX THEN yS := (yS + DY) MOD SVY END;
            UNTIL M.st[yS,xS] > 0.5;
            FOR k := 0 TO NL-1 DO 
              A[k].st[yA,xA] := S[k].st[yS,xS]
            END;
          END
        END
      END;
    END 
  END SprinklePixels;

PROCEDURE ComputeEnergyFromDistrs(
    READONLY AF: FeatureDistr;
    READONLY SF: FeatureDistr;
    cmp: Comparator;
    VAR E: LONGREAL;
    VAR dEdAF: FeatureDistr;
  ) =
  (*
    Computes the energy "E" and its derivative "dEdAF" with respect to
    the feature distribution "AF". *)
  BEGIN
    WITH
      NA = NUMBER(AF),
      NS = NUMBER(SF),
      NF = NUMBER(SF[0]),
      qA = 1.0d0/FLOAT(NA, LONG),
      qS = -1.0d0/FLOAT(NS, LONG),
      junk = NEW(REF Features, NF)^
    DO
      E := 0.0d0;
      FOR i := 0 TO NA-1 DO
        FOR k := 0 TO NF-1 DO dEdAF[i,k] := 0.0d0 END
      END;
      (* Compute attraction energy: *)
      FOR iS := 0 TO NS-1 DO
        FOR k := 0 TO NF-1 DO junk[k] := 0.0d0 END;
        FOR iA := 0 TO NA-1 DO
          cmp.eval(AF[iA], qA, dEdAF[iA], SF[iS], qS, junk, E);
        END;
      END;
      (* Compute repulsion energy: *)
      FOR iA := 0 TO NA-1 DO
        FOR jA := iA+1 TO NA-1 DO
          cmp.eval(AF[iA], qA, dEdAF[iA], AF[jA], qA, dEdAF[jA], E);
        END;
      END
    END
  END ComputeEnergyFromDistrs;

PROCEDURE ComputeSampleFeatureDistr(
    READONLY S: MaskedTx;
    map: Mapper.T;
  ): REF FeatureDistr =
  (*
    Computes the multiset of local features from the 
    sample multilayer texture "S[0..NL-1]", considering only pixels
    where "S[NL]" is 1.
  *)

  BEGIN
    WITH
      NL = NUMBER(S) - 1,
      SUX = S[0].UX,
      SVY = S[0].VY,
      tF = NEW(REF FeatureDistr, SUX*SVY, map.NF)^
    DO

      VAR kF: NAT := 0;

      PROCEDURE CallMapper(<*UNUSED*> x, y:NAT; READONLY n: Nhoods) =
        BEGIN
          IF map.defined(n[NL]) THEN 
            map.eval(SUBARRAY(n, 0, NL), tF[kF]);
            INC(kF)
          END
        END CallMapper;

      BEGIN
        Tx.MultEnumerateNeighborhoods(S, CallMapper, periodic := FALSE);
        WITH
          SF = NEW(REF FeatureDistr, kF, map.NF)
        DO
          SF^ := SUBARRAY(tF, 0, kF);
          RETURN SF
        END
      END
    END
  END ComputeSampleFeatureDistr;
  
PROCEDURE WriteFeatureDistr(
    READONLY F: FeatureDistr;
    name: TEXT;
  ) =
  (*
    Prints the feature samples in "F", in a format suitable
    for gnuplot, etc.  The file name is "name" with ".fdt" appended.
    *)
  <* FATAL Thread.Alerted, Wr.Failure, OSError.E *>

  BEGIN
    IF NUMBER(F) = 0 THEN RETURN END;
    WITH
      N = NUMBER(F),
      M = NUMBER(F[0]),
      wr = FileWr.Open(name & ".ftd")
    DO
      FOR i := 0 TO N-1 DO
        FOR k := 0 TO M-1 DO
          Wr.PutText(wr, FLR(F[i,k], 8));
        END;
        Wr.PutText(wr, "\n")
      END;
      Wr.Close(wr)
    END
  END WriteFeatureDistr;

<*UNUSED*>
PROCEDURE FoldPixelRange(READONLY A: ThickTx): BOOL =
  (*
    Folds pixel values that are outside [-1 _ +1] into that interval.
    returns TRUE if any pixel changed. *)
  VAR changed: BOOL := FALSE;
  BEGIN
    FOR k := 0 TO LAST(A) DO
      WITH 
        t = A[k],
        UX = t.UX,
        VY = t.VY
      DO
        FOR y := 0 TO VY-1 DO
          WITH ty = t.st[y] DO
            FOR x := 0 TO UX-1 DO 
              WITH c = ty[x] DO
                IF c > +1.1 OR c < -1.1 THEN 
                  changed := TRUE;
                  c := 0.9/c
                END;
              END
            END
          END
        END
      END
    END;
    RETURN changed
  END FoldPixelRange;

PROCEDURE ComputeTextureFeatureDistr(
    READONLY A: ThickTx;
    map: Mapper.T;
    VAR AF: FeatureDistr;
  ) =
  (*
    Computes the multiset of local features from the multilayer texture "t". *)

  BEGIN
    WITH
      UX = A[0].UX,
      VY = A[0].VY,
      NN = UX*VY
    DO
      <* ASSERT NUMBER(AF) = NN *>

      VAR kF: NAT := 0;

      PROCEDURE CallMapper(<*UNUSED*> x, y:NAT; READONLY n: Nhoods) =
        BEGIN
          map.eval(n, AF[kF]);
          INC(kF)
        END CallMapper;

      BEGIN
        Tx.MultEnumerateNeighborhoods(A, CallMapper, periodic := TRUE);
        <* ASSERT kF = NN *>
      END
    END
  END ComputeTextureFeatureDistr;

PROCEDURE FeatureDerivsToTextureDerivs(
    READONLY dEdAF: FeatureDistr;
    READONLY AF: FeatureDistr;
    READONLY A: ThickTx;
    map: Mapper.T;
    READONLY dEdA: ThickTx;
    VAR dEdn: Nhoods;  (* Work area *)
  ) =
  (*
    Computes the energy derivatives in texture space "dEdA" from
    the derivatives in feature space "dEdAF".
  *)
  BEGIN
    WITH
      NL = NUMBER(A),
      UX = A[0].UX,
      VY = A[0].VY,
      NN = UX*VY
    DO
      <* ASSERT NUMBER(AF) = UX*VY *>

      VAR kF: NAT := 0;

      PROCEDURE CallMapper(x, y: NAT; READONLY n: Nhoods) =
        BEGIN
          map.diff(dEdAF[kF], AF[kF], n, dEdn);
          FOR k := 0 TO NL-1 DO
            Tx.SplatNeighborhood(dEdA[k], x, y, dEdn[k])
          END;
          INC(kF)
        END CallMapper;

      BEGIN
        FOR k := 0 TO NL-1 DO Tx.Zero(dEdA[k]) END;
        Tx.MultEnumerateNeighborhoods(A, CallMapper, periodic := TRUE);
        <* ASSERT kF = NN *>
        FOR k := 0 TO NL-1 DO Tx.FoldBorder(dEdA[k]) END;
      END
    END
  END FeatureDerivsToTextureDerivs;

PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>
  VAR o: Options;
  BEGIN
    WITH pp = NEW(ParseParams.T).init(stderr) DO
      TRY

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

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

        pp.getKeyword("-size");
        o.UX := pp.getNextInt(1, 4096);
        o.VY := pp.getNextInt(1, 4096);

        IF pp.keywordPresent("-offset") THEN
          o.VX := pp.getNextInt(0, o.UX-1)
        ELSE
          o.VX := 0
        END;

        IF pp.keywordPresent("-maxEvals") THEN
          o.maxEvals := pp.getNextInt(1, 100000)
        ELSE
          o.maxEvals := 50
        END;

        IF pp.keywordPresent("-initStep") THEN
          o.initStep := pp.getNextLongReal(0.00001d0, 10000.0d0)
        ELSE
          o.initStep := 0.01d0
        END;
        
        IF pp.keywordPresent("-random") THEN
          o.random := pp.getNextLongReal(0.0d0, 10.0d0)
        ELSE
          o.random := 0.0d0
        END;
        
        o.fold := pp.keywordPresent("-fold");

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

        pp.getKeyword("-map");
        WITH mapName = pp.getNext() DO
          IF Text.Equal(mapName, "Phi") THEN
            o.map := NEW(PhiFeatures2D.T).init(NL := 1)
          ELSIF Text.Equal(mapName, "Csi") THEN
            o.map := NEW(CsiFeatures2D.T).init(NL := 1)
          ELSIF Text.Equal(mapName, "FunSlope") THEN
            o.map := NEW(FunSlopeFeatures2D.T).init(NL := 1)
          ELSIF Text.Equal(mapName, "FunGrad") THEN
            o.map := NEW(FunGradFeatures2D.T).init(NL := 1)
          ELSIF Text.Equal(mapName, "Cross") THEN
            o.map := NEW(CrossFeatures2D.T).init(NL := 1)
          ELSE
            pp.error("invalid \"-map\" option = " & mapName)
          END
        END;

        pp.getKeyword("-cmp");
        WITH cmpName = pp.getNext() DO
          IF Text.Equal(cmpName, "Coulomb") THEN
            o.cmp := NEW(CoulombComparator,
              NF := 5,
              fuzz := pp.getNextLongReal(0.00001d0, 10000.0d0)
            )
         ELSE
            pp.error("invalid \"-cmp\" option = " & cmpName)
          END
        END;

        pp.finish();
      EXCEPT
      | ParseParams.Error =>
          Wr.PutText(stderr, "Usage: MREvolver2D\\\n");
          Wr.PutText(stderr, "  -inName <name> \\\n");
          Wr.PutText(stderr, "  -outName <name> \\\n");
          Wr.PutText(stderr, "  -size <width> <height> [ -offset <offset> ] \\\n");
          Wr.PutText(stderr, "  -map { Phi | Cross |... } \\\n");
          Wr.PutText(stderr, "  -cmp { Coulomb <fuzz> | ... } \\\n");
          Wr.PutText(stderr, "  [ -maxEvals <num> ] \\\n");
          Wr.PutText(stderr, "  [ -initStep <timestep> ] \\\n");
          Wr.PutText(stderr, "  [ -random <n.nn> ] \\\n");
          Wr.PutText(stderr, "  [ -fold ] [ -norm ] \n");
          Process.Exit (1);
      END;
    END;
    RETURN o
  END GetOptions;
  
PROCEDURE Pow(b: INT; e: NAT): INT =
  VAR r: INT;
  BEGIN
    r := 1;
    WHILE e > 0 DO r := r * b; DEC(e) END;
    RETURN r
  END Pow;

PROCEDURE SQR(x: LONG): LONG =
  BEGIN
    RETURN x*x
  END SQR;

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

PROCEDURE ELR(x: LONG; prec: NAT := 6): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(x, prec := prec-1, style := Fmt.Style.Sci), prec + 7)
  END ELR;

<*UNUSED*>
PROCEDURE FR(x: REAL; prec: NAT := 6): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.Real(x, prec := prec, style := Fmt.Style.Fix), prec + 6)
  END FR;

<*UNUSED*>
PROCEDURE ER(x: REAL; prec: NAT := 6): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.Real(x, prec := prec-1, style := Fmt.Style.Sci), prec + 6)
  END ER;

PROCEDURE FI(i: INT; width: NAT; fill: CHAR := ' '): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.Int(i), width, fill)
  END FI;

PROCEDURE FB(x: BOOL): TEXT =
  BEGIN
    IF x THEN RETURN "TRUE" ELSE RETURN "FALSE" END
  END FB;

PROCEDURE Message(t: TEXT) =
  <* FATAL Thread.Alerted, Wr.Failure *>
  BEGIN
    Wr.PutText(stderr, t);
    Wr.PutChar(stderr, '\n');
    Wr.Flush(stderr);
  END Message;

<*UNUSED*>
PROCEDURE NewLongs(READONLY a: LONGS): REF LONGS =
  (*
    Modula-3 is so silly....
  *)
  BEGIN
    WITH n = NEW(REF LONGS, NUMBER(a)) DO n^ := a; RETURN n END
  END NewLongs;

(* DISTRIBUTION COMPARATORS *)

TYPE
  CoulombComparator = Comparator OBJECT
      fuzz: LONG;
    OVERRIDES
      print := CoulombPrint;
      eval := CoulombEval;
    END;

PROCEDURE CoulombPrint(s: CoulombComparator; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Coulomb: electric potential model\n");
    Wr.PutText(wr, "  fuzz radius = " & FLR(s.fuzz) & "\n");
    Wr.Flush(wr)
  END CoulombPrint;

PROCEDURE CoulombEval(
    s: CoulombComparator;
    READONLY A: Features;
    qA: LONG;
    VAR dEdA: Features;
    READONLY B: Features;
    qB: LONG;
    VAR dEdB: Features;
    VAR E: LONGREAL;
  ) =
  (*
    Adds to "E" the energy due to the charge pair "(A,qA)" and "(b,qB)";
    and to "dEdA", "dEdB" the correponding derivatives. *)
  VAR h: LONG := s.fuzz * s.fuzz;
  BEGIN
    FOR k := 0 TO LAST(A) DO
      WITH rABk = A[k] - B[k] DO
        h := h + rABk*rABk; (* dhdAk := 2*rABk; dhdBk := -2*rABk *)
      END;
    END;
    WITH
      r = Math.sqrt(h), (* drdh := 1/(2r) *)
      e = qA*qB/r,      (* dedr := -qA*qB/r^2 = -e/r *)
      c = -e/h
    DO
      E := E + e;
      FOR k := 0 TO LAST(A) DO
        WITH dedAk = c * (A[k] - B[k]) DO
          dEdA[k] := dEdA[k] + dedAk;  (* dedr*drdh*dhdAk = (-e/h)*(rABk)) *)
          dEdB[k] := dEdB[k] - dedAk
        END;
      END;
    END
  END CoulombEval;

(* TEXTURE EVALUATORS *)

TYPE 
  StatEvaluator = TxEvaluator OBJECT
      fold: BOOL;     (* TRUE to constrain pixels to [-1 _ +1] *)
      norm: BOOL;     (* TRUE to normalize the gradient locally *)

      (* All fields below this point are set by the "init" method. *)
      map: Mapper.T;           (* Maps neighborhoods to local features *)
      cmp: Comparator;         (* Compares two feature distributions *)
      SF: REF FeatureDistr;    (* The reference distribution, from sample. *)
      
       
      (* These fields are used internally by the "eval" and "grad" methods: *)
      gScale: Tx.T;             (* Scale factors for local gradient normaliz. *)
      dEdn: REF Nhoods;         (* Derivatives of "E" term rel nhood. *)

      (* These fields are set by the "eval" method and used by "grad": *)
      AF: REF FeatureDistr;     (* Feature distribution of "A" *)
      dEdAF: REF FeatureDistr;  (* Derivatives of "E" relative to "AF" *)

    METHODS

      init(
        READONLY S: MaskedTx;   (* Reference texture *)
        AUX, AVX, AVY: NAT;     (* Dimensions of target texture *)
        map: Mapper.T;          (* Feature mapper to use *)
        cmp: Comparator;        (* Feature distribution comparator to use *)
        fold: BOOL;             (* TRUE to constrain pixels to [-1 _ +1] *)
        norm: BOOL;             (* TRUE to normalize the gradient locally *)
        fName: TEXT := "";      (* File name for feature file, or "" if none *)
      ): StatEvaluator := StatInit;
      
    OVERRIDES
      print := StatPrint;
      eval := StatEval;
      grad := StatGrad;
    END;

PROCEDURE StatInit(
    s: StatEvaluator;
    READONLY S: MaskedTx;
    AUX, AVX, AVY: NAT;
    map: Mapper.T;
    cmp: Comparator;
    fold: BOOL;
    norm: BOOL;
    fName: TEXT := "";
  ): StatEvaluator =
  BEGIN
    s.SF := ComputeSampleFeatureDistr(S, map);
    IF NOT Text.Equal(fName, "") THEN WriteFeatureDistr(s.SF^, fName) END;
    s.AF := NEW(REF FeatureDistr, AUX*AVY, map.NF);
    WITH
      NS = NUMBER(s.SF^),
      NA = NUMBER(s.AF^),
      NL = map.NL
    DO
      s.NL := NL;
      s.flat := (NA = 0) OR (NS = 0);
      s.dEdAF := NEW(REF FeatureDistr, NA, map.NF);
      s.gScale := Tx.New(AUX, AVX, AVY, ZeroPixel, ZeroPixel);
      s.dEdn := NEW(REF Nhoods, NL);
      s.map := map;
      s.cmp := cmp;
      s.fold := fold;
      s.norm := norm;
      RETURN s
    END
  END StatInit;
  
PROCEDURE StatPrint(s: StatEvaluator; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    <* ASSERT s.SF # NIL *>
    WITH
      NSF = NUMBER(s.SF^)
    DO
      Wr.PutText(wr, "StatEvaluator: evaluates a texture by comparing its\n");
      Wr.PutText(wr, "  feature distribution with that of a sample texture.\n");
      Wr.PutText(wr, "--- feature mapper ----------------------------------\n");
      s.map.print(wr);
      Wr.PutText(wr, "--- distribution comparator -------------------------\n");
      s.cmp.print(wr);
      Wr.PutText(wr, "--- other parameters --------------------------------\n");
      Wr.PutText(wr, "  size of sample distribution " & Fmt.Int(NSF) & "\n");
      Wr.PutText(wr, "  constrain pixel values to [-1 _ +1] = " & FB(s.fold) & "\n");
      Wr.PutText(wr, "  normalize the gradient locally = " & FB(s.norm) & "\n");
      Wr.Flush(wr)
    END
  END StatPrint;

PROCEDURE StatEval(
    s: StatEvaluator;
    READONLY A: ThickTx; 
    VAR E: LONG;
    fName: TEXT := "";
  ) = 
  BEGIN
    ComputeTextureFeatureDistr(A, s.map, s.AF^);
    IF NOT Text.Equal(fName, "") THEN WriteFeatureDistr(s.AF^, fName) END;
    ComputeEnergyFromDistrs(s.AF^, s.SF^, s.cmp, E, s.dEdAF^)
  END StatEval;
  
PROCEDURE StatGrad(
    s: StatEvaluator; 
    READONLY A: ThickTx; 
    READONLY dEdA: ThickTx;
  ) = 
  BEGIN
    FeatureDerivsToTextureDerivs(s.dEdAF^, s.AF^, A, s.map, dEdA, s.dEdn^);
    IF s.fold THEN AdjustGradient(A, dEdA) END;
    IF s.norm THEN NormalizeGradientLocally(dEdA, s.gScale) END;
  END StatGrad;
  
PROCEDURE AdjustGradient(
    READONLY A: ThickTx;
    READONLY dEdA: ThickTx;
  ) =
  (*
    Adjusts gradient so that pixels on or outside the range limits
    have gradients pointing outwards. *)
  BEGIN
    FOR k := 0 TO LAST(A) DO
      WITH 
        t = A[k], dt = dEdA[k],
        UX = t.UX,
        VY = t.VY
      DO
        FOR y := 0 TO VY-1 DO
          WITH ty = t.st[y], dty = dt.st[y] DO
            FOR x := 0 TO UX-1 DO 
              WITH c = ty[x], dc = dty[x] DO
                IF c > +1.0 THEN 
                  dc := c - 1.0
                ELSIF c < -1.0 THEN
                  dc := c + 1.0
                END
              END
            END
          END
        END
      END
    END;
  END AdjustGradient;
  
PROCEDURE NormalizeGradientLocally(
    READONLY dEdA: ThickTx;  (* I/O: Gradient to normalize *)
    READONLY gScale: Tx.T;   (* WORK: scale factors *)
  ) = 
  
  CONST R0 = 0.25d0;
        R1 = 0.125d0;
        R2 = 0.0625d0;
        
  VAR fNL: LONG := FLOAT(NUMBER(dEdA), LONG);
        
  PROCEDURE ComputeLocalScale (* : Tx.MultNeighborhoodVisitProc *) (
      x, y: CARDINAL; 
      READONLY n: Nhoods;
    ) =
    VAR sumSqr: LONG := 0.0d0;
        xx: INT := x; 
        yy: INT := y;
    BEGIN
      FOR k := 0 TO LAST(n) DO 
        WITH 
          nk = n[k],
          nmm = FLOAT(nk[-1,-1], LONG),
          nmo = FLOAT(nk[-1,00], LONG),
          nmp = FLOAT(nk[-1,+1], LONG),

          nom = FLOAT(nk[00,-1], LONG),
          noo = FLOAT(nk[00,00], LONG),
          nop = FLOAT(nk[00,+1], LONG),

          npm = FLOAT(nk[+1,-1], LONG),
          npo = FLOAT(nk[+1,00], LONG),
          npp = FLOAT(nk[+1,+1], LONG)
        DO 
          sumSqr := sumSqr 
            + R0 * ( noo * noo )
            + R1 * ( nom*nom + nop*nop + nmo*nmo + npo*npo )
            + R2 * ( nmm*nmm + nmp*nmp + npm*npm + npp*npp )
        END
      END;
      Tx.Reduce(gScale, xx, yy);
      gScale.st[yy,xx] := FLOAT(fNL/Math.sqrt(sumSqr))
    END ComputeLocalScale;
  
  BEGIN
    Tx.MultEnumerateNeighborhoods(dEdA, ComputeLocalScale, TRUE);
    FOR k := 0 TO LAST(dEdA) DO
      WITH 
        d = dEdA[k],
        UX = d.UX,
        VY = d.VY
      DO
        FOR y := 0 TO VY-1 DO
          WITH dy = d.st[y], sy = gScale.st[y] DO
            FOR x := 0 TO UX-1 DO 
              WITH d = dy[x], s = sy[x] DO d := d * s END
            END
          END
        END
      END
    END;
  END NormalizeGradientLocally;
  
BEGIN
  Main();
END MREvolver2D.