diff --git a/src/kernel_services/analysis/filter/field.ml b/src/kernel_services/analysis/filter/field.ml
new file mode 100644
index 0000000000000000000000000000000000000000..d7b53d09c64e66c0cae29340445818ce99432278
--- /dev/null
+++ b/src/kernel_services/analysis/filter/field.ml
@@ -0,0 +1,52 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+module type S = sig
+
+ type scalar
+ include Datatype.S_with_collections with type t = scalar
+
+ val zero : scalar
+ val one : scalar
+ val infinity : scalar
+
+ val of_int : int -> scalar
+ val of_float : float -> scalar
+ val to_float : scalar -> float
+
+ val neg : scalar -> scalar
+ val abs : scalar -> scalar
+ val max : scalar -> scalar -> scalar
+ val min : scalar -> scalar -> scalar
+
+ val ( = ) : scalar -> scalar -> bool
+ val ( <= ) : scalar -> scalar -> bool
+ val ( < ) : scalar -> scalar -> bool
+ val ( >= ) : scalar -> scalar -> bool
+ val ( > ) : scalar -> scalar -> bool
+
+ val ( + ) : scalar -> scalar -> scalar
+ val ( - ) : scalar -> scalar -> scalar
+ val ( * ) : scalar -> scalar -> scalar
+ val ( / ) : scalar -> scalar -> scalar
+
+end
diff --git a/src/kernel_services/analysis/filter/finite.ml b/src/kernel_services/analysis/filter/finite.ml
new file mode 100644
index 0000000000000000000000000000000000000000..2799fe8141f66bcd6ac703aaf260ba1d92d74245
--- /dev/null
+++ b/src/kernel_services/analysis/filter/finite.ml
@@ -0,0 +1,39 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+open Nat
+
+type 'n finite = int
+
+let first : type n. n succ finite = 0
+let last : type n. n succ nat -> n succ finite = fun n -> Nat.to_int n - 1
+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 option =
+ fun limit n -> if 0 <= n && n < Nat.to_int limit then Some n else None
+
+let for_each (type n) (f : n finite -> 'a -> 'a) (limit : n nat) acc =
+ let acc = ref acc in
+ for i = 0 to Nat.to_int limit - 1 do acc := f i !acc done ;
+ !acc
diff --git a/src/kernel_services/analysis/filter/finite.mli b/src/kernel_services/analysis/filter/finite.mli
new file mode 100644
index 0000000000000000000000000000000000000000..0e6d08866fb1f1f036437c2deb8a0c079562d2f2
--- /dev/null
+++ b/src/kernel_services/analysis/filter/finite.mli
@@ -0,0 +1,48 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+open Nat
+
+
+
+(* The type [n finite] encodes all finite sets of cardinal [n]. It is used by
+ the module Linear to represent accesses to vectors and matrices coefficients,
+ statically ensuring that no out of bounds access can be performed. *)
+type 'n finite
+
+val first : 'n succ finite
+val last : 'n succ nat -> 'n succ finite
+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 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). *)
+val to_int : 'n finite -> int
+
+(* The call [for_each acc limit f] folds over each finite elements of a set of
+ cardinal limit, computing f at each step. The function complexity is O(n). *)
+val for_each : ('n finite -> 'a -> 'a) -> 'n nat -> 'a -> 'a
diff --git a/src/kernel_services/analysis/filter/linear.ml b/src/kernel_services/analysis/filter/linear.ml
new file mode 100644
index 0000000000000000000000000000000000000000..1dd2fc44fd08b4f711c79ff02d00d1a97e57bd37
--- /dev/null
+++ b/src/kernel_services/analysis/filter/linear.ml
@@ -0,0 +1,166 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+open Nat
+open Finite
+
+
+
+module Space (Field : Field.S) = struct
+
+ type scalar = Field.scalar
+
+ type ('n, 'm) m = { data : scalar Parray.t ; rows : 'n nat ; cols : 'm nat }
+ type ('n, 'm) matrix = M : ('n succ, 'm succ) m -> ('n succ, 'm succ) matrix
+ type 'n vector = ('n, zero succ) matrix
+
+
+
+ type 'n row = Format.formatter -> 'n finite -> unit
+ let pretty (type n m) (row : n row) fmt (M m : (n, m) matrix) =
+ let cut () = Format.pp_print_cut fmt () in
+ let first () = Format.fprintf fmt "@[<h>⌈%a⌉@]" row Finite.first in
+ let mid i = Format.fprintf fmt "@[<h>|%a|@]" row i in
+ let last () = Format.fprintf fmt "@[<h>⌋%a⌊@]" row Finite.(last m.rows) in
+ let row i () =
+ if Finite.(i = first) then first ()
+ else if Finite.(i = last m.rows) then (cut () ; last ())
+ else (cut () ; mid i)
+ in
+ Format.pp_open_vbox fmt 0 ;
+ Finite.for_each row m.rows () ;
+ Format.pp_close_box fmt ()
+
+
+
+ module Vector = struct
+
+ let pretty_row (type n) fmt (M { data ; _ } : n vector) =
+ Format.pp_open_hbox fmt () ;
+ Parray.pretty ~sep:"@ " Field.pretty fmt data ;
+ Format.pp_close_box fmt ()
+
+ let init size f =
+ let data = Parray.init (Nat.to_int size) (fun _ -> Field.zero) in
+ let set i data = Parray.set data (Finite.to_int i) (f i) in
+ let data = Finite.for_each set size data in
+ M { data ; rows = size ; cols = Nat.one }
+
+ let size (type n) (M vector : n vector) : n nat = vector.rows
+ let repeat n size = init size (fun _ -> n)
+ let zero size = repeat Field.zero size
+
+ let get (type n) (i : n finite) (M vec : n vector) : scalar =
+ Parray.get vec.data (Finite.to_int i)
+
+ let pretty (type n) fmt (vector : n vector) =
+ let get fmt (i : n finite) = Field.pretty fmt (get i vector) in
+ pretty get fmt vector
+
+ let set (type n) (i : n finite) scalar (M vec : n vector) : n vector =
+ M { vec with data = Parray.set vec.data (Finite.to_int i) scalar }
+
+ let norm (type n) (v : n vector) : scalar =
+ let max i r = Field.(max (abs (get i v)) r) in
+ Finite.for_each max (size v) Field.zero
+
+ let ( * ) (type n) (l : n vector) (r : n vector) =
+ let inner i acc = Field.(acc + get i l * get i r) in
+ Finite.for_each inner (size l) Field.zero
+
+ let base (type n) (i : n succ finite) (dimension : n succ nat) =
+ zero dimension |> set i Field.one
+
+ end
+
+
+
+ module Matrix = struct
+
+ let index cols i j = i * Nat.to_int cols + j
+
+ let get (type n m) (i : n finite) (j : m finite) (M m : (n, m) matrix) =
+ let i = Finite.to_int i and j = Finite.to_int j in
+ Parray.get m.data (index m.cols i j)
+
+ let set (type n m) i j num (M m : (n, m) matrix) : (n, m) matrix =
+ let i = Finite.to_int i and j = Finite.to_int j in
+ let data = Parray.set m.data (index m.cols i j) num in
+ M { m with data }
+
+ let row row (M m) = Vector.init m.cols @@ fun i -> get row i (M m)
+ let col col (M m) = Vector.init m.rows @@ fun i -> get i col (M m)
+
+ let dimensions : type n m. (n, m) matrix -> n nat * m nat =
+ fun (M m) -> m.rows, m.cols
+
+ let pretty (type n m) fmt (M m : (n, m) matrix) =
+ let row fmt i = Vector.pretty_row fmt (row i (M m)) in
+ pretty row fmt (M m)
+
+ let init n m init =
+ let rows = Nat.to_int n and cols = Nat.to_int m in
+ let t = Parray.init (rows * cols) (fun _ -> Field.zero) in
+ let index i j = index m (Finite.to_int i) (Finite.to_int j) in
+ let set i j data = Parray.set data (index i j) (init i j) in
+ let data = Finite.(for_each (fun i t -> for_each (set i) m t) n t) in
+ M { data ; rows = n ; cols = m }
+
+ let zero n m = init n m (fun _ _ -> Field.zero)
+ let id n = Finite.for_each (fun i m -> set i i Field.one m) n (zero n n)
+
+ let transpose : type n m. (n, m) matrix -> (m, n) matrix =
+ fun (M m) -> init m.cols m.rows (fun j i -> get i j (M m))
+
+ type ('n, 'm) add = ('n, 'm) matrix -> ('n, 'm) matrix -> ('n, 'm) matrix
+ let ( + ) : type n m. (n, m) add = fun (M l) (M r) ->
+ let ( + ) i j = Field.(get i j (M l) + get i j (M r)) in
+ init l.rows l.cols ( + )
+
+ type ('n, 'm, 'p) mul = ('n, 'm) matrix -> ('m, 'p) matrix -> ('n, 'p) matrix
+ let ( * ) : type n m p. (n, m, p) mul = fun (M l) (M r) ->
+ let ( * ) i j = Vector.(row i (M l) * col j (M r)) in
+ init l.rows r.cols ( * )
+
+ let norm : type n m. (n, m) matrix -> scalar = fun (M m) ->
+ let add v j r = Field.(abs (Vector.get j v) + r) in
+ let sum v = Finite.for_each (add v) (Vector.size v) Field.zero in
+ let max i res = Field.max res (row i (M m) |> sum) in
+ Finite.for_each max m.rows Field.zero
+
+ let power (type n) (M m : (n, n) matrix) : int -> (n, n) matrix =
+ let n = dimensions (M m) |> fst in
+ let cache = Datatype.Int.Hashtbl.create 17 in
+ let find i = Datatype.Int.Hashtbl.find_opt cache i in
+ let save i v = Datatype.Int.Hashtbl.add cache i v ; v in
+ let rec pow e =
+ if e < 0 then raise (Invalid_argument "negative exponent") ;
+ match find e with
+ | Some r -> r
+ | None when Stdlib.(e = 0) -> id n
+ | None when Stdlib.(e = 1) -> M m
+ | None -> let h = pow (e / 2) in save e (pow (e mod 2) * h * h)
+ in pow
+
+ end
+
+end
diff --git a/src/kernel_services/analysis/filter/linear.mli b/src/kernel_services/analysis/filter/linear.mli
new file mode 100644
index 0000000000000000000000000000000000000000..843b2f969f6b3c8fee5aba3108aad7bcd5e45034
--- /dev/null
+++ b/src/kernel_services/analysis/filter/linear.mli
@@ -0,0 +1,122 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+open Nat
+open Finite
+
+
+
+(* Definition of a linear space over a field. Used by Linear_filter to represent
+ and compute linear filters invariants. *)
+
+module Space (Field : Field.S) : sig
+
+ (* The type of scalars in the field ð•‚. *)
+ type scalar = Field.scalar
+
+ (* The type of matrices in ð•‚â¿Ë£áµ *)
+ type ('n, 'm) matrix
+
+ (* Type representing a column vector of ð•‚â¿. One can use [transpose] if one
+ needs a row vector, for example when printing it. *)
+ type 'n vector = ('n, zero succ) matrix
+
+
+
+ module Vector : sig
+
+ val pretty : Format.formatter -> 'n vector -> unit
+
+ (* The call [zero n] returns the 0 vector in ð•‚â¿. *)
+ val zero : 'n succ nat -> 'n succ vector
+
+ (* The call [repeat x n] returns a vector in ð•‚â¿ which each dimension
+ containing the scalar x. *)
+ val repeat : scalar -> 'n succ nat -> 'n succ vector
+
+ (* The call [base i n] returns the i-th base vector in the orthonormal space
+ of ð•‚â¿. In other words, the returned vector contains zero except for the
+ i-th coordinate, which contains one. *)
+ val base : 'n succ finite -> 'n succ nat -> 'n succ vector
+
+ (* The call [set i x v] returns a new vector of the same linear space as
+ [v], and with the same coordinates, except for the i-th one, which is
+ set to the scalar [x]. *)
+ val set : 'n finite -> scalar -> 'n vector -> 'n vector
+
+ (* The call [size v] for [v] a vector of ð•‚â¿ returns n. *)
+ val size : 'n vector -> 'n nat
+
+ (* The call [norm v] computes the ∞-norm of [v], i.e the maximum of the
+ absolute values of its coordinates. *)
+ val norm : 'n vector -> scalar
+
+ end
+
+
+
+ module Matrix : sig
+
+ val pretty : Format.formatter -> ('n, 'm) matrix -> unit
+
+ (* The call [id n] returns the identity matrix in ð•‚â¿Ë£â¿. *)
+ val id : 'n succ nat -> ('n succ, 'n succ) matrix
+
+ (* The call [zero n m] returns the 0 matrix in ð•‚â¿Ë£áµ. *)
+ val zero : 'n succ nat -> 'm succ nat -> ('n succ, 'm succ) matrix
+
+ (* The call [get i j m] returns the coefficient of the i-th row and
+ the j-th column. *)
+ val get : 'n finite -> 'm finite -> ('n, 'm) matrix -> scalar
+
+ (* The call [set i j x m] returns a new matrix of the same linear space as
+ [m], and with the same coefficients, except for the one of the i-th row
+ and the j-th column, which is set to the scalar [x]. *)
+ val set : 'n finite -> 'm finite -> scalar -> ('n, 'm) matrix -> ('n, 'm) matrix
+
+ (* The call [norm m] computes the ∞-norm of [m], i.e the maximum of the
+ absolute sums of the rows of [m]. *)
+ val norm : ('n, 'm) matrix -> scalar
+
+ (* The call [transpose m] for m in ð•‚â¿Ë£áµ returns a new matrix in ð•‚áµË£â¿ with
+ all the coefficients transposed. *)
+ val transpose : ('n, 'm) matrix -> ('m, 'n) matrix
+
+ (* The call [dimensions m] for m in ð•‚â¿Ë£áµ returns the pair (n, m). *)
+ val dimensions : ('m, 'n) matrix -> 'm nat * 'n nat
+
+ (* Matrices addition. The dimensions compatibility is statically ensured. *)
+ val ( + ) : ('n, 'm) matrix -> ('n, 'm) matrix -> ('n, 'm) matrix
+
+ (* Matrices multiplication. The dimensions compatibility is statically
+ ensured. *)
+ val ( * ) : ('n, 'm) matrix -> ('m, 'p) matrix -> ('n, 'p) matrix
+
+ (* Matrix exponentiation. The call [power m] returns a memoized function.
+ When one needs to compute several exponentiations of the same matrix, one
+ should perform the call [power m] once and used the returned function
+ each times one needs it. *)
+ val power : ('n, 'n) matrix -> (int -> ('n, 'n) matrix)
+
+ end
+
+end
diff --git a/src/kernel_services/analysis/filter/linear_filter.ml b/src/kernel_services/analysis/filter/linear_filter.ml
new file mode 100644
index 0000000000000000000000000000000000000000..d63f8837d6e955a467fbccaabc57e0c42fe8d9c8
--- /dev/null
+++ b/src/kernel_services/analysis/filter/linear_filter.ml
@@ -0,0 +1,98 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+
+
+module Make (Field : Field.S) = struct
+
+ module Linear = Linear.Space (Field)
+ open Linear
+ open Nat
+
+
+
+ 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
+ ; center : 'n vector
+ ; measure : 'm vector
+ }
+
+
+
+ let create ~state ~input ~center ~measure =
+ Filter { state ; input ; center ; measure }
+
+ 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 "@[<v>" ;
+ Format.fprintf fmt "Filter:@ @ " ;
+ Format.fprintf fmt "- State :@ @ @[<v>%a@]@ @ " Matrix.pretty f.state ;
+ Format.fprintf fmt "- Input :@ @ @[<v>%a@]@ @ " Matrix.pretty f.input ;
+ Format.fprintf fmt "@]"
+
+
+
+ type 'n invariant = ('n, zero succ succ) matrix
+
+ let invariant rows = Matrix.zero rows Nat.(zero |> succ |> succ)
+ let lower i inv = Matrix.get i Finite.first inv
+ let upper i inv = Matrix.get i Finite.(next first) inv
+ let bounds i inv = (lower i inv, upper i inv)
+ let set_lower i bound inv = Matrix.set i Finite.first bound inv
+ let set_upper i bound inv = Matrix.set i Finite.(next first) bound inv
+
+
+
+ let check_convergence matrix =
+ let norm = Matrix.norm matrix in
+ if Field.(norm < one) then Some norm else None
+
+ type ('n, 'm) compute = ('n, 'm) filter -> int -> 'n invariant option
+ let invariant : type n m. (n, m) compute = fun (Filter f) e ->
+ let open Option.Operators in
+ let state = Matrix.power f.state in
+ let measure = Vector.norm f.measure in
+ let order, _ = Matrix.dimensions f.input in
+ let base i = Vector.base i order |> Matrix.transpose in
+ let* StrictlyPositive exponent = Nat.of_strictly_positive_int e in
+ let+ spectral = state (Nat.to_int exponent) |> check_convergence in
+ (* Computation of the inputs contribution for the i-th state dimension *)
+ let input base e = Matrix.(base * state e * f.input |> norm) in
+ let add_input base e res = Field.(res + input base (Finite.to_int e)) in
+ let input i = Finite.for_each (base i |> add_input) exponent Field.zero in
+ (* Computation of the center contribution for the i-th state dimension *)
+ let center e = Matrix.(state e * f.center) in
+ let add_center e res = Matrix.(res + center (Finite.to_int e)) in
+ let center = Finite.for_each add_center exponent (Vector.zero order) in
+ let center i = Matrix.(base i * center |> norm) in
+ (* Bounds computation for each state dimension *)
+ let numerator sign i = Field.(center i + sign * measure * input i) in
+ let bound sign i = Field.(numerator sign i / (one - spectral)) in
+ let lower i inv = set_lower i (bound Field.(neg one) i) inv in
+ let upper i inv = set_upper i (bound Field.one i) inv in
+ Finite.(invariant order |> for_each lower order |> for_each upper order)
+
+end
diff --git a/src/kernel_services/analysis/filter/linear_filter.mli b/src/kernel_services/analysis/filter/linear_filter.mli
new file mode 100644
index 0000000000000000000000000000000000000000..352c681ada85a88ada7861f43dc93fdb57ea658f
--- /dev/null
+++ b/src/kernel_services/analysis/filter/linear_filter.mli
@@ -0,0 +1,91 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+open Nat
+open Finite
+
+(* A filter corresponds to the recursive equation
+ X[k + 1] = AX[k] + Bε[k + 1] + C where :
+ - n is the filter's order and m its number of inputs ;
+ - X[k] ∈ â„^n is the filter's state at iteration [k] ;
+ - ε[k] ∈ â„^m is the filters's inputs at iteration [k] ;
+ - A ∈ â„^(n×n) is the filter's state matrix ;
+ - B ∈ â„^(n×m) is the filter's input matrix ;
+ - C ∈ â„^n is the filter's center.
+
+ The goal of this module is to compute filters invariants, i.e bounds for
+ each of the filter's state dimensions when the iterations goes to infinity.
+ To do so, it only suppose that, at each iteration, each input εi is bounded
+ by [-λi .. λi]. Each input is thus supposed centered around zero but each
+ one can have different bounds.
+
+ {!Linear_filter_test} is an example using this module. *)
+module Make (Field : Field.S) : sig
+
+ module Linear : module type of Linear.Space (Field)
+
+ (* A value of type [(n, m) filter] describes a linear filter of order n (i.e
+ with n state variables) and with m inputs. *)
+ type ('n, 'm) filter
+
+ (* Create a filter's representation. The inputs are as following :
+ - state is the filter's state matrix ;
+ - input is the filter's input matrix ;
+ - center is the filter's center ;
+ - measure is a vector representing upper bounds for the filter's inputs. *)
+ val create :
+ state : ('n succ, 'n succ) Linear.matrix ->
+ input : ('n succ, 'm succ) Linear.matrix ->
+ center : 'n succ Linear.vector ->
+ measure : 'm succ Linear.vector ->
+ ('n succ, 'm succ) filter
+
+ val pretty : Format.formatter -> ('n, 'm) filter -> unit
+
+ (* Representation of a filter's invariant. Bounds for each dimension can be
+ accessed using the corresponding functions. *)
+ type 'n invariant
+ val lower : 'n finite -> 'n invariant -> Field.scalar
+ val upper : 'n finite -> 'n invariant -> Field.scalar
+ val bounds : 'n finite -> 'n invariant -> Field.scalar * Field.scalar
+
+ (* Invariant computation. The computation of [invariant filter k] relies on
+ the search of an exponent such as the norm of the state matrix is strictly
+ lower than one. For the filter to converge, there must exist an α such as,
+ for every β greater than α, ||A^β|| < 1 with A the filter's state matrix.
+ As such, the search does not have to find α, but instead any exponent such
+ as the property is satisfied. As the computed invariant will be more
+ precise with a larger exponent, the computation always uses [k], the
+ largest authorized exponent, and thus only check that indeed ||A^k|| < 1.
+ If the property is not verified, the function returns None as it cannot
+ prove that the filter converges.
+
+ The only thing limiting the invariant optimality is [k]. However, for most
+ simple filters, k = 200 will gives exact bounds 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. Indeed, the
+ invariant computation complexity is bounded by O(kn^3 + mn^2) with [n]
+ the filter's order and [m] its number of inputs. It is thus linear in
+ the targeted exponent. *)
+ val invariant : ('n, 'm) filter -> int -> 'n invariant option
+
+end
diff --git a/src/kernel_services/analysis/filter/linear_filter_test.ml b/src/kernel_services/analysis/filter/linear_filter_test.ml
new file mode 100644
index 0000000000000000000000000000000000000000..f239d3adbeefad5443c3761cb5adb613c2516253
--- /dev/null
+++ b/src/kernel_services/analysis/filter/linear_filter_test.ml
@@ -0,0 +1,147 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+module Rational = struct
+
+ type scalar = Q.t
+
+ module Type = struct
+ include Datatype.Serializable_undefined
+ type t = scalar
+ let name = "Linear.Filter.Test.Rational"
+ let reprs = [ Q.zero ]
+ let compare = Q.compare
+ let equal = Q.equal
+ let hash q = Z.hash (Q.num q) + 11 * Z.hash (Q.den q)
+ end
+
+ include Datatype.Make_with_collections (Type)
+
+ let pretty fmt n =
+ let ten = Z.of_int 10 in
+ let sign = if Q.sign n >= 0 then "" else "-" in
+ let num = Q.num n |> Z.abs and den = Q.den n |> Z.abs in
+ let finish n = Z.Compare.(n >= den || n = Z.zero) in
+ let rec f e n = if finish n then (n, e) else f (e + 1) Z.(n * ten) in
+ let num, exponent = f 0 num in
+ let default fmt n = Format.fprintf fmt "%1.7f" (Q.to_float n) in
+ if exponent > 0 then
+ let number = Q.make num den in
+ Format.fprintf fmt "%s%aE-%d" sign default number exponent
+ else default fmt n
+
+ include Q
+ let infinity = Q.inf
+ let ( = ) = Q.equal
+
+end
+
+
+
+module Filter = Linear_filter.Make (Rational)
+module Linear = Filter.Linear
+
+let max_exponent = 200
+let fin size n = Finite.of_int size n |> Option.get
+let set row col i j n = Linear.Matrix.set (fin row i) (fin col j) n
+
+let pretty_bounds invariant fmt i =
+ let l, u = Filter.bounds i invariant in
+ Format.fprintf fmt "@[<h>[%a .. %a]@]" Rational.pretty l Rational.pretty u
+
+let pretty_invariant order fmt = function
+ | None -> Format.fprintf fmt "%s" (Unicode.top_string ())
+ | Some invariant ->
+ let pp f i = pretty_bounds invariant f i in
+ let pp f i = Format.fprintf f "@[<h>* %d : %a@]@," (Finite.to_int i) pp i in
+ let pretty fmt () = Finite.for_each (fun i () -> pp fmt i) order () in
+ Format.fprintf fmt "@[<v>%a@]" pretty ()
+
+
+
+(* Invariant computation for the filter:
+ X = 0.68 * X - 0.68 * Y + E1;
+ Y = 0.68 * X + 0.68 * Y + E2;
+ with E1 ∈ [-1 .. 1] and E2 ∈ [-1 .. 1]. *)
+module Circle = struct
+
+ let order = Nat.(succ one)
+
+ let state =
+ Linear.Matrix.zero order order
+ |> set order order 0 0 Rational.(of_float 0.68)
+ |> set order order 0 1 Rational.(of_float ~-.0.68)
+ |> set order order 1 0 Rational.(of_float 0.68)
+ |> set order order 1 1 Rational.(of_float 0.68)
+
+ let input = Linear.Matrix.id order
+
+ let measure = Linear.Vector.repeat Rational.one order
+
+ let center = Linear.Vector.zero order
+
+ let compute () =
+ let filter = Filter.create ~state ~input ~center ~measure in
+ let invariant = Filter.invariant filter max_exponent in
+ Kernel.result "@[<v>Circle :@,%a@,@]" (pretty_invariant order) invariant
+
+end
+
+
+
+(* Invariant computation for the filter:
+ X = 1.5 * X - 0.7 * Y + E + 1;
+ Y = X + 1;
+ with E ∈ [-0.1 .. 0.1]. *)
+module Simple = struct
+
+ let order = Nat.(succ one)
+ let delay = Nat.one
+
+ let state =
+ Linear.Matrix.zero order order
+ |> set order order 0 0 Rational.(of_float 1.5)
+ |> set order order 0 1 Rational.(of_float ~-.0.7)
+ |> set order order 1 0 Rational.(of_float 1.)
+ |> set order order 1 1 Rational.(of_float 0.)
+
+ let input =
+ Linear.Matrix.zero order delay
+ |> set order delay 0 0 Rational.one
+ |> set order delay 1 0 Rational.zero
+
+ let measure = Linear.Vector.repeat (Rational.of_float 0.1) delay
+
+ let center = Linear.Vector.repeat Rational.one order
+
+ let compute () =
+ let filter = Filter.create ~state ~input ~center ~measure in
+ let invariant = Filter.invariant filter max_exponent in
+ Kernel.result "@[<v>Simple :@,%a@,@]" (pretty_invariant order) invariant
+
+end
+
+
+
+let run () =
+ Circle.compute () ;
+ Simple.compute ()
diff --git a/src/kernel_services/analysis/filter/linear_filter_test.mli b/src/kernel_services/analysis/filter/linear_filter_test.mli
new file mode 100644
index 0000000000000000000000000000000000000000..e973f651b1aef8f89f6fc9c1e2a110da2801685f
--- /dev/null
+++ b/src/kernel_services/analysis/filter/linear_filter_test.mli
@@ -0,0 +1,24 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+(** Tests the Linear_filter module: invariant computation of filters. *)
+val run : unit -> unit
diff --git a/src/kernel_services/analysis/filter/nat.ml b/src/kernel_services/analysis/filter/nat.ml
new file mode 100644
index 0000000000000000000000000000000000000000..28e197482e85b7fad514d2dcec7c6f425506cee1
--- /dev/null
+++ b/src/kernel_services/analysis/filter/nat.ml
@@ -0,0 +1,39 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+type zero = |
+type 'n succ = |
+type 'n nat = int
+type positive_or_null = PositiveOrNull : 'n nat -> positive_or_null
+type strictly_positive = StrictlyPositive : 'n succ nat -> strictly_positive
+
+let zero : zero nat = 0
+let one : zero succ nat = 1
+let succ : type n. n nat -> n succ nat = fun n -> n + 1
+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 =
+ if n < 0 then None else Some (PositiveOrNull n)
+
+let of_strictly_positive_int 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
new file mode 100644
index 0000000000000000000000000000000000000000..2bc2534b7908e7c65568c5fcd4b4090061e04c7b
--- /dev/null
+++ b/src/kernel_services/analysis/filter/nat.mli
@@ -0,0 +1,48 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+(* Encoding of the Peano arithmetic in OCaml type system. A value of type
+ [n nat] contains [n] at the value and the type level, allowing to express
+ properties on objects sizes and accesses for example. It is used by the
+ module Linear to represent vectors and matrices dimensions. *)
+type zero = |
+type 'n succ = |
+type 'n nat
+type positive_or_null = PositiveOrNull : 'n nat -> positive_or_null
+type strictly_positive = StrictlyPositive : 'n succ nat -> strictly_positive
+
+val zero : zero nat
+val one : zero succ nat
+val succ : 'n nat -> 'n succ nat
+val prev : 'n succ nat -> 'n nat
+
+(* The call [to_int n] returns an integer equal to n. This function complexity
+ is O(1). *)
+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 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 option
diff --git a/src/libraries/utils/parray.ml b/src/libraries/utils/parray.ml
new file mode 100644
index 0000000000000000000000000000000000000000..fd022290d664dfb8d8758c0063a60941b75ed2a8
--- /dev/null
+++ b/src/libraries/utils/parray.ml
@@ -0,0 +1,70 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+type 'a t = 'a data ref
+and 'a data = Memory of 'a array | Diff of int * 'a * 'a t
+
+let init n f = ref (Memory (Array.init n f))
+
+let rec reroot t =
+ match !t with
+ | Memory mem -> mem
+ | Diff (i, v, t') ->
+ let mem = reroot t' in
+ let v' = Array.get mem i in
+ Array.set mem i v ;
+ t := Memory mem ;
+ t' := Diff (i, v', t) ;
+ mem
+
+let pretty ?(sep = format_of_string " ") pp fmt t =
+ let mem = reroot t in
+ for i = 0 to Array.length mem - 1 do
+ if i > 0 then Format.fprintf fmt sep ;
+ pp fmt mem.(i)
+ done
+
+let get t i =
+ match !t with
+ | Memory mem -> Array.get mem i
+ | Diff _ -> Array.get (reroot t) i
+
+let set t i v =
+ let mem = reroot t in
+ let old = Array.get mem i in
+ Array.set mem i v ;
+ let res = ref (Memory mem) in
+ t := Diff (i, old, res) ;
+ res
+
+let fold f t acc =
+ let mem = reroot t and index = ref 0 in
+ let f' acc v = let acc = f !index v acc in incr index ; acc in
+ Array.fold_left f' acc mem
+
+let map f t =
+ let mem = reroot t |> Array.copy in
+ let res = ref t in
+ for i = 0 to Array.length mem - 1 do
+ res := set !res i (f mem.(i))
+ done ;
+ !res
diff --git a/src/libraries/utils/parray.mli b/src/libraries/utils/parray.mli
new file mode 100644
index 0000000000000000000000000000000000000000..b7eb12ecc0c03550494978d9b130e1d52977cb21
--- /dev/null
+++ b/src/libraries/utils/parray.mli
@@ -0,0 +1,38 @@
+(**************************************************************************)
+(* *)
+(* This file is part of Frama-C. *)
+(* *)
+(* Copyright (C) 2007-2023 *)
+(* CEA (Commissariat à l'énergie atomique et aux énergies *)
+(* alternatives) *)
+(* *)
+(* you can redistribute it and/or modify it under the terms of the GNU *)
+(* Lesser General Public License as published by the Free Software *)
+(* Foundation, version 2.1. *)
+(* *)
+(* It is distributed in the hope that it will be useful, *)
+(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
+(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *)
+(* GNU Lesser General Public License for more details. *)
+(* *)
+(* See the GNU Lesser General Public License version 2.1 *)
+(* for more details (enclosed in the file licenses/LGPLv2.1). *)
+(* *)
+(**************************************************************************)
+
+(* 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 *)
+
+type 'a t
+
+val init : int -> (int -> 'a) -> 'a t
+val get : 'a t -> int -> 'a
+val set : 'a t -> int -> 'a -> 'a t
+val fold : (int -> 'a -> 'b -> 'b) -> 'a t -> 'b -> 'b
+val map : ('a -> 'a) -> 'a t -> 'a t
+
+val pretty :
+ ?sep : Pretty_utils.sformat ->
+ (Format.formatter -> 'a -> unit) ->
+ Format.formatter -> 'a t -> unit
diff --git a/tests/float/linear_filter_test.i b/tests/float/linear_filter_test.i
new file mode 100644
index 0000000000000000000000000000000000000000..7cb38611d3a0c2ae466b20e6ee7bf17f9b0836d9
--- /dev/null
+++ b/tests/float/linear_filter_test.i
@@ -0,0 +1,7 @@
+/* run.config
+ MODULE: @PTEST_NAME@
+ OPT:
+*/
+/* run.config*
+ DONTRUN:
+*/
diff --git a/tests/float/linear_filter_test.ml b/tests/float/linear_filter_test.ml
new file mode 100644
index 0000000000000000000000000000000000000000..420df28ae1ef75ecd621cacb90d67cdc5b4cf87b
--- /dev/null
+++ b/tests/float/linear_filter_test.ml
@@ -0,0 +1,4 @@
+(* Programmatic tests of the invariant computation of linear filters by the
+ module Linear_filter. Run by linear_filter_test.i. *)
+
+let () = Boot.Main.extend Linear_filter_test.run
diff --git a/tests/float/oracle/linear_filter_test.res.oracle b/tests/float/oracle/linear_filter_test.res.oracle
new file mode 100644
index 0000000000000000000000000000000000000000..f69acd94d653045be3b4fb16c72e12aea07e4b1b
--- /dev/null
+++ b/tests/float/oracle/linear_filter_test.res.oracle
@@ -0,0 +1,7 @@
+[kernel] Parsing linear_filter_test.i (no preprocessing)
+[kernel] Circle :
+ * 0 : [-31.3829787 .. 31.3829787]
+ * 1 : [-31.3829787 .. 31.3829787]
+[kernel] Simple :
+ * 0 : [5.7267253E-1 .. 2.4273275]
+ * 1 : [1.5726725 .. 3.4273275]