MODULE BITSignalsFromCands EXPORTS Main;
(* Last edited on 1999-12-31 09:56:46 by hcgl *)

(* Converts candidates into pairs of signals. *)

(*
  Reads a list of candidates. For each candidates, 
  converts each curve segment into a signal (the shape function),
  and writes it to disk.  Also writes the mean and difference 
  signals.
*)

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

FROM Math IMPORT sqrt;
FROM PZShape IMPORT ComputeSpectrum, ComputePowerSpectrum;
FROM PZTypes IMPORT LONG, LONGS, INT, NAT;

FROM Stdio IMPORT stderr;

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

TYPE
  Options = RECORD
      (* Data identification parameters: *)
      input: TEXT;              (* Input candidate file name. *)
      output: TEXT;             (* Output directory. *)
      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" *)
    END;

TYPE
  ChainData = PZLR3Chain.ReadData;

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

 PROCEDURE Main() =
  BEGIN
    WITH
      o = GetOptions()
    DO
      TRY FS.CreateDirectory(o.output) EXCEPT OSError.E => (*OK*) END;
      WITH    
        candData = PZCandidate.Read(FileRd.Open(o.input & ".can")),
        cand = candData.c^,
        inCmt = candData.cmt,
        needed = PZCandidate.GetCurves(cand)^,
        chData = PZLR3Chain.ReadAll(
          o.chainPrefix, o.band, o.extension, sel := needed, 
          dir := o.chainDir, headerOnly := FALSE
        ),
        ch = chData.chain^
      DO
        FOR iCand := 0 TO LAST(cand) DO
          ProcessOneCandidate(iCand, cand[iCand], ch, o, inCmt)
        END
      END
    END
  END Main;

PROCEDURE ProcessOneCandidate(
    iCand: NAT;
    READONLY cand: PZCandidate.T;
    READONLY ch: ARRAY OF ChainData;
    READONLY o: Options;
    inCmt: TEXT;
  ) =
  (*
    Let "a,b" be the two curves that constitute candidate "cand".
    This procedure extracts the chains and shape functions of the
    two segments of "cand" and writes them to disk.  It also writes
    the mean and difference signals, and the Fourier transforms
    of all the signals.
    
    The files are written to directory "DDD/NNNNNN" where "DDD" is
    "o.output" and "NNNNNN" is the candidate's number.  
    The file names are "a.flc" "b.flc", "d.fsh", etc.
  *)
  BEGIN
    WITH
      aSeg = cand.seg[0],
      cha = ch[aSeg.cvx],

      bSeg = cand.seg[1],
      chb = ch[bSeg.cvx],

      a = PZLR3Chain.ExtractSegment(cha.c^, aSeg)^,
      b = PZLR3Chain.ExtractSegment(chb.c^, bSeg)^,

      dirName = o.output & "/" & Fmt.Pad(Fmt.Int(iCand), 6, '0'),
      candComment = inCmt & "\n" & CandComment(iCand, cand, ch),

      sa = ExtractSegmentShape(a, dirName, "a", candComment)^,
      sb = ExtractSegmentShape(b, dirName, "b", candComment)^,
      
      sm = ComputeMeanSignal(sa, sb)^,
      sd = ComputeDiffSignal(sa, sb)^
    DO 
      ProcessSignal(sa, dirName, "a", candComment);
      ProcessSignal(sb, dirName, "b", candComment);
      ProcessSignal(sd, dirName, "d", candComment);
      ProcessSignal(sm, dirName, "m", candComment);
      Wr.PutText(stderr, "\n");
    END
  END ProcessOneCandidate;

PROCEDURE ExtractSegmentShape(
    READONLY c: PZLR3Chain.T;  (* Mother chain *)
    dirName: TEXT;             (* Output file name prefix *)
    curveName: TEXT;           (* Curve name ("a" or "b"). *)
    candCmt: 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(dirName) EXCEPT OSError.E => (*OK*) END;
    WITH
      rs = PZShape.Signal(c),
      curveCmt = candCmt & "\n" & "side = \"" & curveName & "\""
    DO
      WITH fileName = dirName & "/" & curveName & ".flc", wr = FileWr.Open(fileName) DO
        Wr.PutText(stderr, "# writing " & fileName & "\n");
        PZLR3Chain.Write(wr, curveCmt, c, unit := unitForCoord);
        Wr.Close(wr)
      END;
      RETURN rs
    END
  END ExtractSegmentShape;
  
PROCEDURE ProcessSignal(
    READONLY s: PZLRChain.T;
    dirName: TEXT;             (* Output file name prefix *)
    curveName: TEXT;           (* Curve name ("a" or "b"). *)
    candCmt: TEXT;             (* Candidate's comments *)
  ) =
  (*
    Writes to disk the signal "s", and also its Fourier transform
    and power spectra.
  *)
  BEGIN
    TRY FS.CreateDirectory(dirName) EXCEPT OSError.E => (*OK*) END;
    WITH 
      fft = ComputeSpectrum(s)^,
      N = NUMBER(fft),
      unitFourier = unitForShape / sqrt(FLOAT(N, LONG)) / 10.0d0,
      unitPower = 1000.0d0 * unitFourier*unitFourier,
      pwr = ComputePowerSpectrum(fft)^, 
      comment = candCmt & "\n" & SignalComment(curveName)
    DO
      WITH fileName = dirName & "/" & curveName & ".fsh", wr = FileWr.Open(fileName) DO
        Wr.PutText(stderr, "# writing " & fileName & "\n");
        PZLRChain.Write(wr, comment & "\nshape function", s, 
          stride := 0.0d0, unit := unitForShape/2.0d0
        );
        Wr.Close(wr)
      END;
      WITH fileName = dirName & "/" & 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 := unitFourier
        );
        Wr.Close(wr)
      END;
      WITH fileName = dirName & "/" & 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 := unitPower
        );
        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 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 SignalComment(curveName: TEXT): TEXT = 
  BEGIN
    IF Text.Equal(curveName, "a") THEN
      RETURN "first fragment's shape function"
    ELSIF Text.Equal(curveName, "b") THEN
      RETURN "second fragment's shape function"
    ELSIF 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 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();

        pp.finish();
      EXCEPT
      | ParseParams.Error =>
          Wr.PutText(stderr, "Usage: BITExtractSubCands \\\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");
         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 BITSignalsFromCands.