program MergeHistograms2;

{$mode objfpc}{$H+}

uses
  {$IFDEF UNIX}{$IFDEF UseCThreads}
  cthreads,
  {$ENDIF}{$ENDIF}
  Classes, SysUtils, CustApp, Math
  { you can add units after this };

type

  { TMergeHistograms }

  TMergeHistograms = class(TCustomApplication)
  protected
    procedure DoRun; override;
  public
    constructor Create(TheOwner: TComponent); override;
    destructor Destroy; override;
    procedure WriteHelp; virtual;
  end;

{ TMergeHistograms }

const
  MaxRStr       = 'Max R = ';
  PStr          = 'P = ';
  HistStartStr  = '--- Cluster histogram ---';
  HistStartStr2 = 'Size';
  Delim         = #9;

type
  THistEntry = record
    S: LongInt;
    Num: Extended;//LongInt;
    Prob: Extended;
  end;
  THist = array of THistEntry;
  TPercOut = record
    FN: AnsiString;
    MaxR: LongInt;
    P: Extended;
    Hist: THist;
    HistNum: LongInt;
    Skip: Boolean;
    // ...
  end;

const
  EmptyHistEntry: THistEntry = (S: 0; Num: 0; Prob: 0;);

var
  PercOut: array of TPercOut;
  PercOutNum: LongInt;
  UniqueS: array of LongInt;
  UniqueSNum, UniqueSNumGuess: LongInt;
  BoolAverage: Boolean;

procedure ReadInput;
var
  FN: AnsiString;
  procedure Init;
  begin
    SetLength(PercOut, 10); // Initialize
    PercOutNum := 0;
    UniqueSNumGuess := 0;
  end;
  procedure ReadPercOutput(const AFN: AnsiString); // Works with Percolation2D_v4
  var
    FH: TextFile;
    Str: AnsiString;
    I: LongInt;
    Mode: Byte;
  begin
    if not FileExists(FN)
      then
        begin
          WriteLn('[WARNING] File "', FN, '" does not exist!');
          Exit;
        end;

    PercOutNum := PercOutNum+1;
    if PercOutNum>Length(PercOut)
      then
        SetLength(PercOut, 2*PercOutNum+1);

    with PercOut[PercOutNum-1] do
      begin
        FN := AFN;
        MaxR := -1;
        P := -1;
        SetLength(Hist, 10); // Initialize
        HistNum := 0;
        Skip := False;

        try
          Assign(FH, AFN);
          Reset(FH);

          Mode := 0;
          while not EOF(FH) do
            begin
              ReadLn(FH, Str);

              case Mode of
                0:
                  begin
                    I := Pos(MaxRStr, Str);
                    if I>0
                      then
                        begin
                          Str := Copy(Str, I+Length(MaxRStr));
                          MaxR := StrToInt(Str);
                          Continue;
                        end;

                    I := Pos(PStr, Str);
                    if I>0
                      then
                        begin
                          Str := Copy(Str, I+Length(PStr));
                          P := StrToFloat(Str);
                          Continue;
                        end;

                    I := Pos(HistStartStr, Str);
                    if I>0
                      then
                        begin
                          Mode := 1;
                          Continue;
                        end;
                  end;
                1:
                  if Pos(HistStartStr2, Str)>0
                    then
                      Mode := 2;
                2:
                  begin
                    if Str=''
                      then
                        Break;

                    HistNum := HistNum+1;
                    if HistNum>Length(Hist)
                      then
                        SetLength(Hist, 2*HistNum+1);

                    with Hist[HistNum-1] do
                      begin
                        I := Pos(Delim, Str);
                        S := StrToInt(Copy(Str, 1, I-1));

                        Str := Copy(Str, I+1);
                        I := Pos(Delim, Str);
                        Num := StrToInt(Copy(Str, 1, I-1));

                        Str := Copy(Str, I+1);
                        Prob := StrToFloat(Str);
                      end;
                  end;
              end;
            end;
        finally
          Close(FH);
        end;

        if Length(Hist)<>HistNum
          then
            SetLength(Hist, HistNum);

        UniqueSNumGuess := Max(UniqueSNumGuess, HistNum);
      end;
  end;
  procedure Finish;
  begin
    if Length(PercOut)<>PercOutNum
      then
        SetLength(PercOut, PercOutNum);
  end;
begin
  Init;

  while True do
    begin
      ReadLn(FN);
//WriteLn(FN);
      if (FN='')
        then
          Break;

      ReadPercOutput(FN);
//WriteLn(PercOut[PercOutNum-1].HistNum);
    end;

  Finish;
end;

