[Kernel] Comments about the invariant computation and other stuffs

- The functions `Nat.of_int`, `Nat.of_strictly_positive_int` and
  `Finite.of_int` now returns an option instead of forcing valid
  bounds on their inputs.

- The spectral exponent search is simpler and returns the maximal
  exponent in the given limit that returns a spectral radius lower
  than one.

- A bug has been fixed in the invariant computation. The sum over
  the inputs missed an element.

src/kernel_services/analysis/filter/finite.ml

@@ -33,8 +33,8 @@ let next : type n. n finite -> n succ finite = fun n -> n + 1
 let ( == ) : type n. n finite -> n finite -> bool = fun l r -> l = r
 let to_int : type n. n finite -> int = fun n -> n
 
-let of_int : type n. n succ nat -> int -> n succ finite =
-  fun limit n -> min (max 0 n) (Nat.to_int limit)
+let of_int : type n. n succ nat -> int -> n succ finite option =
+  fun limit n -> if 0 <= n && n < Nat.to_int limit then Some n else None
 
 let for_each (type n) acc (limit : n nat) (f : n finite -> 'a -> 'a) =
   let acc = ref acc in

src/kernel_services/analysis/filter/finite.mli

@@ -30,10 +30,9 @@ val next : 'n finite -> 'n succ finite
 val ( == ) : 'n finite -> 'n finite -> bool
 
 (* The call [of_int limit n] returns a finite value representing the n-nd
-   element of a finite set of cardinal limit. If n is strictly negative then it
-   returns the first element. If n is greater or equal to the limit then it
-   returns the last element. This function complexity is O(1). *)
-val of_int : 'n succ nat -> int -> 'n succ finite
+   element of a finite set of cardinal limit. If n is not in the bounds, none is
+   returned. This function complexity is O(1). *)
+val of_int : 'n succ nat -> int -> 'n succ finite option
 
 (* The call [to_int n] returns an integer equal to n. This function complexity
    is O(1). *)

src/kernel_services/analysis/filter/linear_filter.ml

@@ -28,46 +28,47 @@ module Make (Field : Field.S) = struct
 
   module Linear = Linear.Space (Field)
   open Linear.Types
+  open Field.Types
+  open Linear
 
 
-  (* Invariant search fuel parameters. TODO: let the user change those. *)
-  let fuel = 50
-  let coarse_refinement_fuel = 5
-  let fine_refinement_fuel = 50
-  let refinement_threshold = 200
 
-  (* Find a first exponant such as the spectral radius is lower than one. *)
-  let rec first_exponant m exp fuel =
-    if Field.(Linear.Matrix.norm m <= one) then (exp, m)
-    else if fuel < 0 then raise (Invalid_argument "Divergeant filter")
-    else first_exponant Linear.Matrix.(m * m) (exp * 2) (fuel - 1)
+  let rec first_steps max steps matrix exponent =
+    let steps = (matrix, exponent) :: steps in
+    if Field.(Matrix.norm matrix < one) then
+      if exponent * 2 > max then Some steps
+      else first_steps max steps Matrix.(matrix * matrix) (exponent * 2)
+    else if exponent <= max then
+      first_steps max steps Matrix.(matrix * matrix) (exponent * 2)
+    else None
 
