examples/space_saving.mlw

   1
   2 (** Space-Saving Algorithm
   3
   4   This is an online algorithm to find out frequent elements in a data
   5   stream. Say we want to detect an element occurring more than N/2 times
   6   in a stream of N elements. We maintain two values (x1 and x2) and two
   7   counters (n1 and n2). If the next value is x1 or x2, we increment the
   8   corresponding counter. Otherwise, we replace the value with the smallest
   9   counter with the next value, *and we increment the corresponding counter*.
  10   If the stream contains a value occurring more than N/2 times, it is
  11   guaranteed to be either x1 or x2.
  12
  13   This generalizes to k values being monitored.
  14
  15   The algorithm is described here:
  16
  17     Metwally, A., Agrawal, D., El Abbadi, A.
  18     Efficient Computation of Frequent and Top-k Elements in Data Streams.
  19     ICDT 2005. LNCS vol 3363.
  20
  21   See also mjrty.mlw for a related algorithm.
  22
  23   Author: Jean-Christophe Filliâtre (CNRS)
  24 *)
  25
  26 use int.Int
  27 use int.MinMax
  28 use map.Occ
  29
  30 (** The elements of the stream. The only thing we can do is to test
  31     elements for equality. *)
  32 type elt
  33
  34 val (=) (x y: elt) : bool
  35   ensures { result <-> x = y }
  36
  37 (** We introduce a dummy value in order to initialize variables in the code. *)
  38 val constant dummy: elt
  39
  40 (** We model an online algorithm with a stream `s` of elements
  41     and a function `next` to get the next element from the stream.
  42     The reference `n` is the number of elements retrieved so far. *)
  43 val ghost s: int -> elt
  44   ensures { forall i. result i <> dummy }
  45
  46 val ghost ref n: int
  47
  48 val next () : elt
  49   requires { n >= 0 }
  50   writes   { n }
  51   ensures  { n = old n + 1 }
  52   ensures  { result = s (old n) <> dummy }
  53
  54 (** Let us start gently with k=2 values monitored. *)
  55
  56 let space_saving_2 () : unit
  57   requires { n = 0 }
  58   diverges
  59 = let ref x1 = dummy in
  60   let ref n1 = 0 in
  61   let ref x2 = dummy in
  62   let ref n2 = 0 in
  63   while true do
  64     invariant { n1 + n2 = n >= 0 }
  65     invariant { 0 <= occ x1 s 0 n <= n1 }
  66     invariant { 0 <= occ x2 s 0 n <= n2 }
  67     invariant { if n1 = 0 then x1 = dummy else x1 <> dummy }
  68     invariant { if n2 = 0 then x2 = dummy else x2 <> dummy }
  69     invariant { n1 > 0 -> n2 > 0 -> x1 <> x2 }
  70     invariant { forall y. y <> x1 -> y <> x2 -> occ y s 0 n <= min n1 n2 }
  71     (* 1. We show that the desired property is a consequence of the
  72        invariants above: any frequent value v (i.e. occurring strictly
  73        more than min(n1, n2) times) is either x1 or x2. *)
  74     assert { [@expl:thm] forall v. occ v s 0 n > min n1 n2 -> v = x1 || v = x2 };
  75     (* and beside, we have min(n1, n2) <= n/2 (from the first invariant)
  76        so any value occurring >n/2 times is either x1 or x2. *)
  77     assert { [@expl:thm] forall v. 2 * occ v s 0 n > n -> v = x1 || v = x2 };
  78     (* 2. Read the next value and update the state. *)
  79     let x = next () in
  80     if x = x1 then n1 <- n1 + 1
  81     else if x = x2 then n2 <- n2 + 1
  82     else if n1 <= n2 then (x1 <- x; n1 <- n1 + 1)
  83                      else (x2 <- x; n2 <- n2 + 1)
  84   done
  85
  86 (** Note: for k=2 (i.e. two values monitored), this is less precise
  87     than MJRTY (see mjrty.mlw). Indeed, Space-Saving only tells us
  88     that a value with more than N/2 occurrences, if any, is either x1
  89     or x2, while MJRTY determines what would be this value. *)
  90
  91 (** Now we go for the general case of `k` values being monitored,
  92     for any k >= 2. *)
  93
  94 use array.Array
  95 use array.ArraySum
  96
  97 (** The index for the minimum value of a non-empty array. *)
  98 let function minimum (a: array int) : (m: int)
  99   requires { length a > 0 }
 100   ensures  { 0 <= m < length a }
 101   ensures  { forall i. 0 <= i < length a -> a[m] <= a[i] }
 102 = let ref m = 0 in
 103   for i = 1 to length a - 1 do
 104     invariant { 0 <= m < length a }
 105     invariant { forall j. 0 <= j < i -> a[m] <= a[j] }
 106     if a[i] < a[m] then m <- i
 107   done;
 108   return m
 109
 110 predicate occurs (v: elt) (a: array elt) =
 111   exists i. 0 <= i < length a /\ v = a[i]
 112
 113 (** It is a bit of a pity that we have to split sums like this to help
 114     SMT solvers... *)
 115 let increment (ghost k: int) (c: array int) (i: int) (ghost n: int) : unit
 116   requires { 0 <= i < length c = k }
 117   requires { sum c 0 k = n - 1 }
 118   ensures  { c[i] = old c[i] + 1 }
 119   ensures  { forall j. 0 <= j < k -> j <> i -> c[j] = old c[j] }
 120   ensures  { sum c 0 k = n }
 121 = assert { sum c 0 k = sum c 0 i + sum c i (i+1) + sum c (i+1) k };
 122   c[i] <- c[i] + 1;
 123   assert { sum c 0 k = sum c 0 i + sum c i (i+1) + sum c (i+1) k }
 124
 125 (** Finds x in array e of size k, or returns k if not present. *)
 126 let find (k: int) (x: elt) (e: array elt) : (i: int)
 127   requires { length e = k }
 128   ensures  { 0 <= i <= k }
 129   ensures  { forall j. 0 <= j < i -> e[j] <> x }
 130   ensures  { i < k -> e[i] = x }
 131 = for i = 0 to k-1 do
 132     invariant { forall j. 0 <= j < i -> e[j] <> x }
 133     if e[i] = x then return i
 134   done;
 135   return k
 136
 137 (** Let us help SMT solvers with non-linear arithmetic. *)
 138 let lemma minimum_k (k: int) (c: array int) (n: int)
 139   requires { length c = k >= 2 }
 140   requires { sum c 0 k = n >= 0 }
 141   ensures  { k * c[minimum c] <= n }
 142 = let m = minimum c in
 143   for i = 0 to k - 1 do invariant { i * c[m] <= sum c 0 i } () done
 144
 145 (** Space-Saving Algorithm with `k` values being monitored. *)
 146
 147 let space_saving_k (k: int) : unit
 148   requires { k >= 2 }
 149   requires { n = 0 }
 150   diverges
 151 = let ref e = Array.make k dummy in
 152   let ref c = Array.make k 0 in
 153   while true do
 154     invariant { sum c 0 k = n >= 0 }
 155     invariant { forall i. 0 <= i < k -> 0 <= occ e[i] s 0 n <= c[i] }
 156     invariant { forall i. 0 <= i < k ->
 157                   if c[i] = 0 then e[i] = dummy else e[i] <> dummy }
 158     invariant { forall i. 0 <= i < k -> c[i] > 0 ->
 159                 forall j. 0 <= j < k -> c[j] > 0 -> i <> j -> e[i] <> e[j] }
 160     invariant { forall y. (forall i. 0 <= i < k -> y <> e[i]) ->
 161                           occ y s 0 n <= c[minimum c] }
 162     (* 1. We show that the desired property is a consequence of the
 163        invariants above: any frequent value `v` (i.e. occurring strictly
 164        more than min(c) times) is in `e`. *)
 165     assert { [@expl:thm] forall v. occ v s 0 n > c[minimum c] -> occurs v e };
 166     (* and beside, we have min(c) <= n/k (from the first invariant) *)
 167     minimum_k k c n;
 168     (* so any value occurring >n/k times is in `e`. *)
 169     assert { [@expl:thm] forall v. k * occ v s 0 n > n -> occurs v e
 170                                 by k * occ v s 0 n > k * c[minimum c] };
 171     (* 2. Read the next value and update the state. *)
 172     let x = next () in
 173     let i = find k x e in
 174     if i < k then
 175       increment k c i n
 176     else (
 177       let m = minimum c in
 178       e[m] <- x;
 179       increment k c m n;
 180     )
 181   done