Merge branch 'feature/34-add-lse+softmax-fns' into develop

Fixes #34

Author: delphidabbler

collection/663.dat

@@ -0,0 +1,15 @@
+function LSE(const A: array of Double): Double;
+var
+  MaxElem: Double;
+  Elem: Double;
+  Sum: Double;
+begin
+  if System.Length(A) = 0 then
+    raise SysUtils.EArgumentException.Create('Empty array');
+  // Using the centering "trick": see https://rpubs.com/FJRubio/LSE
+  MaxElem := MaxOfArray(A);
+  Sum := 0.0;
+  for Elem in A do
+    Sum := Sum + System.Exp(Elem - MaxElem);
+  Result := System.Ln(Sum) + MaxElem;
+end;

collection/664.dat

@@ -0,0 +1,10 @@
+function SoftMax(const A: array of Double): Types.TDoubleDynArray;
+var
+  LSEOfA: Double;
+  Idx: Integer;
+begin
+  LSEOfA := LSE(A); // raise EArgumentException if A is empty
+  System.SetLength(Result, System.Length(A));
+  for Idx := 0 to Pred(System.Length(A)) do
+    Result[Idx] := System.Exp(A[Idx] - LSEOfA);
+end;

collection/maths.ini

@@ -1841,3 +1841,33 @@ AdvancedTest.URL="https://github.com/delphidabbler/code-snippets/tree/master/tes
 Snip=662.dat
 DelphiXE=Y
 Delphi12A=Y
+
+[LSE]
+DisplayName=LSE
+DescEx="<p>Returns the logarithm of the sum of the exponentials of the given array of floating pointing point numbers.</p><p>An <var>EArgumentException</var> exception is raised if the array is empty.</p>"
+Extra="<p>The mathematical definition of <em>LSE</em> is <em>LSE(x1,...,xn) = log(exp(x1) + ... + exp(xn))</em>. The version of the algorithm used here uses an <a href="https://gregorygundersen.com/blog/2020/02/09/log-sum-exp/">algebraic trick</a> to minimise the risk of overflow.</p><p>For more information see <a href="https://nhigham.com/2021/01/05/what-is-the-log-sum-exp-function/">What Is the Log-Sum-Exp Function?</a> by Nick Higham.</p>"
+Kind=routine
+Units=SysUtils
+Depends=MaxOfArray_Double
+SeeAlso=Softmax
+TestInfo=advanced
+AdvancedTest.Level=unit-tests
+AdvancedTest.URL="https://github.com/delphidabbler/code-snippets/tree/master/tests/Cat-Maths"
+Snip=663.dat
+DelphiXE=Y
+Delphi12A=Y
+
+[Softmax]
+DisplayName=Softmax
+DescEx="<p>Applies the <a href="https://en.wikipedia.org/wiki/Softmax_function"><em>softmax</em></a> function to each element of floating point array <var>A</var> and normalizes them into a probability distribution proportional to the exponentials of the elements of <var>A</var>. The normalised values are returned as an array of the same size as <var>A</var>.</p><p>An <var>EArgumentException</var> exception is raised if <var>A</var> is empty.</p>"
+Extra="<p>The <em>softmax</em> function is often used in statistics and in machine learning. It is closely related to the <em>log-sum-exp</em> (aka <em>LSE</em>) function since, for vector <strong>x</strong>=(x1,...,xn), <em>softmax(<strong>x</strong>)=exp[<strong>x</strong>&#8722;LSE(<strong>x</strong>)]</em>.</p>"
+Kind=routine
+Units=SysUtils
+Depends=LSE
+SeeAlso=LSE,MaxOfArray_Double
+TestInfo=advanced
+AdvancedTest.Level=unit-tests
+AdvancedTest.URL="https://github.com/delphidabbler/code-snippets/tree/master/tests/Cat-Maths"
+Snip=664.dat
+DelphiXE=Y
+Delphi12A=Y

tests/Cat-Maths/TestUMathsCatSnippets.pas

@@ -23,6 +23,7 @@   TestMathsCatSnippets = class(TTestCase)
     procedure TestPowNZN_EOverflow;
     procedure TestDigitPowerSum_EOverflow;
     procedure TestDigitPowerSum_EArgumentException;
+    procedure TestLSE_EArgumentException;
     function EqualArrays(const Left, Right: TBytes): Boolean;
     function ReverseArray(const A: TBytes): TBytes;
   published