-  (* Refine the exponant to improve the precision. *)
-  let rec refine base m exp coarse fine =
-    if coarse < 0 || fine < 0 || exp > refinement_threshold
-    then exp, m
-    else if (exp * 2) > refinement_threshold
-    then refine base Linear.Matrix.(base * m) (exp + 1) coarse (fine - 1)
-    else refine base Linear.Matrix.(m * m) (exp * 2) (coarse - 1) fine
+  let rec refine max = function
+    | [] -> None
+    | [ (matrix, exponent) ] -> Some (Matrix.norm matrix, exponent)
+    | (matrix, exponent) :: (matrix', exponent') :: previous ->
+      let exponent'' = exponent + exponent' in
+      if exponent'' > max
+      then refine max ((matrix, exponent) :: previous)
+      else refine max ((Matrix.(matrix * matrix'), exponent'') :: previous)
 
-  (* Looking for an exponant n such as norm(A^n) ≤ 1, which implies that the
-     spectral radius of A is lower than 1. The exponant is refined to improve the
-     invariant's precision. *)
-  let find_exponant base =
-    let exponant, matrix = first_exponant base 1 fuel in
-    let coarse = coarse_refinement_fuel and fine = fine_refinement_fuel in
-    let exponant, matrix = refine base matrix exponant coarse fine in
-    (Linear.Matrix.norm matrix, exponant)
+  let find_exponent max_acceptable_exponent base =
+    let open Option.Operators in
+    let* steps = first_steps max_acceptable_exponent [] base 1 in
+    refine max_acceptable_exponent steps
 
 
 
   module Types = struct
-    type ('n, 'm) filter = Filter : ('n succ, 'm succ) data -> ('n, 'm) filter
+
+    type ('n, 'm) filter =
+      | Filter : ('n succ, 'm succ) data -> ('n succ, 'm succ) filter
+
     and ('n, 'm) data =
       { state : ('n, 'n) matrix
       ; input : ('n, 'm) matrix
       ; measure : 'm vector
       }
+
   end
 
   open Types
@@ -76,24 +77,26 @@ module Make (Field : Field.S) = struct
 
   let create state input measure = Filter { state ; input ; measure }
 
-  let pretty fmt (Filter f) =
-    let open Linear in
+  type ('n, 'm) formatter = Format.formatter -> ('n, 'm) filter -> unit
+  let pretty : type n m. (n, m) formatter = fun fmt (Filter f) ->
     Format.fprintf fmt "Filter:@.@." ;
     Format.fprintf fmt "- State :@.@.  @[<v>%a@]@.@." Matrix.pretty f.state ;
     Format.fprintf fmt "- Input :@.@.  @[<v>%a@]@.@." Matrix.pretty f.input
 
-  let sum order f g stop =
-    let rec compute (acc, res) i =
-      if i > 0 then compute (f acc, Field.(res + g acc)) (i - 1) else res
-    in compute (Linear.Matrix.id order, Field.zero) (stop - 1)
-
-  let invariant (Filter f) =
-    let order, _ = Linear.Matrix.dimensions f.input in
-    let spectral, exponant = find_exponant f.state in
-    let power p = Linear.Matrix.(f.state * p) in
-    let norm  p = Linear.Matrix.(p * f.input |> norm) in
+  let sum order p norm stop =
+    let ( + ) res acc = Field.(res + norm acc) in
+    let rec aux (m, r) i = if i >= 0 then aux (p m, r + m) (i - 1) else r in
+    aux (Matrix.id order, Field.zero) (stop - 1)
+
+  type ('n, 'm) invariant = ('n, 'm) filter -> int -> scalar option
+  let invariant : type n m. (n, m) invariant = fun (Filter f) max ->
+    let open Option.Operators in
+    let order, _ = Matrix.dimensions f.input in
+    let+ spectral, exponant = find_exponent max f.state in
+    let power p = Matrix.(f.state * p) in
+    let norm  p = Matrix.(p * f.input |> norm) in
     let sum = sum order power norm exponant in
-    let bound = Field.(Linear.Vector.norm f.measure * sum / (one - spectral)) in
+    let bound = Field.(Vector.norm f.measure * sum / (one - spectral)) in
     let order = Field.of_int (Nat.to_int order) in
     Field.(bound / order)
 

src/kernel_services/analysis/filter/linear_filter.mli

@@ -28,15 +28,30 @@ module Make (Field : Field.S) : sig
   open Linear.Types
   open Field.Types
 
+  (* A value of type [(n, m) filter] describes a linear filter of order n (i.e
+     with n state variables) and with m inputs. *)
   module Types : sig type ('n, 'm) filter end
   open Types
 
   val create :
     ('n succ, 'n succ) matrix ->
     ('n succ, 'm succ) matrix ->
-    'm succ vector -> ('n, 'm) filter
+    'm succ vector -> ('n succ, 'm succ) filter
 
   val pretty : Format.formatter -> ('n, 'm) filter -> unit
-  val invariant : ('n, 'm) filter -> scalar
+
+  (* Invariant computation. The computation of [invariant filter max] relies on
+     the search of an exponent such as the norm of the state matrix is strictly
+     lower than one. This search depth is bounded by [max]. If no exponent is
+     found before this limit is reached, the function returns None. If an
+     exponent [e] is found, the invariant computation complexity is bounded by
+     O(e * (n^3 + n^2 * m)) with [n] the filter's order and [m] its number of
+     inputs. Only the invariant's upper bound [λ] is returned, the filter's
+     invariant is thus bounded by ±λ. The only thing that limit the optimality
+     of this bound is [max], the initial search depth. However, for most simple
+     filters, a depth of 200 will gives an exact upper bound up to at least ten
+     digits, which is more than enough. Moreover, for those simple filters, the
+     computation is immediate, even when using rational numbers. *)
+  val invariant : ('n, 'm) filter -> int -> scalar option
 
 end

src/kernel_services/analysis/filter/nat.ml

@@ -37,11 +37,7 @@ let prev : type n. n succ nat -> n nat = fun n -> n - 1
 let to_int : type n. n nat -> int = fun n -> n
 
 let of_int n =
-  let next (PositiveOrNull n) = PositiveOrNull (succ n) in
-  let rec aux acc n = if n <= 0 then acc else aux (next acc) (n - 1) in
-  aux (PositiveOrNull zero) n
+  if n < 0 then None else Some (PositiveOrNull n)
 
 let of_strictly_positive_int n =
-  let next (StrictlyPositive n) = StrictlyPositive (succ n) in
-  let rec aux acc n = if n <= 1 then acc else aux (next acc) (n - 1) in
-  aux (StrictlyPositive one) n
+  if n < 1 then None else Some (StrictlyPositive n)

src/kernel_services/analysis/filter/nat.mli

@@ -41,8 +41,8 @@ val to_int : 'n nat -> int
 
 (* Returns a positive or null natural. If the given parameter is stricly
    negative then 0 is returned. This function complexity is O(n). *)
-val of_int : int -> positive_or_null
+val of_int : int -> positive_or_null option
 
 (* Returns a strictly positive natural. If the given parameter is less or equal
    than zero, then 1 is returned. This function complexity is O(n). *)
-val of_strictly_positive_int : int -> strictly_positive
+val of_strictly_positive_int : int -> strictly_positive option

src/kernel_services/analysis/filter/parray.mli → src/libraries/utils/parray.mli

@@ -20,6 +20,10 @@
 (*                                                                        *)
 (**************************************************************************)
 
+(* Persistent array, based on "A Persistent Union-Find Data Structure" by
+   Sylvain Conchon and Jean-Chistophe Filliâtre. For further details, see
+   https://www.lri.fr/~filliatr/ftp/publis/puf-wml07.pdf *)
+
 module Types : sig
   type 'a array
 end