From 06b9500dde883283cc12fb794475a926a06cd867 Mon Sep 17 00:00:00 2001
From: Maxime Jacquemin <maxime.jacquemin@cea.fr>
Date: Tue, 23 Jan 2024 16:01:03 +0100
Subject: [PATCH] [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 | 4 +-
.../analysis/filter/finite.mli | 7 +-
.../analysis/filter/linear_filter.ml | 81 ++++++++++---------
.../analysis/filter/linear_filter.mli | 19 ++++-
src/kernel_services/analysis/filter/nat.ml | 8 +-
src/kernel_services/analysis/filter/nat.mli | 4 +-
.../filter => libraries/utils}/parray.ml | 0
.../filter => libraries/utils}/parray.mli | 4 +
8 files changed, 72 insertions(+), 55 deletions(-)
rename src/{kernel_services/analysis/filter => libraries/utils}/parray.ml (100%)
rename src/{kernel_services/analysis/filter => libraries/utils}/parray.mli (90%)
diff --git a/src/kernel_services/analysis/filter/finite.ml b/src/kernel_services/analysis/filter/finite.ml
index f992d181c71..a3851bc9961 100644
--- a/src/kernel_services/analysis/filter/finite.ml
+++ b/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
diff --git a/src/kernel_services/analysis/filter/finite.mli b/src/kernel_services/analysis/filter/finite.mli
index 1f2a6b986a7..36d4f1ebd55 100644
--- a/src/kernel_services/analysis/filter/finite.mli
+++ b/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). *)
diff --git a/src/kernel_services/analysis/filter/linear_filter.ml b/src/kernel_services/analysis/filter/linear_filter.ml
index 20192fe0398..977c7043350 100644
--- a/src/kernel_services/analysis/filter/linear_filter.ml
+++ b/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)
diff --git a/src/kernel_services/analysis/filter/linear_filter.mli b/src/kernel_services/analysis/filter/linear_filter.mli
index 9eb759f98ed..c7e2b8f25b9 100644
--- a/src/kernel_services/analysis/filter/linear_filter.mli
+++ b/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
diff --git a/src/kernel_services/analysis/filter/nat.ml b/src/kernel_services/analysis/filter/nat.ml
index c85a25a5a47..3f9705c9aa3 100644
--- a/src/kernel_services/analysis/filter/nat.ml
+++ b/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)
diff --git a/src/kernel_services/analysis/filter/nat.mli b/src/kernel_services/analysis/filter/nat.mli
index 9942dfbee1a..d3d6cf83fea 100644
--- a/src/kernel_services/analysis/filter/nat.mli
+++ b/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
diff --git a/src/kernel_services/analysis/filter/parray.ml b/src/libraries/utils/parray.ml
similarity index 100%
rename from src/kernel_services/analysis/filter/parray.ml
rename to src/libraries/utils/parray.ml
diff --git a/src/kernel_services/analysis/filter/parray.mli b/src/libraries/utils/parray.mli
similarity index 90%
rename from src/kernel_services/analysis/filter/parray.mli
rename to src/libraries/utils/parray.mli
index 96df4db916f..4d28c548a5d 100644
--- a/src/kernel_services/analysis/filter/parray.mli
+++ b/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
--
GitLab