program StatisticsConsole;

// Atlas Plot Pro quick start: statistics on known data, then a heat map
// and a box plot written to image files without a window. Every result is
// checked against a value worked out by hand; a mismatch raises and the
// program exits with code 1.

{$APPTYPE CONSOLE}

uses
  System.SysUtils, System.StrUtils, System.Math, System.IOUtils, Winapi.Windows, Vcl.Graphics,
  Atlas.PlotAxis, Atlas.PlotStats, Atlas.BoxPlot, Atlas.HeatMap;

procedure Expect(Condition: Boolean; const What: string);
begin
  if not Condition then raise Exception.Create('Unexpected result: ' + What);
  Writeln('  ok: ', What);
end;

function Near(const A, B: Double): Boolean;
begin
  Result := Abs(A - B) < 1E-9;
end;

procedure Quantiles;
var Values: TArray<Double>;
begin
  Writeln('Quantiles of 1..10 in shuffled order');
  Values := [7, 2, 9, 4, 5, 3, 8, 6, 1, 10];
  Writeln(Format('  Q1 %.2f  median %.2f  Q3 %.2f  mean %.2f  sd %.4f',
    [TAtlasStats.Quantile(Values, 0.25), TAtlasStats.Quantile(Values, 0.5),
     TAtlasStats.Quantile(Values, 0.75), TAtlasStats.Mean(Values), TAtlasStats.StdDev(Values)]));
  Expect(Near(TAtlasStats.Quantile(Values, 0.25), 3.25), 'first quartile interpolates between ranks 3 and 4');
  Expect(Near(TAtlasStats.Quantile(Values, 0.5), 5.5), 'median of an even count is the middle pair average');
  Expect(Near(TAtlasStats.Mean(Values), 5.5) and Near(TAtlasStats.StdDev(Values), Sqrt(82.5 / 9)), 'mean and sample standard deviation');
end;

procedure BoxStatistics;
var Box: TAtlasBoxStats;
begin
  Writeln('Box statistics of 1..9 and one slow value, 40');
  Box := TAtlasStats.Box([1, 2, 3, 4, 5, 6, 7, 8, 9, 40]);
  Writeln(Format('  box %.2f..%.2f  whiskers %.0f..%.0f  %d outlier(s)',
    [Box.Q1, Box.Q3, Box.LowWhisker, Box.HighWhisker, Length(Box.Outliers)]));
  Expect(Near(Box.Q1, 3.25) and Near(Box.Median, 5.5) and Near(Box.Q3, 7.75), 'quartiles');
  Expect(Near(Box.LowWhisker, 1) and Near(Box.HighWhisker, 9), 'whiskers stop at the last value within 1.5 IQR');
  Expect((Length(Box.Outliers) = 1) and Near(Box.Outliers[0], 40), '40 is the only outlier');
end;

procedure Histogram;
var Bars: TAtlasHistogram; I: Integer;
begin
  Writeln('Histogram of 1 2 2 3 3 3 4 4 5 in four bins');
  Bars := TAtlasStats.Histogram([1, 2, 2, 3, 3, 3, 4, 4, 5], 4);
  for I := 0 to High(Bars.Counts) do
    Writeln(Format('  [%.0f, %.0f%s: %d', [Bars.Edges[I], Bars.Edges[I + 1], IfThen(I = High(Bars.Counts), ']', ')'), Bars.Counts[I]]));
  Expect((Bars.Counts[0] = 1) and (Bars.Counts[1] = 2) and (Bars.Counts[2] = 3) and (Bars.Counts[3] = 3), 'bins are half open, the last one closed');
  Expect((Bars.BinOf(2) = 1) and (Bars.BinOf(5) = 3) and (Bars.BinOf(6) = -1), 'BinOf finds the bin; outside the range is -1');
  Expect(Length(TAtlasStats.Histogram([1, 2, 2, 3, 3, 3, 4, 4, 5]).Counts) = 5, 'Bins = 0 applies Sturges'' rule: ceil(log2 9) + 1 = 5');
end;