@@ -52,7 +53,7 @@   TestMathsCatSnippets = class(TTestCase)
     procedure TestMaxOfArray_Integer;
     procedure TestMaxOfArray_Int64;
     procedure TestMaxOfArray_Single;
-    procedure TestMaxOfArray_Double;
+    procedure TestMaxOfArray_Double;  // required by LSE
     procedure TestMaxOfArray_Extended;
     procedure TestPowNZN;   // required by DigitPowerSum
     procedure TestPowNZZ;
@@ -83,6 +84,8 @@   TestMathsCatSnippets = class(TTestCase)
     procedure TestDigitPowerSum;  // required by IsNarcissistic
     procedure TestIsPalindromic;
     procedure TestIsNarcissistic;
+    procedure TestLSE;    // required by SoftMax
+    procedure TestSoftMax;
   end;
 
 implementation
@@ -753,6 +756,31 @@ procedure TestMathsCatSnippets.TestLCD;
   CheckEquals(9, LCD(-9, -9), 'LCD(-9, -9)');
 end;
 
+procedure TestMathsCatSnippets.TestLSE;
+const
+  Fudge = 0.000001;
+  A1: array [1..7] of Double = (-35.0, 20.78, 42.56, -27.8, 41.576, 0.0, 57.945);
+  A2: array [1..7] of Double = (-35.0, 20.78, 42.56, -27.8, 41.576, 0.0, 20.78);
+  A5: array [1..3] of Double = (-430.0, -399.83, -300.00);
+  A6: array [1..10] of Double = (-12.0, 4.0, -6.0, 11.0, 10.0, 3.0, -3.0, 9.0, -8.0, 7.0);
+begin
+  // Hand calculated
+  CheckTrue(SameValue(57.945000285961067157769252279369, LSE(A1)), '#1');
+  // Calculated using http://mycalcsolutions.com/calculator?mathematics;stat_prob;softmax
+  CheckTrue(SameValue(42.87759, LSE(A2), Fudge), '#2');
+  CheckTrue(SameValue(-35.0, LSE([-35.0]), Fudge), '#3');
+  CheckTrue(SameValue(0.0, LSE([0.0]), Fudge), '#4');
+  CheckTrue(SameValue(-300.0, LSE(A5), Fudge), '#5');
+  CheckTrue(SameValue(11.420537, LSE(A6), Fudge), '#6');
+  // Check empty array exception
+  CheckException(TestLSE_EArgumentException, EArgumentException, 'EArgumentException');
+end;
+
+procedure TestMathsCatSnippets.TestLSE_EArgumentException;
+begin
+  LSE([]);
+end;
+
 procedure TestMathsCatSnippets.TestMaxOfArray_Double;
 var
   A: TDoubleDynArray;
@@ -1203,6 +1231,42 @@ procedure TestMathsCatSnippets.TestResizeRect_B;
   CheckEquals(-4, RectHeight(R), '3: RectHeight');
 end;
 
+procedure TestMathsCatSnippets.TestSoftMax;
+
+  function ArraysEqual(const Left, Right: array of Double): Boolean;
+  const
+    Fudge = 0.000001;
+  var
+    Idx: Integer;
+  begin
+    Result := True;
+    if Length(Left) <> Length(Right) then
+      Exit(False);
+    for Idx := Low(Left) to High(Left) do
+      if not SameValue(Left[Idx], Right[Idx], Fudge) then
+        Exit(False);
+  end;
+const
+  A1: array [1..7] of Double = (-35.0, 20.78, 42.56, -27.8, 41.576, 0.0, 57.945);
+  E1: array [1..7] of Double = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0);
+  A2: array [1..7] of Double = (-35.0, 20.78, 42.56, -27.8, 41.576, 0.0, 20.78);
+  E2: array [1..7] of Double = (0.0, 0.0, 0.727901, 0.0, 0.272099, 0.0, 0.0);
+  A5: array [1..3] of Double = (-430.0, -399.83, -300.0);
+  E5: array [1..3] of Double = (0.0, 0.0, 1.0);
+  A6: array [1..10] of Double = (-12.0, 4.0, -6.0, 11.0, 10.0, 3.0, -3.0, 9.0, -8.0, 7.0);
+  E6: array [1..10] of Double = (0.0, 0.000599, 0.0, 0.656694, 0.241584, 0.00022, 0.000001, 0.088874, 0, 0.012028);
+  A7: array [1..3] of Double = (1430.0, 1430.83, 1440.47);
+  E7: array [1..3] of Double = (0.000028, 0.000065, 0.999907);
+begin
+  CheckTrue(ArraysEqual(E1, SoftMax(A1)), '#1');
+  CheckTrue(ArraysEqual(E2, SoftMax(A2)), '#2');
+  CheckTrue(ArraysEqual([1.0], SoftMax([-35.0])), '#3');
+  CheckTrue(ArraysEqual([1.0], SoftMax([0.0])), '#4');
+  CheckTrue(ArraysEqual(E5, SoftMax(A5)), '#6');
+  CheckTrue(ArraysEqual(E6, SoftMax(A6)), '#6');
+  CheckTrue(ArraysEqual(E7, SoftMax(A7)), '#7');
+end;
+
 procedure TestMathsCatSnippets.TestStretchRect_A;
 var
   R0, R1, R2: TRect;

