src/HOL/Library/sum_of_squares.ML
 author chaieb Tue Jun 09 11:10:33 2009 +0200 (2009-06-09) changeset 31512 27118561c2e0 parent 31131 d9752181691a child 32150 4ed2865f3a56 permissions -rw-r--r--
Tuned sos tactic to reject non SOS goals
 chaieb@31119 ` 1` ```structure Sos = ``` chaieb@31119 ` 2` ```struct ``` chaieb@31119 ` 3` chaieb@31119 ` 4` ```val rat_0 = Rat.zero; ``` chaieb@31119 ` 5` ```val rat_1 = Rat.one; ``` chaieb@31119 ` 6` ```val rat_2 = Rat.two; ``` chaieb@31119 ` 7` ```val rat_10 = Rat.rat_of_int 10; ``` chaieb@31119 ` 8` ```val rat_1_2 = rat_1 // rat_2; ``` chaieb@31119 ` 9` ```val max = curry IntInf.max; ``` chaieb@31119 ` 10` ```val min = curry IntInf.min; ``` chaieb@31119 ` 11` chaieb@31119 ` 12` ```val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int; ``` chaieb@31119 ` 13` ```val numerator_rat = Rat.quotient_of_rat #> fst #> Rat.rat_of_int; ``` chaieb@31119 ` 14` ```fun int_of_rat a = ``` chaieb@31119 ` 15` ``` case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int"; ``` chaieb@31119 ` 16` ```fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y)); ``` chaieb@31119 ` 17` chaieb@31119 ` 18` ```fun rat_pow r i = ``` chaieb@31119 ` 19` ``` let fun pow r i = ``` chaieb@31119 ` 20` ``` if i = 0 then rat_1 else ``` chaieb@31119 ` 21` ``` let val d = pow r (i div 2) ``` chaieb@31119 ` 22` ``` in d */ d */ (if i mod 2 = 0 then rat_1 else r) ``` chaieb@31119 ` 23` ``` end ``` chaieb@31119 ` 24` ``` in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end; ``` chaieb@31119 ` 25` chaieb@31119 ` 26` ```fun round_rat r = ``` chaieb@31119 ` 27` ``` let val (a,b) = Rat.quotient_of_rat (Rat.abs r) ``` chaieb@31119 ` 28` ``` val d = a div b ``` chaieb@31119 ` 29` ``` val s = if r = b then d + 1 else d) end ``` chaieb@31119 ` 32` chaieb@31119 ` 33` ```val abs_rat = Rat.abs; ``` chaieb@31119 ` 34` ```val pow2 = rat_pow rat_2; ``` chaieb@31119 ` 35` ```val pow10 = rat_pow rat_10; ``` chaieb@31119 ` 36` chaieb@31119 ` 37` ```val debugging = ref false; ``` chaieb@31119 ` 38` chaieb@31119 ` 39` ```exception Sanity; ``` chaieb@31119 ` 40` chaieb@31119 ` 41` ```exception Unsolvable; ``` chaieb@31119 ` 42` chaieb@31119 ` 43` ```(* Turn a rational into a decimal string with d sig digits. *) ``` chaieb@31119 ` 44` chaieb@31119 ` 45` ```local ``` chaieb@31119 ` 46` ```fun normalize y = ``` chaieb@31119 ` 47` ``` if abs_rat y =/ rat_1 then normalize (y // rat_10) + 1 ``` chaieb@31119 ` 49` ``` else 0 ``` chaieb@31119 ` 50` ``` in ``` chaieb@31119 ` 51` ```fun decimalize d x = ``` chaieb@31119 ` 52` ``` if x =/ rat_0 then "0.0" else ``` chaieb@31119 ` 53` ``` let ``` chaieb@31119 ` 54` ``` val y = Rat.abs x ``` chaieb@31119 ` 55` ``` val e = normalize y ``` chaieb@31119 ` 56` ``` val z = pow10(~ e) */ y +/ rat_1 ``` chaieb@31119 ` 57` ``` val k = int_of_rat (round_rat(pow10 d */ z)) ``` chaieb@31119 ` 58` ``` in (if x a ``` chaieb@31119 ` 69` ``` | h::t => itern (k + 1) t f (f h k a); ``` chaieb@31119 ` 70` chaieb@31119 ` 71` ```fun iter (m,n) f a = ``` chaieb@31119 ` 72` ``` if n < m then a ``` chaieb@31119 ` 73` ``` else iter (m+1,n) f (f m a); ``` chaieb@31119 ` 74` chaieb@31119 ` 75` ```(* The main types. *) ``` chaieb@31119 ` 76` chaieb@31119 ` 77` ```fun strict_ord ord (x,y) = case ord (x,y) of LESS => LESS | _ => GREATER ``` chaieb@31119 ` 78` chaieb@31119 ` 79` ```structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord); ``` chaieb@31119 ` 80` chaieb@31119 ` 81` ```type vector = int* Rat.rat Intfunc.T; ``` chaieb@31119 ` 82` chaieb@31119 ` 83` ```type matrix = (int*int)*(Rat.rat Intpairfunc.T); ``` chaieb@31119 ` 84` chaieb@31119 ` 85` ```type monomial = int Ctermfunc.T; ``` chaieb@31119 ` 86` chaieb@31119 ` 87` ```val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t)) ``` chaieb@31119 ` 88` ``` fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2) ``` chaieb@31119 ` 89` ```structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord) ``` chaieb@31119 ` 90` chaieb@31119 ` 91` ```type poly = Rat.rat Monomialfunc.T; ``` chaieb@31119 ` 92` chaieb@31119 ` 93` ``` fun iszero (k,r) = r =/ rat_0; ``` chaieb@31119 ` 94` chaieb@31119 ` 95` ```fun fold_rev2 f l1 l2 b = ``` chaieb@31119 ` 96` ``` case (l1,l2) of ``` chaieb@31119 ` 97` ``` ([],[]) => b ``` chaieb@31119 ` 98` ``` | (h1::t1,h2::t2) => f h1 h2 (fold_rev2 f t1 t2 b) ``` chaieb@31119 ` 99` ``` | _ => error "fold_rev2"; ``` chaieb@31119 ` 100` ``` ``` chaieb@31119 ` 101` ```(* Vectors. Conventionally indexed 1..n. *) ``` chaieb@31119 ` 102` chaieb@31119 ` 103` ```fun vector_0 n = (n,Intfunc.undefined):vector; ``` chaieb@31119 ` 104` chaieb@31119 ` 105` ```fun dim (v:vector) = fst v; ``` chaieb@31119 ` 106` chaieb@31119 ` 107` ```fun vector_const c n = ``` chaieb@31119 ` 108` ``` if c =/ rat_0 then vector_0 n ``` chaieb@31119 ` 109` ``` else (n,fold_rev (fn k => Intfunc.update (k,c)) (1 upto n) Intfunc.undefined) :vector; ``` chaieb@31119 ` 110` chaieb@31119 ` 111` ```val vector_1 = vector_const rat_1; ``` chaieb@31119 ` 112` chaieb@31119 ` 113` ```fun vector_cmul c (v:vector) = ``` chaieb@31119 ` 114` ``` let val n = dim v ``` chaieb@31119 ` 115` ``` in if c =/ rat_0 then vector_0 n ``` chaieb@31119 ` 116` ``` else (n,Intfunc.mapf (fn x => c */ x) (snd v)) ``` chaieb@31119 ` 117` ``` end; ``` chaieb@31119 ` 118` chaieb@31119 ` 119` ```fun vector_neg (v:vector) = (fst v,Intfunc.mapf Rat.neg (snd v)) :vector; ``` chaieb@31119 ` 120` chaieb@31119 ` 121` ```fun vector_add (v1:vector) (v2:vector) = ``` chaieb@31119 ` 122` ``` let val m = dim v1 ``` chaieb@31119 ` 123` ``` val n = dim v2 ``` chaieb@31119 ` 124` ``` in if m <> n then error "vector_add: incompatible dimensions" ``` chaieb@31119 ` 125` ``` else (n,Intfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd v1) (snd v2)) :vector ``` chaieb@31119 ` 126` ``` end; ``` chaieb@31119 ` 127` chaieb@31119 ` 128` ```fun vector_sub v1 v2 = vector_add v1 (vector_neg v2); ``` chaieb@31119 ` 129` chaieb@31119 ` 130` ```fun vector_dot (v1:vector) (v2:vector) = ``` chaieb@31119 ` 131` ``` let val m = dim v1 ``` chaieb@31119 ` 132` ``` val n = dim v2 ``` chaieb@31119 ` 133` ``` in if m <> n then error "vector_dot: incompatible dimensions" ``` chaieb@31119 ` 134` ``` else Intfunc.fold (fn (i,x) => fn a => x +/ a) ``` chaieb@31119 ` 135` ``` (Intfunc.combine (curry op */) (fn x => x =/ rat_0) (snd v1) (snd v2)) rat_0 ``` chaieb@31119 ` 136` ``` end; ``` chaieb@31119 ` 137` chaieb@31119 ` 138` ```fun vector_of_list l = ``` chaieb@31119 ` 139` ``` let val n = length l ``` chaieb@31119 ` 140` ``` in (n,fold_rev2 (curry Intfunc.update) (1 upto n) l Intfunc.undefined) :vector ``` chaieb@31119 ` 141` ``` end; ``` chaieb@31119 ` 142` chaieb@31119 ` 143` ```(* Matrices; again rows and columns indexed from 1. *) ``` chaieb@31119 ` 144` chaieb@31119 ` 145` ```fun matrix_0 (m,n) = ((m,n),Intpairfunc.undefined):matrix; ``` chaieb@31119 ` 146` chaieb@31119 ` 147` ```fun dimensions (m:matrix) = fst m; ``` chaieb@31119 ` 148` chaieb@31119 ` 149` ```fun matrix_const c (mn as (m,n)) = ``` chaieb@31119 ` 150` ``` if m <> n then error "matrix_const: needs to be square" ``` chaieb@31119 ` 151` ``` else if c =/ rat_0 then matrix_0 mn ``` chaieb@31119 ` 152` ``` else (mn,fold_rev (fn k => Intpairfunc.update ((k,k), c)) (1 upto n) Intpairfunc.undefined) :matrix;; ``` chaieb@31119 ` 153` chaieb@31119 ` 154` ```val matrix_1 = matrix_const rat_1; ``` chaieb@31119 ` 155` chaieb@31119 ` 156` ```fun matrix_cmul c (m:matrix) = ``` chaieb@31119 ` 157` ``` let val (i,j) = dimensions m ``` chaieb@31119 ` 158` ``` in if c =/ rat_0 then matrix_0 (i,j) ``` chaieb@31119 ` 159` ``` else ((i,j),Intpairfunc.mapf (fn x => c */ x) (snd m)) ``` chaieb@31119 ` 160` ``` end; ``` chaieb@31119 ` 161` chaieb@31119 ` 162` ```fun matrix_neg (m:matrix) = ``` chaieb@31119 ` 163` ``` (dimensions m, Intpairfunc.mapf Rat.neg (snd m)) :matrix; ``` chaieb@31119 ` 164` chaieb@31119 ` 165` ```fun matrix_add (m1:matrix) (m2:matrix) = ``` chaieb@31119 ` 166` ``` let val d1 = dimensions m1 ``` chaieb@31119 ` 167` ``` val d2 = dimensions m2 ``` chaieb@31119 ` 168` ``` in if d1 <> d2 ``` chaieb@31119 ` 169` ``` then error "matrix_add: incompatible dimensions" ``` chaieb@31119 ` 170` ``` else (d1,Intpairfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd m1) (snd m2)) :matrix ``` chaieb@31119 ` 171` ``` end;; ``` chaieb@31119 ` 172` chaieb@31119 ` 173` ```fun matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2); ``` chaieb@31119 ` 174` chaieb@31119 ` 175` ```fun row k (m:matrix) = ``` chaieb@31119 ` 176` ``` let val (i,j) = dimensions m ``` chaieb@31119 ` 177` ``` in (j, ``` chaieb@31119 ` 178` ``` Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then Intfunc.update (j,c) a else a) (snd m) Intfunc.undefined ) : vector ``` chaieb@31119 ` 179` ``` end; ``` chaieb@31119 ` 180` chaieb@31119 ` 181` ```fun column k (m:matrix) = ``` chaieb@31119 ` 182` ``` let val (i,j) = dimensions m ``` chaieb@31119 ` 183` ``` in (i, ``` chaieb@31119 ` 184` ``` Intpairfunc.fold (fn ((i,j), c) => fn a => if j = k then Intfunc.update (i,c) a else a) (snd m) Intfunc.undefined) ``` chaieb@31119 ` 185` ``` : vector ``` chaieb@31119 ` 186` ``` end; ``` chaieb@31119 ` 187` chaieb@31119 ` 188` ```fun transp (m:matrix) = ``` chaieb@31119 ` 189` ``` let val (i,j) = dimensions m ``` chaieb@31119 ` 190` ``` in ``` chaieb@31119 ` 191` ``` ((j,i),Intpairfunc.fold (fn ((i,j), c) => fn a => Intpairfunc.update ((j,i), c) a) (snd m) Intpairfunc.undefined) :matrix ``` chaieb@31119 ` 192` ``` end; ``` chaieb@31119 ` 193` chaieb@31119 ` 194` ```fun diagonal (v:vector) = ``` chaieb@31119 ` 195` ``` let val n = dim v ``` chaieb@31119 ` 196` ``` in ((n,n),Intfunc.fold (fn (i, c) => fn a => Intpairfunc.update ((i,i), c) a) (snd v) Intpairfunc.undefined) : matrix ``` chaieb@31119 ` 197` ``` end; ``` chaieb@31119 ` 198` chaieb@31119 ` 199` ```fun matrix_of_list l = ``` chaieb@31119 ` 200` ``` let val m = length l ``` chaieb@31119 ` 201` ``` in if m = 0 then matrix_0 (0,0) else ``` chaieb@31119 ` 202` ``` let val n = length (hd l) ``` chaieb@31119 ` 203` ``` in ((m,n),itern 1 l (fn v => fn i => itern 1 v (fn c => fn j => Intpairfunc.update ((i,j), c))) Intpairfunc.undefined) ``` chaieb@31119 ` 204` ``` end ``` chaieb@31119 ` 205` ``` end; ``` chaieb@31119 ` 206` chaieb@31119 ` 207` ```(* Monomials. *) ``` chaieb@31119 ` 208` chaieb@31119 ` 209` ```fun monomial_eval assig (m:monomial) = ``` chaieb@31119 ` 210` ``` Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (Ctermfunc.apply assig x) k) ``` chaieb@31119 ` 211` ``` m rat_1; ``` chaieb@31119 ` 212` ```val monomial_1 = (Ctermfunc.undefined:monomial); ``` chaieb@31119 ` 213` chaieb@31119 ` 214` ```fun monomial_var x = Ctermfunc.onefunc (x, 1) :monomial; ``` chaieb@31119 ` 215` chaieb@31119 ` 216` ```val (monomial_mul:monomial->monomial->monomial) = ``` chaieb@31119 ` 217` ``` Ctermfunc.combine (curry op +) (K false); ``` chaieb@31119 ` 218` chaieb@31119 ` 219` ```fun monomial_pow (m:monomial) k = ``` chaieb@31119 ` 220` ``` if k = 0 then monomial_1 ``` chaieb@31119 ` 221` ``` else Ctermfunc.mapf (fn x => k * x) m; ``` chaieb@31119 ` 222` chaieb@31119 ` 223` ```fun monomial_divides (m1:monomial) (m2:monomial) = ``` chaieb@31119 ` 224` ``` Ctermfunc.fold (fn (x, k) => fn a => Ctermfunc.tryapplyd m2 x 0 >= k andalso a) m1 true;; ``` chaieb@31119 ` 225` chaieb@31119 ` 226` ```fun monomial_div (m1:monomial) (m2:monomial) = ``` chaieb@31119 ` 227` ``` let val m = Ctermfunc.combine (curry op +) ``` chaieb@31119 ` 228` ``` (fn x => x = 0) m1 (Ctermfunc.mapf (fn x => ~ x) m2) ``` chaieb@31119 ` 229` ``` in if Ctermfunc.fold (fn (x, k) => fn a => k >= 0 andalso a) m true then m ``` chaieb@31119 ` 230` ``` else error "monomial_div: non-divisible" ``` chaieb@31119 ` 231` ``` end; ``` chaieb@31119 ` 232` chaieb@31119 ` 233` ```fun monomial_degree x (m:monomial) = ``` chaieb@31119 ` 234` ``` Ctermfunc.tryapplyd m x 0;; ``` chaieb@31119 ` 235` chaieb@31119 ` 236` ```fun monomial_lcm (m1:monomial) (m2:monomial) = ``` chaieb@31119 ` 237` ``` fold_rev (fn x => Ctermfunc.update (x, max (monomial_degree x m1) (monomial_degree x m2))) ``` chaieb@31119 ` 238` ``` (gen_union (is_equal o cterm_ord) (Ctermfunc.dom m1, Ctermfunc.dom m2)) (Ctermfunc.undefined :monomial); ``` chaieb@31119 ` 239` chaieb@31119 ` 240` ```fun monomial_multidegree (m:monomial) = ``` chaieb@31119 ` 241` ``` Ctermfunc.fold (fn (x, k) => fn a => k + a) m 0;; ``` chaieb@31119 ` 242` chaieb@31119 ` 243` ```fun monomial_variables m = Ctermfunc.dom m;; ``` chaieb@31119 ` 244` chaieb@31119 ` 245` ```(* Polynomials. *) ``` chaieb@31119 ` 246` chaieb@31119 ` 247` ```fun eval assig (p:poly) = ``` chaieb@31119 ` 248` ``` Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0; ``` chaieb@31119 ` 249` chaieb@31119 ` 250` ```val poly_0 = (Monomialfunc.undefined:poly); ``` chaieb@31119 ` 251` chaieb@31119 ` 252` ```fun poly_isconst (p:poly) = ``` chaieb@31119 ` 253` ``` Monomialfunc.fold (fn (m, c) => fn a => Ctermfunc.is_undefined m andalso a) p true; ``` chaieb@31119 ` 254` chaieb@31119 ` 255` ```fun poly_var x = Monomialfunc.onefunc (monomial_var x,rat_1) :poly; ``` chaieb@31119 ` 256` chaieb@31119 ` 257` ```fun poly_const c = ``` chaieb@31119 ` 258` ``` if c =/ rat_0 then poly_0 else Monomialfunc.onefunc(monomial_1, c); ``` chaieb@31119 ` 259` chaieb@31119 ` 260` ```fun poly_cmul c (p:poly) = ``` chaieb@31119 ` 261` ``` if c =/ rat_0 then poly_0 ``` chaieb@31119 ` 262` ``` else Monomialfunc.mapf (fn x => c */ x) p; ``` chaieb@31119 ` 263` chaieb@31119 ` 264` ```fun poly_neg (p:poly) = (Monomialfunc.mapf Rat.neg p :poly);; ``` chaieb@31119 ` 265` chaieb@31119 ` 266` ```fun poly_add (p1:poly) (p2:poly) = ``` chaieb@31119 ` 267` ``` (Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2 :poly); ``` chaieb@31119 ` 268` chaieb@31119 ` 269` ```fun poly_sub p1 p2 = poly_add p1 (poly_neg p2); ``` chaieb@31119 ` 270` chaieb@31119 ` 271` ```fun poly_cmmul (c,m) (p:poly) = ``` chaieb@31119 ` 272` ``` if c =/ rat_0 then poly_0 ``` chaieb@31119 ` 273` ``` else if Ctermfunc.is_undefined m ``` chaieb@31119 ` 274` ``` then Monomialfunc.mapf (fn d => c */ d) p ``` chaieb@31119 ` 275` ``` else Monomialfunc.fold (fn (m', d) => fn a => (Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0; ``` chaieb@31119 ` 276` chaieb@31119 ` 277` ```fun poly_mul (p1:poly) (p2:poly) = ``` chaieb@31119 ` 278` ``` Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0; ``` chaieb@31119 ` 279` chaieb@31119 ` 280` ```fun poly_div (p1:poly) (p2:poly) = ``` chaieb@31119 ` 281` ``` if not(poly_isconst p2) ``` chaieb@31119 ` 282` ``` then error "poly_div: non-constant" else ``` chaieb@31119 ` 283` ``` let val c = eval Ctermfunc.undefined p2 ``` chaieb@31119 ` 284` ``` in if c =/ rat_0 then error "poly_div: division by zero" ``` chaieb@31119 ` 285` ``` else poly_cmul (Rat.inv c) p1 ``` chaieb@31119 ` 286` ``` end; ``` chaieb@31119 ` 287` chaieb@31119 ` 288` ```fun poly_square p = poly_mul p p; ``` chaieb@31119 ` 289` chaieb@31119 ` 290` ```fun poly_pow p k = ``` chaieb@31119 ` 291` ``` if k = 0 then poly_const rat_1 ``` chaieb@31119 ` 292` ``` else if k = 1 then p ``` chaieb@31119 ` 293` ``` else let val q = poly_square(poly_pow p (k div 2)) in ``` chaieb@31119 ` 294` ``` if k mod 2 = 1 then poly_mul p q else q end; ``` chaieb@31119 ` 295` chaieb@31119 ` 296` ```fun poly_exp p1 p2 = ``` chaieb@31119 ` 297` ``` if not(poly_isconst p2) ``` chaieb@31119 ` 298` ``` then error "poly_exp: not a constant" ``` chaieb@31119 ` 299` ``` else poly_pow p1 (int_of_rat (eval Ctermfunc.undefined p2)); ``` chaieb@31119 ` 300` chaieb@31119 ` 301` ```fun degree x (p:poly) = ``` chaieb@31119 ` 302` ``` Monomialfunc.fold (fn (m,c) => fn a => max (monomial_degree x m) a) p 0; ``` chaieb@31119 ` 303` chaieb@31119 ` 304` ```fun multidegree (p:poly) = ``` chaieb@31119 ` 305` ``` Monomialfunc.fold (fn (m, c) => fn a => max (monomial_multidegree m) a) p 0; ``` chaieb@31119 ` 306` chaieb@31119 ` 307` ```fun poly_variables (p:poly) = ``` chaieb@31119 ` 308` ``` sort cterm_ord (Monomialfunc.fold_rev (fn (m, c) => curry (gen_union (is_equal o cterm_ord)) (monomial_variables m)) p []);; ``` chaieb@31119 ` 309` chaieb@31119 ` 310` ```(* Order monomials for human presentation. *) ``` chaieb@31119 ` 311` chaieb@31119 ` 312` ```fun cterm_ord (t,t') = TermOrd.fast_term_ord (term_of t, term_of t'); ``` chaieb@31119 ` 313` chaieb@31119 ` 314` ```val humanorder_varpow = prod_ord cterm_ord (rev_order o int_ord); ``` chaieb@31119 ` 315` chaieb@31119 ` 316` ```local ``` chaieb@31119 ` 317` ``` fun ord (l1,l2) = case (l1,l2) of ``` chaieb@31119 ` 318` ``` (_,[]) => LESS ``` chaieb@31119 ` 319` ``` | ([],_) => GREATER ``` chaieb@31119 ` 320` ``` | (h1::t1,h2::t2) => ``` chaieb@31119 ` 321` ``` (case humanorder_varpow (h1, h2) of ``` chaieb@31119 ` 322` ``` LESS => LESS ``` chaieb@31119 ` 323` ``` | EQUAL => ord (t1,t2) ``` chaieb@31119 ` 324` ``` | GREATER => GREATER) ``` chaieb@31119 ` 325` ```in fun humanorder_monomial m1 m2 = ``` chaieb@31119 ` 326` ``` ord (sort humanorder_varpow (Ctermfunc.graph m1), ``` chaieb@31119 ` 327` ``` sort humanorder_varpow (Ctermfunc.graph m2)) ``` chaieb@31119 ` 328` ```end; ``` chaieb@31119 ` 329` chaieb@31119 ` 330` ```fun fold1 f l = case l of ``` chaieb@31119 ` 331` ``` [] => error "fold1" ``` chaieb@31119 ` 332` ``` | [x] => x ``` chaieb@31119 ` 333` ``` | (h::t) => f h (fold1 f t); ``` chaieb@31119 ` 334` chaieb@31119 ` 335` ```(* Conversions to strings. *) ``` chaieb@31119 ` 336` chaieb@31119 ` 337` ```fun string_of_vector min_size max_size (v:vector) = ``` chaieb@31119 ` 338` ``` let val n_raw = dim v ``` chaieb@31119 ` 339` ``` in if n_raw = 0 then "[]" else ``` chaieb@31119 ` 340` ``` let ``` chaieb@31119 ` 341` ``` val n = max min_size (min n_raw max_size) ``` chaieb@31119 ` 342` ``` val xs = map (Rat.string_of_rat o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) ``` chaieb@31119 ` 343` ``` in "[" ^ fold1 (fn s => fn t => s ^ ", " ^ t) xs ^ ``` chaieb@31119 ` 344` ``` (if n_raw > max_size then ", ...]" else "]") ``` chaieb@31119 ` 345` ``` end ``` chaieb@31119 ` 346` ``` end; ``` chaieb@31119 ` 347` chaieb@31119 ` 348` ```fun string_of_matrix max_size (m:matrix) = ``` chaieb@31119 ` 349` ``` let ``` chaieb@31119 ` 350` ``` val (i_raw,j_raw) = dimensions m ``` chaieb@31119 ` 351` ``` val i = min max_size i_raw ``` chaieb@31119 ` 352` ``` val j = min max_size j_raw ``` chaieb@31119 ` 353` ``` val rstr = map (fn k => string_of_vector j j (row k m)) (1 upto i) ``` chaieb@31119 ` 354` ``` in "["^ fold1 (fn s => fn t => s^";\n "^t) rstr ^ ``` chaieb@31119 ` 355` ``` (if j > max_size then "\n ...]" else "]") ``` chaieb@31119 ` 356` ``` end; ``` chaieb@31119 ` 357` chaieb@31119 ` 358` ```fun string_of_term t = ``` chaieb@31119 ` 359` ``` case t of ``` chaieb@31119 ` 360` ``` a\$b => "("^(string_of_term a)^" "^(string_of_term b)^")" ``` chaieb@31119 ` 361` ``` | Abs x => ``` chaieb@31119 ` 362` ``` let val (xn, b) = Term.dest_abs x ``` chaieb@31119 ` 363` ``` in "(\\"^xn^"."^(string_of_term b)^")" ``` chaieb@31119 ` 364` ``` end ``` chaieb@31119 ` 365` ``` | Const(s,_) => s ``` chaieb@31119 ` 366` ``` | Free (s,_) => s ``` chaieb@31119 ` 367` ``` | Var((s,_),_) => s ``` chaieb@31119 ` 368` ``` | _ => error "string_of_term"; ``` chaieb@31119 ` 369` chaieb@31119 ` 370` ```val string_of_cterm = string_of_term o term_of; ``` chaieb@31119 ` 371` chaieb@31119 ` 372` ```fun string_of_varpow x k = ``` chaieb@31119 ` 373` ``` if k = 1 then string_of_cterm x ``` chaieb@31119 ` 374` ``` else string_of_cterm x^"^"^string_of_int k; ``` chaieb@31119 ` 375` chaieb@31119 ` 376` ```fun string_of_monomial m = ``` chaieb@31119 ` 377` ``` if Ctermfunc.is_undefined m then "1" else ``` chaieb@31119 ` 378` ``` let val vps = fold_rev (fn (x,k) => fn a => string_of_varpow x k :: a) ``` chaieb@31119 ` 379` ``` (sort humanorder_varpow (Ctermfunc.graph m)) [] ``` chaieb@31119 ` 380` ``` in fold1 (fn s => fn t => s^"*"^t) vps ``` chaieb@31119 ` 381` ``` end; ``` chaieb@31119 ` 382` chaieb@31119 ` 383` ```fun string_of_cmonomial (c,m) = ``` chaieb@31119 ` 384` ``` if Ctermfunc.is_undefined m then Rat.string_of_rat c ``` chaieb@31119 ` 385` ``` else if c =/ rat_1 then string_of_monomial m ``` chaieb@31119 ` 386` ``` else Rat.string_of_rat c ^ "*" ^ string_of_monomial m;; ``` chaieb@31119 ` 387` chaieb@31119 ` 388` ```fun string_of_poly (p:poly) = ``` chaieb@31119 ` 389` ``` if Monomialfunc.is_undefined p then "<<0>>" else ``` chaieb@31119 ` 390` ``` let ``` chaieb@31119 ` 391` ``` val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1 m2) (Monomialfunc.graph p) ``` chaieb@31119 ` 392` ``` val s = fold (fn (m,c) => fn a => ``` chaieb@31119 ` 393` ``` if c >" ``` chaieb@31119 ` 399` ``` end; ``` chaieb@31119 ` 400` chaieb@31119 ` 401` ```(* Conversion from HOL term. *) ``` chaieb@31119 ` 402` chaieb@31119 ` 403` ```local ``` chaieb@31119 ` 404` ``` val neg_tm = @{cterm "uminus :: real => _"} ``` chaieb@31119 ` 405` ``` val add_tm = @{cterm "op + :: real => _"} ``` chaieb@31119 ` 406` ``` val sub_tm = @{cterm "op - :: real => _"} ``` chaieb@31119 ` 407` ``` val mul_tm = @{cterm "op * :: real => _"} ``` chaieb@31119 ` 408` ``` val inv_tm = @{cterm "inverse :: real => _"} ``` chaieb@31119 ` 409` ``` val div_tm = @{cterm "op / :: real => _"} ``` chaieb@31119 ` 410` ``` val pow_tm = @{cterm "op ^ :: real => _"} ``` chaieb@31119 ` 411` ``` val zero_tm = @{cterm "0:: real"} ``` chaieb@31119 ` 412` ``` val is_numeral = can (HOLogic.dest_number o term_of) ``` chaieb@31119 ` 413` ``` fun is_comb t = case t of _\$_ => true | _ => false ``` chaieb@31119 ` 414` ``` fun poly_of_term tm = ``` chaieb@31119 ` 415` ``` if tm aconvc zero_tm then poly_0 ``` chaieb@31119 ` 416` ``` else if RealArith.is_ratconst tm ``` chaieb@31119 ` 417` ``` then poly_const(RealArith.dest_ratconst tm) ``` chaieb@31119 ` 418` ``` else ``` chaieb@31119 ` 419` ``` (let val (lop,r) = Thm.dest_comb tm ``` chaieb@31119 ` 420` ``` in if lop aconvc neg_tm then poly_neg(poly_of_term r) ``` chaieb@31119 ` 421` ``` else if lop aconvc inv_tm then ``` chaieb@31119 ` 422` ``` let val p = poly_of_term r ``` chaieb@31119 ` 423` ``` in if poly_isconst p ``` chaieb@31119 ` 424` ``` then poly_const(Rat.inv (eval Ctermfunc.undefined p)) ``` chaieb@31119 ` 425` ``` else error "poly_of_term: inverse of non-constant polyomial" ``` chaieb@31119 ` 426` ``` end ``` chaieb@31119 ` 427` ``` else (let val (opr,l) = Thm.dest_comb lop ``` chaieb@31119 ` 428` ``` in ``` chaieb@31119 ` 429` ``` if opr aconvc pow_tm andalso is_numeral r ``` chaieb@31119 ` 430` ``` then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r) ``` chaieb@31119 ` 431` ``` else if opr aconvc add_tm ``` chaieb@31119 ` 432` ``` then poly_add (poly_of_term l) (poly_of_term r) ``` chaieb@31119 ` 433` ``` else if opr aconvc sub_tm ``` chaieb@31119 ` 434` ``` then poly_sub (poly_of_term l) (poly_of_term r) ``` chaieb@31119 ` 435` ``` else if opr aconvc mul_tm ``` chaieb@31119 ` 436` ``` then poly_mul (poly_of_term l) (poly_of_term r) ``` chaieb@31119 ` 437` ``` else if opr aconvc div_tm ``` chaieb@31119 ` 438` ``` then let ``` chaieb@31119 ` 439` ``` val p = poly_of_term l ``` chaieb@31119 ` 440` ``` val q = poly_of_term r ``` chaieb@31119 ` 441` ``` in if poly_isconst q then poly_cmul (Rat.inv (eval Ctermfunc.undefined q)) p ``` chaieb@31119 ` 442` ``` else error "poly_of_term: division by non-constant polynomial" ``` chaieb@31119 ` 443` ``` end ``` chaieb@31119 ` 444` ``` else poly_var tm ``` chaieb@31119 ` 445` ``` ``` chaieb@31119 ` 446` ``` end ``` chaieb@31119 ` 447` ``` handle CTERM ("dest_comb",_) => poly_var tm) ``` chaieb@31119 ` 448` ``` end ``` chaieb@31119 ` 449` ``` handle CTERM ("dest_comb",_) => poly_var tm) ``` chaieb@31119 ` 450` ```in ``` chaieb@31119 ` 451` ```val poly_of_term = fn tm => ``` chaieb@31119 ` 452` ``` if type_of (term_of tm) = @{typ real} then poly_of_term tm ``` chaieb@31119 ` 453` ``` else error "poly_of_term: term does not have real type" ``` chaieb@31119 ` 454` ```end; ``` chaieb@31119 ` 455` chaieb@31119 ` 456` ```(* String of vector (just a list of space-separated numbers). *) ``` chaieb@31119 ` 457` chaieb@31119 ` 458` ```fun sdpa_of_vector (v:vector) = ``` chaieb@31119 ` 459` ``` let ``` chaieb@31119 ` 460` ``` val n = dim v ``` chaieb@31119 ` 461` ``` val strs = map (decimalize 20 o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) ``` chaieb@31119 ` 462` ``` in fold1 (fn x => fn y => x ^ " " ^ y) strs ^ "\n" ``` chaieb@31119 ` 463` ``` end; ``` chaieb@31119 ` 464` chaieb@31119 ` 465` ```fun increasing f ord (x,y) = ord (f x, f y); ``` chaieb@31119 ` 466` ```fun triple_int_ord ((a,b,c),(a',b',c')) = ``` chaieb@31119 ` 467` ``` prod_ord int_ord (prod_ord int_ord int_ord) ``` chaieb@31119 ` 468` ``` ((a,(b,c)),(a',(b',c'))); ``` chaieb@31119 ` 469` ```structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord); ``` chaieb@31119 ` 470` chaieb@31119 ` 471` ```(* String for block diagonal matrix numbered k. *) ``` chaieb@31119 ` 472` chaieb@31119 ` 473` ```fun sdpa_of_blockdiagonal k m = ``` chaieb@31119 ` 474` ``` let ``` chaieb@31119 ` 475` ``` val pfx = string_of_int k ^" " ``` chaieb@31119 ` 476` ``` val ents = ``` chaieb@31119 ` 477` ``` Inttriplefunc.fold (fn ((b,i,j), c) => fn a => if i > j then a else ((b,i,j),c)::a) m [] ``` chaieb@31119 ` 478` ``` val entss = sort (increasing fst triple_int_ord ) ents ``` chaieb@31119 ` 479` ``` in fold_rev (fn ((b,i,j),c) => fn a => ``` chaieb@31119 ` 480` ``` pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^ ``` chaieb@31119 ` 481` ``` " " ^ decimalize 20 c ^ "\n" ^ a) entss "" ``` chaieb@31119 ` 482` ``` end; ``` chaieb@31119 ` 483` chaieb@31119 ` 484` ```(* String for a matrix numbered k, in SDPA sparse format. *) ``` chaieb@31119 ` 485` chaieb@31119 ` 486` ```fun sdpa_of_matrix k (m:matrix) = ``` chaieb@31119 ` 487` ``` let ``` chaieb@31119 ` 488` ``` val pfx = string_of_int k ^ " 1 " ``` chaieb@31119 ` 489` ``` val ms = Intpairfunc.fold (fn ((i,j), c) => fn a => if i > j then a else ((i,j),c)::a) (snd m) [] ``` chaieb@31119 ` 490` ``` val mss = sort (increasing fst (prod_ord int_ord int_ord)) ms ``` chaieb@31119 ` 491` ``` in fold_rev (fn ((i,j),c) => fn a => ``` chaieb@31119 ` 492` ``` pfx ^ string_of_int i ^ " " ^ string_of_int j ^ ``` chaieb@31119 ` 493` ``` " " ^ decimalize 20 c ^ "\n" ^ a) mss "" ``` chaieb@31119 ` 494` ``` end;; ``` chaieb@31119 ` 495` chaieb@31119 ` 496` ```(* ------------------------------------------------------------------------- *) ``` chaieb@31119 ` 497` ```(* String in SDPA sparse format for standard SDP problem: *) ``` chaieb@31119 ` 498` ```(* *) ``` chaieb@31119 ` 499` ```(* X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD *) ``` chaieb@31119 ` 500` ```(* Minimize obj_1 * v_1 + ... obj_m * v_m *) ``` chaieb@31119 ` 501` ```(* ------------------------------------------------------------------------- *) ``` chaieb@31119 ` 502` chaieb@31119 ` 503` ```fun sdpa_of_problem comment obj mats = ``` chaieb@31119 ` 504` ``` let ``` chaieb@31119 ` 505` ``` val m = length mats - 1 ``` chaieb@31119 ` 506` ``` val (n,_) = dimensions (hd mats) ``` chaieb@31119 ` 507` ``` in "\"" ^ comment ^ "\"\n" ^ ``` chaieb@31119 ` 508` ``` string_of_int m ^ "\n" ^ ``` chaieb@31119 ` 509` ``` "1\n" ^ ``` chaieb@31119 ` 510` ``` string_of_int n ^ "\n" ^ ``` chaieb@31119 ` 511` ``` sdpa_of_vector obj ^ ``` chaieb@31119 ` 512` ``` fold_rev2 (fn k => fn m => fn a => sdpa_of_matrix (k - 1) m ^ a) (1 upto length mats) mats "" ``` chaieb@31119 ` 513` ``` end; ``` chaieb@31119 ` 514` chaieb@31119 ` 515` ```fun index_char str chr pos = ``` chaieb@31119 ` 516` ``` if pos >= String.size str then ~1 ``` chaieb@31119 ` 517` ``` else if String.sub(str,pos) = chr then pos ``` chaieb@31119 ` 518` ``` else index_char str chr (pos + 1); ``` chaieb@31119 ` 519` ```fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b); ``` chaieb@31119 ` 520` ```fun rat_of_string s = ``` chaieb@31119 ` 521` ``` let val n = index_char s #"/" 0 in ``` chaieb@31119 ` 522` ``` if n = ~1 then s |> IntInf.fromString |> valOf |> Rat.rat_of_int ``` chaieb@31119 ` 523` ``` else ``` chaieb@31119 ` 524` ``` let val SOME numer = IntInf.fromString(String.substring(s,0,n)) ``` chaieb@31119 ` 525` ``` val SOME den = IntInf.fromString (String.substring(s,n+1,String.size s - n - 1)) ``` chaieb@31119 ` 526` ``` in rat_of_quotient(numer, den) ``` chaieb@31119 ` 527` ``` end ``` chaieb@31119 ` 528` ``` end; ``` chaieb@31119 ` 529` chaieb@31119 ` 530` ```fun isspace x = x = " " ; ``` chaieb@31119 ` 531` ```fun isnum x = x mem_string ["0","1","2","3","4","5","6","7","8","9"] ``` chaieb@31119 ` 532` chaieb@31119 ` 533` ```(* More parser basics. *) ``` chaieb@31119 ` 534` chaieb@31119 ` 535` ```local ``` chaieb@31119 ` 536` ``` open Scan ``` chaieb@31119 ` 537` ```in ``` chaieb@31119 ` 538` ``` val word = this_string ``` chaieb@31119 ` 539` ``` fun token s = ``` chaieb@31119 ` 540` ``` repeat (\$\$ " ") |-- word s --| repeat (\$\$ " ") ``` chaieb@31119 ` 541` ``` val numeral = one isnum ``` chaieb@31119 ` 542` ``` val decimalint = bulk numeral >> (rat_of_string o implode) ``` chaieb@31119 ` 543` ``` val decimalfrac = bulk numeral ``` chaieb@31119 ` 544` ``` >> (fn s => rat_of_string(implode s) // pow10 (length s)) ``` chaieb@31119 ` 545` ``` val decimalsig = ``` chaieb@31119 ` 546` ``` decimalint -- option (Scan.\$\$ "." |-- decimalfrac) ``` chaieb@31119 ` 547` ``` >> (fn (h,NONE) => h | (h,SOME x) => h +/ x) ``` chaieb@31119 ` 548` ``` fun signed prs = ``` chaieb@31119 ` 549` ``` \$\$ "-" |-- prs >> Rat.neg ``` chaieb@31119 ` 550` ``` || \$\$ "+" |-- prs ``` chaieb@31119 ` 551` ``` || prs; ``` chaieb@31119 ` 552` chaieb@31119 ` 553` ```fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs ``` chaieb@31119 ` 554` chaieb@31119 ` 555` ``` val exponent = (\$\$ "e" || \$\$ "E") |-- signed decimalint; ``` chaieb@31119 ` 556` chaieb@31119 ` 557` ``` val decimal = signed decimalsig -- (emptyin rat_0|| exponent) ``` chaieb@31119 ` 558` ``` >> (fn (h, x) => h */ pow10 (int_of_rat x)); ``` chaieb@31119 ` 559` ```end; ``` chaieb@31119 ` 560` chaieb@31119 ` 561` ``` fun mkparser p s = ``` chaieb@31119 ` 562` ``` let val (x,rst) = p (explode s) ``` chaieb@31119 ` 563` ``` in if null rst then x ``` chaieb@31119 ` 564` ``` else error "mkparser: unparsed input" ``` chaieb@31119 ` 565` ``` end;; ``` chaieb@31119 ` 566` ```val parse_decimal = mkparser decimal; ``` chaieb@31119 ` 567` chaieb@31119 ` 568` ```fun fix err prs = ``` chaieb@31119 ` 569` ``` prs || (fn x=> error err); ``` chaieb@31119 ` 570` chaieb@31119 ` 571` ```fun listof prs sep err = ``` chaieb@31119 ` 572` ``` prs -- Scan.bulk (sep |-- fix err prs) >> uncurry cons; ``` chaieb@31119 ` 573` chaieb@31119 ` 574` ```(* Parse back a vector. *) ``` chaieb@31119 ` 575` chaieb@31119 ` 576` ``` val vector = ``` chaieb@31119 ` 577` ``` token "{" |-- listof decimal (token ",") "decimal" --| token "}" ``` chaieb@31119 ` 578` ``` >> vector_of_list ``` chaieb@31119 ` 579` ``` val parse_vector = mkparser vector ``` chaieb@31119 ` 580` ``` fun skipupto dscr prs inp = ``` chaieb@31119 ` 581` ``` (dscr |-- prs ``` chaieb@31119 ` 582` ``` || Scan.one (K true) |-- skipupto dscr prs) inp ``` chaieb@31119 ` 583` ``` fun ignore inp = ((),[]) ``` chaieb@31119 ` 584` ``` fun sdpaoutput inp = skipupto (word "xVec" -- token "=") ``` chaieb@31119 ` 585` ``` (vector --| ignore) inp ``` chaieb@31119 ` 586` ``` fun csdpoutput inp = ((decimal -- Scan.bulk (Scan.\$\$ " " |-- Scan.option decimal) >> (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp ``` chaieb@31119 ` 587` ``` val parse_sdpaoutput = mkparser sdpaoutput ``` chaieb@31119 ` 588` ``` val parse_csdpoutput = mkparser csdpoutput ``` chaieb@31119 ` 589` chaieb@31119 ` 590` ```(* Also parse the SDPA output to test success (CSDP yields a return code). *) ``` chaieb@31119 ` 591` chaieb@31119 ` 592` ```local ``` chaieb@31119 ` 593` ``` val prs = ``` chaieb@31119 ` 594` ``` skipupto (word "phase.value" -- token "=") ``` chaieb@31119 ` 595` ``` (Scan.option (Scan.\$\$ "p") -- Scan.option (Scan.\$\$ "d") ``` chaieb@31119 ` 596` ``` -- (word "OPT" || word "FEAS")) ``` chaieb@31119 ` 597` ```in ``` chaieb@31119 ` 598` ``` fun sdpa_run_succeeded s = ``` chaieb@31119 ` 599` ``` (prs (explode s); true) handle _ => false ``` chaieb@31119 ` 600` ```end; ``` chaieb@31119 ` 601` chaieb@31119 ` 602` ```(* The default parameters. Unfortunately this goes to a fixed file. *) ``` chaieb@31119 ` 603` chaieb@31119 ` 604` ```val sdpa_default_parameters = ``` chaieb@31119 ` 605` ```"100 unsigned int maxIteration; \n1.0E-7 double 0.0 < epsilonStar;\n1.0E2 double 0.0 < lambdaStar;\n2.0 double 1.0 < omegaStar;\n-1.0E5 double lowerBound;\n1.0E5 double upperBound;\n0.1 double 0.0 <= betaStar < 1.0;\n0.2 double 0.0 <= betaBar < 1.0, betaStar <= betaBar;\n0.9 double 0.0 < gammaStar < 1.0;\n1.0E-7 double 0.0 < epsilonDash;\n";; ``` chaieb@31119 ` 606` chaieb@31119 ` 607` ```(* These were suggested by Makoto Yamashita for problems where we are *) ``` chaieb@31119 ` 608` ```(* right at the edge of the semidefinite cone, as sometimes happens. *) ``` chaieb@31119 ` 609` chaieb@31119 ` 610` ```val sdpa_alt_parameters = ``` chaieb@31119 ` 611` ```"1000 unsigned int maxIteration;\n1.0E-7 double 0.0 < epsilonStar;\n1.0E4 double 0.0 < lambdaStar;\n2.0 double 1.0 < omegaStar;\n-1.0E5 double lowerBound;\n1.0E5 double upperBound;\n0.1 double 0.0 <= betaStar < 1.0;\n0.2 double 0.0 <= betaBar < 1.0, betaStar <= betaBar;\n0.9 double 0.0 < gammaStar < 1.0;\n1.0E-7 double 0.0 < epsilonDash;\n";; ``` chaieb@31119 ` 612` chaieb@31119 ` 613` ```val sdpa_params = sdpa_alt_parameters;; ``` chaieb@31119 ` 614` chaieb@31119 ` 615` ```(* CSDP parameters; so far I'm sticking with the defaults. *) ``` chaieb@31119 ` 616` chaieb@31119 ` 617` ```val csdp_default_parameters = ``` chaieb@31119 ` 618` ```"axtol=1.0e-8\natytol=1.0e-8\nobjtol=1.0e-8\npinftol=1.0e8\ndinftol=1.0e8\nmaxiter=100\nminstepfrac=0.9\nmaxstepfrac=0.97\nminstepp=1.0e-8\nminstepd=1.0e-8\nusexzgap=1\ntweakgap=0\naffine=0\nprintlevel=1\n";; ``` chaieb@31119 ` 619` chaieb@31119 ` 620` ```val csdp_params = csdp_default_parameters;; ``` chaieb@31119 ` 621` chaieb@31119 ` 622` ```fun tmp_file pre suf = ``` chaieb@31119 ` 623` ``` let val i = string_of_int (round (random())) ``` chaieb@31119 ` 624` ``` val name = Path.append (Path.variable "ISABELLE_TMP") (Path.explode (pre ^ i ^ suf)) ``` chaieb@31119 ` 625` ``` in ``` chaieb@31119 ` 626` ``` if File.exists name then tmp_file pre suf ``` chaieb@31119 ` 627` ``` else name ``` chaieb@31119 ` 628` ``` end; ``` chaieb@31119 ` 629` chaieb@31119 ` 630` ```(* Now call SDPA on a problem and parse back the output. *) ``` chaieb@31119 ` 631` chaieb@31119 ` 632` ```fun run_sdpa dbg obj mats = ``` chaieb@31119 ` 633` ``` let ``` chaieb@31119 ` 634` ``` val input_file = tmp_file "sos" ".dat-s" ``` chaieb@31119 ` 635` ``` val output_file = tmp_file "sos" ".out" ``` chaieb@31119 ` 636` ``` val params_file = tmp_file "param" ".sdpa" ``` chaieb@31119 ` 637` ``` val current_dir = File.pwd() ``` chaieb@31119 ` 638` ``` val _ = File.write input_file ``` chaieb@31119 ` 639` ``` (sdpa_of_problem "" obj mats) ``` chaieb@31119 ` 640` ``` val _ = File.write params_file sdpa_params ``` chaieb@31119 ` 641` ``` val _ = File.cd (Path.variable "ISABELLE_TMP") ``` chaieb@31119 ` 642` ``` val _ = File.system_command ("sdpa "^ (Path.implode input_file) ^ " " ^ ``` chaieb@31119 ` 643` ``` (Path.implode output_file) ^ ``` chaieb@31119 ` 644` ``` (if dbg then "" else "> /dev/null")) ``` chaieb@31119 ` 645` ``` val opr = File.read output_file ``` chaieb@31119 ` 646` ``` in if not(sdpa_run_succeeded opr) then error "sdpa: call failed" ``` chaieb@31119 ` 647` ``` else ``` chaieb@31119 ` 648` ``` let val res = parse_sdpaoutput opr ``` chaieb@31119 ` 649` ``` in ((if dbg then () ``` chaieb@31119 ` 650` ``` else (File.rm input_file; File.rm output_file ; File.cd current_dir)); ``` chaieb@31119 ` 651` ``` res) ``` chaieb@31119 ` 652` ``` end ``` chaieb@31119 ` 653` ``` end; ``` chaieb@31119 ` 654` chaieb@31119 ` 655` ```fun sdpa obj mats = run_sdpa (!debugging) obj mats; ``` chaieb@31119 ` 656` chaieb@31119 ` 657` ```(* The same thing with CSDP. *) ``` chaieb@31119 ` 658` chaieb@31119 ` 659` ```fun run_csdp dbg obj mats = ``` chaieb@31119 ` 660` ``` let ``` chaieb@31119 ` 661` ``` val input_file = tmp_file "sos" ".dat-s" ``` chaieb@31119 ` 662` ``` val output_file = tmp_file "sos" ".out" ``` chaieb@31119 ` 663` ``` val params_file = tmp_file "param" ".csdp" ``` chaieb@31119 ` 664` ``` val current_dir = File.pwd() ``` chaieb@31119 ` 665` ``` val _ = File.write input_file (sdpa_of_problem "" obj mats) ``` chaieb@31119 ` 666` ``` val _ = File.write params_file csdp_params ``` chaieb@31119 ` 667` ``` val _ = File.cd (Path.variable "ISABELLE_TMP") ``` chaieb@31119 ` 668` ``` val rv = system ("csdp "^(Path.implode input_file) ^ " " ``` chaieb@31119 ` 669` ``` ^ (Path.implode output_file) ^ ``` chaieb@31119 ` 670` ``` (if dbg then "" else "> /dev/null")) ``` chaieb@31119 ` 671` ``` val opr = File.read output_file ``` chaieb@31119 ` 672` ``` val res = parse_csdpoutput opr ``` chaieb@31119 ` 673` ``` in ``` chaieb@31119 ` 674` ``` ((if dbg then () ``` chaieb@31119 ` 675` ``` else (File.rm input_file; File.rm output_file ; File.cd current_dir)); ``` chaieb@31119 ` 676` ``` (rv,res)) ``` chaieb@31119 ` 677` ``` end; ``` chaieb@31119 ` 678` chaieb@31119 ` 679` ```fun csdp obj mats = ``` chaieb@31119 ` 680` ``` let ``` chaieb@31119 ` 681` ``` val (rv,res) = run_csdp (!debugging) obj mats ``` chaieb@31119 ` 682` ``` in ``` chaieb@31119 ` 683` ``` ((if rv = 1 orelse rv = 2 then error "csdp: Problem is infeasible" ``` chaieb@31119 ` 684` ``` else if rv = 3 then writeln "csdp warning: Reduced accuracy" ``` chaieb@31119 ` 685` ``` else if rv <> 0 then error ("csdp: error "^string_of_int rv) ``` chaieb@31119 ` 686` ``` else ()); ``` chaieb@31119 ` 687` ``` res) ``` chaieb@31119 ` 688` ``` end; ``` chaieb@31119 ` 689` chaieb@31119 ` 690` ```(* Try some apparently sensible scaling first. Note that this is purely to *) ``` chaieb@31119 ` 691` ```(* get a cleaner translation to floating-point, and doesn't affect any of *) ``` chaieb@31119 ` 692` ```(* the results, in principle. In practice it seems a lot better when there *) ``` chaieb@31119 ` 693` ```(* are extreme numbers in the original problem. *) ``` chaieb@31119 ` 694` chaieb@31119 ` 695` ``` (* Version for (int*int) keys *) ``` chaieb@31119 ` 696` ```local ``` chaieb@31119 ` 697` ``` fun max_rat x y = if x fn a => lcm_rat (denominator_rat c) a) amat acc ``` chaieb@31119 ` 700` ``` fun maximal_element fld amat acc = ``` chaieb@31119 ` 701` ``` fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc ``` chaieb@31119 ` 702` ```fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x ``` chaieb@31119 ` 703` ``` in Real.fromLargeInt a / Real.fromLargeInt b end; ``` chaieb@31119 ` 704` ```in ``` chaieb@31119 ` 705` chaieb@31119 ` 706` ```fun pi_scale_then solver (obj:vector) mats = ``` chaieb@31119 ` 707` ``` let ``` chaieb@31119 ` 708` ``` val cd1 = fold_rev (common_denominator Intpairfunc.fold) mats (rat_1) ``` chaieb@31119 ` 709` ``` val cd2 = common_denominator Intfunc.fold (snd obj) (rat_1) ``` chaieb@31119 ` 710` ``` val mats' = map (Intpairfunc.mapf (fn x => cd1 */ x)) mats ``` chaieb@31119 ` 711` ``` val obj' = vector_cmul cd2 obj ``` chaieb@31119 ` 712` ``` val max1 = fold_rev (maximal_element Intpairfunc.fold) mats' (rat_0) ``` chaieb@31119 ` 713` ``` val max2 = maximal_element Intfunc.fold (snd obj') (rat_0) ``` chaieb@31119 ` 714` ``` val scal1 = pow2 (20 - trunc(Math.ln (float_of_rat max1) / Math.ln 2.0)) ``` chaieb@31119 ` 715` ``` val scal2 = pow2 (20 - trunc(Math.ln (float_of_rat max2) / Math.ln 2.0)) ``` chaieb@31119 ` 716` ``` val mats'' = map (Intpairfunc.mapf (fn x => x */ scal1)) mats' ``` chaieb@31119 ` 717` ``` val obj'' = vector_cmul scal2 obj' ``` chaieb@31119 ` 718` ``` in solver obj'' mats'' ``` chaieb@31119 ` 719` ``` end ``` chaieb@31119 ` 720` ```end; ``` chaieb@31119 ` 721` chaieb@31119 ` 722` ```(* Try some apparently sensible scaling first. Note that this is purely to *) ``` chaieb@31119 ` 723` ```(* get a cleaner translation to floating-point, and doesn't affect any of *) ``` chaieb@31119 ` 724` ```(* the results, in principle. In practice it seems a lot better when there *) ``` chaieb@31119 ` 725` ```(* are extreme numbers in the original problem. *) ``` chaieb@31119 ` 726` chaieb@31119 ` 727` ``` (* Version for (int*int*int) keys *) ``` chaieb@31119 ` 728` ```local ``` chaieb@31119 ` 729` ``` fun max_rat x y = if x fn a => lcm_rat (denominator_rat c) a) amat acc ``` chaieb@31119 ` 732` ``` fun maximal_element fld amat acc = ``` chaieb@31119 ` 733` ``` fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc ``` chaieb@31119 ` 734` ```fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x ``` chaieb@31119 ` 735` ``` in Real.fromLargeInt a / Real.fromLargeInt b end; ``` chaieb@31119 ` 736` ```fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0) ``` chaieb@31119 ` 737` ```in ``` chaieb@31119 ` 738` chaieb@31119 ` 739` ```fun tri_scale_then solver (obj:vector) mats = ``` chaieb@31119 ` 740` ``` let ``` chaieb@31119 ` 741` ``` val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1) ``` chaieb@31119 ` 742` ``` val cd2 = common_denominator Intfunc.fold (snd obj) (rat_1) ``` chaieb@31119 ` 743` ``` val mats' = map (Inttriplefunc.mapf (fn x => cd1 */ x)) mats ``` chaieb@31119 ` 744` ``` val obj' = vector_cmul cd2 obj ``` chaieb@31119 ` 745` ``` val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0) ``` chaieb@31119 ` 746` ``` val max2 = maximal_element Intfunc.fold (snd obj') (rat_0) ``` chaieb@31119 ` 747` ``` val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0)) ``` chaieb@31119 ` 748` ``` val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0)) ``` chaieb@31119 ` 749` ``` val mats'' = map (Inttriplefunc.mapf (fn x => x */ scal1)) mats' ``` chaieb@31119 ` 750` ``` val obj'' = vector_cmul scal2 obj' ``` chaieb@31119 ` 751` ``` in solver obj'' mats'' ``` chaieb@31119 ` 752` ``` end ``` chaieb@31119 ` 753` ```end; ``` chaieb@31119 ` 754` chaieb@31119 ` 755` ```(* Round a vector to "nice" rationals. *) ``` chaieb@31119 ` 756` chaieb@31119 ` 757` ```fun nice_rational n x = round_rat (n */ x) // n;; ``` chaieb@31119 ` 758` ```fun nice_vector n ((d,v) : vector) = ``` chaieb@31119 ` 759` ``` (d, Intfunc.fold (fn (i,c) => fn a => ``` chaieb@31119 ` 760` ``` let val y = nice_rational n c ``` chaieb@31119 ` 761` ``` in if c =/ rat_0 then a ``` chaieb@31119 ` 762` ``` else Intfunc.update (i,y) a end) v Intfunc.undefined):vector ``` chaieb@31119 ` 763` chaieb@31119 ` 764` chaieb@31119 ` 765` ```(* Reduce linear program to SDP (diagonal matrices) and test with CSDP. This *) ``` chaieb@31119 ` 766` ```(* one tests A [-1;x1;..;xn] >= 0 (i.e. left column is negated constants). *) ``` chaieb@31119 ` 767` chaieb@31119 ` 768` ```fun linear_program_basic a = ``` chaieb@31119 ` 769` ``` let ``` chaieb@31119 ` 770` ``` val (m,n) = dimensions a ``` chaieb@31119 ` 771` ``` val mats = map (fn j => diagonal (column j a)) (1 upto n) ``` chaieb@31119 ` 772` ``` val obj = vector_const rat_1 m ``` chaieb@31119 ` 773` ``` val (rv,res) = run_csdp false obj mats ``` chaieb@31119 ` 774` ``` in if rv = 1 orelse rv = 2 then false ``` chaieb@31119 ` 775` ``` else if rv = 0 then true ``` chaieb@31119 ` 776` ``` else error "linear_program: An error occurred in the SDP solver" ``` chaieb@31119 ` 777` ``` end; ``` chaieb@31119 ` 778` chaieb@31119 ` 779` ```(* Alternative interface testing A x >= b for matrix A, vector b. *) ``` chaieb@31119 ` 780` chaieb@31119 ` 781` ```fun linear_program a b = ``` chaieb@31119 ` 782` ``` let val (m,n) = dimensions a ``` chaieb@31119 ` 783` ``` in if dim b <> m then error "linear_program: incompatible dimensions" ``` chaieb@31119 ` 784` ``` else ``` chaieb@31119 ` 785` ``` let ``` chaieb@31119 ` 786` ``` val mats = diagonal b :: map (fn j => diagonal (column j a)) (1 upto n) ``` chaieb@31119 ` 787` ``` val obj = vector_const rat_1 m ``` chaieb@31119 ` 788` ``` val (rv,res) = run_csdp false obj mats ``` chaieb@31119 ` 789` ``` in if rv = 1 orelse rv = 2 then false ``` chaieb@31119 ` 790` ``` else if rv = 0 then true ``` chaieb@31119 ` 791` ``` else error "linear_program: An error occurred in the SDP solver" ``` chaieb@31119 ` 792` ``` end ``` chaieb@31119 ` 793` ``` end; ``` chaieb@31119 ` 794` chaieb@31119 ` 795` ```(* Test whether a point is in the convex hull of others. Rather than use *) ``` chaieb@31119 ` 796` ```(* computational geometry, express as linear inequalities and call CSDP. *) ``` chaieb@31119 ` 797` ```(* This is a bit lazy of me, but it's easy and not such a bottleneck so far. *) ``` chaieb@31119 ` 798` chaieb@31119 ` 799` ```fun in_convex_hull pts pt = ``` chaieb@31119 ` 800` ``` let ``` chaieb@31119 ` 801` ``` val pts1 = (1::pt) :: map (fn x => 1::x) pts ``` chaieb@31119 ` 802` ``` val pts2 = map (fn p => map (fn x => ~x) p @ p) pts1 ``` chaieb@31119 ` 803` ``` val n = length pts + 1 ``` chaieb@31119 ` 804` ``` val v = 2 * (length pt + 1) ``` chaieb@31119 ` 805` ``` val m = v + n - 1 ``` chaieb@31119 ` 806` ``` val mat = ((m,n), ``` chaieb@31119 ` 807` ``` itern 1 pts2 (fn pts => fn j => itern 1 pts ``` chaieb@31119 ` 808` ``` (fn x => fn i => Intpairfunc.update ((i,j), Rat.rat_of_int x))) ``` chaieb@31119 ` 809` ``` (iter (1,n) (fn i => Intpairfunc.update((v + i,i+1), rat_1)) ``` chaieb@31119 ` 810` ``` Intpairfunc.undefined)) ``` chaieb@31119 ` 811` ``` in linear_program_basic mat ``` chaieb@31119 ` 812` ``` end; ``` chaieb@31119 ` 813` chaieb@31119 ` 814` ```(* Filter down a set of points to a minimal set with the same convex hull. *) ``` chaieb@31119 ` 815` chaieb@31119 ` 816` ```local ``` chaieb@31119 ` 817` ``` fun augment1 (m::ms) = if in_convex_hull ms m then ms else ms@[m] ``` chaieb@31119 ` 818` ``` fun augment m ms = funpow 3 augment1 (m::ms) ``` chaieb@31119 ` 819` ```in ``` chaieb@31119 ` 820` ```fun minimal_convex_hull mons = ``` chaieb@31119 ` 821` ``` let val mons' = fold_rev augment (tl mons) [hd mons] ``` chaieb@31119 ` 822` ``` in funpow (length mons') augment1 mons' ``` chaieb@31119 ` 823` ``` end ``` chaieb@31119 ` 824` ```end; ``` chaieb@31119 ` 825` chaieb@31119 ` 826` ```fun dest_ord f x = is_equal (f x); ``` chaieb@31119 ` 827` chaieb@31119 ` 828` ```(* Stuff for "equations" ((int*int*int)->num functions). *) ``` chaieb@31119 ` 829` chaieb@31119 ` 830` ```fun tri_equation_cmul c eq = ``` chaieb@31119 ` 831` ``` if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq; ``` chaieb@31119 ` 832` chaieb@31119 ` 833` ```fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2; ``` chaieb@31119 ` 834` chaieb@31119 ` 835` ```fun tri_equation_eval assig eq = ``` chaieb@31119 ` 836` ``` let fun value v = Inttriplefunc.apply assig v ``` chaieb@31119 ` 837` ``` in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0 ``` chaieb@31119 ` 838` ``` end; ``` chaieb@31119 ` 839` chaieb@31119 ` 840` ```(* Eliminate among linear equations: return unconstrained variables and *) ``` chaieb@31119 ` 841` ```(* assignments for the others in terms of them. We give one pseudo-variable *) ``` chaieb@31119 ` 842` ```(* "one" that's used for a constant term. *) ``` chaieb@31119 ` 843` chaieb@31119 ` 844` ```local ``` chaieb@31119 ` 845` ``` fun extract_first p l = case l of (* FIXME : use find_first instead *) ``` chaieb@31119 ` 846` ``` [] => error "extract_first" ``` chaieb@31119 ` 847` ``` | h::t => if p h then (h,t) else ``` chaieb@31119 ` 848` ``` let val (k,s) = extract_first p t in (k,h::s) end ``` chaieb@31119 ` 849` ```fun eliminate vars dun eqs = case vars of ``` chaieb@31119 ` 850` ``` [] => if forall Inttriplefunc.is_undefined eqs then dun ``` chaieb@31119 ` 851` ``` else raise Unsolvable ``` chaieb@31119 ` 852` ``` | v::vs => ``` chaieb@31119 ` 853` ``` ((let ``` chaieb@31119 ` 854` ``` val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs ``` chaieb@31119 ` 855` ``` val a = Inttriplefunc.apply eq v ``` chaieb@31119 ` 856` ``` val eq' = tri_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq) ``` chaieb@31119 ` 857` ``` fun elim e = ``` chaieb@31119 ` 858` ``` let val b = Inttriplefunc.tryapplyd e v rat_0 ``` chaieb@31119 ` 859` ``` in if b =/ rat_0 then e else ``` chaieb@31119 ` 860` ``` tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq) ``` chaieb@31119 ` 861` ``` end ``` chaieb@31119 ` 862` ``` in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs) ``` chaieb@31119 ` 863` ``` end) ``` chaieb@31119 ` 864` ``` handle ERROR _ => eliminate vs dun eqs) ``` chaieb@31119 ` 865` ```in ``` chaieb@31119 ` 866` ```fun tri_eliminate_equations one vars eqs = ``` chaieb@31119 ` 867` ``` let ``` chaieb@31119 ` 868` ``` val assig = eliminate vars Inttriplefunc.undefined eqs ``` chaieb@31119 ` 869` ``` val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] ``` chaieb@31119 ` 870` ``` in (distinct (dest_ord triple_int_ord) vs, assig) ``` chaieb@31119 ` 871` ``` end ``` chaieb@31119 ` 872` ```end; ``` chaieb@31119 ` 873` chaieb@31119 ` 874` ```(* Eliminate all variables, in an essentially arbitrary order. *) ``` chaieb@31119 ` 875` chaieb@31119 ` 876` ```fun tri_eliminate_all_equations one = ``` chaieb@31119 ` 877` ``` let ``` chaieb@31119 ` 878` ``` fun choose_variable eq = ``` chaieb@31119 ` 879` ``` let val (v,_) = Inttriplefunc.choose eq ``` chaieb@31119 ` 880` ``` in if is_equal (triple_int_ord(v,one)) then ``` chaieb@31119 ` 881` ``` let val eq' = Inttriplefunc.undefine v eq ``` chaieb@31119 ` 882` ``` in if Inttriplefunc.is_undefined eq' then error "choose_variable" ``` chaieb@31119 ` 883` ``` else fst (Inttriplefunc.choose eq') ``` chaieb@31119 ` 884` ``` end ``` chaieb@31119 ` 885` ``` else v ``` chaieb@31119 ` 886` ``` end ``` chaieb@31119 ` 887` ``` fun eliminate dun eqs = case eqs of ``` chaieb@31119 ` 888` ``` [] => dun ``` chaieb@31119 ` 889` ``` | eq::oeqs => ``` chaieb@31119 ` 890` ``` if Inttriplefunc.is_undefined eq then eliminate dun oeqs else ``` chaieb@31119 ` 891` ``` let val v = choose_variable eq ``` chaieb@31119 ` 892` ``` val a = Inttriplefunc.apply eq v ``` chaieb@31119 ` 893` ``` val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a) ``` chaieb@31119 ` 894` ``` (Inttriplefunc.undefine v eq) ``` chaieb@31119 ` 895` ``` fun elim e = ``` chaieb@31119 ` 896` ``` let val b = Inttriplefunc.tryapplyd e v rat_0 ``` chaieb@31119 ` 897` ``` in if b =/ rat_0 then e ``` chaieb@31119 ` 898` ``` else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq) ``` chaieb@31119 ` 899` ``` end ``` chaieb@31119 ` 900` ``` in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun)) ``` chaieb@31119 ` 901` ``` (map elim oeqs) ``` chaieb@31119 ` 902` ``` end ``` chaieb@31119 ` 903` ```in fn eqs => ``` chaieb@31119 ` 904` ``` let ``` chaieb@31119 ` 905` ``` val assig = eliminate Inttriplefunc.undefined eqs ``` chaieb@31119 ` 906` ``` val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] ``` chaieb@31119 ` 907` ``` in (distinct (dest_ord triple_int_ord) vs,assig) ``` chaieb@31119 ` 908` ``` end ``` chaieb@31119 ` 909` ```end; ``` chaieb@31119 ` 910` ``` ``` chaieb@31119 ` 911` ```(* Solve equations by assigning arbitrary numbers. *) ``` chaieb@31119 ` 912` chaieb@31119 ` 913` ```fun tri_solve_equations one eqs = ``` chaieb@31119 ` 914` ``` let ``` chaieb@31119 ` 915` ``` val (vars,assigs) = tri_eliminate_all_equations one eqs ``` chaieb@31119 ` 916` ``` val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars ``` chaieb@31119 ` 917` ``` (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1)) ``` chaieb@31119 ` 918` ``` val ass = ``` chaieb@31119 ` 919` ``` Inttriplefunc.combine (curry op +/) (K false) ``` chaieb@31119 ` 920` ``` (Inttriplefunc.mapf (tri_equation_eval vfn) assigs) vfn ``` chaieb@31119 ` 921` ``` in if forall (fn e => tri_equation_eval ass e =/ rat_0) eqs ``` chaieb@31119 ` 922` ``` then Inttriplefunc.undefine one ass else raise Sanity ``` chaieb@31119 ` 923` ``` end; ``` chaieb@31119 ` 924` chaieb@31119 ` 925` ```(* Multiply equation-parametrized poly by regular poly and add accumulator. *) ``` chaieb@31119 ` 926` chaieb@31119 ` 927` ```fun tri_epoly_pmul p q acc = ``` chaieb@31119 ` 928` ``` Monomialfunc.fold (fn (m1, c) => fn a => ``` chaieb@31119 ` 929` ``` Monomialfunc.fold (fn (m2,e) => fn b => ``` chaieb@31119 ` 930` ``` let val m = monomial_mul m1 m2 ``` chaieb@31119 ` 931` ``` val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined ``` chaieb@31119 ` 932` ``` in Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b ``` chaieb@31119 ` 933` ``` end) q a) p acc ; ``` chaieb@31119 ` 934` chaieb@31119 ` 935` ```(* Usual operations on equation-parametrized poly. *) ``` chaieb@31119 ` 936` chaieb@31119 ` 937` ```fun tri_epoly_cmul c l = ``` chaieb@31119 ` 938` ``` if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (tri_equation_cmul c) l;; ``` chaieb@31119 ` 939` chaieb@31119 ` 940` ```val tri_epoly_neg = tri_epoly_cmul (Rat.rat_of_int ~1); ``` chaieb@31119 ` 941` chaieb@31119 ` 942` ```val tri_epoly_add = Inttriplefunc.combine tri_equation_add Inttriplefunc.is_undefined; ``` chaieb@31119 ` 943` chaieb@31119 ` 944` ```fun tri_epoly_sub p q = tri_epoly_add p (tri_epoly_neg q);; ``` chaieb@31119 ` 945` chaieb@31119 ` 946` ```(* Stuff for "equations" ((int*int)->num functions). *) ``` chaieb@31119 ` 947` chaieb@31119 ` 948` ```fun pi_equation_cmul c eq = ``` chaieb@31119 ` 949` ``` if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq; ``` chaieb@31119 ` 950` chaieb@31119 ` 951` ```fun pi_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2; ``` chaieb@31119 ` 952` chaieb@31119 ` 953` ```fun pi_equation_eval assig eq = ``` chaieb@31119 ` 954` ``` let fun value v = Inttriplefunc.apply assig v ``` chaieb@31119 ` 955` ``` in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0 ``` chaieb@31119 ` 956` ``` end; ``` chaieb@31119 ` 957` chaieb@31119 ` 958` ```(* Eliminate among linear equations: return unconstrained variables and *) ``` chaieb@31119 ` 959` ```(* assignments for the others in terms of them. We give one pseudo-variable *) ``` chaieb@31119 ` 960` ```(* "one" that's used for a constant term. *) ``` chaieb@31119 ` 961` chaieb@31119 ` 962` ```local ``` chaieb@31119 ` 963` ```fun extract_first p l = case l of ``` chaieb@31119 ` 964` ``` [] => error "extract_first" ``` chaieb@31119 ` 965` ``` | h::t => if p h then (h,t) else ``` chaieb@31119 ` 966` ``` let val (k,s) = extract_first p t in (k,h::s) end ``` chaieb@31119 ` 967` ```fun eliminate vars dun eqs = case vars of ``` chaieb@31119 ` 968` ``` [] => if forall Inttriplefunc.is_undefined eqs then dun ``` chaieb@31119 ` 969` ``` else raise Unsolvable ``` chaieb@31119 ` 970` ``` | v::vs => ``` chaieb@31119 ` 971` ``` let ``` chaieb@31119 ` 972` ``` val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs ``` chaieb@31119 ` 973` ``` val a = Inttriplefunc.apply eq v ``` chaieb@31119 ` 974` ``` val eq' = pi_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq) ``` chaieb@31119 ` 975` ``` fun elim e = ``` chaieb@31119 ` 976` ``` let val b = Inttriplefunc.tryapplyd e v rat_0 ``` chaieb@31119 ` 977` ``` in if b =/ rat_0 then e else ``` chaieb@31119 ` 978` ``` pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq) ``` chaieb@31119 ` 979` ``` end ``` chaieb@31119 ` 980` ``` in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs) ``` chaieb@31119 ` 981` ``` end ``` chaieb@31119 ` 982` ``` handle ERROR _ => eliminate vs dun eqs ``` chaieb@31119 ` 983` ```in ``` chaieb@31119 ` 984` ```fun pi_eliminate_equations one vars eqs = ``` chaieb@31119 ` 985` ``` let ``` chaieb@31119 ` 986` ``` val assig = eliminate vars Inttriplefunc.undefined eqs ``` chaieb@31119 ` 987` ``` val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] ``` chaieb@31119 ` 988` ``` in (distinct (dest_ord triple_int_ord) vs, assig) ``` chaieb@31119 ` 989` ``` end ``` chaieb@31119 ` 990` ```end; ``` chaieb@31119 ` 991` chaieb@31119 ` 992` ```(* Eliminate all variables, in an essentially arbitrary order. *) ``` chaieb@31119 ` 993` chaieb@31119 ` 994` ```fun pi_eliminate_all_equations one = ``` chaieb@31119 ` 995` ``` let ``` chaieb@31119 ` 996` ``` fun choose_variable eq = ``` chaieb@31119 ` 997` ``` let val (v,_) = Inttriplefunc.choose eq ``` chaieb@31119 ` 998` ``` in if is_equal (triple_int_ord(v,one)) then ``` chaieb@31119 ` 999` ``` let val eq' = Inttriplefunc.undefine v eq ``` chaieb@31119 ` 1000` ``` in if Inttriplefunc.is_undefined eq' then error "choose_variable" ``` chaieb@31119 ` 1001` ``` else fst (Inttriplefunc.choose eq') ``` chaieb@31119 ` 1002` ``` end ``` chaieb@31119 ` 1003` ``` else v ``` chaieb@31119 ` 1004` ``` end ``` chaieb@31119 ` 1005` ``` fun eliminate dun eqs = case eqs of ``` chaieb@31119 ` 1006` ``` [] => dun ``` chaieb@31119 ` 1007` ``` | eq::oeqs => ``` chaieb@31119 ` 1008` ``` if Inttriplefunc.is_undefined eq then eliminate dun oeqs else ``` chaieb@31119 ` 1009` ``` let val v = choose_variable eq ``` chaieb@31119 ` 1010` ``` val a = Inttriplefunc.apply eq v ``` chaieb@31119 ` 1011` ``` val eq' = pi_equation_cmul ((Rat.rat_of_int ~1) // a) ``` chaieb@31119 ` 1012` ``` (Inttriplefunc.undefine v eq) ``` chaieb@31119 ` 1013` ``` fun elim e = ``` chaieb@31119 ` 1014` ``` let val b = Inttriplefunc.tryapplyd e v rat_0 ``` chaieb@31119 ` 1015` ``` in if b =/ rat_0 then e ``` chaieb@31119 ` 1016` ``` else pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq) ``` chaieb@31119 ` 1017` ``` end ``` chaieb@31119 ` 1018` ``` in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun)) ``` chaieb@31119 ` 1019` ``` (map elim oeqs) ``` chaieb@31119 ` 1020` ``` end ``` chaieb@31119 ` 1021` ```in fn eqs => ``` chaieb@31119 ` 1022` ``` let ``` chaieb@31119 ` 1023` ``` val assig = eliminate Inttriplefunc.undefined eqs ``` chaieb@31119 ` 1024` ``` val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] ``` chaieb@31119 ` 1025` ``` in (distinct (dest_ord triple_int_ord) vs,assig) ``` chaieb@31119 ` 1026` ``` end ``` chaieb@31119 ` 1027` ```end; ``` chaieb@31119 ` 1028` ``` ``` chaieb@31119 ` 1029` ```(* Solve equations by assigning arbitrary numbers. *) ``` chaieb@31119 ` 1030` chaieb@31119 ` 1031` ```fun pi_solve_equations one eqs = ``` chaieb@31119 ` 1032` ``` let ``` chaieb@31119 ` 1033` ``` val (vars,assigs) = pi_eliminate_all_equations one eqs ``` chaieb@31119 ` 1034` ``` val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars ``` chaieb@31119 ` 1035` ``` (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1)) ``` chaieb@31119 ` 1036` ``` val ass = ``` chaieb@31119 ` 1037` ``` Inttriplefunc.combine (curry op +/) (K false) ``` chaieb@31119 ` 1038` ``` (Inttriplefunc.mapf (pi_equation_eval vfn) assigs) vfn ``` chaieb@31119 ` 1039` ``` in if forall (fn e => pi_equation_eval ass e =/ rat_0) eqs ``` chaieb@31119 ` 1040` ``` then Inttriplefunc.undefine one ass else raise Sanity ``` chaieb@31119 ` 1041` ``` end; ``` chaieb@31119 ` 1042` chaieb@31119 ` 1043` ```(* Multiply equation-parametrized poly by regular poly and add accumulator. *) ``` chaieb@31119 ` 1044` chaieb@31119 ` 1045` ```fun pi_epoly_pmul p q acc = ``` chaieb@31119 ` 1046` ``` Monomialfunc.fold (fn (m1, c) => fn a => ``` chaieb@31119 ` 1047` ``` Monomialfunc.fold (fn (m2,e) => fn b => ``` chaieb@31119 ` 1048` ``` let val m = monomial_mul m1 m2 ``` chaieb@31119 ` 1049` ``` val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined ``` chaieb@31119 ` 1050` ``` in Monomialfunc.update (m,pi_equation_add (pi_equation_cmul c e) es) b ``` chaieb@31119 ` 1051` ``` end) q a) p acc ; ``` chaieb@31119 ` 1052` chaieb@31119 ` 1053` ```(* Usual operations on equation-parametrized poly. *) ``` chaieb@31119 ` 1054` chaieb@31119 ` 1055` ```fun pi_epoly_cmul c l = ``` chaieb@31119 ` 1056` ``` if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (pi_equation_cmul c) l;; ``` chaieb@31119 ` 1057` chaieb@31119 ` 1058` ```val pi_epoly_neg = pi_epoly_cmul (Rat.rat_of_int ~1); ``` chaieb@31119 ` 1059` chaieb@31119 ` 1060` ```val pi_epoly_add = Inttriplefunc.combine pi_equation_add Inttriplefunc.is_undefined; ``` chaieb@31119 ` 1061` chaieb@31119 ` 1062` ```fun pi_epoly_sub p q = pi_epoly_add p (pi_epoly_neg q);; ``` chaieb@31119 ` 1063` chaieb@31119 ` 1064` ```fun allpairs f l1 l2 = fold_rev (fn x => (curry (op @)) (map (f x) l2)) l1 []; ``` chaieb@31119 ` 1065` chaieb@31119 ` 1066` ```(* Hence produce the "relevant" monomials: those whose squares lie in the *) ``` chaieb@31119 ` 1067` ```(* Newton polytope of the monomials in the input. (This is enough according *) ``` chaieb@31119 ` 1068` ```(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal, *) ``` chaieb@31119 ` 1069` ```(* vol 45, pp. 363--374, 1978. *) ``` chaieb@31119 ` 1070` ```(* *) ``` chaieb@31119 ` 1071` ```(* These are ordered in sort of decreasing degree. In particular the *) ``` chaieb@31119 ` 1072` ```(* constant monomial is last; this gives an order in diagonalization of the *) ``` chaieb@31119 ` 1073` ```(* quadratic form that will tend to display constants. *) ``` chaieb@31119 ` 1074` chaieb@31119 ` 1075` ```fun newton_polytope pol = ``` chaieb@31119 ` 1076` ``` let ``` chaieb@31119 ` 1077` ``` val vars = poly_variables pol ``` chaieb@31119 ` 1078` ``` val mons = map (fn m => map (fn x => monomial_degree x m) vars) ``` chaieb@31119 ` 1079` ``` (Monomialfunc.dom pol) ``` chaieb@31119 ` 1080` ``` val ds = map (fn x => (degree x pol + 1) div 2) vars ``` chaieb@31119 ` 1081` ``` val all = fold_rev (fn n => allpairs cons (0 upto n)) ds [[]] ``` chaieb@31119 ` 1082` ``` val mons' = minimal_convex_hull mons ``` chaieb@31119 ` 1083` ``` val all' = ``` chaieb@31119 ` 1084` ``` filter (fn m => in_convex_hull mons' (map (fn x => 2 * x) m)) all ``` chaieb@31119 ` 1085` ``` in map (fn m => fold_rev2 (fn v => fn i => fn a => if i = 0 then a else Ctermfunc.update (v,i) a) ``` chaieb@31119 ` 1086` ``` vars m monomial_1) (rev all') ``` chaieb@31119 ` 1087` ``` end; ``` chaieb@31119 ` 1088` chaieb@31119 ` 1089` ```(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form. *) ``` chaieb@31119 ` 1090` chaieb@31119 ` 1091` ```local ``` chaieb@31119 ` 1092` ```fun diagonalize n i m = ``` chaieb@31119 ` 1093` ``` if Intpairfunc.is_undefined (snd m) then [] ``` chaieb@31119 ` 1094` ``` else ``` chaieb@31119 ` 1095` ``` let val a11 = Intpairfunc.tryapplyd (snd m) (i,i) rat_0 ``` chaieb@31119 ` 1096` ``` in if a11 fn a => ``` chaieb@31119 ` 1104` ``` let val y = c // a11 ``` chaieb@31119 ` 1105` ``` in if y = rat_0 then a else Intfunc.update (i,y) a ``` chaieb@31119 ` 1106` ``` end) (snd v) Intfunc.undefined) ``` chaieb@31119 ` 1107` ``` fun upt0 x y a = if y = rat_0 then a else Intpairfunc.update (x,y) a ``` chaieb@31119 ` 1108` ``` val m' = ``` chaieb@31119 ` 1109` ``` ((n,n), ``` chaieb@31119 ` 1110` ``` iter (i+1,n) (fn j => ``` chaieb@31119 ` 1111` ``` iter (i+1,n) (fn k => ``` chaieb@31119 ` 1112` ``` (upt0 (j,k) (Intpairfunc.tryapplyd (snd m) (j,k) rat_0 -/ Intfunc.tryapplyd (snd v) j rat_0 */ Intfunc.tryapplyd (snd v') k rat_0)))) ``` chaieb@31119 ` 1113` ``` Intpairfunc.undefined) ``` chaieb@31119 ` 1114` ``` in (a11,v')::diagonalize n (i + 1) m' ``` chaieb@31119 ` 1115` ``` end ``` chaieb@31119 ` 1116` ``` end ``` chaieb@31119 ` 1117` ```in ``` chaieb@31119 ` 1118` ```fun diag m = ``` chaieb@31119 ` 1119` ``` let ``` chaieb@31119 ` 1120` ``` val nn = dimensions m ``` chaieb@31119 ` 1121` ``` val n = fst nn ``` chaieb@31119 ` 1122` ``` in if snd nn <> n then error "diagonalize: non-square matrix" ``` chaieb@31119 ` 1123` ``` else diagonalize n 1 m ``` chaieb@31119 ` 1124` ``` end ``` chaieb@31119 ` 1125` ```end; ``` chaieb@31119 ` 1126` chaieb@31119 ` 1127` ```fun gcd_rat a b = Rat.rat_of_int (Integer.gcd (int_of_rat a) (int_of_rat b)); ``` chaieb@31119 ` 1128` chaieb@31119 ` 1129` ```(* Adjust a diagonalization to collect rationals at the start. *) ``` chaieb@31119 ` 1130` ``` (* FIXME : Potentially polymorphic keys, but here only: integers!! *) ``` chaieb@31119 ` 1131` ```local ``` chaieb@31119 ` 1132` ``` fun upd0 x y a = if y =/ rat_0 then a else Intfunc.update(x,y) a; ``` chaieb@31119 ` 1133` ``` fun mapa f (d,v) = ``` chaieb@31119 ` 1134` ``` (d, Intfunc.fold (fn (i,c) => fn a => upd0 i (f c) a) v Intfunc.undefined) ``` chaieb@31119 ` 1135` ``` fun adj (c,l) = ``` chaieb@31119 ` 1136` ``` let val a = ``` chaieb@31119 ` 1137` ``` Intfunc.fold (fn (i,c) => fn a => lcm_rat a (denominator_rat c)) ``` chaieb@31119 ` 1138` ``` (snd l) rat_1 // ``` chaieb@31119 ` 1139` ``` Intfunc.fold (fn (i,c) => fn a => gcd_rat a (numerator_rat c)) ``` chaieb@31119 ` 1140` ``` (snd l) rat_0 ``` chaieb@31119 ` 1141` ``` in ((c // (a */ a)),mapa (fn x => a */ x) l) ``` chaieb@31119 ` 1142` ``` end ``` chaieb@31119 ` 1143` ```in ``` chaieb@31119 ` 1144` ```fun deration d = if null d then (rat_0,d) else ``` chaieb@31119 ` 1145` ``` let val d' = map adj d ``` chaieb@31119 ` 1146` ``` val a = fold (lcm_rat o denominator_rat o fst) d' rat_1 // ``` chaieb@31119 ` 1147` ``` fold (gcd_rat o numerator_rat o fst) d' rat_0 ``` chaieb@31119 ` 1148` ``` in ((rat_1 // a),map (fn (c,l) => (a */ c,l)) d') ``` chaieb@31119 ` 1149` ``` end ``` chaieb@31119 ` 1150` ```end; ``` chaieb@31119 ` 1151` ``` ``` chaieb@31119 ` 1152` ```(* Enumeration of monomials with given multidegree bound. *) ``` chaieb@31119 ` 1153` chaieb@31119 ` 1154` ```fun enumerate_monomials d vars = ``` chaieb@31119 ` 1155` ``` if d < 0 then [] ``` chaieb@31119 ` 1156` ``` else if d = 0 then [Ctermfunc.undefined] ``` chaieb@31119 ` 1157` ``` else if null vars then [monomial_1] else ``` chaieb@31119 ` 1158` ``` let val alts = ``` chaieb@31119 ` 1159` ``` map (fn k => let val oths = enumerate_monomials (d - k) (tl vars) ``` chaieb@31119 ` 1160` ``` in map (fn ks => if k = 0 then ks else Ctermfunc.update (hd vars, k) ks) oths end) (0 upto d) ``` chaieb@31119 ` 1161` ``` in fold1 (curry op @) alts ``` chaieb@31119 ` 1162` ``` end; ``` chaieb@31119 ` 1163` chaieb@31119 ` 1164` ```(* Enumerate products of distinct input polys with degree <= d. *) ``` chaieb@31119 ` 1165` ```(* We ignore any constant input polynomials. *) ``` chaieb@31119 ` 1166` ```(* Give the output polynomial and a record of how it was derived. *) ``` chaieb@31119 ` 1167` chaieb@31119 ` 1168` ```local ``` chaieb@31119 ` 1169` ``` open RealArith ``` chaieb@31119 ` 1170` ```in ``` chaieb@31119 ` 1171` ```fun enumerate_products d pols = ``` chaieb@31119 ` 1172` ```if d = 0 then [(poly_const rat_1,Rational_lt rat_1)] ``` chaieb@31119 ` 1173` ```else if d < 0 then [] else ``` chaieb@31119 ` 1174` ```case pols of ``` chaieb@31119 ` 1175` ``` [] => [(poly_const rat_1,Rational_lt rat_1)] ``` chaieb@31119 ` 1176` ``` | (p,b)::ps => ``` chaieb@31119 ` 1177` ``` let val e = multidegree p ``` chaieb@31119 ` 1178` ``` in if e = 0 then enumerate_products d ps else ``` chaieb@31119 ` 1179` ``` enumerate_products d ps @ ``` chaieb@31119 ` 1180` ``` map (fn (q,c) => (poly_mul p q,Product(b,c))) ``` chaieb@31119 ` 1181` ``` (enumerate_products (d - e) ps) ``` chaieb@31119 ` 1182` ``` end ``` chaieb@31119 ` 1183` ```end; ``` chaieb@31119 ` 1184` chaieb@31119 ` 1185` ```(* Convert regular polynomial. Note that we treat (0,0,0) as -1. *) ``` chaieb@31119 ` 1186` chaieb@31119 ` 1187` ```fun epoly_of_poly p = ``` chaieb@31119 ` 1188` ``` Monomialfunc.fold (fn (m,c) => fn a => Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p Monomialfunc.undefined; ``` chaieb@31119 ` 1189` chaieb@31119 ` 1190` ```(* String for block diagonal matrix numbered k. *) ``` chaieb@31119 ` 1191` chaieb@31119 ` 1192` ```fun sdpa_of_blockdiagonal k m = ``` chaieb@31119 ` 1193` ``` let ``` chaieb@31119 ` 1194` ``` val pfx = string_of_int k ^" " ``` chaieb@31119 ` 1195` ``` val ents = ``` chaieb@31119 ` 1196` ``` Inttriplefunc.fold ``` chaieb@31119 ` 1197` ``` (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a) ``` chaieb@31119 ` 1198` ``` m [] ``` chaieb@31119 ` 1199` ``` val entss = sort (increasing fst triple_int_ord) ents ``` chaieb@31119 ` 1200` ``` in fold_rev (fn ((b,i,j),c) => fn a => ``` chaieb@31119 ` 1201` ``` pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^ ``` chaieb@31119 ` 1202` ``` " " ^ decimalize 20 c ^ "\n" ^ a) entss "" ``` chaieb@31119 ` 1203` ``` end; ``` chaieb@31119 ` 1204` chaieb@31119 ` 1205` ```(* SDPA for problem using block diagonal (i.e. multiple SDPs) *) ``` chaieb@31119 ` 1206` chaieb@31119 ` 1207` ```fun sdpa_of_blockproblem comment nblocks blocksizes obj mats = ``` chaieb@31119 ` 1208` ``` let val m = length mats - 1 ``` chaieb@31119 ` 1209` ``` in "\"" ^ comment ^ "\"\n" ^ ``` chaieb@31119 ` 1210` ``` string_of_int m ^ "\n" ^ ``` chaieb@31119 ` 1211` ``` string_of_int nblocks ^ "\n" ^ ``` chaieb@31119 ` 1212` ``` (fold1 (fn s => fn t => s^" "^t) (map string_of_int blocksizes)) ^ ``` chaieb@31119 ` 1213` ``` "\n" ^ ``` chaieb@31119 ` 1214` ``` sdpa_of_vector obj ^ ``` chaieb@31119 ` 1215` ``` fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a) ``` chaieb@31119 ` 1216` ``` (1 upto length mats) mats "" ``` chaieb@31119 ` 1217` ``` end; ``` chaieb@31119 ` 1218` chaieb@31119 ` 1219` ```(* Hence run CSDP on a problem in block diagonal form. *) ``` chaieb@31119 ` 1220` chaieb@31119 ` 1221` ```fun run_csdp dbg nblocks blocksizes obj mats = ``` chaieb@31119 ` 1222` ``` let ``` chaieb@31119 ` 1223` ``` val input_file = tmp_file "sos" ".dat-s" ``` chaieb@31119 ` 1224` ``` val output_file = tmp_file "sos" ".out" ``` chaieb@31119 ` 1225` ``` val params_file = tmp_file "param" ".csdp" ``` chaieb@31119 ` 1226` ``` val _ = File.write input_file ``` chaieb@31119 ` 1227` ``` (sdpa_of_blockproblem "" nblocks blocksizes obj mats) ``` chaieb@31119 ` 1228` ``` val _ = File.write params_file csdp_params ``` chaieb@31119 ` 1229` ``` val current_dir = File.pwd() ``` chaieb@31119 ` 1230` ``` val _ = File.cd (Path.variable "ISABELLE_TMP") ``` chaieb@31119 ` 1231` ``` val rv = system ("csdp "^(Path.implode input_file) ^ " " ``` chaieb@31119 ` 1232` ``` ^ (Path.implode output_file) ^ ``` chaieb@31119 ` 1233` ``` (if dbg then "" else "> /dev/null")) ``` chaieb@31119 ` 1234` ``` val opr = File.read output_file ``` chaieb@31119 ` 1235` ``` val res = parse_csdpoutput opr ``` chaieb@31119 ` 1236` ``` in ``` chaieb@31119 ` 1237` ``` ((if dbg then () ``` chaieb@31119 ` 1238` ``` else (File.rm input_file ; File.rm output_file ; File.cd current_dir)); ``` chaieb@31119 ` 1239` ``` (rv,res)) ``` chaieb@31119 ` 1240` ``` end; ``` chaieb@31119 ` 1241` chaieb@31119 ` 1242` ```fun csdp nblocks blocksizes obj mats = ``` chaieb@31119 ` 1243` ``` let ``` chaieb@31119 ` 1244` ``` val (rv,res) = run_csdp (!debugging) nblocks blocksizes obj mats ``` chaieb@31119 ` 1245` ``` in ((if rv = 1 orelse rv = 2 then error "csdp: Problem is infeasible" ``` chaieb@31119 ` 1246` ``` else if rv = 3 then writeln "csdp warning: Reduced accuracy" ``` chaieb@31119 ` 1247` ``` else if rv <> 0 then error ("csdp: error "^string_of_int rv) ``` chaieb@31119 ` 1248` ``` else ()); ``` chaieb@31119 ` 1249` ``` res) ``` chaieb@31119 ` 1250` ``` end; ``` chaieb@31119 ` 1251` chaieb@31119 ` 1252` ```(* 3D versions of matrix operations to consider blocks separately. *) ``` chaieb@31119 ` 1253` chaieb@31119 ` 1254` ```val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0); ``` chaieb@31119 ` 1255` ```fun bmatrix_cmul c bm = ``` chaieb@31119 ` 1256` ``` if c =/ rat_0 then Inttriplefunc.undefined ``` chaieb@31119 ` 1257` ``` else Inttriplefunc.mapf (fn x => c */ x) bm; ``` chaieb@31119 ` 1258` chaieb@31119 ` 1259` ```val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1); ``` chaieb@31119 ` 1260` ```fun bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);; ``` chaieb@31119 ` 1261` chaieb@31119 ` 1262` ```(* Smash a block matrix into components. *) ``` chaieb@31119 ` 1263` chaieb@31119 ` 1264` ```fun blocks blocksizes bm = ``` chaieb@31119 ` 1265` ``` map (fn (bs,b0) => ``` chaieb@31119 ` 1266` ``` let val m = Inttriplefunc.fold ``` chaieb@31119 ` 1267` ``` (fn ((b,i,j),c) => fn a => if b = b0 then Intpairfunc.update ((i,j),c) a else a) bm Intpairfunc.undefined ``` chaieb@31119 ` 1268` ``` val d = Intpairfunc.fold (fn ((i,j),c) => fn a => max a (max i j)) m 0 ``` chaieb@31119 ` 1269` ``` in (((bs,bs),m):matrix) end) ``` chaieb@31119 ` 1270` ``` (blocksizes ~~ (1 upto length blocksizes));; ``` chaieb@31119 ` 1271` chaieb@31119 ` 1272` ```(* FIXME : Get rid of this !!!*) ``` chaieb@31119 ` 1273` ```fun tryfind f [] = error "tryfind" ``` chaieb@31119 ` 1274` ``` | tryfind f (x::xs) = (f x handle ERROR _ => tryfind f xs); ``` chaieb@31119 ` 1275` chaieb@31119 ` 1276` chaieb@31119 ` 1277` ```(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *) ``` chaieb@31119 ` 1278` chaieb@31119 ` 1279` chaieb@31119 ` 1280` ```local ``` chaieb@31119 ` 1281` ``` open RealArith ``` chaieb@31119 ` 1282` ```in ``` chaieb@31119 ` 1283` ```fun real_positivnullstellensatz_general linf d eqs leqs pol = ``` chaieb@31119 ` 1284` ```let ``` chaieb@31119 ` 1285` ``` val vars = fold_rev (curry (gen_union (op aconvc)) o poly_variables) ``` chaieb@31119 ` 1286` ``` (pol::eqs @ map fst leqs) [] ``` chaieb@31119 ` 1287` ``` val monoid = if linf then ``` chaieb@31119 ` 1288` ``` (poly_const rat_1,Rational_lt rat_1):: ``` chaieb@31119 ` 1289` ``` (filter (fn (p,c) => multidegree p <= d) leqs) ``` chaieb@31119 ` 1290` ``` else enumerate_products d leqs ``` chaieb@31119 ` 1291` ``` val nblocks = length monoid ``` chaieb@31119 ` 1292` ``` fun mk_idmultiplier k p = ``` chaieb@31119 ` 1293` ``` let ``` chaieb@31119 ` 1294` ``` val e = d - multidegree p ``` chaieb@31119 ` 1295` ``` val mons = enumerate_monomials e vars ``` chaieb@31119 ` 1296` ``` val nons = mons ~~ (1 upto length mons) ``` chaieb@31119 ` 1297` ``` in (mons, ``` chaieb@31119 ` 1298` ``` fold_rev (fn (m,n) => Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons Monomialfunc.undefined) ``` chaieb@31119 ` 1299` ``` end ``` chaieb@31119 ` 1300` chaieb@31119 ` 1301` ``` fun mk_sqmultiplier k (p,c) = ``` chaieb@31119 ` 1302` ``` let ``` chaieb@31119 ` 1303` ``` val e = (d - multidegree p) div 2 ``` chaieb@31119 ` 1304` ``` val mons = enumerate_monomials e vars ``` chaieb@31119 ` 1305` ``` val nons = mons ~~ (1 upto length mons) ``` chaieb@31119 ` 1306` ``` in (mons, ``` chaieb@31119 ` 1307` ``` fold_rev (fn (m1,n1) => ``` chaieb@31119 ` 1308` ``` fold_rev (fn (m2,n2) => fn a => ``` chaieb@31119 ` 1309` ``` let val m = monomial_mul m1 m2 ``` chaieb@31119 ` 1310` ``` in if n1 > n2 then a else ``` chaieb@31119 ` 1311` ``` let val c = if n1 = n2 then rat_1 else rat_2 ``` chaieb@31119 ` 1312` ``` val e = Monomialfunc.tryapplyd a m Inttriplefunc.undefined ``` chaieb@31119 ` 1313` ``` in Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a ``` chaieb@31119 ` 1314` ``` end ``` chaieb@31119 ` 1315` ``` end) nons) ``` chaieb@31119 ` 1316` ``` nons Monomialfunc.undefined) ``` chaieb@31119 ` 1317` ``` end ``` chaieb@31119 ` 1318` chaieb@31119 ` 1319` ``` val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid) ``` chaieb@31119 ` 1320` ``` val (idmonlist,ids) = split_list(map2 mk_idmultiplier (1 upto length eqs) eqs) ``` chaieb@31119 ` 1321` ``` val blocksizes = map length sqmonlist ``` chaieb@31119 ` 1322` ``` val bigsum = ``` chaieb@31119 ` 1323` ``` fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids ``` chaieb@31119 ` 1324` ``` (fold_rev2 (fn (p,c) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs ``` chaieb@31119 ` 1325` ``` (epoly_of_poly(poly_neg pol))) ``` chaieb@31119 ` 1326` ``` val eqns = Monomialfunc.fold (fn (m,e) => fn a => e::a) bigsum [] ``` chaieb@31119 ` 1327` ``` val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns ``` chaieb@31119 ` 1328` ``` val qvars = (0,0,0)::pvs ``` chaieb@31119 ` 1329` ``` val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig ``` chaieb@31119 ` 1330` ``` fun mk_matrix v = ``` chaieb@31119 ` 1331` ``` Inttriplefunc.fold (fn ((b,i,j), ass) => fn m => ``` chaieb@31119 ` 1332` ``` if b < 0 then m else ``` chaieb@31119 ` 1333` ``` let val c = Inttriplefunc.tryapplyd ass v rat_0 ``` chaieb@31119 ` 1334` ``` in if c = rat_0 then m else ``` chaieb@31119 ` 1335` ``` Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m) ``` chaieb@31119 ` 1336` ``` end) ``` chaieb@31119 ` 1337` ``` allassig Inttriplefunc.undefined ``` chaieb@31119 ` 1338` ``` val diagents = Inttriplefunc.fold ``` chaieb@31119 ` 1339` ``` (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a) ``` chaieb@31119 ` 1340` ``` allassig Inttriplefunc.undefined ``` chaieb@31119 ` 1341` chaieb@31119 ` 1342` ``` val mats = map mk_matrix qvars ``` chaieb@31119 ` 1343` ``` val obj = (length pvs, ``` chaieb@31119 ` 1344` ``` itern 1 pvs (fn v => fn i => Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0)) ``` chaieb@31119 ` 1345` ``` Intfunc.undefined) ``` chaieb@31119 ` 1346` ``` val raw_vec = if null pvs then vector_0 0 ``` chaieb@31119 ` 1347` ``` else tri_scale_then (csdp nblocks blocksizes) obj mats ``` chaieb@31119 ` 1348` ``` fun int_element (d,v) i = Intfunc.tryapplyd v i rat_0 ``` chaieb@31119 ` 1349` ``` fun cterm_element (d,v) i = Ctermfunc.tryapplyd v i rat_0 ``` chaieb@31119 ` 1350` chaieb@31119 ` 1351` ``` fun find_rounding d = ``` chaieb@31119 ` 1352` ``` let ``` chaieb@31119 ` 1353` ``` val _ = if !debugging ``` chaieb@31119 ` 1354` ``` then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n") ``` chaieb@31119 ` 1355` ``` else () ``` chaieb@31119 ` 1356` ``` val vec = nice_vector d raw_vec ``` chaieb@31119 ` 1357` ``` val blockmat = iter (1,dim vec) ``` chaieb@31119 ` 1358` ``` (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a) ``` chaieb@31119 ` 1359` ``` (bmatrix_neg (nth mats 0)) ``` chaieb@31119 ` 1360` ``` val allmats = blocks blocksizes blockmat ``` chaieb@31119 ` 1361` ``` in (vec,map diag allmats) ``` chaieb@31119 ` 1362` ``` end ``` chaieb@31119 ` 1363` ``` val (vec,ratdias) = ``` chaieb@31119 ` 1364` ``` if null pvs then find_rounding rat_1 ``` chaieb@31119 ` 1365` ``` else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @ ``` chaieb@31119 ` 1366` ``` map pow2 (5 upto 66)) ``` chaieb@31119 ` 1367` ``` val newassigs = ``` chaieb@31119 ` 1368` ``` fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k)) ``` chaieb@31119 ` 1369` ``` (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1)) ``` chaieb@31119 ` 1370` ``` val finalassigs = ``` chaieb@31119 ` 1371` ``` Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs ``` chaieb@31119 ` 1372` ``` fun poly_of_epoly p = ``` chaieb@31119 ` 1373` ``` Monomialfunc.fold (fn (v,e) => fn a => Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a) ``` chaieb@31119 ` 1374` ``` p Monomialfunc.undefined ``` chaieb@31119 ` 1375` ``` fun mk_sos mons = ``` chaieb@31119 ` 1376` ``` let fun mk_sq (c,m) = ``` chaieb@31119 ` 1377` ``` (c,fold_rev (fn k=> fn a => Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a) ``` chaieb@31119 ` 1378` ``` (1 upto length mons) Monomialfunc.undefined) ``` chaieb@31119 ` 1379` ``` in map mk_sq ``` chaieb@31119 ` 1380` ``` end ``` chaieb@31119 ` 1381` ``` val sqs = map2 mk_sos sqmonlist ratdias ``` chaieb@31119 ` 1382` ``` val cfs = map poly_of_epoly ids ``` chaieb@31119 ` 1383` ``` val msq = filter (fn (a,b) => not (null b)) (map2 pair monoid sqs) ``` chaieb@31119 ` 1384` ``` fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0 ``` chaieb@31119 ` 1385` ``` val sanity = ``` chaieb@31119 ` 1386` ``` fold_rev (fn ((p,c),s) => poly_add (poly_mul p (eval_sq s))) msq ``` chaieb@31119 ` 1387` ``` (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs ``` chaieb@31119 ` 1388` ``` (poly_neg pol)) ``` chaieb@31119 ` 1389` chaieb@31119 ` 1390` ```in if not(Monomialfunc.is_undefined sanity) then raise Sanity else ``` chaieb@31119 ` 1391` ``` (cfs,map (fn (a,b) => (snd a,b)) msq) ``` chaieb@31119 ` 1392` ``` end ``` chaieb@31119 ` 1393` chaieb@31119 ` 1394` chaieb@31119 ` 1395` ```end; ``` chaieb@31119 ` 1396` chaieb@31119 ` 1397` ```(* Iterative deepening. *) ``` chaieb@31119 ` 1398` chaieb@31119 ` 1399` ```fun deepen f n = ``` chaieb@31119 ` 1400` ``` (writeln ("Searching with depth limit " ^ string_of_int n) ; (f n handle ERROR s => (writeln ("failed with message: " ^ s) ; deepen f (n+1)))) ``` chaieb@31119 ` 1401` chaieb@31119 ` 1402` ```(* The ordering so we can create canonical HOL polynomials. *) ``` chaieb@31119 ` 1403` chaieb@31119 ` 1404` ```fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon); ``` chaieb@31119 ` 1405` chaieb@31119 ` 1406` ```fun monomial_order (m1,m2) = ``` chaieb@31119 ` 1407` ``` if Ctermfunc.is_undefined m2 then LESS ``` chaieb@31119 ` 1408` ``` else if Ctermfunc.is_undefined m1 then GREATER ``` chaieb@31119 ` 1409` ``` else ``` chaieb@31119 ` 1410` ``` let val mon1 = dest_monomial m1 ``` chaieb@31119 ` 1411` ``` val mon2 = dest_monomial m2 ``` chaieb@31119 ` 1412` ``` val deg1 = fold (curry op + o snd) mon1 0 ``` chaieb@31119 ` 1413` ``` val deg2 = fold (curry op + o snd) mon2 0 ``` chaieb@31119 ` 1414` ``` in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS ``` chaieb@31119 ` 1415` ``` else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2) ``` chaieb@31119 ` 1416` ``` end; ``` chaieb@31119 ` 1417` chaieb@31119 ` 1418` ```fun dest_poly p = ``` chaieb@31119 ` 1419` ``` map (fn (m,c) => (c,dest_monomial m)) ``` chaieb@31119 ` 1420` ``` (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p)); ``` chaieb@31119 ` 1421` chaieb@31119 ` 1422` ```(* Map back polynomials and their composites to HOL. *) ``` chaieb@31119 ` 1423` chaieb@31119 ` 1424` ```local ``` chaieb@31119 ` 1425` ``` open Thm Numeral RealArith ``` chaieb@31119 ` 1426` ```in ``` chaieb@31119 ` 1427` chaieb@31119 ` 1428` ```fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x) ``` chaieb@31119 ` 1429` ``` (mk_cnumber @{ctyp nat} k) ``` chaieb@31119 ` 1430` chaieb@31119 ` 1431` ```fun cterm_of_monomial m = ``` chaieb@31119 ` 1432` ``` if Ctermfunc.is_undefined m then @{cterm "1::real"} ``` chaieb@31119 ` 1433` ``` else ``` chaieb@31119 ` 1434` ``` let ``` chaieb@31119 ` 1435` ``` val m' = dest_monomial m ``` chaieb@31119 ` 1436` ``` val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' [] ``` chaieb@31119 ` 1437` ``` in fold1 (fn s => fn t => capply (capply @{cterm "op * :: real => _"} s) t) vps ``` chaieb@31119 ` 1438` ``` end ``` chaieb@31119 ` 1439` chaieb@31119 ` 1440` ```fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c ``` chaieb@31119 ` 1441` ``` else if c = Rat.one then cterm_of_monomial m ``` chaieb@31119 ` 1442` ``` else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m); ``` chaieb@31119 ` 1443` chaieb@31119 ` 1444` ```fun cterm_of_poly p = ``` chaieb@31119 ` 1445` ``` if Monomialfunc.is_undefined p then @{cterm "0::real"} ``` chaieb@31119 ` 1446` ``` else ``` chaieb@31119 ` 1447` ``` let ``` chaieb@31119 ` 1448` ``` val cms = map cterm_of_cmonomial ``` chaieb@31119 ` 1449` ``` (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p)) ``` chaieb@31119 ` 1450` ``` in fold1 (fn t1 => fn t2 => capply(capply @{cterm "op + :: real => _"} t1) t2) cms ``` chaieb@31119 ` 1451` ``` end; ``` chaieb@31119 ` 1452` chaieb@31119 ` 1453` ```fun cterm_of_sqterm (c,p) = Product(Rational_lt c,Square(cterm_of_poly p)); ``` chaieb@31119 ` 1454` chaieb@31119 ` 1455` ```fun cterm_of_sos (pr,sqs) = if null sqs then pr ``` chaieb@31119 ` 1456` ``` else Product(pr,fold1 (fn a => fn b => Sum(a,b)) (map cterm_of_sqterm sqs)); ``` chaieb@31119 ` 1457` chaieb@31119 ` 1458` ```end ``` chaieb@31119 ` 1459` chaieb@31119 ` 1460` ```(* Interface to HOL. *) ``` chaieb@31119 ` 1461` ```local ``` chaieb@31119 ` 1462` ``` open Thm Conv RealArith ``` chaieb@31119 ` 1463` ``` val concl = dest_arg o cprop_of ``` chaieb@31119 ` 1464` ``` fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS ``` chaieb@31119 ` 1465` ```in ``` chaieb@31119 ` 1466` ``` (* FIXME: Replace tryfind by get_first !! *) ``` chaieb@31119 ` 1467` ```fun real_nonlinear_prover ctxt = ``` chaieb@31119 ` 1468` ``` let ``` chaieb@31119 ` 1469` ``` val {add,mul,neg,pow,sub,main} = Normalizer.semiring_normalizers_ord_wrapper ctxt ``` chaieb@31119 ` 1470` ``` (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"})) ``` chaieb@31119 ` 1471` ``` simple_cterm_ord ``` chaieb@31119 ` 1472` ``` val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv, ``` chaieb@31119 ` 1473` ``` real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main) ``` chaieb@31119 ` 1474` ``` fun mainf translator (eqs,les,lts) = ``` chaieb@31119 ` 1475` ``` let ``` chaieb@31119 ` 1476` ``` val eq0 = map (poly_of_term o dest_arg1 o concl) eqs ``` chaieb@31119 ` 1477` ``` val le0 = map (poly_of_term o dest_arg o concl) les ``` chaieb@31119 ` 1478` ``` val lt0 = map (poly_of_term o dest_arg o concl) lts ``` chaieb@31119 ` 1479` ``` val eqp0 = map (fn (t,i) => (t,Axiom_eq i)) (eq0 ~~ (0 upto (length eq0 - 1))) ``` chaieb@31119 ` 1480` ``` val lep0 = map (fn (t,i) => (t,Axiom_le i)) (le0 ~~ (0 upto (length le0 - 1))) ``` chaieb@31119 ` 1481` ``` val ltp0 = map (fn (t,i) => (t,Axiom_lt i)) (lt0 ~~ (0 upto (length lt0 - 1))) ``` chaieb@31119 ` 1482` ``` val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0 ``` chaieb@31119 ` 1483` ``` val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0 ``` chaieb@31119 ` 1484` ``` val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0 ``` chaieb@31119 ` 1485` ``` fun trivial_axiom (p,ax) = ``` chaieb@31119 ` 1486` ``` case ax of ``` chaieb@31119 ` 1487` ``` Axiom_eq n => if eval Ctermfunc.undefined p <>/ Rat.zero then nth eqs n ``` chaieb@31119 ` 1488` ``` else error "trivial_axiom: Not a trivial axiom" ``` chaieb@31119 ` 1489` ``` | Axiom_le n => if eval Ctermfunc.undefined p if eval Ctermfunc.undefined p <=/ Rat.zero then nth lts n ``` chaieb@31119 ` 1492` ``` else error "trivial_axiom: Not a trivial axiom" ``` chaieb@31119 ` 1493` ``` | _ => error "trivial_axiom: Not a trivial axiom" ``` chaieb@31119 ` 1494` ``` in ``` chaieb@31119 ` 1495` ``` ((let val th = tryfind trivial_axiom (keq @ klep @ kltp) ``` chaieb@31119 ` 1496` ``` in fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv field_comp_conv) th end) ``` chaieb@31119 ` 1497` ``` handle ERROR _ => ( ``` chaieb@31119 ` 1498` ``` let ``` chaieb@31119 ` 1499` ``` val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one) ``` chaieb@31119 ` 1500` ``` val leq = lep @ ltp ``` chaieb@31119 ` 1501` ``` fun tryall d = ``` chaieb@31119 ` 1502` ``` let val e = multidegree pol ``` chaieb@31119 ` 1503` ``` val k = if e = 0 then 0 else d div e ``` chaieb@31119 ` 1504` ``` val eq' = map fst eq ``` chaieb@31119 ` 1505` ``` in tryfind (fn i => (d,i,real_positivnullstellensatz_general false d eq' leq ``` chaieb@31119 ` 1506` ``` (poly_neg(poly_pow pol i)))) ``` chaieb@31119 ` 1507` ``` (0 upto k) ``` chaieb@31119 ` 1508` ``` end ``` chaieb@31119 ` 1509` ``` val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0 ``` chaieb@31119 ` 1510` ``` val proofs_ideal = ``` chaieb@31119 ` 1511` ``` map2 (fn q => fn (p,ax) => Eqmul(cterm_of_poly q,ax)) cert_ideal eq ``` chaieb@31119 ` 1512` ``` val proofs_cone = map cterm_of_sos cert_cone ``` chaieb@31119 ` 1513` ``` val proof_ne = if null ltp then Rational_lt Rat.one else ``` chaieb@31119 ` 1514` ``` let val p = fold1 (fn s => fn t => Product(s,t)) (map snd ltp) ``` chaieb@31119 ` 1515` ``` in funpow i (fn q => Product(p,q)) (Rational_lt Rat.one) ``` chaieb@31119 ` 1516` ``` end ``` chaieb@31119 ` 1517` ``` val proof = fold1 (fn s => fn t => Sum(s,t)) ``` chaieb@31119 ` 1518` ``` (proof_ne :: proofs_ideal @ proofs_cone) ``` chaieb@31119 ` 1519` ``` in writeln "Translating proof certificate to HOL"; ``` chaieb@31119 ` 1520` ``` translator (eqs,les,lts) proof ``` chaieb@31119 ` 1521` ``` end)) ``` chaieb@31119 ` 1522` ``` end ``` chaieb@31119 ` 1523` ``` in mainf end ``` chaieb@31119 ` 1524` ```end ``` chaieb@31119 ` 1525` chaieb@31119 ` 1526` ```fun C f x y = f y x; ``` chaieb@31119 ` 1527` ``` (* FIXME : This is very bad!!!*) ``` chaieb@31119 ` 1528` ```fun subst_conv eqs t = ``` chaieb@31119 ` 1529` ``` let ``` chaieb@31119 ` 1530` ``` val t' = fold (Thm.cabs o Thm.lhs_of) eqs t ``` chaieb@31119 ` 1531` ``` in Conv.fconv_rule (Thm.beta_conversion true) (fold (C combination) eqs (reflexive t')) ``` chaieb@31119 ` 1532` ``` end ``` chaieb@31119 ` 1533` chaieb@31119 ` 1534` ```(* A wrapper that tries to substitute away variables first. *) ``` chaieb@31119 ` 1535` chaieb@31119 ` 1536` ```local ``` chaieb@31119 ` 1537` ``` open Thm Conv RealArith ``` chaieb@31119 ` 1538` ``` fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS ``` chaieb@31119 ` 1539` ``` val concl = dest_arg o cprop_of ``` chaieb@31119 ` 1540` ``` val shuffle1 = ``` chaieb@31119 ` 1541` ``` fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps) }) ``` chaieb@31119 ` 1542` ``` val shuffle2 = ``` chaieb@31119 ` 1543` ``` fconv_rule (rewr_conv @{lemma "(x + a == y) == (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps)}) ``` chaieb@31119 ` 1544` ``` fun substitutable_monomial fvs tm = case term_of tm of ``` chaieb@31119 ` 1545` ``` Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm) ``` chaieb@31119 ` 1546` ``` else error "substitutable_monomial" ``` chaieb@31119 ` 1547` ``` | @{term "op * :: real => _"}\$c\$(t as Free _ ) => ``` chaieb@31119 ` 1548` ``` if is_ratconst (dest_arg1 tm) andalso not (member (op aconvc) fvs (dest_arg tm)) ``` chaieb@31119 ` 1549` ``` then (dest_ratconst (dest_arg1 tm),dest_arg tm) else error "substitutable_monomial" ``` chaieb@31119 ` 1550` ``` | @{term "op + :: real => _"}\$s\$t => ``` chaieb@31119 ` 1551` ``` (substitutable_monomial (add_cterm_frees (dest_arg tm) fvs) (dest_arg1 tm) ``` chaieb@31119 ` 1552` ``` handle ERROR _ => substitutable_monomial (add_cterm_frees (dest_arg1 tm) fvs) (dest_arg tm)) ``` chaieb@31119 ` 1553` ``` | _ => error "substitutable_monomial" ``` chaieb@31119 ` 1554` chaieb@31119 ` 1555` ``` fun isolate_variable v th = ``` chaieb@31119 ` 1556` ``` let val w = dest_arg1 (cprop_of th) ``` chaieb@31119 ` 1557` ``` in if v aconvc w then th ``` chaieb@31119 ` 1558` ``` else case term_of w of ``` chaieb@31119 ` 1559` ``` @{term "op + :: real => _"}\$s\$t => ``` chaieb@31119 ` 1560` ``` if dest_arg1 w aconvc v then shuffle2 th ``` chaieb@31119 ` 1561` ``` else isolate_variable v (shuffle1 th) ``` chaieb@31119 ` 1562` ``` | _ => error "isolate variable : This should not happen?" ``` chaieb@31119 ` 1563` ``` end ``` chaieb@31119 ` 1564` ```in ``` chaieb@31119 ` 1565` chaieb@31119 ` 1566` ```fun real_nonlinear_subst_prover ctxt = ``` chaieb@31119 ` 1567` ``` let ``` chaieb@31119 ` 1568` ``` val {add,mul,neg,pow,sub,main} = Normalizer.semiring_normalizers_ord_wrapper ctxt ``` chaieb@31119 ` 1569` ``` (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"})) ``` chaieb@31119 ` 1570` ``` simple_cterm_ord ``` chaieb@31119 ` 1571` chaieb@31119 ` 1572` ``` val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv, ``` chaieb@31119 ` 1573` ``` real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main) ``` chaieb@31119 ` 1574` chaieb@31119 ` 1575` ``` fun make_substitution th = ``` chaieb@31119 ` 1576` ``` let ``` chaieb@31119 ` 1577` ``` val (c,v) = substitutable_monomial [] (dest_arg1(concl th)) ``` chaieb@31119 ` 1578` ``` val th1 = Drule.arg_cong_rule (capply @{cterm "op * :: real => _"} (cterm_of_rat (Rat.inv c))) (mk_meta_eq th) ``` chaieb@31119 ` 1579` ``` val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1 ``` chaieb@31119 ` 1580` ``` in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2) ``` chaieb@31119 ` 1581` ``` end ``` chaieb@31119 ` 1582` ``` fun oprconv cv ct = ``` chaieb@31119 ` 1583` ``` let val g = Thm.dest_fun2 ct ``` chaieb@31119 ` 1584` ``` in if g aconvc @{cterm "op <= :: real => _"} ``` chaieb@31119 ` 1585` ``` orelse g aconvc @{cterm "op < :: real => _"} ``` chaieb@31119 ` 1586` ``` then arg_conv cv ct else arg1_conv cv ct ``` chaieb@31119 ` 1587` ``` end ``` chaieb@31119 ` 1588` ``` fun mainf translator = ``` chaieb@31119 ` 1589` ``` let ``` chaieb@31119 ` 1590` ``` fun substfirst(eqs,les,lts) = ``` chaieb@31119 ` 1591` ``` ((let ``` chaieb@31119 ` 1592` ``` val eth = tryfind make_substitution eqs ``` chaieb@31119 ` 1593` ``` val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv))) ``` chaieb@31119 ` 1594` ``` in substfirst ``` chaieb@31119 ` 1595` ``` (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t ``` chaieb@31119 ` 1596` ``` aconvc @{cterm "0::real"}) (map modify eqs), ``` chaieb@31119 ` 1597` ``` map modify les,map modify lts) ``` chaieb@31119 ` 1598` ``` end) ``` chaieb@31119 ` 1599` ``` handle ERROR _ => real_nonlinear_prover ctxt translator (rev eqs, rev les, rev lts)) ``` chaieb@31119 ` 1600` ``` in substfirst ``` chaieb@31119 ` 1601` ``` end ``` chaieb@31119 ` 1602` chaieb@31119 ` 1603` chaieb@31119 ` 1604` ``` in mainf ``` chaieb@31119 ` 1605` ``` end ``` chaieb@31119 ` 1606` chaieb@31119 ` 1607` ```(* Overall function. *) ``` chaieb@31119 ` 1608` chaieb@31119 ` 1609` ```fun real_sos ctxt t = gen_prover_real_arith ctxt (real_nonlinear_subst_prover ctxt) t; ``` chaieb@31119 ` 1610` ```end; ``` chaieb@31119 ` 1611` chaieb@31131 ` 1612` ```(* A tactic *) ``` chaieb@31131 ` 1613` ```fun strip_all ct = ``` chaieb@31131 ` 1614` ``` case term_of ct of ``` chaieb@31131 ` 1615` ``` Const("all",_) \$ Abs (xn,xT,p) => ``` chaieb@31131 ` 1616` ``` let val (a,(v,t')) = (apsnd (Thm.dest_abs (SOME xn)) o Thm.dest_comb) ct ``` chaieb@31131 ` 1617` ``` in apfst (cons v) (strip_all t') ``` chaieb@31131 ` 1618` ``` end ``` chaieb@31131 ` 1619` ```| _ => ([],ct) ``` chaieb@31131 ` 1620` chaieb@31131 ` 1621` ```fun core_sos_conv ctxt t = Drule.arg_cong_rule @{cterm Trueprop} (real_sos ctxt (Thm.dest_arg t) RS @{thm Eq_TrueI}) ``` chaieb@31512 ` 1622` chaieb@31512 ` 1623` ```val known_sos_constants = ``` chaieb@31512 ` 1624` ``` [@{term "op ==>"}, @{term "Trueprop"}, ``` chaieb@31512 ` 1625` ``` @{term "op -->"}, @{term "op &"}, @{term "op |"}, ``` chaieb@31512 ` 1626` ``` @{term "Not"}, @{term "op = :: bool => _"}, ``` chaieb@31512 ` 1627` ``` @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"}, ``` chaieb@31512 ` 1628` ``` @{term "op = :: real => _"}, @{term "op < :: real => _"}, ``` chaieb@31512 ` 1629` ``` @{term "op <= :: real => _"}, ``` chaieb@31512 ` 1630` ``` @{term "op + :: real => _"}, @{term "op - :: real => _"}, ``` chaieb@31512 ` 1631` ``` @{term "op * :: real => _"}, @{term "uminus :: real => _"}, ``` chaieb@31512 ` 1632` ``` @{term "op / :: real => _"}, @{term "inverse :: real => _"}, ``` chaieb@31512 ` 1633` ``` @{term "op ^ :: real => _"}, @{term "abs :: real => _"}, ``` chaieb@31512 ` 1634` ``` @{term "min :: real => _"}, @{term "max :: real => _"}, ``` chaieb@31512 ` 1635` ``` @{term "0::real"}, @{term "1::real"}, @{term "number_of :: int => real"}, ``` chaieb@31512 ` 1636` ``` @{term "number_of :: int => nat"}, ``` chaieb@31512 ` 1637` ``` @{term "Int.Bit0"}, @{term "Int.Bit1"}, ``` chaieb@31512 ` 1638` ``` @{term "Int.Pls"}, @{term "Int.Min"}]; ``` chaieb@31512 ` 1639` chaieb@31512 ` 1640` ```fun check_sos kcts ct = ``` chaieb@31512 ` 1641` ``` let ``` chaieb@31512 ` 1642` ``` val t = term_of ct ``` chaieb@31512 ` 1643` ``` val _ = if not (null (Term.add_tfrees t []) ``` chaieb@31512 ` 1644` ``` andalso null (Term.add_tvars t [])) ``` chaieb@31512 ` 1645` ``` then error "SOS: not sos. Additional type varables" else () ``` chaieb@31512 ` 1646` ``` val fs = Term.add_frees t [] ``` chaieb@31512 ` 1647` ``` val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs ``` chaieb@31512 ` 1648` ``` then error "SOS: not sos. Variables with type not real" else () ``` chaieb@31512 ` 1649` ``` val vs = Term.add_vars t [] ``` chaieb@31512 ` 1650` ``` val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs ``` chaieb@31512 ` 1651` ``` then error "SOS: not sos. Variables with type not real" else () ``` chaieb@31512 ` 1652` ``` val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t []) ``` chaieb@31512 ` 1653` ``` val _ = if null ukcs then () ``` chaieb@31512 ` 1654` ``` else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs)) ``` chaieb@31512 ` 1655` ```in () end ``` chaieb@31512 ` 1656` chaieb@31131 ` 1657` ```fun core_sos_tac ctxt = CSUBGOAL (fn (ct, i) => ``` chaieb@31512 ` 1658` ``` let val _ = check_sos known_sos_constants ct ``` chaieb@31512 ` 1659` ``` val (avs, p) = strip_all ct ``` chaieb@31131 ` 1660` ``` val th = standard (fold_rev forall_intr avs (real_sos ctxt (Thm.dest_arg p))) ``` chaieb@31131 ` 1661` ``` in rtac th i end); ``` chaieb@31131 ` 1662` chaieb@31131 ` 1663` ```fun default_SOME f NONE v = SOME v ``` chaieb@31131 ` 1664` ``` | default_SOME f (SOME v) _ = SOME v; ``` chaieb@31131 ` 1665` chaieb@31131 ` 1666` ```fun lift_SOME f NONE a = f a ``` chaieb@31131 ` 1667` ``` | lift_SOME f (SOME a) _ = SOME a; ``` chaieb@31131 ` 1668` chaieb@31131 ` 1669` chaieb@31131 ` 1670` ```local ``` chaieb@31131 ` 1671` ``` val is_numeral = can (HOLogic.dest_number o term_of) ``` chaieb@31131 ` 1672` ```in ``` chaieb@31131 ` 1673` ```fun get_denom b ct = case term_of ct of ``` chaieb@31131 ` 1674` ``` @{term "op / :: real => _"} \$ _ \$ _ => ``` chaieb@31131 ` 1675` ``` if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct) ``` chaieb@31131 ` 1676` ``` else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct)) (Thm.dest_arg ct, b) ``` chaieb@31131 ` 1677` ``` | @{term "op < :: real => _"} \$ _ \$ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct) ``` chaieb@31131 ` 1678` ``` | @{term "op <= :: real => _"} \$ _ \$ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct) ``` chaieb@31131 ` 1679` ``` | _ \$ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct) ``` chaieb@31131 ` 1680` ``` | _ => NONE ``` chaieb@31131 ` 1681` ```end; ``` chaieb@31131 ` 1682` chaieb@31131 ` 1683` ```fun elim_one_denom_tac ctxt = ``` chaieb@31131 ` 1684` ```CSUBGOAL (fn (P,i) => ``` chaieb@31131 ` 1685` ``` case get_denom false P of ``` chaieb@31131 ` 1686` ``` NONE => no_tac ``` chaieb@31131 ` 1687` ``` | SOME (d,ord) => ``` chaieb@31131 ` 1688` ``` let ``` chaieb@31131 ` 1689` ``` val ss = simpset_of (ProofContext.theory_of ctxt) addsimps @{thms field_simps} ``` chaieb@31131 ` 1690` ``` addsimps [@{thm nonzero_power_divide}, @{thm power_divide}] ``` chaieb@31131 ` 1691` ``` val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)] ``` chaieb@31131 ` 1692` ``` (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto} ``` chaieb@31131 ` 1693` ``` else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast}) ``` chaieb@31131 ` 1694` ``` in (rtac th i THEN Simplifier.asm_full_simp_tac ss i) end); ``` chaieb@31131 ` 1695` chaieb@31131 ` 1696` ```fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i); ``` chaieb@31131 ` 1697` chaieb@31131 ` 1698` ```fun sos_tac ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac ctxt ``` chaieb@31131 ` 1699` chaieb@31131 ` 1700` chaieb@31512 ` 1701` ```end; ```