procedure CorrelationAndFit;
var X, Y: TArray<Double>; I: Integer; Fit: TAtlasPolynomial;
begin
  Writeln('Correlation and least-squares fits');
  Writeln(Format('  r = %.2f', [TAtlasStats.Correlation([1, 2, 3, 4, 5], [2, 1, 4, 3, 5])]));
  Expect(Near(TAtlasStats.Correlation([1, 2, 3, 4, 5], [2, 1, 4, 3, 5]), 0.8), 'Pearson r of a partly ordered pair is 0.8');
  SetLength(X, 9);
  SetLength(Y, 9);
  for I := 0 to 8 do
  begin
    X[I] := I - 4;
    Y[I] := 2 * X[I] * X[I] - 3 * X[I] + 1;   // an exact parabola
  end;
  Fit := TAtlasStats.PolynomialFit(X, Y, 2);
  Writeln(Format('  quadratic: degree %d, R squared %.4f, value at 10 = %.4f', [Fit.Degree, Fit.RSquared, Fit.Evaluate(10)]));
  Expect((Fit.Degree = 2) and Near(Fit.RSquared, 1) and (Abs(Fit.Evaluate(10) - 171) < 1E-6), 'the parabola is recovered, also outside the data');
  Fit := TAtlasStats.PolynomialFit(X, Y, 1);
  Writeln(Format('  straight line: R squared %.4f', [Fit.RSquared]));
  Expect(Abs(Fit.RSquared - 0.30474) < 1E-4, 'a line explains less than a third of a parabola');
end;

procedure Smoothing;
var Smooth: TArray<Double>;
begin
  Writeln('Exponential smoothing of 10 20 20 0 with alpha 0.5');
  Smooth := TAtlasStats.ExponentialSmoothing([10, 20, 20, 0], 0.5);
  Writeln(Format('  %.2f %.2f %.2f %.2f', [Smooth[0], Smooth[1], Smooth[2], Smooth[3]]));
  Expect(Near(Smooth[0], 10) and Near(Smooth[1], 15) and Near(Smooth[2], 17.5) and Near(Smooth[3], 8.75), 'each value moves halfway towards the input');
end;

procedure Rejections;
var Rejected: Boolean;
begin
  Writeln('Invalid input raises EAtlasPlotError');
  Rejected := False;
  try
    TAtlasStats.PolynomialFit([1, 2], [3, 4], 2);   // two points cannot fix a parabola
  except
    on EAtlasPlotError do Rejected := True;
  end;
  Expect(Rejected, 'a fit with too few distinct arguments is rejected');
end;

procedure ExportCharts;
var Map: TAtlasHeatMap; Box: TAtlasBoxPlot; Svg: string;
begin
  Writeln('Charts rendered without a window');
  Map := TAtlasHeatMap.Create(nil);
  try
    Map.Title := 'Cache hit rate (%) by node and hour';
    Map.RowLabels.CommaText := 'node-a,node-b,node-c';
    Map.ColumnLabels.CommaText := '00,06,12,18';
    // Percentages, because FormatFloat prints a % sign as it is and does not multiply by 100.
    Map.LoadMatrix([[91, 88, 72, 80], [93, NaN, 75, 83], [89, 86, 70, 79]]);
    Map.LowColor := clWhite;
    Map.HighColor := clBlack;
    Map.SetScale(50, 100);
    Map.ShowValues := True;
    Map.ValueFormat := '0"%"';
    Map.SaveToFile('heatmap.png', 480, 300);
    Expect(TFile.GetSize('heatmap.png') > 1000, 'heatmap.png written');
    Expect(Map.IsEmpty(1, 1) and (Map.ColorOf(75) = TColor(RGB(128, 128, 128))), 'NaN is an empty cell; 75 is halfway on the 50..100 scale');
  finally
    Map.Free;
  end;
  Box := TAtlasBoxPlot.Create(nil);
  try
    Box.Title := 'Response time by service';
    Box.ValueAxis.Title := 'ms';
    Box.AddCategory('Gateway', [12, 15, 15, 17, 18, 18, 18, 20, 21, 60]);
    Box.AddCategory('Search', [30, 31, 35, 36, 38, 40, 41, 45]);
    Box.SaveToFile('boxplot.svg', 480, 300);
    Svg := TFile.ReadAllText('boxplot.svg', TEncoding.UTF8);
    Expect((Pos('>Gateway<', Svg) > 0) and (Pos('>Search<', Svg) > 0), 'boxplot.svg names both categories');
    Expect(Near(Box.Stats(0).Median, 18) and (Length(Box.Stats(0).Outliers) = 1), 'the chart keeps the statistics of each sample');
  finally
    Box.Free;
  end;
end;

begin
  try
    FormatSettings := TFormatSettings.Invariant;
    Quantiles;
    BoxStatistics;
    Histogram;
    CorrelationAndFit;
    Smoothing;
    Rejections;
    ExportCharts;
    Writeln('All results verified.');
  except
    on E: Exception do
    begin
      Writeln(E.ClassName, ': ', E.Message);
      ExitCode := 1;
    end;
  end;
end.