tests/Cat-Maths/UMathsCatSnippets.pas

@@ -6,7 +6,7 @@
  * The unit is copyright © 2005-2024 by Peter Johnson & Contributors and is
  * licensed under the MIT License (https://opensource.org/licenses/MIT).
  * 
- * Generated on : Thu, 09 Jan 2025 15:04:31 GMT.
+ * Generated on : Fri, 10 Jan 2025 09:52:09 GMT.
  * Generated by : DelphiDabbler CodeSnip Release 4.24.0.
  * 
  * The latest version of CodeSnip is available from the CodeSnip GitHub project
@@ -227,6 +227,13 @@ function IsRectNormal(const R: Windows.TRect): Boolean;
 }
 function LCD(A, B: Integer): Integer;
 
+{
+  Returns the logarithm of the sum of the exponentials of the given array of
+  floating pointing point numbers.
+  An EArgumentException exception is raised if the array is empty.
+}
+function LSE(const A: array of Double): Double;
+
 {
   Returns the maximum value contained in the given array of double precision
   floating point values.
@@ -488,6 +495,15 @@ function SignOfInt(const Value: Int64): Integer;
 }
 procedure SimplifyFraction(var Num, Denom: Int64);
 
+{
+  Applies the softmax function to each element of floating point array A and
+  normalizes them into a probability distribution proportional to the
+  exponentials of the elements of A. The normalised values are returned as an
+  array of the same size as A.
+  An EArgumentException exception is raised if A is empty.
+}
+function SoftMax(const A: array of Double): Types.TDoubleDynArray;
+
 {
   Stretches rectangle R by the given scaling factors and returns the result.
   The rectangle's width is scaled by ScalingX and its height by ScalingY.
@@ -1198,6 +1214,27 @@ function LCD(A, B: Integer): Integer;
   Result := Abs((A * B)) div GCD(A, B);
 end;
 
+{
+  Returns the logarithm of the sum of the exponentials of the given array of
+  floating pointing point numbers.
+  An EArgumentException exception is raised if the array is empty.
+}
+function LSE(const A: array of Double): Double;
+var
+  MaxElem: Double;
+  Elem: Double;
+  Sum: Double;
+begin
+  if System.Length(A) = 0 then
+    raise SysUtils.EArgumentException.Create('Empty array');
+  // Using the centering "trick": see https://rpubs.com/FJRubio/LSE
+  MaxElem := MaxOfArray(A);
+  Sum := 0.0;
+  for Elem in A do
+    Sum := Sum + System.Exp(Elem - MaxElem);
+  Result := System.Ln(Sum) + MaxElem;
+end;
+
 {
   Returns the maximum value contained in the given array of double precision
   floating point values.
@@ -1950,6 +1987,24 @@ procedure SimplifyFraction(var Num, Denom: Int64);
   Denom := Denom div CommonFactor;
 end;
 
+{
+  Applies the softmax function to each element of floating point array A and
+  normalizes them into a probability distribution proportional to the
+  exponentials of the elements of A. The normalised values are returned as an
+  array of the same size as A.
+  An EArgumentException exception is raised if A is empty.
+}
+function SoftMax(const A: array of Double): Types.TDoubleDynArray;
+var
+  LSEOfA: Double;
+  Idx: Integer;
+begin
+  LSEOfA := LSE(A); // raise EArgumentException if A is empty
+  System.SetLength(Result, System.Length(A));
+  for Idx := 0 to Pred(System.Length(A)) do
+    Result[Idx] := System.Exp(A[Idx] - LSEOfA);
+end;
+
 {
   Stretches rectangle R by the given scaling factors and returns the result.
   The rectangle's width is scaled by ScalingX and its height by ScalingY.