From 4ae36de4def3c07b3a6df21632405ba951312e16 Mon Sep 17 00:00:00 2001 From: lukstafi Date: Tue, 17 Feb 2015 03:16:01 +0100 Subject: [PATCH] A slightly modified variant of the Gauss algorithm example, passes. --- examples/liquid_gauss.gadt | 23 +- examples/liquid_gauss.gadti.target | 293 ++++++---------------- examples/liquid_gauss_harder.gadt | 114 +++++++++ examples/liquid_gauss_harder.gadti.target | 74 ++++++ src/InvarGenTTest.ml | 6 +- 5 files changed, 278 insertions(+), 232 deletions(-) create mode 100644 examples/liquid_gauss_harder.gadt create mode 100644 examples/liquid_gauss_harder.gadti.target diff --git a/examples/liquid_gauss.gadt b/examples/liquid_gauss.gadt index 906f269..26b61f3 100644 --- a/examples/liquid_gauss.gadt +++ b/examples/liquid_gauss.gadt @@ -87,25 +87,24 @@ let rowElim r s n i = let gauss data = let n = matrix_dim1 data in + let m = matrix_dim2 data in - let rec rowMax i = - let x = fabs (matrix_get data i i) in - let rec loop j x mx = - eif j + 1 <= n then - let y = fabs (matrix_get data j i) in - eif (less x y) then loop (j+1) y j - else loop (j+1) x mx - else mx in - loop (i+1) x i in + let rec rowMax i j x mx = + eif j + 1 <= n then + let y = fabs (matrix_get data j i) in + eif (less x y) then rowMax i (j+1) y j + else rowMax i (j+1) x mx + else mx in let rec loop1 i = if i + 1 <= n then - let mx = rowMax i in - norm (getRow data mx) (n+1) i; + let x = fabs (matrix_get data i i) in + let mx = rowMax i (i+1) x i in + norm (getRow data mx) m i; rowSwap data i mx; let rec loop2 j = if j + 1 <= n then ( - rowElim (getRow data i) (getRow data j) (n+1) i; + rowElim (getRow data i) (getRow data j) m i; loop2 (j+1)) else () in loop2 (i+1); diff --git a/examples/liquid_gauss.gadti.target b/examples/liquid_gauss.gadti.target index 3f81408..44fd43a 100644 --- a/examples/liquid_gauss.gadti.target +++ b/examples/liquid_gauss.gadti.target @@ -1,219 +1,74 @@ -(* -(* Here is a test *) - -let main () = - let - val data = alloc(3, alloc(4, 0.0)) - - val _ = update2 (data, 0, 0, 1.0) - val _ = update2 (data, 0, 1, 1.0) - val _ = update2 (data, 0, 2, 1.0) - val _ = update2 (data, 0, 3, 0.0 -. 4.0) - - val _ = update (data, 1, alloc (4, 0.0)) - val _ = update2 (data, 1, 0, 0.0 -. 2.0) - val _ = update2 (data, 1, 1, 3.0) - val _ = update2 (data, 1, 2, 1.0) - val _ = update2 (data, 1, 3, 7.0) - - val _ = update (data, 2, alloc (4, 0.0)) - val _ = update2 (data, 2, 0, 3.0) - val _ = update2 (data, 2, 1, 0.0 -. 1.0) - val _ = update2 (data, 2, 2, 2.0) - val _ = update2 (data, 2, 3, 7.0) - - val mat = Matrix (3, 4, data) - - val _ = print_matrix (mat) - val _ = gauss (mat) - val _ = print_matrix (mat) - in - () - end -withtype unit -> unit -*) -(* -(* - * An implementation of the Gaussian elimination approach in DML - *) - -fun fabs (f) = - if (f >. 0.0) then f else (0.0 -. f) -withtype <> => float -> float - -fun('a) rowSwap (data, i, j) = - let - val temp = sub (data, i) - val _ = update (data, i, sub (data, j)) - in - update (data, j, temp) - end -withtype {m:nat,n:nat,i:nat,j:nat | i < m, j < m} <> => - ('a array(n)) array(m) * int(i) * int(j) -> unit - -fun norm (r, n, i) = - let - val x = sub (r, i) - val _ = update (r, i, 1.0) - fun loop (k) = - if k < n then - let - val _ = update (r, k, sub (r, k) /. x) - in - loop (k+1) - end - else () - withtype {k:nat | k <= n} => int(k) -> unit - in - loop (i+1) - end -withtype {n:nat, i:nat | i < n} <> => - float array(n) * int(n) * int(i) -> unit - -fun rowElim (r, s, n, i) = - let - val x = sub (s, i) - val _ = update (s, i, 0.0) - fun loop (k) = - if k < n then - let - val _ = update (s, k, sub (s, k) -. x *. sub (r, k)) - in - loop (k+1) - end - else () - withtype {k:nat | k <= n} => int(k) -> unit - in - loop (i+1) - end -withtype {n:nat, i:nat | i < n} <> => - float array(n) * float array(n) * int(n) * int(i) -> unit - -fun rowMax (data, m, i) = - let - val x = fabs (sub2 (data, i, i)) - fun loop (j, x, max) = - if j < m then - let - val y = fabs (sub2 (data, j, i)) - in - if (y >. x) then loop (j+1, y, j) - else loop (j+1, x, max) - end - else max - withtype {j:nat, max:nat | max < j <= m} => - int(j) * float * int(max) -> [a:nat | a < m] int(a) - in - loop (i+1, x, i) - end -withtype {m:nat,n:nat,i:nat | i < m, i < n} <> => - (float array(n)) array(m) * int(m) * int(i) -> [a:nat | a < m] int(a) - -datatype 'a matrix with (nat, nat) = - {m:nat, n:nat} - Matrix(m, n) of int(m) * int(n) * ('a array(n)) array(m) - -fun gauss (mat) = - let - val Matrix (n, _, data) = mat - fun loop1 (i) = - if i < n then - let - val max = rowMax(data, n, i) - val _ = norm (sub (data, max), n+1, i) - val _ = rowSwap (data, i, max) - fun loop2 (j) = - if j < n then - let - val _ = rowElim (sub (data, i), - sub (data, j), n+1, i) - in - loop2 (j+1) - end - else () - withtype {j:nat | j <= n} => int(j) -> unit - val _ = loop2 (i+1) - in - loop1 (i+1) - end - else () - withtype {i:nat | i <= n} => int(i) -> unit - in - loop1 (0) - end -withtype {n:pos} <> => float matrix(n,n+1) -> unit - -fun print_array (data, i, j) = - let - fun loop (k) = - if k < j then - let - val _ = print_string ("\t") - val _ = print_float (sub (data, k)) - in - loop (k+1) - end - else print_string "\n" - withtype {k:int | i < k <= j} => int(k) -> unit - in - if i < j then - let - val _ = print_float (sub (data, i)) - in - loop (i+1) - end - else print_string "\n" - end -withtype {n:nat, i:int, j:int | 0 <= i <= j <= n} => - float array(n) * int(i) * int(j) -> unit - -fun print_matrix (mat) = - let - val Matrix (m, n, data) = mat - fun loop (i) = - if i < m then - let - val _ = print_array (sub (data, i), 0, n) - in - loop (i+1) - end - else print_string ("\n") - withtype {i:nat | i <= m} => int(i) -> unit - in - loop (0) - end -withtype {m:nat,n:nat} float matrix(m,n) -> unit - -(* Here is a test *) - -fun main () = - let - val data = alloc(3, alloc(4, 0.0)) - - val _ = update2 (data, 0, 0, 1.0) - val _ = update2 (data, 0, 1, 1.0) - val _ = update2 (data, 0, 2, 1.0) - val _ = update2 (data, 0, 3, 0.0 -. 4.0) - - val _ = update (data, 1, alloc (4, 0.0)) - val _ = update2 (data, 1, 0, 0.0 -. 2.0) - val _ = update2 (data, 1, 1, 3.0) - val _ = update2 (data, 1, 2, 1.0) - val _ = update2 (data, 1, 3, 7.0) - - val _ = update (data, 2, alloc (4, 0.0)) - val _ = update2 (data, 2, 0, 3.0) - val _ = update2 (data, 2, 1, 0.0 -. 1.0) - val _ = update2 (data, 2, 2, 2.0) - val _ = update2 (data, 2, 3, 7.0) - - val mat = Matrix (3, 4, data) - - val _ = print_matrix (mat) - val _ = gauss (mat) - val _ = print_matrix (mat) - in - () - end -withtype unit -> unit -*) \ No newline at end of file +external type Matrix : num * num + +external val matrix_make : + ∀n, k[0 ≤ n ∧ 0 ≤ k]. Num n → Num k → Matrix (n, k) + +external val matrix_get : + ∀n, k, i, j[0 ≤ i ∧ i + 1 ≤ n ∧ 0 ≤ j ∧ j + 1 ≤ k]. + Matrix (n, k) → Num i → Num j → Float + +external val matrix_set : + ∀n, k, i, j[0 ≤ i ∧ i + 1 ≤ n ∧ 0 ≤ j ∧ j + 1 ≤ k]. + Matrix (n, k) → Num i → Num j → Float → () + +external val matrix_dim1 : + ∀n, k[0 ≤ n ∧ 0 ≤ k]. Matrix (n, k) → Num n + +external val matrix_dim2 : + ∀n, k[0 ≤ n ∧ 0 ≤ k]. Matrix (n, k) → Num k + +datatype Array : type * num + +external val array_make : ∀n, a[0 ≤ n]. Num n → a → Array (a, n) + +external val array_get : + ∀n, k, a[0 ≤ k ∧ k + 1 ≤ n]. Array (a, n) → Num k → a + +external val array_set : + ∀n, k, a[0 ≤ k ∧ k + 1 ≤ n]. Array (a, n) → Num k → a → () + +external val array_length : ∀n, a[0 ≤ n]. Array (a, n) → Num n + +external val n2f : ∀n. Num n → Float + +external val equal : ∀a. a → a → Bool + +external val leq : ∀a. a → a → Bool + +external val less : ∀a. a → a → Bool + +external val minus : Float → Float → Float + +external val plus : Float → Float → Float + +external val mult : Float → Float → Float + +external val div : Float → Float → Float + +external val fabs : Float → Float + +external val fl0 : Float + +external val fl1 : Float + +val getRow : + ∀i, k, n[0 ≤ n ∧ 0 ≤ k ∧ k + 1 ≤ i]. + Matrix (i, n) → Num k → Array (Float, n) + +val putRow : + ∀i, j, k, n[0 ≤ n ∧ 0 ≤ k ∧ k + 1 ≤ i ∧ n ≤ j]. + Matrix (i, j) → Num k → Array (Float, n) → () + +val rowSwap : + ∀i, j, k, n[0 ≤ n ∧ 0 ≤ k ∧ k + 1 ≤ i ∧ n + 1 ≤ i ∧ + 0 ≤ j]. Matrix (i, j) → Num k → Num n → () + +val norm : + ∀i, k, n[0 ≤ n ∧ n + 1 ≤ i ∧ k ≤ i]. + Array (Float, i) → Num k → Num n → () + +val rowElim : + ∀i, j, k, n[0 ≤ n ∧ n + 1 ≤ i ∧ k ≤ i ∧ k ≤ j]. + Array (Float, j) → Array (Float, i) → Num k → Num n → () + +val gauss : ∀k, n[0 ≤ n ∧ 1 ≤ k ∧ n ≤ k]. Matrix (n, k) → () diff --git a/examples/liquid_gauss_harder.gadt b/examples/liquid_gauss_harder.gadt new file mode 100644 index 0000000..906f269 --- /dev/null +++ b/examples/liquid_gauss_harder.gadt @@ -0,0 +1,114 @@ +external type Matrix : num * num = + "(float, Bigarray.float64_elt, Bigarray.c_layout) Bigarray.Array2.t" +external let matrix_make : + ∀n, k [0≤n ∧ 0≤k]. Num n → Num k → Matrix (n, k) = + "fun a b -> Bigarray.Array2.create Bigarray.float64 Bigarray.c_layout a b" +external let matrix_get : + ∀n, k, i, j [0≤i ∧ i+1≤n ∧ 0≤j ∧ j+1≤k]. + Matrix (n, k) → Num i → Num j → Float = "Bigarray.Array2.get" +external let matrix_set : + ∀n, k, i, j [0≤i ∧ i+1≤n ∧ 0≤j ∧ j+1≤k]. + Matrix (n, k) → Num i → Num j → Float → () = "Bigarray.Array2.set" +external let matrix_dim1 : + ∀n, k [0≤n ∧ 0≤k]. Matrix (n, k) → Num n = "Bigarray.Array2.dim1" +external let matrix_dim2 : + ∀n, k [0≤n ∧ 0≤k]. Matrix (n, k) → Num k = "Bigarray.Array2.dim2" + +datatype Array : type * num +external let array_make : + ∀n, a [0≤n]. Num n → a → Array (a, n) = "fun a b -> Array.make a b" +external let array_get : + ∀n, k, a [0≤k ∧ k+1≤n]. Array (a, n) → Num k → a = "fun a b -> Array.get a b" +external let array_set : + ∀n, k, a [0≤k ∧ k+1≤n]. Array (a, n) → Num k → a → () = + "fun a b c -> Array.set a b c" +external let array_length : + ∀n, a [0≤n]. Array (a, n) → Num n = "fun a -> Array.length a" + +external let n2f : ∀n. Num n → Float = "float_of_int" +external let equal : ∀a. a → a → Bool = "fun x y -> x = y" +external let leq : ∀a. a → a → Bool = "fun x y -> x <= y" +external let less : ∀a. a → a → Bool = "fun x y -> x < y" + +external let minus : Float → Float → Float = "(-.)" +external let plus : Float → Float → Float = "(+.)" +external let mult : Float → Float → Float = "( *. )" +external let div : Float → Float → Float = "( /. )" +external let fabs : Float → Float = "abs_float" +external let fl0 : Float = "0.0" +external let fl1 : Float = "1.0" + +let getRow data i = + let stride = matrix_dim2 data in + let rowData = array_make stride fl0 in + let rec extract j = + if j + 1 <= stride then ( + array_set rowData j (matrix_get data i j); + (* lukstafi: the call below missing in the original source? *) + extract (j + 1)) + else () in + extract 0; + rowData + +let putRow data i row = + let stride = array_length row in + let rec put j = + if j + 1 <= stride then ( + matrix_set data i j (array_get row j); + (* lukstafi: the call below missing in the original source? *) + put (j + 1)) + else () in + put 0 + +let rowSwap data i j = + let temp = getRow data i in + putRow data i (getRow data j); + putRow data j temp + +let norm r n i = + let x = array_get r i in + array_set r i fl1; + let rec loop k = + if k + 1 <= n then ( + array_set r k (div (array_get r k) x); + loop (k+1)) + else () in + loop (i+1) + +let rowElim r s n i = + let x = array_get s i in + array_set s i fl0; + let rec loop k = + if k + 1 <= n then ( + array_set s k (minus (array_get s k) (mult x (array_get r k))); + loop (k+1)) + else () in + loop (i+1) + +let gauss data = + let n = matrix_dim1 data in + + let rec rowMax i = + let x = fabs (matrix_get data i i) in + let rec loop j x mx = + eif j + 1 <= n then + let y = fabs (matrix_get data j i) in + eif (less x y) then loop (j+1) y j + else loop (j+1) x mx + else mx in + loop (i+1) x i in + + let rec loop1 i = + if i + 1 <= n then + let mx = rowMax i in + norm (getRow data mx) (n+1) i; + rowSwap data i mx; + let rec loop2 j = + if j + 1 <= n then ( + rowElim (getRow data i) (getRow data j) (n+1) i; + loop2 (j+1)) + else () in + loop2 (i+1); + loop1 (i+1) + else () in + loop1 0 diff --git a/examples/liquid_gauss_harder.gadti.target b/examples/liquid_gauss_harder.gadti.target new file mode 100644 index 0000000..020d6c8 --- /dev/null +++ b/examples/liquid_gauss_harder.gadti.target @@ -0,0 +1,74 @@ +external type Matrix : num * num + +external val matrix_make : + ∀n, k[0 ≤ n ∧ 0 ≤ k]. Num n → Num k → Matrix (n, k) + +external val matrix_get : + ∀n, k, i, j[0 ≤ i ∧ i + 1 ≤ n ∧ 0 ≤ j ∧ j + 1 ≤ k]. + Matrix (n, k) → Num i → Num j → Float + +external val matrix_set : + ∀n, k, i, j[0 ≤ i ∧ i + 1 ≤ n ∧ 0 ≤ j ∧ j + 1 ≤ k]. + Matrix (n, k) → Num i → Num j → Float → () + +external val matrix_dim1 : + ∀n, k[0 ≤ n ∧ 0 ≤ k]. Matrix (n, k) → Num n + +external val matrix_dim2 : + ∀n, k[0 ≤ n ∧ 0 ≤ k]. Matrix (n, k) → Num k + +datatype Array : type * num + +external val array_make : ∀n, a[0 ≤ n]. Num n → a → Array (a, n) + +external val array_get : + ∀n, k, a[0 ≤ k ∧ k + 1 ≤ n]. Array (a, n) → Num k → a + +external val array_set : + ∀n, k, a[0 ≤ k ∧ k + 1 ≤ n]. Array (a, n) → Num k → a → () + +external val array_length : ∀n, a[0 ≤ n]. Array (a, n) → Num n + +external val n2f : ∀n. Num n → Float + +external val equal : ∀a. a → a → Bool + +external val leq : ∀a. a → a → Bool + +external val less : ∀a. a → a → Bool + +external val minus : Float → Float → Float + +external val plus : Float → Float → Float + +external val mult : Float → Float → Float + +external val div : Float → Float → Float + +external val fabs : Float → Float + +external val fl0 : Float + +external val fl1 : Float + +val getRow : + ∀i, k, n[0 ≤ n ∧ 0 ≤ k ∧ k + 1 ≤ i]. + Matrix (i, n) → Num k → Array (Float, n) + +val putRow : + ∀i, j, k, n[0 ≤ n ∧ 0 ≤ k ∧ k + 1 ≤ i ∧ n ≤ j]. + Matrix (i, j) → Num k → Array (Float, n) → () + +val rowSwap : + ∀i, j, k, n[0 ≤ n ∧ 0 ≤ k ∧ k + 1 ≤ i ∧ n + 1 ≤ i ∧ + 0 ≤ j]. Matrix (i, j) → Num k → Num n → () + +val norm : + ∀i, k, n[0 ≤ n ∧ n + 1 ≤ i ∧ k ≤ i]. + Array (Float, i) → Num k → Num n → () + +val rowElim : + ∀i, j, k, n[0 ≤ n ∧ n + 1 ≤ i ∧ k ≤ i ∧ k ≤ j]. + Array (Float, j) → Array (Float, i) → Num k → Num n → () + +val gauss : ∀n[0 ≤ n]. Matrix (n, n + 1) → () diff --git a/src/InvarGenTTest.ml b/src/InvarGenTTest.ml index e55e6fc..c3a8e05 100644 --- a/src/InvarGenTTest.ml +++ b/src/InvarGenTTest.ml @@ -640,9 +640,13 @@ let tests = "InvarGenT" >::: [ test_case ~prefer_bound_to_outer:true "liquid_gauss_simpler" ()); "liquid_gauss" >:: (fun () -> - todo "FIXME"; skip_if !debug "debug"; test_case "liquid_gauss" ()); + "liquid_gauss_harder" >:: + (fun () -> + todo "too hard for current numerical abudction"; + skip_if !debug "debug"; + test_case "liquid_gauss_harder" ()); "liquid_fft_simpler" >:: (fun () -> todo "Analysis postponed after InvarGenT 2.0";