Implementation, testing, and examples of Fisher-Yates

dfyconfig.toml

@@ -2,4 +2,4 @@ includes = ["src/**/*.dfy"]
 excludes = []
 
 [options]
-enforce-determinism = true
+enforce-determinism = false

docs/dafny/ExamplesExternUniformPowerOfTwo.dfy

@@ -6,6 +6,7 @@
 module RandomExamples.ExternUniform {
   import Rationals
   import DafnyVMC
+  import Helper
 
   method Main()
     decreases *
@@ -124,5 +125,15 @@ module RandomExamples.ExternUniform {
     // https://www.wolframalpha.com/input?i=ReplaceAll%5BE%5E%28-x%5E2+%2F+%282+*%CF%83%5E2%29%29+%2F+Sum%5BE%5E%28-y%5E2%2F%282+%CF%83%5E2%29%29%2C+%7By%2C+-Infinity%2C+Infinity%7D%5D%2C+%7Bx+-%3E+0%2C+%CF%83+-%3E+1.4%7D%5D
     // https://www.wolframalpha.com/input?i=ReplaceAll%5BE%5E%28-x%5E2+%2F+%282+*%CF%83%5E2%29%29+%2F+Sum%5BE%5E%28-y%5E2%2F%282+%CF%83%5E2%29%29%2C+%7By%2C+-Infinity%2C+Infinity%7D%5D%2C+%7Bx+-%3E+1%2C+%CF%83+-%3E+1.4%7D%5D
     print "Estimated probabilities for DiscreteGaussianSample(1.4): ", count0 as real / n as real, " (should be around 0.284959) and ", count1 as real / n as real, ", ", countneg1 as real / n as real, " (should both be around 0.220797)\n";
+
+    // Fisher-Yates Example
+    print "Ten permutations of 012: ";
+    var arr: array<nat> := new nat[3](i => i); // [0, 1, 2]
+    for i := 0 to 10 {
+      var arrCopy := arr;
+      r.Shuffle(arrCopy);
+      print Helper.SeqToString(Helper.ArrayToSeq(arrCopy), Helper.NatToString), ", ";
+    }
+    print "\n";
   }
 }

docs/dafny/ExamplesFoundational.dfy

@@ -6,6 +6,7 @@
 module RandomExamples.Foundational {
   import Rationals
   import DafnyVMC
+  import Helper
 
   method Main()
     decreases *
@@ -124,5 +125,15 @@ module RandomExamples.Foundational {
     // https://www.wolframalpha.com/input?i=ReplaceAll%5BE%5E%28-x%5E2+%2F+%282+*%CF%83%5E2%29%29+%2F+Sum%5BE%5E%28-y%5E2%2F%282+%CF%83%5E2%29%29%2C+%7By%2C+-Infinity%2C+Infinity%7D%5D%2C+%7Bx+-%3E+0%2C+%CF%83+-%3E+1.4%7D%5D
     // https://www.wolframalpha.com/input?i=ReplaceAll%5BE%5E%28-x%5E2+%2F+%282+*%CF%83%5E2%29%29+%2F+Sum%5BE%5E%28-y%5E2%2F%282+%CF%83%5E2%29%29%2C+%7By%2C+-Infinity%2C+Infinity%7D%5D%2C+%7Bx+-%3E+1%2C+%CF%83+-%3E+1.4%7D%5D
     print "Estimated probabilities for DiscreteGaussianSample(1.4): ", count0 as real / n as real, " (should be around 0.284959) and ", count1 as real / n as real, ", ", countneg1 as real / n as real, " (should both be around 0.220797)\n";
+
+    // Fisher-Yates Example
+    print "Ten permutations of 012: ";
+    var arr: array<nat> := new nat[3](i => i); // [0, 1, 2]
+    for i := 0 to 10 {
+      var arrCopy := arr;
+      r.Shuffle(arrCopy);
+      print Helper.SeqToString(Helper.ArrayToSeq(arrCopy), Helper.NatToString), ", ";
+    }
+    print "\n";
   }
 }

src/Dafny-VMC.dfy

