Skip to content

Commit

Permalink
Add model and state correctness property of BernoulliExpNeg (#87)
Browse files Browse the repository at this point in the history
<small>By submitting this pull request, I confirm that my contribution
is made under the terms of the [MIT
license](https://github.com/dafny-lang/dafny/blob/master/LICENSE.txt).</small>
  • Loading branch information
fzaiser authored Oct 6, 2023
1 parent 69d62ac commit 3007adf
Show file tree
Hide file tree
Showing 7 changed files with 194 additions and 20 deletions.
10 changes: 10 additions & 0 deletions audit.log
Original file line number Diff line number Diff line change
@@ -1,7 +1,17 @@
src/Distributions/BernoulliExpNeg/Correctness.dfy(13,17): Correctness: Declaration has explicit `{:axiom}` attribute.
src/Distributions/BernoulliExpNeg/Correctness.dfy(17,17): SampleIsIndepFn: Declaration has explicit `{:axiom}` attribute.
src/Distributions/BernoulliExpNeg/Implementation.dfy(34,6): BernoulliExpNegSample: Definition has `assume {:axiom}` statement in body.
src/Distributions/BernoulliExpNeg/Implementation.dfy(52,6): BernoulliExpNegSampleCaseLe1: Definition has `assume {:axiom}` statement in body.
src/Distributions/BernoulliExpNeg/Model.dfy(30,4): GammaReductionLoop: Definition has `assume {:axiom}` statement in body.
src/Distributions/BernoulliExpNeg/Model.dfy(60,4): GammaLe1Loop: Definition has `assume {:axiom}` statement in body.
src/Distributions/Coin/Interface.dfy(20,6): CoinSample: Definition has `assume {:axiom}` statement in body.
src/Distributions/Uniform/Correctness.dfy(157,17): SampleIsIndepFn: Declaration has explicit `{:axiom}` attribute.
src/Distributions/Uniform/Implementation.dfy(43,6): UniformSample: Definition has `assume {:axiom}` statement in body.
src/Distributions/Uniform/Model.dfy(34,17): SampleTerminates: Declaration has explicit `{:axiom}` attribute.
src/Math/Exponential.dfy(11,17): EvalOne: Declaration has explicit `{:axiom}` attribute.
src/Math/Exponential.dfy(2,26): Exp: Declaration has explicit `{:axiom}` attribute.
src/Math/Exponential.dfy(4,17): FunctionalEquation: Declaration has explicit `{:axiom}` attribute.
src/Math/Exponential.dfy(7,17): Increasing: Declaration has explicit `{:axiom}` attribute.
src/Math/MeasureTheory.dfy(150,17): CountableSumOfZeroesIsZero: Declaration has explicit `{:axiom}` attribute.
src/Math/MeasureTheory.dfy(25,4): CountableSum: Definition has `assume {:axiom}` statement in body.
src/ProbabilisticProgramming/Independence.dfy(28,27): IsIndepFn: Declaration has explicit `{:axiom}` attribute.
Expand Down
20 changes: 20 additions & 0 deletions src/Distributions/BernoulliExpNeg/Correctness.dfy
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
/*******************************************************************************
* Copyright by the contributors to the Dafny Project
* SPDX-License-Identifier: MIT
*******************************************************************************/

module BernoulliExpNeg.Correctness {
import Rationals
import Exponential
import RandomNumberGenerator
import Independence
import Model

lemma {:axiom} Correctness(gamma: Rationals.Rational)
requires 0 <= gamma.numer
ensures RandomNumberGenerator.mu(iset s | Model.Sample(gamma)(s).0) == Exponential.Exp(-Rationals.ToReal(gamma))

lemma {:axiom} SampleIsIndepFn(gamma: Rationals.Rational)
requires 0 <= gamma.numer
ensures Independence.IsIndepFn(Model.Sample(gamma))
}
52 changes: 33 additions & 19 deletions src/Distributions/BernoulliExpNeg/Implementation.dfy
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
module BernoulliExpNeg.Implementation {
import Rationals
import Interface
import Model

trait {:termination false} Trait extends Interface.Trait {

Expand All @@ -14,28 +15,41 @@ module BernoulliExpNeg.Implementation {
modifies this
requires gamma.numer >= 0
decreases *
ensures (c, s) == Model.Sample(gamma)(old(s))
{
if gamma.numer <= gamma.denom {
var k := 1;
var a := BernoulliSample(Rationals.Rational(gamma.numer, k * gamma.denom));
while a
decreases *
{
k := k + 1;
a := BernoulliSample(Rationals.Rational(gamma.numer, k * gamma.denom));
}
c := k % 2 == 1;
var gamma' := gamma;
var b := true;
while b && gamma'.numer >= gamma'.denom
decreases gamma'.numer
invariant gamma'.numer >= 0
{
b := BernoulliExpNegSampleCaseLe1(Rationals.Int(1));
gamma' := Rationals.Rational(gamma'.numer - gamma'.denom, gamma'.denom);
}
if b {
c:= BernoulliExpNegSampleCaseLe1(gamma');
} else {
var k := 1;
while k <= Rationals.Floor(gamma) {
var b := BernoulliExpNegSample(Rationals.Int(1));
if !b {
return false;
}
k := k + 1;
}
c:= BernoulliExpNegSample(Rationals.FractionalPart(gamma));
c := false;
}
assume {:axiom} (c, s) == Model.Sample(gamma)(old(s));
}

method BernoulliExpNegSampleCaseLe1(gamma: Rationals.Rational) returns (c: bool)
modifies this
requires 0 <= gamma.numer <= gamma.denom
decreases *
ensures (c, s) == Model.SampleGammaLe1(gamma)(old(s))
{
var k := 0;
var a := true;
while a
decreases *
{
k := k + 1;
a := BernoulliSample(Rationals.Rational(gamma.numer, k * gamma.denom));
}
c := k % 2 == 1;
assume {:axiom} (c, s) == Model.SampleGammaLe1(gamma)(old(s));
}

}
Expand Down
2 changes: 2 additions & 0 deletions src/Distributions/BernoulliExpNeg/Interface.dfy
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,15 @@
module BernoulliExpNeg.Interface {
import Rationals
import Bernoulli
import Model

trait {:termination false} Trait extends Bernoulli.Interface.Trait {

method BernoulliExpNegSample(gamma: Rationals.Rational) returns (c: bool)
modifies this
decreases *
requires gamma.numer >= 0
ensures (c, s) == Model.Sample(gamma)(old(s))

}
}
76 changes: 76 additions & 0 deletions src/Distributions/BernoulliExpNeg/Model.dfy
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
/*******************************************************************************
* Copyright by the contributors to the Dafny Project
* SPDX-License-Identifier: MIT
*******************************************************************************/

module BernoulliExpNeg.Model {
import Rationals
import Uniform
import Monad
import WhileAndUntil
import BernoulliModel = Bernoulli.Model

function Sample(gamma: Rationals.Rational): Monad.Hurd<bool>
requires gamma.denom != 0
requires gamma.numer >= 0
{
Monad.Bind(
GammaReductionLoop(gamma),
(bgamma: (bool, Rationals.Rational)) =>
if bgamma.0 then
SampleGammaLe1(bgamma.1)
else
Monad.Return(false)
)
}

function GammaReductionLoop(gamma: Rationals.Rational): Monad.Hurd<(bool, Rationals.Rational)>
requires gamma.numer >= 0
{
assume {:axiom} false; // assume termination
WhileAndUntil.ProbWhile(
(bgamma: (bool, Rationals.Rational)) => bgamma.0 && bgamma.1.denom <= bgamma.1.numer,
GammaReductionLoopIter,
(true, gamma)
)
}

function GammaReductionLoopIter(bgamma: (bool, Rationals.Rational)): Monad.Hurd<(bool, Rationals.Rational)>
requires bgamma.1.numer >= 0
{
Monad.Bind(
SampleGammaLe1(Rationals.Int(1)),
b' => Monad.Return((b', Rationals.Rational(bgamma.1.numer - bgamma.1.denom, bgamma.1.denom)))
)
}

function SampleGammaLe1(gamma: Rationals.Rational): Monad.Hurd<bool>
{
if 0 <= gamma.numer <= gamma.denom
then Monad.Bind(
GammaLe1Loop(gamma, (true, 0)),
(ak: (bool, nat)) => Monad.Return(ak.1 % 2 == 1)
)
else Monad.Return(false) // to keep this function total, we return a dummy value here
}

function GammaLe1Loop(gamma: Rationals.Rational, ak: (bool, nat)): Monad.Hurd<(bool, nat)>
requires 0 <= gamma.numer <= gamma.denom
{
assume {:axiom} false; // assume termination
WhileAndUntil.ProbWhile(
(ak: (bool, nat)) => ak.0,
(ak: (bool, nat)) => GammaLe1LoopIter(gamma, ak),
ak
)
}

function GammaLe1LoopIter(gamma: Rationals.Rational, ak: (bool, nat)): Monad.Hurd<(bool, nat)>
requires 0 <= gamma.numer <= gamma.denom
{
Monad.Bind(
BernoulliModel.Sample(gamma.numer, (ak.1 + 1) * gamma.denom),
(a': bool) => Monad.Return((a', ak.1 + 1))
)
}
}
2 changes: 1 addition & 1 deletion src/Distributions/Uniform/Implementation.dfy
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ module Uniform.Implementation {
while u >= n
decreases *
invariant Model.Sample(n)(old(s)) == Model.Sample(n)(prev_s)
invariant (u, s) == UniformPowerOfTwo.Model.Sample(2 * n)(prev_s)
invariant (u, s) == Model.Proposal(n)(prev_s)
{
prev_s := s;
u := UniformPowerOfTwoSample(2 * n);
Expand Down
52 changes: 52 additions & 0 deletions src/Math/Exponential.dfy
Original file line number Diff line number Diff line change
@@ -0,0 +1,52 @@
module Exponential {
ghost function {:axiom} Exp(x: real): real

lemma {:axiom} FunctionalEquation(x: real, y: real)
ensures Exp(x + y) == Exp(x) * Exp(y)

lemma {:axiom} Increasing(x: real, y: real)
requires x < y
ensures Exp(x) < Exp(y)

lemma {:axiom} EvalOne()
ensures 2.718281828 <= Exp(1.0) <= 2.718281829

lemma Positive(x: real)
ensures Exp(x) > 0.0
{
assert Exp(x) >= 0.0 by {
var sqrt := Exp(x / 2.0);
calc {
Exp(x);
{ FunctionalEquation(x / 2.0, x / 2.0); }
sqrt * sqrt;
>=
0.0;
}
}
if Exp(x) == 0.0 {
calc {
0.0;
Exp(x);
< { Increasing(x, x + 1.0); }
Exp(x + 1.0);
{ FunctionalEquation(x, 1.0); }
Exp(x) * Exp(1.0);
==
0.0;
}
}
}

lemma EvalZero()
ensures Exp(0.0) == 1.0
{
var one := Exp(0.0);
assert one > 0.0 by {
Positive(0.0);
}
assert one * one == one by {
FunctionalEquation(0.0, 0.0);
}
}
}

0 comments on commit 3007adf

Please sign in to comment.