procedure FindSolution; // Merge histograms
  procedure FindUniqueS;
  var
    I: LongInt;
    J: array of LongInt;
    MinS: LongInt;
  begin
    UniqueSNum := 0;
    SetLength(UniqueS, UniqueSNumGuess);

    SetLength(J, PercOutNum);
    for I:=Low(J) to High(J) do
      J[I] := 0;

    while True do
      begin
        MinS := 0;
        for I:=Low(J) to High(J) do
          with PercOut[I] do
            if (J[I]<HistNum) and ((MinS<=0) or (Hist[J[I]].S<MinS))
              then
                MinS := Hist[J[I]].S;
        if MinS<=0
          then
            Break;

        UniqueSNum := UniqueSNum+1;
        if UniqueSNum>Length(UniqueS)
          then
            SetLength(UniqueS, 2*UniqueSNum+1);

        UniqueS[UniqueSNum-1] := MinS;

        for I:=Low(J) to High(J) do
          with PercOut[I] do
            if (J[I]<HistNum) and (Hist[J[I]].S<=MinS)
              then
                J[I] := J[I]+1;
      end;

    if Length(UniqueS)<>UniqueSNum
      then
        SetLength(UniqueS, UniqueSNum);
  end;
  procedure ExpandHist;
  var
    I, J, K: LongInt;
  begin
    for I:=Low(PercOut) to High(PercOut) do
      with PercOut[I] do
        begin
          J := HistNum-1;

          HistNum := UniqueSNum;
          SetLength(Hist, HistNum); // Expand

          for K:=High(UniqueS) downto Low(UniqueS) do
            if Hist[J].S=UniqueS[K]
              then
                begin
                  Hist[K] := Hist[J];
                  J := J-1;
                end
              else
                begin
                  Hist[K] := EmptyHistEntry;
                  Hist[K].S := UniqueS[K];
                end;
        end;
  end;
  procedure AverageHist;
  var
    I, J, K, N: LongInt;
    Alpha: Extended;
  begin
    for I:=Low(PercOut) to High(PercOut) do
      with PercOut[I] do
        if not Skip
          then
            begin
              N := 1;
              for J:=I+1 to High(PercOut) do
                if Abs(PercOut[J].P-P)<1e-14 //PercOut[J].P=P // BEWARE: AD HOC TOLERANCE!
                  then
                    begin
                      // Average the histograms // DOES NOT WORK FOR DIFFERENT MaxR!
                      N := N+1;
                      Alpha := 1/N;
                      for K:=Low(UniqueS) to High(UniqueS) do
                        with Hist[K] do
                          begin
                            Num  := (1-Alpha)*Num  + Alpha*PercOut[J].Hist[K].Num;
                            Prob := (1-Alpha)*Prob + Alpha*PercOut[J].Hist[K].Prob;
                          end;
                      PercOut[J].Skip := True;
                    end;
            end;
  end;
begin
  FindUniqueS;
  ExpandHist;
  if BoolAverage
    then
      AverageHist;
end;

procedure WriteOutput;
var
  I, J: LongInt;
  First: Boolean;
begin
  WriteLn('Cluster sizes S');
  First := True;
  for J:=Low(UniqueS) to High(UniqueS) do
    begin
      if First
        then
          First := False
        else
          Write(Delim);
      Write(UniqueS[J]);
    end;
  WriteLn;

  WriteLn;

  WriteLn('P', Delim, 'Number[S]');
  for I:=Low(PercOut) to High(PercOut) do
    with PercOut[I] do
      if not Skip
        then
          begin
            Write(P:0:10);
            for J:=Low(Hist) to High(Hist) do
              if BoolAverage //Write(Delim, Hist[J].Num);
                then
                  Write(Delim, Hist[J].Num:0:5) //
                else
                  Write(Delim, Round(Hist[J].Num));
            WriteLn;
          end;

  WriteLn;

  WriteLn('P', Delim, 'Probability[S]');
  for I:=Low(PercOut) to High(PercOut) do
    with PercOut[I] do
      if not Skip
        then
          begin
            Write(P:0:10);
            for J:=Low(Hist) to High(Hist) do
              Write(Delim, Hist[J].Prob:0:10);
            WriteLn;
          end;

  WriteLn;

  WriteLn('P', Delim, 'P*Probability[S]');
  for I:=Low(PercOut) to High(PercOut) do
    with PercOut[I] do
      if not Skip
        then
          begin
            Write(P:0:10);
            for J:=Low(Hist) to High(Hist) do
              Write(Delim, (P * Hist[J].Prob):0:10);
            WriteLn;
          end;
end;

procedure CleanUp;
var
  I: LongInt;
begin
  for I:=Low(PercOut) to High(PercOut) do
    SetLength(PercOut[I].Hist, 0);
  SetLength(PercOut, 0);
  SetLength(UniqueS, 0);
end;

procedure TMergeHistograms.DoRun;
var
  ErrorMsg: String;
begin
  // quick check parameters
  ErrorMsg:=CheckOptions('ha', 'help avg');
  if ErrorMsg<>'' then begin
    ShowException(Exception.Create(ErrorMsg));
    Terminate;
    Exit;
  end;

  // parse parameters
  if HasOption('h', 'help') then begin
    WriteHelp;
    Terminate;
    Exit;
  end;
  BoolAverage := HasOption('a', 'avg');

  ReadInput;
  FindSolution;
  WriteOutput;
  CleanUp;

  // stop program loop
  Terminate;
end;

constructor TMergeHistograms.Create(TheOwner: TComponent);
begin
  inherited Create(TheOwner);
  StopOnException:=True;
end;

destructor TMergeHistograms.Destroy;
begin
  inherited Destroy;
end;

procedure TMergeHistograms.WriteHelp;
begin
  { add your help code here }
  writeln('Usage: ', ExeName, ' -h');
  writeln('Usage: ', ExeName, ' -avg');
end;

var
  Application: TMergeHistograms;
begin
  Application:=TMergeHistograms.Create(nil);
  Application.Title:='MergeHistograms';
  Application.Run;
  Application.Free;
end.