@@ -11,12 +11,13 @@ module DafnyVMC {
   import BernoulliExpNeg
   import DiscreteGaussian
   import DiscreteLaplace
+  import FisherYates
 
-  class DRandomFoundational extends Coin.Interface.Trait, UniformPowerOfTwo.Implementation.Trait, Bernoulli.Implementation.Trait, Uniform.Implementation.TraitFoundational, BernoulliExpNeg.Implementation.Trait, DiscreteGaussian.Implementation.Trait, DiscreteLaplace.Implementation.Trait {
+  class DRandomFoundational extends Coin.Interface.Trait, UniformPowerOfTwo.Implementation.Trait, Bernoulli.Implementation.Trait, Uniform.Implementation.TraitFoundational, BernoulliExpNeg.Implementation.Trait, DiscreteGaussian.Implementation.Trait, DiscreteLaplace.Implementation.Trait, FisherYates.Implementation.Trait  {
     constructor {:extern} ()
   }
 
-  class DRandomExternUniformPowerOfTwo extends Coin.Interface.Trait, UniformPowerOfTwo.Implementation.TraitExtern, Bernoulli.Implementation.Trait, Uniform.Implementation.TraitFoundational, BernoulliExpNeg.Implementation.Trait, DiscreteGaussian.Implementation.Trait, DiscreteLaplace.Implementation.Trait {
+  class DRandomExternUniformPowerOfTwo extends Coin.Interface.Trait, UniformPowerOfTwo.Implementation.TraitExtern, Bernoulli.Implementation.Trait, Uniform.Implementation.TraitFoundational, BernoulliExpNeg.Implementation.Trait, DiscreteGaussian.Implementation.Trait, DiscreteLaplace.Implementation.Trait, FisherYates.Implementation.Trait  {
     constructor {:extern} ()
   }
 }

src/Math/Helper.dfy

@@ -8,6 +8,27 @@ module Helper {
    Definitions
   ************/
 
+  function ArrayToSeq<T>(a: array<T>): (xs: seq<T>)
+    reads a
+    ensures a.Length == |xs|
+  {
+    a[..]
+  }
+
+  function NatToString(n: nat): string {
+    match n
+    case 0 => "0" case 1 => "1" case 2 => "2" case 3 => "3" case 4 => "4"
+    case 5 => "5" case 6 => "6" case 7 => "7" case 8 => "8" case 9 => "9"
+    case _ => NatToString(n / 10) + NatToString(n % 10)
+  }
+
+  function SeqToString<T>(xs: seq<T>, printer: T -> string): string {
+    if |xs| == 0 then
+      ""
+    else
+      printer(xs[0]) + SeqToString(xs[1..], printer)
+  }
+
   /*******
    Lemmas
   *******/

src/Math/Permutations.dfy

@@ -0,0 +1,50 @@
+/*******************************************************************************
+ *  Copyright by the contributors to the Dafny Project
+ *  SPDX-License-Identifier: MIT
+ *******************************************************************************/
+
+module Permutations {
+  import NatArith
+
+  /************
+   Definitions
+  ************/
+
+  function NumberOfPermutationsOf<T(==)>(s: seq<T>): nat {
+    |CalculateAllPermutationsOf(s)|
+  }
+
+  function CalculateAllPermutationsOf<T(==)>(s: seq<T>): (x: set<seq<T>>)
+    ensures forall p | p in x :: |p| == |s|
+  {
+    if |s| == 0 then
+      {s}
+    else
+      set p, i | p in CalculateAllPermutationsOf(s[1..]) && 0 <= i <= |s|-1 :: InsertAt(p, s[0], i)
+  }
+
+  function InsertAt<T>(s: seq<T>, x: T, i: nat): seq<T>
+    requires i <= |s|
+  {
+    s[..i] + [x] + s[i..]
+  }
+
+  /*******
+   Lemmas
+  *******/
+
+  lemma CalculateAllPermutationsOfIsNonEmpty<T>(s: seq<T>)
+    ensures s in CalculateAllPermutationsOf(s)
+    ensures NumberOfPermutationsOf(s) > 0
+  {
+    assert s in CalculateAllPermutationsOf(s) by {
+      if |s| == 0 {
+      } else {
+        assert s[1..] in CalculateAllPermutationsOf(s[1..]) by {
+          CalculateAllPermutationsOfIsNonEmpty(s[1..]);
+        }
+        assert InsertAt(s[1..], s[0], 0) == s;
+      }
+    }
+  }
+}

src/Util/FisherYates/Implementation.dfy

@@ -0,0 +1,33 @@
+/*******************************************************************************
+ *  Copyright by the contributors to the Dafny Project
+ *  SPDX-License-Identifier: MIT
+ *******************************************************************************/
+
+module FisherYates.Implementation {
+  import Interface
+
+  trait {:termination false} Trait extends Interface.Trait {
+
+    method Shuffle<T(!new)>(a: array<T>)
+      decreases *
+      modifies this, a
+    {
+      var i := 0;
+      while i < a.Length {
+        var j := UniformIntervalSample(i, a.Length);
+        Swap(a, i, j);
+        i := i + 1;
+      }
+    }
+
+    method Swap<T>(a: array<T>, i: nat, j: nat)
+      modifies a
+      requires i <= j
+      requires 0 <= i < a.Length
+      requires 0 <= j < a.Length
+    {
+      a[i], a[j] := a[j], a[i];
+    }
+
+  }
+}

src/Util/FisherYates/Interface.dfy

@@ -0,0 +1,16 @@
+/*******************************************************************************
+ *  Copyright by the contributors to the Dafny Project
+ *  SPDX-License-Identifier: MIT
+ *******************************************************************************/
+
+module FisherYates.Interface {
+  import Uniform
+
+  trait {:termination false} Trait extends Uniform.Interface.Trait {
+
+    method Shuffle<T(!new)>(a: array<T>)
+      decreases *
+      modifies this, a
+
+  }
+}

tests/Tests.dfy

@@ -13,17 +13,10 @@ module Tests {
   import DiscreteLaplace
   import DiscreteGaussian
   import NatArith
-
-  function Abs(x: real): real {
-    if x < 0.0 then -x else x
-  }
-
-  function NatToString(n: nat): string {
-    match n
-    case 0 => "0" case 1 => "1" case 2 => "2" case 3 => "3" case 4 => "4"
-    case 5 => "5" case 6 => "6" case 7 => "7" case 8 => "8" case 9 => "9"
-    case _ => NatToString(n / 10) + NatToString(n % 10)
-  }
+  import RealArith
+  import Permutations
+  import FisherYates
+  import Helper
 
   method TestBernoulliIsWithin4SigmaOfTrueMean(
     n: nat,
@@ -46,17 +39,17 @@ module Tests {
     requires n > 0
   {
     var empiricalMean := empiricalSum / n as real;
-    var diff := Abs(empiricalMean - trueMean);
+    var diff := RealArith.Abs(empiricalMean - trueMean);
     var threshold := 4.0 * 4.0 * trueVariance / n as real;
-    if diff * diff >= threshold {
+    if diff * diff > threshold {
       print "Test failed: ", description, "\n";
       print "True mean: ", trueMean, "\n";
       print "Empirical mean: ", empiricalMean, "\n";
       print "Difference between empirical and true mean: ", diff, "\n";
       print "squared difference: ", diff * diff, "\n";
       print "sigma squared:      ", trueVariance / n as real, "\n";
     }
-    expect diff * diff < threshold, "Empirical mean should be within 3 sigma of true mean. This individual test may fail with probability of about 6.3e-5.";
+    expect diff * diff <= threshold, "Empirical mean should be within 3 sigma of true mean. This individual test may fail with probability of about 6.3e-5.";
   }
 
   method TestCoin(n: nat, r: Coin.Interface.Trait)
@@ -90,9 +83,9 @@ module Tests {
       sum := sum + l;
     }
     for i := 0 to NatArith.Power(2, k) {
-      TestBernoulliIsWithin4SigmaOfTrueMean(n, a[i] as real, 1.0 / (m as real), "p(" + NatToString(i) + ")");
+      TestBernoulliIsWithin4SigmaOfTrueMean(n, a[i] as real, 1.0 / (m as real), "p(" + Helper.NatToString(i) + ")");
     }
-    TestEmpiricalIsWithin4SigmaOfTrueMean(n, sum as real, (m - 1) as real / 2.0, (m * m - 1) as real / 12.0, "mean of UniformPowerOfTwo(" + NatToString(u) + ")");
+    TestEmpiricalIsWithin4SigmaOfTrueMean(n, sum as real, (m - 1) as real / 2.0, (m * m - 1) as real / 12.0, "mean of UniformPowerOfTwo(" + Helper.NatToString(u) + ")");
   }
 
   method TestUniformPowerOfTwoMean(n: nat, u: nat, r: UniformPowerOfTwo.Interface.Trait)
@@ -109,7 +102,7 @@ module Tests {
       expect 0 <= l < m, "sample not in the right bound";
       sum := sum + l;
     }
-    TestEmpiricalIsWithin4SigmaOfTrueMean(n, sum as real, (m - 1) as real / 2.0, (m * m - 1) as real / 12.0, "mean of UniformPowerOfTwo(" + NatToString(u) + ")");
+    TestEmpiricalIsWithin4SigmaOfTrueMean(n, sum as real, (m - 1) as real / 2.0, (m * m - 1) as real / 12.0, "mean of UniformPowerOfTwo(" + Helper.NatToString(u) + ")");
   }
 
   method TestUniform(n: nat, u: nat, r: Uniform.Interface.Trait)
@@ -128,9 +121,9 @@ module Tests {
       sum := sum + l;
     }
     for i := 0 to u {
-      TestBernoulliIsWithin4SigmaOfTrueMean(n, a[i] as real, 1.0 / (u as real), "p(" + NatToString(i) + ")");
+      TestBernoulliIsWithin4SigmaOfTrueMean(n, a[i] as real, 1.0 / (u as real), "p(" + Helper.NatToString(i) + ")");
     }
-    TestEmpiricalIsWithin4SigmaOfTrueMean(n, sum as real, (u - 1) as real / 2.0, (u * u - 1) as real / 12.0, "mean of Uniform(" + NatToString(u) + ")");
+    TestEmpiricalIsWithin4SigmaOfTrueMean(n, sum as real, (u - 1) as real / 2.0, (u * u - 1) as real / 12.0, "mean of Uniform(" + Helper.NatToString(u) + ")");
   }
 
   method TestUniformMean(n: nat, u: nat, r: Uniform.Interface.Trait)
@@ -146,7 +139,7 @@ module Tests {
       expect 0 <= l < u, "sample not in the right bound";
       sum := sum + l;
     }
-    TestEmpiricalIsWithin4SigmaOfTrueMean(n, sum as real, (u - 1) as real / 2.0, (u * u - 1) as real / 12.0, "mean of Uniform(" + NatToString(u) + ")");
+    TestEmpiricalIsWithin4SigmaOfTrueMean(n, sum as real, (u - 1) as real / 2.0, (u * u - 1) as real / 12.0, "mean of Uniform(" + Helper.NatToString(u) + ")");
   }
 
   method TestUniformInterval(n: nat, r: Uniform.Interface.Trait)
@@ -291,4 +284,38 @@ module Tests {
     var varianceBound := 1.4 * 1.4; // variance of DiscreteGaussian(1.4) is < 1.4^2
     TestEmpiricalIsWithin4SigmaOfTrueMean(n, sum as real, 0.0, varianceBound, "mean");
   }
+
+  // Shuffles an array `a` n-times and verifies that for each permutation that occurs x-times, roughly x/n == 1/|Permutations.NumberOfPermutationsOf(a[..])|
+  method TestFisherYates<T(==, !new)>(n: nat, a: array<T>, r: FisherYates.Interface.Trait, printer: T -> string) 
+    decreases *
+    modifies r
+    modifies a
+    requires n > 0
+  {
+    var aAsSeq: seq<T> := Helper.ArrayToSeq(a);
+    var numberOfPermutations: nat := Permutations.NumberOfPermutationsOf(aAsSeq);
+    var numberOfObservedPermutations: map<seq<T>, nat> := map[];
+    Permutations.CalculateAllPermutationsOfIsNonEmpty(aAsSeq);
+
+    for i := 0 to n {
+      var aCopy := a;
+      r.Shuffle(aCopy);
+      var aCopyAsSeq := Helper.ArrayToSeq(aCopy);
+      if aCopyAsSeq in numberOfObservedPermutations {
+        numberOfObservedPermutations := numberOfObservedPermutations[aCopyAsSeq := numberOfObservedPermutations[aCopyAsSeq] + 1];
+      } else {
+        numberOfObservedPermutations := numberOfObservedPermutations[aCopyAsSeq := 1];
+      }
+    }
+
+    var items := numberOfObservedPermutations.Items;
+
+    while items != {}
+      decreases |items| 
+     {
+      var item :| item in items;
+      items := items - {item};
+      TestBernoulliIsWithin4SigmaOfTrueMean(n, item.1 as real, 1.0 / (numberOfPermutations as real), "p(" + Helper.SeqToString(item.0, printer) + ")");
+    }
+  }
 }

tests/TestsExternUniformPowerOfTwo.dfy

@@ -7,6 +7,7 @@ module TestsExternUniform {
   import NatArith
   import DafnyVMC
   import Tests
+  import Helper
 
   method {:test} TestCoin() {
     var r := new DafnyVMC.DRandomExternUniformPowerOfTwo();
@@ -106,4 +107,28 @@ module TestsExternUniform {
     var r := new DafnyVMC.DRandomExternUniformPowerOfTwo();
     Tests.TestDiscreteGaussian(1_000_000, r);
   }
+
+  method {:test} TestFisherYatesIdentity()
+    decreases *
+  {
+    var a: array<nat> := new nat[4](i => i); // [0, 1, 2, 3]
+    var r := new DafnyVMC.DRandomExternUniformPowerOfTwo();
+    Tests.TestFisherYates(1_000_000, a, r, Helper.NatToString);
+  }
+
+  method {:test} TestFisherYatesConstant()
+    decreases *
+  {
+    var a: array<nat> := new nat[4](i => 0); // [0, 0, 0, 0]
+    var r := new DafnyVMC.DRandomExternUniformPowerOfTwo();
+    Tests.TestFisherYates(1_000_000, a, r, Helper.NatToString);
+  }
+
+  method {:test} TestFisherYatesMixed()
+    decreases *
+  {
+    var a: array<nat> := new nat[] [0, 1, 1, 2]; // [0, 1, 1, 2]
+    var r := new DafnyVMC.DRandomExternUniformPowerOfTwo();
+    Tests.TestFisherYates(1_000_000, a, r, Helper.NatToString);
+  }
 }