Examples/sos.ml

   1 (* ========================================================================= *)
   2 (* Nonlinear universal reals procedure using SOS decomposition.              *)
   3 (* ========================================================================= *)
   4
   5 prioritize_real();;
   6
   7 let debugging = ref false;;
   8
   9 exception Sanity;;
  10
  11 exception Unsolvable;;
  12
  13 (* ------------------------------------------------------------------------- *)
  14 (* Turn a rational into a decimal string with d sig digits.                  *)
  15 (* ------------------------------------------------------------------------- *)
  16
  17 let decimalize =
  18   let rec normalize y =
  19     if abs_num y </ Int 1 // Int 10 then normalize (Int 10 */ y) - 1
  20     else if abs_num y >=/ Int 1 then normalize (y // Int 10) + 1
  21     else 0 in
  22   fun d x ->
  23     if x =/ Int 0 then "0.0" else
  24     let y = abs_num x in
  25     let e = normalize y in
  26     let z = pow10(-e) */ y +/ Int 1 in
  27     let k = round_num(pow10 d */ z) in
  28     (if x </ Int 0 then "-0." else "0.") ^
  29     implode(tl(explode(string_of_num k))) ^
  30     (if e = 0 then "" else "e"^string_of_int e);;
  31
  32 (* ------------------------------------------------------------------------- *)
  33 (* Iterations over numbers, and lists indexed by numbers.                    *)
  34 (* ------------------------------------------------------------------------- *)
  35
  36 let rec itern k l f a =
  37   match l with
  38     [] -> a
  39   | h::t -> itern (k + 1) t f (f h k a);;
  40
  41 let rec iter (m,n) f a =
  42   if n < m then a
  43   else iter (m+1,n) f (f m a);;
  44
  45 (* ------------------------------------------------------------------------- *)
  46 (* The main types.                                                           *)
  47 (* ------------------------------------------------------------------------- *)
  48
  49 type vector = int*(int,num)func;;
  50
  51 type matrix = (int*int)*(int*int,num)func;;
  52
  53 type monomial = (term,int)func;;
  54
  55 type poly = (monomial,num)func;;
  56
  57 (* ------------------------------------------------------------------------- *)
  58 (* Assignment avoiding zeros.                                                *)
  59 (* ------------------------------------------------------------------------- *)
  60
  61 let (|-->) x y a = if y =/ Int 0 then a else (x |-> y) a;;
  62
  63 (* ------------------------------------------------------------------------- *)
  64 (* This can be generic.                                                      *)
  65 (* ------------------------------------------------------------------------- *)
  66
  67 let element (d,v) i = tryapplyd v i (Int 0);;
  68
  69 let mapa f (d,v) =
  70   d,foldl (fun a i c -> (i |--> f(c)) a) undefined v;;
  71
  72 let is_zero (d,v) = is_undefined v;;
  73
  74 (* ------------------------------------------------------------------------- *)
  75 (* Vectors. Conventionally indexed 1..n.                                     *)
  76 (* ------------------------------------------------------------------------- *)
  77
  78 let vec_0 n = (n,undefined:vector);;
  79
  80 let vec_dim (v:vector) = fst v;;
  81
  82 let vec_const c n =
  83   if c =/ Int 0 then vec_0 n
  84   else (n,itlist (fun k -> k |-> c) (1--n) undefined :vector);;
  85
  86 let vec_1 = vec_const (Int 1);;
  87
  88 let vec_cmul c (v:vector) =
  89   let n = vec_dim v in
  90   if c =/ Int 0 then vec_0 n
  91   else n,mapf (fun x -> c */ x) (snd v)
  92
  93 let vec_neg (v:vector) = (fst v,mapf minus_num (snd v) :vector);;
  94
  95 let vec_add (v1:vector) (v2:vector) =
  96   let m = vec_dim v1 and n = vec_dim v2 in
  97   if m <> n then failwith "vec_add: incompatible dimensions" else
  98   (n,combine (+/) (fun x -> x =/ Int 0) (snd v1) (snd v2) :vector);;
  99
 100 let vec_sub v1 v2 = vec_add v1 (vec_neg v2);;
 101
 102 let vec_dot (v1:vector) (v2:vector) =
 103   let m = vec_dim v1 and n = vec_dim v2 in
 104   if m <> n then failwith "vec_add: incompatible dimensions" else
 105   foldl (fun a i x -> x +/ a) (Int 0)
 106         (combine ( */ ) (fun x -> x =/ Int 0) (snd v1) (snd v2));;
 107
 108 let vec_of_list l =
 109   let n = length l in
 110   (n,itlist2 (|->) (1--n) l undefined :vector);;
 111
 112 (* ------------------------------------------------------------------------- *)
 113 (* Matrices; again rows and columns indexed from 1.                          *)
 114 (* ------------------------------------------------------------------------- *)
 115
 116 let matrix_0 (m,n) = ((m,n),undefined:matrix);;
 117
 118 let dimensions (m:matrix) = fst m;;
 119
 120 let matrix_const c (m,n as mn) =
 121   if m <> n then failwith "matrix_const: needs to be square"
 122   else if c =/ Int 0 then matrix_0 mn
 123   else (mn,itlist (fun k -> (k,k) |-> c) (1--n) undefined :matrix);;
 124
 125 let matrix_1 = matrix_const (Int 1);;
 126
 127 let matrix_cmul c (m:matrix) =
 128   let (i,j) = dimensions m in
 129   if c =/ Int 0 then matrix_0 (i,j)
 130   else (i,j),mapf (fun x -> c */ x) (snd m);;
 131
 132 let matrix_neg (m:matrix) = (dimensions m,mapf minus_num (snd m) :matrix);;
 133
 134 let matrix_add (m1:matrix) (m2:matrix) =
 135   let d1 = dimensions m1 and d2 = dimensions m2 in
 136   if d1 <> d2 then failwith "matrix_add: incompatible dimensions"
 137   else (d1,combine (+/) (fun x -> x =/ Int 0) (snd m1) (snd m2) :matrix);;
 138
 139 let matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2);;
 140
 141 let row k (m:matrix) =
 142   let i,j = dimensions m in
 143   (j,
 144    foldl (fun a (i,j) c -> if i = k then (j |-> c) a else a) undefined (snd m)
 145    : vector);;
 146
 147 let column k (m:matrix) =
 148   let i,j = dimensions m in
 149   (i,
 150    foldl (fun a (i,j) c -> if j = k then (i |-> c) a else a) undefined (snd m)
 151    : vector);;
 152
 153 let transp (m:matrix) =
 154   let i,j = dimensions m in
 155   ((j,i),foldl (fun a (i,j) c -> ((j,i) |-> c) a) undefined (snd m) :matrix);;
 156
 157 let diagonal (v:vector) =
 158   let n = vec_dim v in
 159   ((n,n),foldl (fun a i c -> ((i,i) |-> c) a) undefined (snd v) : matrix);;
 160
 161 let matrix_of_list l =
 162   let m = length l in
 163   if m = 0 then matrix_0 (0,0) else
 164   let n = length (hd l) in
 165   (m,n),itern 1 l (fun v i -> itern 1 v (fun c j -> (i,j) |-> c)) undefined;;
 166
 167 (* ------------------------------------------------------------------------- *)
 168 (* Monomials.                                                                *)
 169 (* ------------------------------------------------------------------------- *)
 170
 171 let monomial_eval assig (m:monomial) =
 172   foldl (fun a x k -> a */ power_num (apply assig x) (Int k))
 173         (Int 1) m;;
 174
 175 let monomial_1 = (undefined:monomial);;
 176
 177 let monomial_var x = (x |=> 1 :monomial);;
 178
 179 let (monomial_mul:monomial->monomial->monomial) =
 180   combine (+) (fun x -> false);;
 181
 182 let monomial_pow (m:monomial) k =
 183   if k = 0 then monomial_1
 184   else mapf (fun x -> k * x) m;;
 185
 186 let monomial_divides (m1:monomial) (m2:monomial) =
 187   foldl (fun a x k -> tryapplyd m2 x 0 >= k & a) true m1;;
 188
 189 let monomial_div (m1:monomial) (m2:monomial) =
 190   let m = combine (+) (fun x -> x = 0) m1 (mapf (fun x -> -x) m2) in
 191   if foldl (fun a x k -> k >= 0 & a) true m then m
 192   else failwith "monomial_div: non-divisible";;
 193
 194 let monomial_degree x (m:monomial) = tryapplyd m x 0;;
 195
 196 let monomial_lcm (m1:monomial) (m2:monomial) =
 197   (itlist (fun x -> x |-> max (monomial_degree x m1) (monomial_degree x m2))
 198           (union (dom m1) (dom m2)) undefined :monomial);;
 199
 200 let monomial_multidegree (m:monomial) = foldl (fun a x k -> k + a) 0 m;;
 201
 202 let monomial_variables m = dom m;;
 203
 204 (* ------------------------------------------------------------------------- *)
 205 (* Polynomials.                                                              *)
 206 (* ------------------------------------------------------------------------- *)
 207
 208 let eval assig (p:poly) =
 209   foldl (fun a m c -> a +/ c */ monomial_eval assig m) (Int 0) p;;
 210
 211 let poly_0 = (undefined:poly);;
 212
 213 let poly_isconst (p:poly) = foldl (fun a m c -> m = monomial_1 & a) true p;;
 214
 215 let poly_var x = ((monomial_var x) |=> Int 1 :poly);;
 216
 217 let poly_const c =
 218   if c =/ Int 0 then poly_0 else (monomial_1 |=> c);;
 219
 220 let poly_cmul c (p:poly) =
 221   if c =/ Int 0 then poly_0
 222   else mapf (fun x -> c */ x) p;;
 223
 224 let poly_neg (p:poly) = (mapf minus_num p :poly);;
 225
 226 let poly_add (p1:poly) (p2:poly) =
 227   (combine (+/) (fun x -> x =/ Int 0) p1 p2 :poly);;
 228
 229 let poly_sub p1 p2 = poly_add p1 (poly_neg p2);;
 230
 231 let poly_cmmul (c,m) (p:poly) =
 232   if c =/ Int 0 then poly_0
 233   else if m = monomial_1 then mapf (fun d -> c */ d) p
 234   else foldl (fun a m' d -> (monomial_mul m m' |-> c */ d) a) poly_0 p;;
 235
 236 let poly_mul (p1:poly) (p2:poly) =
 237   foldl (fun a m c -> poly_add (poly_cmmul (c,m) p2) a) poly_0 p1;;
 238
 239 let poly_div (p1:poly) (p2:poly) =
 240   if not(poly_isconst p2) then failwith "poly_div: non-constant" else
 241   let c = eval undefined p2 in
 242   if c =/ Int 0 then failwith "poly_div: division by zero"
 243   else poly_cmul (Int 1 // c) p1;;
 244
 245 let poly_square p = poly_mul p p;;
 246
 247 let rec poly_pow p k =
 248   if k = 0 then poly_const (Int 1)
 249   else if k = 1 then p
 250   else let q = poly_square(poly_pow p (k / 2)) in
 251        if k mod 2 = 1 then poly_mul p q else q;;
 252
 253 let poly_exp p1 p2 =
 254   if not(poly_isconst p2) then failwith "poly_exp: not a constant" else
 255   poly_pow p1 (Num.int_of_num (eval undefined p2));;
 256
 257 let degree x (p:poly) = foldl (fun a m c -> max (monomial_degree x m) a) 0 p;;
 258
 259 let multidegree (p:poly) =
 260   foldl (fun a m c -> max (monomial_multidegree m) a) 0 p;;
 261
 262 let poly_variables (p:poly) =
 263   foldr (fun m c -> union (monomial_variables m)) p [];;
 264
 265 (* ------------------------------------------------------------------------- *)
 266 (* Order monomials for human presentation.                                   *)
 267 (* ------------------------------------------------------------------------- *)
 268
 269 let humanorder_varpow (x1,k1) (x2,k2) = x1 < x2 or x1 = x2 & k1 > k2;;
 270
 271 let humanorder_monomial =
 272   let rec ord l1 l2 = match (l1,l2) with
 273     _,[] -> true
 274   | [],_ -> false
 275   | h1::t1,h2::t2 -> humanorder_varpow h1 h2 or h1 = h2 & ord t1 t2 in
 276   fun m1 m2 -> m1 = m2 or
 277                ord (sort humanorder_varpow (graph m1))
 278                    (sort humanorder_varpow (graph m2));;
 279
 280 (* ------------------------------------------------------------------------- *)
 281 (* Conversions to strings.                                                   *)
 282 (* ------------------------------------------------------------------------- *)
 283
 284 let string_of_vector min_size max_size (v:vector) =
 285   let n_raw = vec_dim v in
 286   if n_raw = 0 then "[]" else
 287   let n = max min_size (min n_raw max_size) in
 288   let xs = map (string_of_num o element v) (1--n) in
 289   "[" ^ end_itlist (fun s t -> s ^ ", " ^ t) xs ^
 290   (if n_raw > max_size then ", ...]" else "]");;
 291
 292 let string_of_matrix max_size (m:matrix) =
 293   let i_raw,j_raw = dimensions m in
 294   let i = min max_size i_raw and j = min max_size j_raw in
 295   let rstr = map (fun k -> string_of_vector j j (row k m)) (1--i) in
 296   "["^end_itlist(fun s t -> s^";\n "^t) rstr ^
 297   (if j > max_size then "\n ...]" else "]");;
 298
 299 let rec string_of_term t =
 300         if (is_comb t) then
 301                 let (a,b) = (dest_comb t) in
 302                 "("^(string_of_term a)^" "^(string_of_term b)^")"
 303         else if (is_abs t) then
 304                 let (a,b) = (dest_abs t) in
 305                 "(\\"^(string_of_term a)^"."^(string_of_term b)^")"
 306         else if (is_const t) then
 307                 let (a,_) = (dest_const t) in a
 308         else if (is_var t) then
 309                 let (a,_) = (dest_var t) in a
 310         else failwith "string_of_term";;
 311
 312 let string_of_varpow x k =
 313   if k = 1 then string_of_term x else string_of_term x^"^"^string_of_int k;;
 314
 315 let string_of_monomial m =
 316   if m = monomial_1 then "1" else
 317   let vps = List.fold_right (fun (x,k) a -> string_of_varpow x k :: a)
 318                             (sort humanorder_varpow (graph m)) [] in
 319   end_itlist (fun s t -> s^"*"^t) vps;;
 320
 321 let string_of_cmonomial (c,m) =
 322   if m = monomial_1 then string_of_num c
 323   else if c =/ Int 1 then string_of_monomial m
 324   else string_of_num c ^ "*" ^ string_of_monomial m;;
 325
 326 let string_of_poly (p:poly) =
 327   if p = poly_0 then "<<0>>" else
 328   let cms = sort (fun (m1,_) (m2,_) -> humanorder_monomial m1 m2) (graph p) in
 329   let s =
 330     List.fold_left (fun a (m,c) ->
 331              if c </ Int 0 then a ^ " - " ^ string_of_cmonomial(minus_num c,m)
 332              else a ^ " + " ^ string_of_cmonomial(c,m))
 333           "" cms in
 334   let s1 = String.sub s 0 3
 335   and s2 = String.sub s 3 (String.length s - 3) in
 336   "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>";;
 337
 338 (* ------------------------------------------------------------------------- *)
 339 (* Printers.                                                                 *)
 340 (* ------------------------------------------------------------------------- *)
 341
 342 let print_vector v = Format.print_string(string_of_vector 0 20 v);;
 343
 344 let print_matrix m = Format.print_string(string_of_matrix 20 m);;
 345
 346 let print_monomial m = Format.print_string(string_of_monomial m);;
 347
 348 let print_poly m = Format.print_string(string_of_poly m);;
 349
 350 #install_printer print_vector;;
 351 #install_printer print_matrix;;
 352 #install_printer print_monomial;;
 353 #install_printer print_poly;;
 354
 355 (* ------------------------------------------------------------------------- *)
 356 (* Conversion from HOL term.                                                 *)
 357 (* ------------------------------------------------------------------------- *)
 358
 359 let poly_of_term =
 360   let neg_tm = `(--):real->real`
 361   and add_tm = `(+):real->real->real`
 362   and sub_tm = `(-):real->real->real`
 363   and mul_tm = `(*):real->real->real`
 364   and inv_tm = `(inv):real->real`
 365   and div_tm = `(/):real->real->real`
 366   and pow_tm = `(pow):real->num->real`
 367   and zero_tm = `&0:real`
 368   and real_ty = `:real` in
 369   let rec poly_of_term tm =
 370     if tm = zero_tm then poly_0
 371     else if is_ratconst tm then poly_const(rat_of_term tm)
 372     else if not(is_comb tm) then poly_var tm else
 373     let lop,r = dest_comb tm in
 374     if lop = neg_tm then poly_neg(poly_of_term r)
 375     else if lop = inv_tm then
 376       let p = poly_of_term r in
 377       if poly_isconst p then poly_const(Int 1 // eval undefined p)
 378       else failwith "poly_of_term: inverse of non-constant polyomial"
 379     else if not(is_comb lop) then poly_var tm else
 380     let op,l = dest_comb lop in
 381     if op = pow_tm & is_numeral r then
 382         poly_pow (poly_of_term l) (dest_small_numeral r)
 383     else if op = add_tm then poly_add (poly_of_term l) (poly_of_term r)
 384     else if op = sub_tm then poly_sub (poly_of_term l) (poly_of_term r)
 385     else if op = mul_tm then poly_mul (poly_of_term l) (poly_of_term r)
 386     else if op = div_tm then
 387       let p = poly_of_term l and q = poly_of_term r in
 388       if poly_isconst q then poly_cmul (Int 1 // eval undefined q) p
 389       else failwith "poly_of_term: division by non-constant polynomial"
 390     else poly_var tm in
 391   fun tm -> if type_of tm = real_ty then poly_of_term tm
 392             else failwith "poly_of_term: term does not have real type";;
 393
 394 (* ------------------------------------------------------------------------- *)
 395 (* String of vector (just a list of space-separated numbers).                *)
 396 (* ------------------------------------------------------------------------- *)
 397
 398 let sdpa_of_vector (v:vector) =
 399   let n = vec_dim v in
 400   let strs = map (decimalize 20 o element v) (1--n) in
 401   end_itlist (fun x y -> x ^ " " ^ y) strs ^ "\n";;
 402
 403 (* ------------------------------------------------------------------------- *)
 404 (* String for block diagonal matrix numbered k.                              *)
 405 (* ------------------------------------------------------------------------- *)
 406
 407 let sdpa_of_blockdiagonal k m =
 408   let pfx = string_of_int k ^" " in
 409   let ents =
 410     foldl (fun a (b,i,j) c -> if i > j then a else ((b,i,j),c)::a) [] m in
 411   let entss = sort (increasing fst) ents in
 412   itlist (fun ((b,i,j),c) a ->
 413      pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
 414      " " ^ decimalize 20 c ^ "\n" ^ a) entss "";;
 415
 416 (* ------------------------------------------------------------------------- *)
 417 (* String for a matrix numbered k, in SDPA sparse format.                    *)
 418 (* ------------------------------------------------------------------------- *)
 419
 420 let sdpa_of_matrix k (m:matrix) =
 421   let pfx = string_of_int k ^ " 1 " in
 422   let ms = foldr (fun (i,j) c a -> if i > j then a else ((i,j),c)::a)
 423                  (snd m) [] in
 424   let mss = sort (increasing fst) ms in
 425   itlist (fun ((i,j),c) a ->
 426      pfx ^ string_of_int i ^ " " ^ string_of_int j ^
 427      " " ^ decimalize 20 c ^ "\n" ^ a) mss "";;
 428
 429 (* ------------------------------------------------------------------------- *)
 430 (* String in SDPA sparse format for standard SDP problem:                    *)
 431 (*                                                                           *)
 432 (*    X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD                *)
 433 (*    Minimize obj_1 * v_1 + ... obj_m * v_m                                 *)
 434 (* ------------------------------------------------------------------------- *)
 435
 436 let sdpa_of_problem comment obj mats =
 437   let m = length mats - 1
 438   and n,_ = dimensions (hd mats) in
 439   "\"" ^ comment ^ "\"\n" ^
 440   string_of_int m ^ "\n" ^
 441   "1\n" ^
 442   string_of_int n ^ "\n" ^
 443   sdpa_of_vector obj ^
 444   itlist2 (fun k m a -> sdpa_of_matrix (k - 1) m ^ a)
 445           (1--length mats) mats "";;
 446
 447 (* ------------------------------------------------------------------------- *)
 448 (* More parser basics.                                                       *)
 449 (* ------------------------------------------------------------------------- *)
 450
 451 let word s =
 452    end_itlist (fun p1 p2 -> (p1 ++ p2) >> (fun (s,t) -> s^t))
 453               (map a (explode s));;
 454 let token s =
 455   many (some isspace) ++ word s ++ many (some isspace)
 456   >> (fun ((_,t),_) -> t);;
 457
 458 let decimal =
 459   let numeral = some isnum in
 460   let decimalint = atleast 1 numeral >> (Num.num_of_string o implode) in
 461   let decimalfrac = atleast 1 numeral
 462     >> (fun s -> Num.num_of_string(implode s) // pow10 (length s)) in
 463   let decimalsig =
 464     decimalint ++ possibly (a "." ++ decimalfrac >> snd)
 465     >> (function (h,[]) -> h | (h,[x]) -> h +/ x) in
 466   let signed prs =
 467        a "-" ++ prs >> (minus_num o snd)
 468     || a "+" ++ prs >> snd
 469     || prs in
 470   let exponent = (a "e" || a "E") ++ signed decimalint >> snd in
 471     signed decimalsig ++ possibly exponent
 472     >> (function (h,[]) -> h | (h,[x]) -> h */ power_num (Int 10) x);;
 473
 474 let mkparser p s =
 475   let x,rst = p(explode s) in
 476   if rst = [] then x else failwith "mkparser: unparsed input";;
 477
 478 let parse_decimal = mkparser decimal;;
 479
 480 (* ------------------------------------------------------------------------- *)
 481 (* Parse back a vector.                                                      *)
 482 (* ------------------------------------------------------------------------- *)
 483
 484 let parse_sdpaoutput,parse_csdpoutput =
 485   let vector =
 486     token "{" ++ listof decimal (token ",") "decimal" ++ token "}"
 487                >> (fun ((_,v),_) -> vec_of_list v) in
 488   let parse_vector = mkparser vector in
 489   let rec skipupto dscr prs inp =
 490       (dscr ++ prs >> snd
 491     || some (fun c -> true) ++ skipupto dscr prs >> snd) inp in
 492   let ignore inp = (),[] in
 493   let sdpaoutput =
 494     skipupto (word "xVec" ++ token "=")
 495              (vector ++ ignore >> fst) in
 496   let csdpoutput =
 497     (decimal ++ many(a " " ++ decimal >> snd) >> (fun (h,t) -> h::t)) ++
 498     (a " " ++ a "\n" ++ ignore) >> (vec_of_list o fst) in
 499   mkparser sdpaoutput,mkparser csdpoutput;;
 500
 501 (* ------------------------------------------------------------------------- *)
 502 (* Also parse the SDPA output to test success (CSDP yields a return code).   *)
 503 (* ------------------------------------------------------------------------- *)
 504
 505 let sdpa_run_succeeded =
 506    let rec skipupto dscr prs inp =
 507       (dscr ++ prs >> snd
 508     || some (fun c -> true) ++ skipupto dscr prs >> snd) inp in
 509   let prs = skipupto (word "phase.value" ++ token "=")
 510                      (possibly (a "p") ++ possibly (a "d") ++
 511                       (word "OPT" || word "FEAS")) in
 512   fun s -> try prs (explode s); true with Noparse -> false;;
 513
 514 (* ------------------------------------------------------------------------- *)
 515 (* The default parameters. Unfortunately this goes to a fixed file.          *)
 516 (* ------------------------------------------------------------------------- *)
 517
 518 let sdpa_default_parameters =
 519 "100     unsigned int maxIteration;
 520 1.0E-7  double 0.0 < epsilonStar;
 521 1.0E2   double 0.0 < lambdaStar;
 522 2.0     double 1.0 < omegaStar;
 523 -1.0E5  double lowerBound;
 524 1.0E5   double upperBound;
 525 0.1     double 0.0 <= betaStar <  1.0;
 526 0.2     double 0.0 <= betaBar  <  1.0, betaStar <= betaBar;
 527 0.9     double 0.0 < gammaStar  <  1.0;
 528 1.0E-7  double 0.0 < epsilonDash;
 529 ";;
 530
 531 (* ------------------------------------------------------------------------- *)
 532 (* These were suggested by Makoto Yamashita for problems where we are        *)
 533 (* right at the edge of the semidefinite cone, as sometimes happens.         *)
 534 (* ------------------------------------------------------------------------- *)
 535
 536 let sdpa_alt_parameters =
 537 "1000    unsigned int maxIteration;
 538 1.0E-7  double 0.0 < epsilonStar;
 539 1.0E4   double 0.0 < lambdaStar;
 540 2.0     double 1.0 < omegaStar;
 541 -1.0E5  double lowerBound;
 542 1.0E5   double upperBound;
 543 0.1     double 0.0 <= betaStar <  1.0;
 544 0.2     double 0.0 <= betaBar  <  1.0, betaStar <= betaBar;
 545 0.9     double 0.0 < gammaStar  <  1.0;
 546 1.0E-7  double 0.0 < epsilonDash;
 547 ";;
 548
 549 let sdpa_params = sdpa_alt_parameters;;
 550
 551 (* ------------------------------------------------------------------------- *)
 552 (* CSDP parameters; so far I'm sticking with the defaults.                   *)
 553 (* ------------------------------------------------------------------------- *)
 554
 555 let csdp_default_parameters =
 556 "axtol=1.0e-8
 557 atytol=1.0e-8
 558 objtol=1.0e-8
 559 pinftol=1.0e8
 560 dinftol=1.0e8
 561 maxiter=100
 562 minstepfrac=0.9
 563 maxstepfrac=0.97
 564 minstepp=1.0e-8
 565 minstepd=1.0e-8
 566 usexzgap=1
 567 tweakgap=0
 568 affine=0
 569 printlevel=1
 570 ";;
 571
 572 let csdp_params = csdp_default_parameters;;
 573
 574 (* ------------------------------------------------------------------------- *)
 575 (* Now call SDPA on a problem and parse back the output.                     *)
 576 (* ------------------------------------------------------------------------- *)
 577
 578 let run_sdpa dbg obj mats =
 579   let input_file = Filename.temp_file "sos" ".dat-s" in
 580   let output_file =
 581     String.sub input_file 0 (String.length input_file - 6) ^ ".out"
 582   and params_file = Filename.concat (!temp_path) "param.sdpa" in
 583   file_of_string input_file (sdpa_of_problem "" obj mats);
 584   file_of_string params_file sdpa_params;
 585   Sys.command("cd "^ !temp_path ^
 586               "; sdpa "^input_file ^ " " ^ output_file ^
 587               (if dbg then "" else "> /dev/null"));
 588   let op = string_of_file output_file in
 589   if not(sdpa_run_succeeded op) then failwith "sdpa: call failed" else
 590   let res = parse_sdpaoutput op in
 591   ((if dbg then ()
 592    else (Sys.remove input_file; Sys.remove output_file));
 593    res);;
 594
 595 let sdpa obj mats = run_sdpa (!debugging) obj mats;;
 596
 597 (* ------------------------------------------------------------------------- *)
 598 (* The same thing with CSDP.                                                 *)
 599 (* ------------------------------------------------------------------------- *)
 600
 601 let run_csdp dbg obj mats =
 602   let input_file = Filename.temp_file "sos" ".dat-s" in
 603   let output_file =
 604     String.sub input_file 0 (String.length input_file - 6) ^ ".out"
 605   and params_file = Filename.concat (!temp_path) "param.csdp" in
 606   file_of_string input_file (sdpa_of_problem "" obj mats);
 607   file_of_string params_file csdp_params;
 608   let rv = Sys.command("cd "^(!temp_path)^"; csdp "^input_file ^
 609                         " " ^ output_file ^
 610                        (if dbg then "" else "> /dev/null")) in
 611   let op = string_of_file output_file in
 612   let res = parse_csdpoutput op in
 613   ((if dbg then ()
 614     else (Sys.remove input_file; Sys.remove output_file));
 615     rv,res);;
 616
 617 let csdp obj mats =
 618   let rv,res = run_csdp (!debugging) obj mats in
 619   (if rv = 1 or rv = 2 then failwith "csdp: Problem is infeasible"
 620    else if rv = 3 then
 621     (Format.print_string "csdp warning: Reduced accuracy";
 622      Format.print_newline())
 623    else if rv <> 0 then failwith("csdp: error "^string_of_int rv)
 624    else ());
 625   res;;
 626
 627 (* ------------------------------------------------------------------------- *)
 628 (* Try some apparently sensible scaling first. Note that this is purely to   *)
 629 (* get a cleaner translation to floating-point, and doesn't affect any of    *)
 630 (* the results, in principle. In practice it seems a lot better when there   *)
 631 (* are extreme numbers in the original problem.                              *)
 632 (* ------------------------------------------------------------------------- *)
 633
 634 let scale_then =
 635   let common_denominator amat acc =
 636     foldl (fun a m c -> lcm_num (denominator c) a) acc amat
 637   and maximal_element amat acc =
 638     foldl (fun maxa m c -> max_num maxa (abs_num c)) acc amat in
 639   fun solver obj mats ->
 640     let cd1 = itlist common_denominator mats (Int 1)
 641     and cd2 = common_denominator (snd obj)  (Int 1) in
 642     let mats' = map (mapf (fun x -> cd1 */ x)) mats
 643     and obj' = vec_cmul cd2 obj in
 644     let max1 = itlist maximal_element mats' (Int 0)
 645     and max2 = maximal_element (snd obj') (Int 0) in
 646     let scal1 = pow2 (20-int_of_float(log(float_of_num max1) /. log 2.0))
 647     and scal2 = pow2 (20-int_of_float(log(float_of_num max2) /. log 2.0)) in
 648     let mats'' = map (mapf (fun x -> x */ scal1)) mats'
 649     and obj'' = vec_cmul scal2 obj' in
 650     solver obj'' mats'';;
 651
 652 (* ------------------------------------------------------------------------- *)
 653 (* Round a vector to "nice" rationals.                                       *)
 654 (* ------------------------------------------------------------------------- *)
 655
 656 let nice_rational n x = round_num (n */ x) // n;;
 657
 658 let nice_vector n = mapa (nice_rational n);;
 659
 660 (* ------------------------------------------------------------------------- *)
 661 (* Reduce linear program to SDP (diagonal matrices) and test with CSDP. This *)
 662 (* one tests A [-1;x1;..;xn] >= 0 (i.e. left column is negated constants).   *)
 663 (* ------------------------------------------------------------------------- *)
 664
 665 let linear_program_basic a =
 666   let m,n = dimensions a in
 667   let mats =  map (fun j -> diagonal (column j a)) (1--n)
 668   and obj = vec_const (Int 1) m in
 669   let rv,res = run_csdp false obj mats in
 670   if rv = 1 or rv = 2 then false
 671   else if rv = 0 then true
 672   else failwith "linear_program: An error occurred in the SDP solver";;
 673
 674 (* ------------------------------------------------------------------------- *)
 675 (* Alternative interface testing A x >= b for matrix A, vector b.            *)
 676 (* ------------------------------------------------------------------------- *)
 677
 678 let linear_program a b =
 679   let m,n = dimensions a in
 680   if vec_dim b <> m then failwith "linear_program: incompatible dimensions" else
 681   let mats = diagonal b :: map (fun j -> diagonal (column j a)) (1--n)
 682   and obj = vec_const (Int 1) m in
 683   let rv,res = run_csdp false obj mats in
 684   if rv = 1 or rv = 2 then false
 685   else if rv = 0 then true
 686   else failwith "linear_program: An error occurred in the SDP solver";;
 687
 688 (* ------------------------------------------------------------------------- *)
 689 (* Test whether a point is in the convex hull of others. Rather than use     *)
 690 (* computational geometry, express as linear inequalities and call CSDP.     *)
 691 (* This is a bit lazy of me, but it's easy and not such a bottleneck so far. *)
 692 (* ------------------------------------------------------------------------- *)
 693
 694 let in_convex_hull pts pt =
 695   let pts1 = (1::pt) :: map (fun x -> 1::x) pts in
 696   let pts2 = map (fun p -> map (fun x -> -x) p @ p) pts1 in
 697   let n = length pts + 1
 698   and v = 2 * (length pt + 1) in
 699   let m = v + n - 1 in
 700   let mat =
 701     (m,n),
 702     itern 1 pts2 (fun pts j -> itern 1 pts (fun x i -> (i,j) |-> Int x))
 703                  (iter (1,n) (fun i -> (v + i,i+1) |-> Int 1) undefined) in
 704   linear_program_basic mat;;
 705
 706 (* ------------------------------------------------------------------------- *)
 707 (* Filter down a set of points to a minimal set with the same convex hull.   *)
 708 (* ------------------------------------------------------------------------- *)
 709
 710 let minimal_convex_hull =
 711   let augment1 (m::ms) = if in_convex_hull ms m then ms else ms@[m] in
 712   let augment m ms = funpow 3 augment1 (m::ms) in
 713   fun mons ->
 714     let mons' = itlist augment (tl mons) [hd mons] in
 715     funpow (length mons') augment1 mons';;
 716
 717 (* ------------------------------------------------------------------------- *)
 718 (* Stuff for "equations" (generic A->num functions).                         *)
 719 (* ------------------------------------------------------------------------- *)
 720
 721 let equation_cmul c eq =
 722   if c =/ Int 0 then undefined else mapf (fun d -> c */ d) eq;;
 723
 724 let equation_add eq1 eq2 = combine (+/) (fun x -> x =/ Int 0) eq1 eq2;;
 725
 726 let equation_eval assig eq =
 727   let value v = apply assig v in
 728   foldl (fun a v c -> a +/ value(v) */ c) (Int 0) eq;;
 729
 730 (* ------------------------------------------------------------------------- *)
 731 (* Eliminate among linear equations: return unconstrained variables and      *)
 732 (* assignments for the others in terms of them. We give one pseudo-variable  *)
 733 (* "one" that's used for a constant term.                                    *)
 734 (* ------------------------------------------------------------------------- *)
 735
 736
 737 let eliminate_equations =
 738   let rec extract_first p l =
 739     match l with
 740       [] -> failwith "extract_first"
 741     | h::t -> if p(h) then h,t else
 742               let k,s = extract_first p t in
 743               k,h::s in
 744   let rec eliminate vars dun eqs =
 745     match vars with
 746       [] -> if forall is_undefined eqs then dun
 747 else raise Unsolvable
 748     | v::vs ->
 749             try let eq,oeqs = extract_first (fun e -> defined e v) eqs in
 750                 let a = apply eq v in
 751                 let eq' = equation_cmul (Int(-1) // a) (undefine v eq) in
 752                 let elim e =
 753                   let b = tryapplyd e v (Int 0) in
 754                   if b =/ Int 0 then e else
 755                   equation_add e (equation_cmul (minus_num b // a) eq) in
 756                 eliminate vs ((v |-> eq') (mapf elim dun)) (map elim oeqs)
 757             with Failure _ -> eliminate vs dun eqs in
 758   fun one vars eqs ->
 759     let assig = eliminate vars undefined eqs in
 760     let vs = foldl (fun a x f -> subtract (dom f) [one] @ a) [] assig in
 761     setify vs,assig;;
 762
 763 (* ------------------------------------------------------------------------- *)
 764 (* Eliminate all variables, in an essentially arbitrary order.               *)
 765 (* ------------------------------------------------------------------------- *)
 766
 767 let eliminate_all_equations one =
 768   let choose_variable eq =
 769     let (v,_) = choose eq in
 770     if v = one then
 771       let eq' = undefine v eq in
 772       if is_undefined eq' then failwith "choose_variable" else
 773       let (w,_) = choose eq' in w
 774     else v in
 775   let rec eliminate dun eqs =
 776     match eqs with
 777       [] -> dun
 778     | eq::oeqs ->
 779         if is_undefined eq then eliminate dun oeqs else
 780         let v = choose_variable eq in
 781         let a = apply eq v in
 782         let eq' = equation_cmul (Int(-1) // a) (undefine v eq) in
 783         let elim e =
 784           let b = tryapplyd e v (Int 0) in
 785           if b =/ Int 0 then e else
 786           equation_add e (equation_cmul (minus_num b // a) eq) in
 787         eliminate ((v |-> eq') (mapf elim dun)) (map elim oeqs) in
 788   fun eqs ->
 789     let assig = eliminate undefined eqs in
 790     let vs = foldl (fun a x f -> subtract (dom f) [one] @ a) [] assig in
 791     setify vs,assig;;
 792
 793 (* ------------------------------------------------------------------------- *)
 794 (* Solve equations by assigning arbitrary numbers.                           *)
 795 (* ------------------------------------------------------------------------- *)
 796
 797 let solve_equations one eqs =
 798   let vars,assigs = eliminate_all_equations one eqs in
 799   let vfn = itlist (fun v -> (v |-> Int 0)) vars (one |=> Int(-1)) in
 800   let ass =
 801     combine (+/) (fun c -> false) (mapf (equation_eval vfn) assigs) vfn in
 802   if forall (fun e -> equation_eval ass e =/ Int 0) eqs
 803   then undefine one ass else raise Sanity;;
 804
 805 (* ------------------------------------------------------------------------- *)
 806 (* Hence produce the "relevant" monomials: those whose squares lie in the    *)
 807 (* Newton polytope of the monomials in the input. (This is enough according  *)
 808 (* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal,       *)
 809 (* vol 45, pp. 363--374, 1978.                                               *)
 810 (*                                                                           *)
 811 (* These are ordered in sort of decreasing degree. In particular the         *)
 812 (* constant monomial is last; this gives an order in diagonalization of the  *)
 813 (* quadratic form that will tend to display constants.                       *)
 814 (* ------------------------------------------------------------------------- *)
 815
 816 let newton_polytope pol =
 817   let vars = poly_variables pol in
 818   let mons = map (fun m -> map (fun x -> monomial_degree x m) vars) (dom pol)
 819   and ds = map (fun x -> (degree x pol + 1) / 2) vars in
 820   let all = itlist (fun n -> allpairs (fun h t -> h::t) (0--n)) ds [[]]
 821   and mons' = minimal_convex_hull mons in
 822   let all' =
 823     filter (fun m -> in_convex_hull mons' (map (fun x -> 2 * x) m)) all in
 824   map (fun m -> itlist2 (fun v i a -> if i = 0 then a else (v |-> i) a)
 825                         vars m monomial_1) (rev all');;
 826
 827 (* ------------------------------------------------------------------------- *)
 828 (* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form.  *)
 829 (* ------------------------------------------------------------------------- *)
 830
 831 let diag m =
 832   let nn = dimensions m in
 833   let n = fst nn in
 834   if snd nn <> n then failwith "diagonalize: non-square matrix" else
 835   let rec diagonalize i m =
 836     if is_zero m then [] else
 837     let a11 = element m (i,i) in
 838     if a11 </ Int 0 then failwith "diagonalize: not PSD"
 839     else if a11 =/ Int 0 then
 840       if is_zero(row i m) then diagonalize (i + 1) m
 841       else failwith "diagonalize: not PSD"
 842     else
 843       let v = row i m in
 844       let v' = mapa (fun a1k -> a1k // a11) v in
 845       let m' =
 846       (n,n),
 847       iter (i+1,n) (fun j ->
 848           iter (i+1,n) (fun k ->
 849               ((j,k) |--> (element m (j,k) -/ element v j */ element v' k))))
 850           undefined in
 851       (a11,v')::diagonalize (i + 1) m' in
 852   diagonalize 1 m;;
 853
 854 (* ------------------------------------------------------------------------- *)
 855 (* Adjust a diagonalization to collect rationals at the start.               *)
 856 (* ------------------------------------------------------------------------- *)
 857
 858 let deration d =
 859   if d = [] then Int 0,d else
 860   let adj(c,l) =
 861     let a = foldl (fun a i c -> lcm_num a (denominator c)) (Int 1) (snd l) //
 862             foldl (fun a i c -> gcd_num a (numerator c)) (Int 0) (snd l) in
 863     (c // (a */ a)),mapa (fun x -> a */ x) l in
 864   let d' = map adj d in
 865   let a = itlist (lcm_num o denominator o fst) d' (Int 1) //
 866           itlist (gcd_num o numerator o fst) d' (Int 0)  in
 867   (Int 1 // a),map (fun (c,l) -> (a */ c,l)) d';;
 868
 869 (* ------------------------------------------------------------------------- *)
 870 (* Enumeration of monomials with given multidegree bound.                    *)
 871 (* ------------------------------------------------------------------------- *)
 872
 873 let rec enumerate_monomials d vars =
 874   if d < 0 then []
 875   else if d = 0 then [undefined]
 876   else if vars = [] then [monomial_1] else
 877   let alts =
 878     map (fun k -> let oths = enumerate_monomials (d - k) (tl vars) in
 879                   map (fun ks -> if k = 0 then ks else (hd vars |-> k) ks) oths)
 880         (0--d) in
 881   end_itlist (@) alts;;
 882
 883 (* ------------------------------------------------------------------------- *)
 884 (* Enumerate products of distinct input polys with degree <= d.              *)
 885 (* We ignore any constant input polynomials.                                 *)
 886 (* Give the output polynomial and a record of how it was derived.            *)
 887 (* ------------------------------------------------------------------------- *)
 888
 889 let rec enumerate_products d pols =
 890   if d = 0 then [poly_const num_1,Rational_lt num_1] else if d < 0 then [] else
 891   match pols with
 892     [] -> [poly_const num_1,Rational_lt num_1]
 893   | (p,b)::ps -> let e = multidegree p in
 894                  if e = 0 then enumerate_products d ps else
 895                  enumerate_products d ps @
 896                  map (fun (q,c) -> poly_mul p q,Product(b,c))
 897                      (enumerate_products (d - e) ps);;
 898
 899 (* ------------------------------------------------------------------------- *)
 900 (* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
 901 (* ------------------------------------------------------------------------- *)
 902
 903 let epoly_pmul p q acc =
 904   foldl (fun a m1 c ->
 905            foldl (fun b m2 e ->
 906                         let m =  monomial_mul m1 m2 in
 907                         let es = tryapplyd b m undefined in
 908                         (m |-> equation_add (equation_cmul c e) es) b)
 909                  a q) acc p;;
 910
 911 (* ------------------------------------------------------------------------- *)
 912 (* Usual operations on equation-parametrized poly.                           *)
 913 (* ------------------------------------------------------------------------- *)
 914
 915 let epoly_cmul c l =
 916   if c =/ Int 0 then undefined else mapf (equation_cmul c) l;;
 917
 918
 919
 920
 921 (* ------------------------------------------------------------------------- *)
 922 (* Convert regular polynomial. Note that we treat (0,0,0) as -1.             *)
 923 (* ------------------------------------------------------------------------- *)
 924
 925 let epoly_of_poly p =
 926   foldl (fun a m c -> (m |-> ((0,0,0) |=> minus_num c)) a) undefined p;;
 927
 928 (* ------------------------------------------------------------------------- *)
 929 (* String for block diagonal matrix numbered k.                              *)
 930 (* ------------------------------------------------------------------------- *)
 931
 932 let sdpa_of_blockdiagonal k m =
 933   let pfx = string_of_int k ^" " in
 934   let ents =
 935     foldl (fun a (b,i,j) c -> if i > j then a else ((b,i,j),c)::a) [] m in
 936   let entss = sort (increasing fst) ents in
 937   itlist (fun ((b,i,j),c) a ->
 938      pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
 939      " " ^ decimalize 20 c ^ "\n" ^ a) entss "";;
 940
 941 (* ------------------------------------------------------------------------- *)
 942 (* SDPA for problem using block diagonal (i.e. multiple SDPs)                *)
 943 (* ------------------------------------------------------------------------- *)
 944
 945 let sdpa_of_blockproblem comment nblocks blocksizes obj mats =
 946   let m = length mats - 1 in
 947   "\"" ^ comment ^ "\"\n" ^
 948   string_of_int m ^ "\n" ^
 949   string_of_int nblocks ^ "\n" ^
 950   (end_itlist (fun s t -> s^" "^t) (map string_of_int blocksizes)) ^
 951   "\n" ^
 952   sdpa_of_vector obj ^
 953   itlist2 (fun k m a -> sdpa_of_blockdiagonal (k - 1) m ^ a)
 954           (1--length mats) mats "";;
 955
 956 (* ------------------------------------------------------------------------- *)
 957 (* Hence run CSDP on a problem in block diagonal form.                       *)
 958 (* ------------------------------------------------------------------------- *)
 959
 960 let run_csdp dbg nblocks blocksizes obj mats =
 961   let input_file = Filename.temp_file "sos" ".dat-s" in
 962   let output_file =
 963     String.sub input_file 0 (String.length input_file - 6) ^ ".out"
 964   and params_file = Filename.concat (!temp_path) "param.csdp" in
 965   file_of_string input_file
 966    (sdpa_of_blockproblem "" nblocks blocksizes obj mats);
 967   file_of_string params_file csdp_params;
 968   let rv = Sys.command("cd "^(!temp_path)^"; csdp "^input_file ^
 969                         " " ^ output_file ^
 970                        (if dbg then "" else "> /dev/null")) in
 971   let op = string_of_file output_file in
 972   let res = parse_csdpoutput op in
 973   ((if dbg then ()
 974     else (Sys.remove input_file; Sys.remove output_file));
 975     rv,res);;
 976
 977 let csdp nblocks blocksizes obj mats =
 978   let rv,res = run_csdp (!debugging) nblocks blocksizes obj mats in
 979   (if rv = 1 or rv = 2 then failwith "csdp: Problem is infeasible"
 980    else if rv = 3 then
 981     (Format.print_string "csdp warning: Reduced accuracy";
 982      Format.print_newline())
 983    else if rv <> 0 then failwith("csdp: error "^string_of_int rv)
 984    else ());
 985   res;;
 986
 987 (* ------------------------------------------------------------------------- *)
 988 (* 3D versions of matrix operations to consider blocks separately.           *)
 989 (* ------------------------------------------------------------------------- *)
 990
 991 let bmatrix_add = combine (+/) (fun x -> x =/ Int 0);;
 992
 993 let bmatrix_cmul c bm =
 994   if c =/ Int 0 then undefined
 995   else mapf (fun x -> c */ x) bm;;
 996
 997 let bmatrix_neg = bmatrix_cmul (Int(-1));;
 998
 999 let bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);;
1000
1001 (* ------------------------------------------------------------------------- *)
1002 (* Smash a block matrix into components.                                     *)
1003 (* ------------------------------------------------------------------------- *)
1004
1005 let blocks blocksizes bm =
1006   map (fun (bs,b0) ->
1007         let m = foldl
1008           (fun a (b,i,j) c -> if b = b0 then ((i,j) |-> c) a else a)
1009           undefined bm in
1010         let d = foldl (fun a (i,j) c -> max a (max i j)) 0 m in
1011         (((bs,bs),m):matrix))
1012       (zip blocksizes (1--length blocksizes));;
1013
1014 (* ------------------------------------------------------------------------- *)
1015 (* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
1016 (* ------------------------------------------------------------------------- *)
1017
1018 let real_positivnullstellensatz_general linf d eqs leqs pol =
1019   let vars = itlist (union o poly_variables) (pol::eqs @ map fst leqs) [] in
1020   let monoid =
1021     if linf then
1022       (poly_const num_1,Rational_lt num_1)::
1023       (filter (fun (p,c) -> multidegree p <= d) leqs)
1024     else enumerate_products d leqs in
1025   let nblocks = length monoid in
1026   let mk_idmultiplier k p =
1027     let e = d - multidegree p in
1028     let mons = enumerate_monomials e vars in
1029     let nons = zip mons (1--length mons) in
1030     mons,
1031     itlist (fun (m,n) -> (m |-> ((-k,-n,n) |=> Int 1))) nons undefined in
1032   let mk_sqmultiplier k (p,c) =
1033     let e = (d - multidegree p) / 2 in
1034     let mons = enumerate_monomials e vars in
1035     let nons = zip mons (1--length mons) in
1036     mons,
1037     itlist (fun (m1,n1) ->
1038       itlist (fun (m2,n2) a ->
1039           let m = monomial_mul m1 m2 in
1040           if n1 > n2 then a else
1041           let c = if n1 = n2 then Int 1 else Int 2 in
1042           let e = tryapplyd a m undefined in
1043           (m |-> equation_add ((k,n1,n2) |=> c) e) a)
1044          nons)
1045        nons undefined in
1046   let sqmonlist,sqs = unzip(map2 mk_sqmultiplier (1--length monoid) monoid)
1047   and idmonlist,ids =  unzip(map2 mk_idmultiplier (1--length eqs) eqs) in
1048   let blocksizes = map length sqmonlist in
1049   let bigsum =
1050     itlist2 (fun p q a -> epoly_pmul p q a) eqs ids
1051             (itlist2 (fun (p,c) s a -> epoly_pmul p s a) monoid sqs
1052                      (epoly_of_poly(poly_neg pol))) in
1053   let eqns = foldl (fun a m e -> e::a) [] bigsum in
1054   let pvs,assig = eliminate_all_equations (0,0,0) eqns in
1055   let qvars = (0,0,0)::pvs in
1056   let allassig = itlist (fun v -> (v |-> (v |=> Int 1))) pvs assig in
1057   let mk_matrix v =
1058     foldl (fun m (b,i,j) ass -> if b < 0 then m else
1059                                 let c = tryapplyd ass v (Int 0) in
1060                                 if c =/ Int 0 then m else
1061                                 ((b,j,i) |-> c) (((b,i,j) |-> c) m))
1062           undefined allassig in
1063   let diagents = foldl
1064     (fun a (b,i,j) e -> if b > 0 & i = j then equation_add e a else a)
1065     undefined allassig in
1066   let mats = map mk_matrix qvars
1067   and obj = length pvs,
1068             itern 1 pvs (fun v i -> (i |--> tryapplyd diagents v (Int 0)))
1069                         undefined in
1070   let raw_vec = if pvs = [] then vec_0 0
1071                 else scale_then (csdp nblocks blocksizes) obj mats in
1072   let find_rounding d =
1073    (if !debugging then
1074      (Format.print_string("Trying rounding with limit "^string_of_num d);
1075       Format.print_newline())
1076     else ());
1077     let vec = nice_vector d raw_vec in
1078     let blockmat = iter (1,vec_dim vec)
1079      (fun i a -> bmatrix_add (bmatrix_cmul (element vec i) (el i mats)) a)
1080      (bmatrix_neg (el 0 mats)) in
1081     let allmats = blocks blocksizes blockmat in
1082     vec,map diag allmats in
1083   let vec,ratdias =
1084     if pvs = [] then find_rounding num_1
1085     else tryfind find_rounding (map Num.num_of_int (1--31) @
1086                                 map pow2 (5--66)) in
1087   let newassigs =
1088     itlist (fun k -> el (k - 1) pvs |-> element vec k)
1089            (1--vec_dim vec) ((0,0,0) |=> Int(-1)) in
1090   let finalassigs =
1091     foldl (fun a v e -> (v |-> equation_eval newassigs e) a) newassigs
1092           allassig in
1093   let poly_of_epoly p =
1094     foldl (fun a v e -> (v |--> equation_eval finalassigs e) a)
1095           undefined p in
1096   let mk_sos mons =
1097     let mk_sq (c,m) =
1098         c,itlist (fun k a -> (el (k - 1) mons |--> element m k) a)
1099                  (1--length mons) undefined in
1100     map mk_sq in
1101   let sqs = map2 mk_sos sqmonlist ratdias
1102   and cfs = map poly_of_epoly ids in
1103   let msq = filter (fun (a,b) -> b <> []) (map2 (fun a b -> a,b) monoid sqs) in
1104   let eval_sq sqs = itlist
1105    (fun (c,q) -> poly_add (poly_cmul c (poly_mul q q))) sqs poly_0 in
1106   let sanity =
1107     itlist (fun ((p,c),s) -> poly_add (poly_mul p (eval_sq s))) msq
1108            (itlist2 (fun p q -> poly_add (poly_mul p q)) cfs eqs
1109                     (poly_neg pol)) in
1110   if not(is_undefined sanity) then raise Sanity else
1111   cfs,map (fun (a,b) -> snd a,b) msq;;
1112
1113 (* ------------------------------------------------------------------------- *)
1114 (* Iterative deepening.                                                      *)
1115 (* ------------------------------------------------------------------------- *)
1116
1117 let rec deepen f n =
1118   try print_string "Searching with depth limit ";
1119       print_int n; print_newline(); f n
1120   with Failure _ -> deepen f (n + 1);;
1121
1122 (* ------------------------------------------------------------------------- *)
1123 (* The ordering so we can create canonical HOL polynomials.                  *)
1124 (* ------------------------------------------------------------------------- *)
1125
1126 let dest_monomial mon = sort (increasing fst) (graph mon);;
1127
1128 let monomial_order =
1129   let rec lexorder l1 l2 =
1130     match (l1,l2) with
1131       [],[] -> true
1132     | vps,[] -> false
1133     | [],vps -> true
1134     | ((x1,n1)::vs1),((x2,n2)::vs2) ->
1135           if x1 < x2 then true
1136           else if x2 < x1 then false
1137           else if n1 < n2 then false
1138           else if n2 < n1 then true
1139           else lexorder vs1 vs2 in
1140   fun m1 m2 ->
1141     if m2 = monomial_1 then true else if m1 = monomial_1 then false else
1142     let mon1 = dest_monomial m1 and mon2 = dest_monomial m2 in
1143     let deg1 = itlist ((+) o snd) mon1 0
1144     and deg2 = itlist ((+) o snd) mon2 0 in
1145     if deg1 < deg2 then false else if deg1 > deg2 then true
1146     else lexorder mon1 mon2;;
1147
1148 let dest_poly p =
1149   map (fun (m,c) -> c,dest_monomial m)
1150       (sort (fun (m1,_) (m2,_) -> monomial_order m1 m2) (graph p));;
1151
1152 (* ------------------------------------------------------------------------- *)
1153 (* Map back polynomials and their composites to HOL.                         *)
1154 (* ------------------------------------------------------------------------- *)
1155
1156 let term_of_varpow =
1157   let pow_tm = `(pow):real->num->real` in
1158   fun x k ->
1159     if k = 1 then x else mk_comb(mk_comb(pow_tm,x),mk_small_numeral k);;
1160
1161 let term_of_monomial =
1162   let one_tm = `&1:real`
1163   and mul_tm = `(*):real->real->real` in
1164   fun m -> if m = monomial_1 then one_tm else
1165            let m' = dest_monomial m in
1166            let vps = itlist (fun (x,k) a -> term_of_varpow x k :: a) m' [] in
1167            end_itlist (fun s t -> mk_comb(mk_comb(mul_tm,s),t)) vps;;
1168
1169 let term_of_cmonomial =
1170   let mul_tm = `(*):real->real->real` in
1171   fun (m,c) ->
1172     if m = monomial_1 then term_of_rat c
1173     else if c =/ num_1 then term_of_monomial m
1174     else mk_comb(mk_comb(mul_tm,term_of_rat c),term_of_monomial m);;
1175
1176 let term_of_poly =
1177   let zero_tm = `&0:real`
1178   and add_tm = `(+):real->real->real` in
1179   fun p ->
1180     if p = poly_0 then zero_tm else
1181     let cms = map term_of_cmonomial
1182      (sort (fun (m1,_) (m2,_) -> monomial_order m1 m2) (graph p)) in
1183     end_itlist (fun t1 t2 -> mk_comb(mk_comb(add_tm,t1),t2)) cms;;
1184
1185 let term_of_sqterm (c,p) =
1186   Product(Rational_lt c,Square(term_of_poly p));;
1187
1188 let term_of_sos (pr,sqs) =
1189   if sqs = [] then pr
1190   else Product(pr,end_itlist (fun a b -> Sum(a,b)) (map term_of_sqterm sqs));;
1191
1192 (* ------------------------------------------------------------------------- *)
1193 (* Interface to HOL.                                                         *)
1194 (* ------------------------------------------------------------------------- *)
1195
1196 let REAL_NONLINEAR_PROVER translator (eqs,les,lts) =
1197   let eq0 = map (poly_of_term o lhand o concl) eqs
1198   and le0 = map (poly_of_term o lhand o concl) les
1199   and lt0 = map (poly_of_term o lhand o concl) lts in
1200   let eqp0 = map (fun (t,i) -> t,Axiom_eq i) (zip eq0 (0--(length eq0 - 1)))
1201   and lep0 = map (fun (t,i) -> t,Axiom_le i) (zip le0 (0--(length le0 - 1)))
1202   and ltp0 = map (fun (t,i) -> t,Axiom_lt i) (zip lt0 (0--(length lt0 - 1))) in
1203   let keq,eq = partition (fun (p,_) -> multidegree p = 0) eqp0
1204   and klep,lep = partition (fun (p,_) -> multidegree p = 0) lep0
1205   and kltp,ltp = partition (fun (p,_) -> multidegree p = 0) ltp0 in
1206   let trivial_axiom (p,ax) =
1207     match ax with
1208       Axiom_eq n when eval undefined p <>/ num_0 -> el n eqs
1209     | Axiom_le n when eval undefined p </ num_0 -> el n les
1210     | Axiom_lt n when eval undefined p <=/ num_0 -> el n lts
1211     | _ -> failwith "not a trivial axiom" in
1212   try let th = tryfind trivial_axiom (keq @ klep @ kltp) in
1213       CONV_RULE (LAND_CONV REAL_POLY_CONV THENC REAL_RAT_RED_CONV) th
1214   with Failure _ ->
1215   let pol = itlist poly_mul (map fst ltp) (poly_const num_1) in
1216   let leq = lep @ ltp in
1217   let tryall d =
1218     let e = multidegree pol in
1219     let k = if e = 0 then 0 else d / e in
1220     let eq' = map fst eq in
1221     tryfind (fun i -> d,i,real_positivnullstellensatz_general false d eq' leq
1222                           (poly_neg(poly_pow pol i)))
1223             (0--k) in
1224   let d,i,(cert_ideal,cert_cone) = deepen tryall 0 in
1225   let proofs_ideal =
1226     map2 (fun q (p,ax) -> Eqmul(term_of_poly q,ax)) cert_ideal eq
1227   and proofs_cone = map term_of_sos cert_cone
1228   and proof_ne =
1229     if ltp = [] then Rational_lt num_1 else
1230     let p = end_itlist (fun s t -> Product(s,t)) (map snd ltp) in
1231     funpow i (fun q -> Product(p,q)) (Rational_lt num_1) in
1232   let proof = end_itlist (fun s t -> Sum(s,t))
1233                          (proof_ne :: proofs_ideal @ proofs_cone) in
1234   print_string("Translating proof certificate to HOL");
1235   print_newline();
1236   translator (eqs,les,lts) proof;;
1237
1238 (* ------------------------------------------------------------------------- *)
1239 (* A wrapper that tries to substitute away variables first.                  *)
1240 (* ------------------------------------------------------------------------- *)
1241
1242 let REAL_NONLINEAR_SUBST_PROVER =
1243   let zero = `&0:real`
1244   and mul_tm = `( * ):real->real->real`
1245   and shuffle1 =
1246     CONV_RULE(REWR_CONV(REAL_ARITH `a + x = (y:real) <=> x = y - a`))
1247   and shuffle2 =
1248     CONV_RULE(REWR_CONV(REAL_ARITH `x + a = (y:real) <=> x = y - a`)) in
1249   let rec substitutable_monomial fvs tm =
1250     match tm with
1251       Var(_,Tyapp("real",[])) when not (mem tm fvs) -> Int 1,tm
1252     | Comb(Comb(Const("real_mul",_),c),(Var(_,_) as t))
1253          when is_ratconst c & not (mem t fvs)
1254           -> rat_of_term c,t
1255     | Comb(Comb(Const("real_add",_),s),t) ->
1256          (try substitutable_monomial (union (frees t) fvs) s
1257           with Failure _ -> substitutable_monomial (union (frees s) fvs) t)
1258     | _ -> failwith "substitutable_monomial"
1259   and isolate_variable v th =
1260     match lhs(concl th) with
1261       x when x = v -> th
1262     | Comb(Comb(Const("real_add",_),(Var(_,Tyapp("real",[])) as x)),t)
1263         when x = v -> shuffle2 th
1264     | Comb(Comb(Const("real_add",_),s),t) ->
1265         isolate_variable v(shuffle1 th) in
1266   let make_substitution th =
1267     let (c,v) = substitutable_monomial [] (lhs(concl th)) in
1268     let th1 = AP_TERM (mk_comb(mul_tm,term_of_rat(Int 1 // c))) th in
1269     let th2 = CONV_RULE(BINOP_CONV REAL_POLY_MUL_CONV) th1 in
1270     CONV_RULE (RAND_CONV REAL_POLY_CONV) (isolate_variable v th2) in
1271   fun translator ->
1272     let rec substfirst(eqs,les,lts) =
1273       try let eth = tryfind make_substitution eqs in
1274           let modify =
1275             CONV_RULE(LAND_CONV(SUBS_CONV[eth] THENC REAL_POLY_CONV)) in
1276           substfirst(filter (fun t -> lhand(concl t) <> zero) (map modify eqs),
1277                      map modify les,map modify lts)
1278       with Failure _ -> REAL_NONLINEAR_PROVER translator (eqs,les,lts) in
1279     substfirst;;
1280
1281 (* ------------------------------------------------------------------------- *)
1282 (* Overall function.                                                         *)
1283 (* ------------------------------------------------------------------------- *)
1284
1285 let REAL_SOS =
1286   let init = GEN_REWRITE_CONV ONCE_DEPTH_CONV [DECIMAL]
1287   and pure = GEN_REAL_ARITH REAL_NONLINEAR_SUBST_PROVER in
1288   fun tm -> let th = init tm in EQ_MP (SYM th) (pure(rand(concl th)));;
1289
1290 (* ------------------------------------------------------------------------- *)
1291 (* Add hacks for division.                                                   *)
1292 (* ------------------------------------------------------------------------- *)
1293
1294 let REAL_SOSFIELD =
1295   let inv_tm = `inv:real->real` in
1296   let prenex_conv =
1297     TOP_DEPTH_CONV BETA_CONV THENC
1298     PURE_REWRITE_CONV[FORALL_SIMP; EXISTS_SIMP; real_div;
1299                       REAL_INV_INV; REAL_INV_MUL; GSYM REAL_POW_INV] THENC
1300     NNFC_CONV THENC DEPTH_BINOP_CONV `(/\)` CONDS_CELIM_CONV THENC
1301     PRENEX_CONV
1302   and setup_conv = NNF_CONV THENC WEAK_CNF_CONV THENC CONJ_CANON_CONV
1303   and core_rule t =
1304     try REAL_ARITH t
1305     with Failure _ -> try REAL_RING t
1306     with Failure _ -> REAL_SOS t
1307   and is_inv =
1308     let is_div = is_binop `(/):real->real->real` in
1309     fun tm -> (is_div tm or (is_comb tm & rator tm = inv_tm)) &
1310               not(is_ratconst(rand tm)) in
1311   let BASIC_REAL_FIELD tm =
1312     let is_freeinv t = is_inv t & free_in t tm in
1313     let itms = setify(map rand (find_terms is_freeinv tm)) in
1314     let hyps = map (fun t -> SPEC t REAL_MUL_RINV) itms in
1315     let tm' = itlist (fun th t -> mk_imp(concl th,t)) hyps tm in
1316     let itms' = map (curry mk_comb inv_tm) itms in
1317     let gvs = map (genvar o type_of) itms' in
1318     let tm'' = subst (zip gvs itms') tm' in
1319     let th1 = setup_conv tm'' in
1320     let cjs = conjuncts(rand(concl th1)) in
1321     let ths = map core_rule cjs in
1322     let th2 = EQ_MP (SYM th1) (end_itlist CONJ ths) in
1323     rev_itlist (C MP) hyps (INST (zip itms' gvs) th2) in
1324   fun tm ->
1325     let th0 = prenex_conv tm in
1326     let tm0 = rand(concl th0) in
1327     let avs,bod = strip_forall tm0 in
1328     let th1 = setup_conv bod in
1329     let ths = map BASIC_REAL_FIELD (conjuncts(rand(concl th1))) in
1330     EQ_MP (SYM th0) (GENL avs (EQ_MP (SYM th1) (end_itlist CONJ ths)));;
1331
1332 (* ------------------------------------------------------------------------- *)
1333 (* Integer version.                                                          *)
1334 (* ------------------------------------------------------------------------- *)
1335
1336 let INT_SOS =
1337   let atom_CONV =
1338     let pth = prove
1339       (`(~(x <= y) <=> y + &1 <= x:int) /\
1340         (~(x < y) <=> y <= x) /\
1341         (~(x = y) <=> x + &1 <= y \/ y + &1 <= x) /\
1342         (x < y <=> x + &1 <= y)`,
1343        REWRITE_TAC[INT_NOT_LE; INT_NOT_LT; INT_NOT_EQ; INT_LT_DISCRETE]) in
1344     GEN_REWRITE_CONV I [pth]
1345   and bub_CONV = GEN_REWRITE_CONV TOP_SWEEP_CONV
1346    [int_eq; int_le; int_lt; int_ge; int_gt;
1347     int_of_num_th; int_neg_th; int_add_th; int_mul_th;
1348     int_sub_th; int_pow_th; int_abs_th; int_max_th; int_min_th] in
1349   let base_CONV = TRY_CONV atom_CONV THENC bub_CONV in
1350   let NNF_NORM_CONV = GEN_NNF_CONV false
1351    (base_CONV,fun t -> base_CONV t,base_CONV(mk_neg t)) in
1352   let init_CONV =
1353     GEN_REWRITE_CONV DEPTH_CONV [FORALL_SIMP; EXISTS_SIMP] THENC
1354     GEN_REWRITE_CONV DEPTH_CONV [INT_GT; INT_GE] THENC
1355     CONDS_ELIM_CONV THENC NNF_NORM_CONV in
1356   let p_tm = `p:bool`
1357   and not_tm = `(~)` in
1358   let pth = TAUT(mk_eq(mk_neg(mk_neg p_tm),p_tm)) in
1359   fun tm ->
1360     let th0 = INST [tm,p_tm] pth
1361     and th1 = NNF_NORM_CONV(mk_neg tm) in
1362     let th2 = REAL_SOS(mk_neg(rand(concl th1))) in
1363     EQ_MP th0 (EQ_MP (AP_TERM not_tm (SYM th1)) th2);;
1364
1365 (* ------------------------------------------------------------------------- *)
1366 (* Natural number version.                                                   *)
1367 (* ------------------------------------------------------------------------- *)
1368
1369 let SOS_RULE tm =
1370   let avs = frees tm in
1371   let tm' = list_mk_forall(avs,tm) in
1372   let th1 = NUM_TO_INT_CONV tm' in
1373   let th2 = INT_SOS (rand(concl th1)) in
1374   SPECL avs (EQ_MP (SYM th1) th2);;
1375
1376 (* ------------------------------------------------------------------------- *)
1377 (* Now pure SOS stuff.                                                       *)
1378 (* ------------------------------------------------------------------------- *)
1379
1380 prioritize_real();;
1381
1382 (* ------------------------------------------------------------------------- *)
1383 (* Some combinatorial helper functions.                                      *)
1384 (* ------------------------------------------------------------------------- *)
1385
1386 let rec allpermutations l =
1387   if l = [] then [[]] else
1388   itlist (fun h acc -> map (fun t -> h::t)
1389                 (allpermutations (subtract l [h])) @ acc) l [];;
1390
1391 let allvarorders l =
1392   map (fun vlis x -> index x vlis) (allpermutations l);;
1393
1394 let changevariables_monomial zoln (m:monomial) =
1395   foldl (fun a x k -> (assoc x zoln |-> k) a) monomial_1 m;;
1396
1397 let changevariables zoln pol =
1398   foldl (fun a m c -> (changevariables_monomial zoln m |-> c) a)
1399         poly_0 pol;;
1400
1401 (* ------------------------------------------------------------------------- *)
1402 (* Return to original non-block matrices.                                    *)
1403 (* ------------------------------------------------------------------------- *)
1404
1405 let sdpa_of_vector (v:vector) =
1406   let n = vec_dim v in
1407   let strs = map (decimalize 20 o element v) (1--n) in
1408   end_itlist (fun x y -> x ^ " " ^ y) strs ^ "\n";;
1409
1410 let sdpa_of_blockdiagonal k m =
1411   let pfx = string_of_int k ^" " in
1412   let ents =
1413     foldl (fun a (b,i,j) c -> if i > j then a else ((b,i,j),c)::a) [] m in
1414   let entss = sort (increasing fst) ents in
1415   itlist (fun ((b,i,j),c) a ->
1416      pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
1417      " " ^ decimalize 20 c ^ "\n" ^ a) entss "";;
1418
1419 let sdpa_of_matrix k (m:matrix) =
1420   let pfx = string_of_int k ^ " 1 " in
1421   let ms = foldr (fun (i,j) c a -> if i > j then a else ((i,j),c)::a)
1422                  (snd m) [] in
1423   let mss = sort (increasing fst) ms in
1424   itlist (fun ((i,j),c) a ->
1425      pfx ^ string_of_int i ^ " " ^ string_of_int j ^
1426      " " ^ decimalize 20 c ^ "\n" ^ a) mss "";;
1427
1428 let sdpa_of_problem comment obj mats =
1429   let m = length mats - 1
1430   and n,_ = dimensions (hd mats) in
1431   "\"" ^ comment ^ "\"\n" ^
1432   string_of_int m ^ "\n" ^
1433   "1\n" ^
1434   string_of_int n ^ "\n" ^
1435   sdpa_of_vector obj ^
1436   itlist2 (fun k m a -> sdpa_of_matrix (k - 1) m ^ a)
1437           (1--length mats) mats "";;
1438
1439 let run_sdpa dbg obj mats =
1440   let input_file = Filename.temp_file "sos" ".dat-s" in
1441   let output_file =
1442     String.sub input_file 0 (String.length input_file - 6) ^ ".out"
1443   and params_file = Filename.concat (!temp_path) "param.sdpa" in
1444   file_of_string input_file (sdpa_of_problem "" obj mats);
1445   file_of_string params_file sdpa_params;
1446   Sys.command("cd "^(!temp_path)^"; sdpa "^input_file ^ " " ^ output_file ^
1447               (if dbg then "" else "> /dev/null"));
1448   let op = string_of_file output_file in
1449   if not(sdpa_run_succeeded op) then failwith "sdpa: call failed" else
1450   let res = parse_sdpaoutput op in
1451   ((if dbg then ()
1452    else (Sys.remove input_file; Sys.remove output_file));
1453    res);;
1454
1455 let sdpa obj mats = run_sdpa (!debugging) obj mats;;
1456
1457 let run_csdp dbg obj mats =
1458   let input_file = Filename.temp_file "sos" ".dat-s" in
1459   let output_file =
1460     String.sub input_file 0 (String.length input_file - 6) ^ ".out"
1461   and params_file = Filename.concat (!temp_path) "param.csdp" in
1462   file_of_string input_file (sdpa_of_problem "" obj mats);
1463   file_of_string params_file csdp_params;
1464   let rv = Sys.command("cd "^(!temp_path)^"; csdp "^input_file ^
1465                        " " ^ output_file ^
1466                        (if dbg then "" else "> /dev/null")) in
1467   let op = string_of_file output_file in
1468   let res = parse_csdpoutput op in
1469   ((if dbg then ()
1470     else (Sys.remove input_file; Sys.remove output_file));
1471     rv,res);;
1472
1473 let csdp obj mats =
1474   let rv,res = run_csdp (!debugging) obj mats in
1475   (if rv = 1 or rv = 2 then failwith "csdp: Problem is infeasible"
1476    else if rv = 3 then
1477     (Format.print_string "csdp warning: Reduced accuracy";
1478      Format.print_newline())
1479    else if rv <> 0 then failwith("csdp: error "^string_of_int rv)
1480    else ());
1481   res;;
1482
1483 (* ------------------------------------------------------------------------- *)
1484 (* Sum-of-squares function with some lowbrow symmetry reductions.            *)
1485 (* ------------------------------------------------------------------------- *)
1486
1487 let sumofsquares_general_symmetry tool pol =
1488   let vars = poly_variables pol
1489   and lpps = newton_polytope pol in
1490   let n = length lpps in
1491   let sym_eqs =
1492     let invariants = filter
1493      (fun vars' ->
1494         is_undefined(poly_sub pol (changevariables (zip vars vars') pol)))
1495      (allpermutations vars) in
1496     let lpps2 = allpairs monomial_mul lpps lpps in
1497     let lpp2_classes =
1498        setify(map (fun m ->
1499          setify(map (fun vars' -> changevariables_monomial (zip vars vars') m)
1500                    invariants)) lpps2) in
1501     let lpns = zip lpps (1--length lpps) in
1502     let lppcs =
1503       filter (fun (m,(n1,n2)) -> n1 <= n2)
1504              (allpairs
1505                (fun (m1,n1) (m2,n2) -> (m1,m2),(n1,n2)) lpns lpns) in
1506     let clppcs = end_itlist (@)
1507        (map (fun ((m1,m2),(n1,n2)) ->
1508                 map (fun vars' ->
1509                         (changevariables_monomial (zip vars vars') m1,
1510                          changevariables_monomial (zip vars vars') m2),(n1,n2))
1511                     invariants)
1512             lppcs) in
1513     let clppcs_dom = setify(map fst clppcs) in
1514     let clppcs_cls = map (fun d -> filter (fun (e,_) -> e = d) clppcs)
1515                          clppcs_dom in
1516     let eqvcls = map (setify o map snd) clppcs_cls in
1517     let mk_eq cls acc =
1518       match cls with
1519         [] -> raise Sanity
1520       | [h] -> acc
1521       | h::t -> map (fun k -> (k |-> Int(-1)) (h |=> Int 1)) t @ acc in
1522     itlist mk_eq eqvcls [] in
1523   let eqs = foldl (fun a x y -> y::a) []
1524    (itern 1 lpps (fun m1 n1 ->
1525         itern 1 lpps (fun m2 n2 f ->
1526                 let m = monomial_mul m1 m2 in
1527                 if n1 > n2 then f else
1528                 let c = if n1 = n2 then Int 1 else Int 2 in
1529                 (m |-> ((n1,n2) |-> c) (tryapplyd f m undefined)) f))
1530        (foldl (fun a m c -> (m |-> ((0,0)|=>c)) a)
1531               undefined pol)) @
1532     sym_eqs in
1533   let pvs,assig = eliminate_all_equations (0,0) eqs in
1534   let allassig = itlist (fun v -> (v |-> (v |=> Int 1))) pvs assig in
1535   let qvars = (0,0)::pvs in
1536   let diagents =
1537     end_itlist equation_add (map (fun i -> apply allassig (i,i)) (1--n)) in
1538   let mk_matrix v =
1539    ((n,n),
1540     foldl (fun m (i,j) ass -> let c = tryapplyd ass v (Int 0) in
1541                               if c =/ Int 0 then m else
1542                               ((j,i) |-> c) (((i,j) |-> c) m))
1543           undefined allassig :matrix) in
1544   let mats = map mk_matrix qvars
1545   and obj = length pvs,
1546             itern 1 pvs (fun v i -> (i |--> tryapplyd diagents v (Int 0)))
1547                 undefined in
1548   let raw_vec = if pvs = [] then vec_0 0 else tool obj mats in
1549   let find_rounding d =
1550    (if !debugging then
1551      (Format.print_string("Trying rounding with limit "^string_of_num d);
1552       Format.print_newline())
1553     else ());
1554     let vec = nice_vector d raw_vec in
1555     let mat = iter (1,vec_dim vec)
1556      (fun i a -> matrix_add (matrix_cmul (element vec i) (el i mats)) a)
1557      (matrix_neg (el 0 mats)) in
1558     deration(diag mat) in
1559   let rat,dia =
1560     if pvs = [] then
1561        let mat = matrix_neg (el 0 mats) in
1562        deration(diag mat)
1563     else
1564        tryfind find_rounding (map Num.num_of_int (1--31) @
1565                               map pow2 (5--66)) in
1566   let poly_of_lin(d,v) =
1567     d,foldl(fun a i c -> (el (i - 1) lpps |-> c) a) undefined (snd v) in
1568   let lins = map poly_of_lin dia in
1569   let sqs = map (fun (d,l) -> poly_mul (poly_const d) (poly_pow l 2)) lins in
1570   let sos = poly_cmul rat (end_itlist poly_add sqs) in
1571   if is_undefined(poly_sub sos pol) then rat,lins else raise Sanity;;
1572
1573 let sumofsquares = sumofsquares_general_symmetry csdp;;
1574
1575 (* ------------------------------------------------------------------------- *)
1576 (* Pure HOL SOS conversion.                                                  *)
1577 (* ------------------------------------------------------------------------- *)
1578
1579 let SOS_CONV =
1580   let mk_square =
1581     let pow_tm = `(pow)` and two_tm = `2` in
1582     fun tm -> mk_comb(mk_comb(pow_tm,tm),two_tm)
1583   and mk_prod = mk_binop `(*)`
1584   and mk_sum = mk_binop `(+)` in
1585   fun tm ->
1586     let k,sos = sumofsquares(poly_of_term tm) in
1587     let mk_sqtm(c,p) =
1588       mk_prod (term_of_rat(k */ c)) (mk_square(term_of_poly p)) in
1589     let tm' = end_itlist mk_sum (map mk_sqtm sos) in
1590     let th = REAL_POLY_CONV tm and th' = REAL_POLY_CONV tm' in
1591     TRANS th (SYM th');;
1592
1593 (* ------------------------------------------------------------------------- *)
1594 (* Attempt to prove &0 <= x by direct SOS decomposition.                     *)
1595 (* ------------------------------------------------------------------------- *)
1596
1597 let PURE_SOS_TAC =
1598   let tac =
1599     MATCH_ACCEPT_TAC(REWRITE_RULE[GSYM REAL_POW_2] REAL_LE_SQUARE) ORELSE
1600     MATCH_ACCEPT_TAC REAL_LE_SQUARE ORELSE
1601     (MATCH_MP_TAC REAL_LE_ADD THEN CONJ_TAC) ORELSE
1602     (MATCH_MP_TAC REAL_LE_MUL THEN CONJ_TAC) ORELSE
1603     CONV_TAC(RAND_CONV REAL_RAT_REDUCE_CONV THENC REAL_RAT_LE_CONV) in
1604   REPEAT GEN_TAC THEN REWRITE_TAC[real_ge] THEN
1605   GEN_REWRITE_TAC I [GSYM REAL_SUB_LE] THEN
1606   CONV_TAC(RAND_CONV SOS_CONV) THEN
1607   REPEAT tac THEN NO_TAC;;
1608
1609 let PURE_SOS tm = prove(tm,PURE_SOS_TAC);;
1610
1611 (* ------------------------------------------------------------------------- *)
1612 (* Examples.                                                                 *)
1613 (* ------------------------------------------------------------------------- *)
1614
1615 (*****
1616
1617 time REAL_SOS
1618   `a1 >= &0 /\ a2 >= &0 /\
1619    (a1 * a1 + a2 * a2 = b1 * b1 + b2 * b2 + &2) /\
1620    (a1 * b1 + a2 * b2 = &0)
1621    ==> a1 * a2 - b1 * b2 >= &0`;;
1622
1623 time REAL_SOS `&3 * x + &7 * a < &4 /\ &3 < &2 * x ==> a < &0`;;
1624
1625 time REAL_SOS
1626   `b pow 2 < &4 * a * c ==> ~(a * x pow 2 + b * x + c = &0)`;;
1627
1628 time REAL_SOS
1629   `(a * x pow 2 + b * x + c = &0) ==> b pow 2 >= &4 * a * c`;;
1630
1631 time REAL_SOS
1632  `&0 <= x /\ x <= &1 /\ &0 <= y /\ y <= &1
1633   ==> x pow 2 + y pow 2 < &1 \/
1634       (x - &1) pow 2 + y pow 2 < &1 \/
1635       x pow 2 + (y - &1) pow 2 < &1 \/
1636       (x - &1) pow 2 + (y - &1) pow 2 < &1`;;
1637
1638 time REAL_SOS
1639  `&0 <= b /\ &0 <= c /\ &0 <= x /\ &0 <= y /\
1640   (x pow 2 = c) /\ (y pow 2 = a pow 2 * c + b)
1641   ==> a * c <= y * x`;;
1642
1643 time REAL_SOS
1644  `&0 <= x /\ &0 <= y /\ &0 <= z /\ x + y + z <= &3
1645   ==> x * y + x * z + y * z >= &3 * x * y * z`;;
1646
1647 time REAL_SOS
1648  `(x pow 2 + y pow 2 + z pow 2 = &1) ==> (x + y + z) pow 2 <= &3`;;
1649
1650 time REAL_SOS
1651  `(w pow 2 + x pow 2 + y pow 2 + z pow 2 = &1)
1652   ==> (w + x + y + z) pow 2 <= &4`;;
1653
1654 time REAL_SOS
1655  `x >= &1 /\ y >= &1 ==> x * y >= x + y - &1`;;
1656
1657 time REAL_SOS
1658  `x > &1 /\ y > &1 ==> x * y > x + y - &1`;;
1659
1660 time REAL_SOS
1661  `abs(x) <= &1
1662   ==> abs(&64 * x pow 7 - &112 * x pow 5 + &56 * x pow 3 - &7 * x) <= &1`;;
1663
1664 time REAL_SOS
1665  `abs(x - z) <= e /\ abs(y - z) <= e /\ &0 <= u /\ &0 <= v /\ (u + v = &1)
1666   ==> abs((u * x + v * y) - z) <= e`;;
1667
1668 (* ------------------------------------------------------------------------- *)
1669 (* One component of denominator in dodecahedral example.                     *)
1670 (* ------------------------------------------------------------------------- *)
1671
1672 time REAL_SOS
1673   `&2 <= x /\ x <= &125841 / &50000 /\
1674    &2 <= y /\ y <= &125841 / &50000 /\
1675    &2 <= z /\ z <= &125841 / &50000
1676    ==> &2 * (x * z + x * y + y * z) - (x * x + y * y + z * z) >= &0`;;
1677
1678 (* ------------------------------------------------------------------------- *)
1679 (* Over a larger but simpler interval.                                       *)
1680 (* ------------------------------------------------------------------------- *)
1681
1682 time REAL_SOS
1683  `&2 <= x /\ x <= &4 /\ &2 <= y /\ y <= &4 /\ &2 <= z /\ z <= &4
1684   ==> &0 <= &2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)`;;
1685
1686 (* ------------------------------------------------------------------------- *)
1687 (* We can do 12. I think 12 is a sharp bound; see PP's certificate.          *)
1688 (* ------------------------------------------------------------------------- *)
1689
1690 time REAL_SOS
1691  `&2 <= x /\ x <= &4 /\ &2 <= y /\ y <= &4 /\ &2 <= z /\ z <= &4
1692   ==> &12 <= &2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)`;;
1693
1694 (* ------------------------------------------------------------------------- *)
1695 (* Gloptipoly example.                                                       *)
1696 (* ------------------------------------------------------------------------- *)
1697
1698 (*** This works but normalization takes minutes
1699
1700 time REAL_SOS
1701  `(x - y - &2 * x pow 4 = &0) /\ &0 <= x /\ x <= &2 /\ &0 <= y /\ y <= &3
1702   ==> y pow 2 - &7 * y - &12 * x + &17 >= &0`;;
1703
1704  ***)
1705
1706 (* ------------------------------------------------------------------------- *)
1707 (* Inequality from sci.math (see "Leon-Sotelo, por favor").                  *)
1708 (* ------------------------------------------------------------------------- *)
1709
1710 time REAL_SOS
1711   `&0 <= x /\ &0 <= y /\ (x * y = &1)
1712    ==> x + y <= x pow 2 + y pow 2`;;
1713
1714 time REAL_SOS
1715   `&0 <= x /\ &0 <= y /\ (x * y = &1)
1716    ==> x * y * (x + y) <= x pow 2 + y pow 2`;;
1717
1718 time REAL_SOS
1719   `&0 <= x /\ &0 <= y ==> x * y * (x + y) pow 2 <= (x pow 2 + y pow 2) pow 2`;;
1720
1721 (* ------------------------------------------------------------------------- *)
1722 (* Some examples over integers and natural numbers.                          *)
1723 (* ------------------------------------------------------------------------- *)
1724
1725 time SOS_RULE `!m n. 2 * m + n = (n + m) + m`;;
1726 time SOS_RULE `!n. ~(n = 0) ==> (0 MOD n = 0)`;;
1727 time SOS_RULE  `!m n. m < n ==> (m DIV n = 0)`;;
1728 time SOS_RULE `!n:num. n <= n * n`;;
1729 time SOS_RULE  `!m n. n * (m DIV n) <= m`;;
1730 time SOS_RULE `!n. ~(n = 0) ==> (0 DIV n = 0)`;;
1731 time SOS_RULE `!m n p. ~(p = 0) /\ m <= n ==> m DIV p <= n DIV p`;;
1732 time SOS_RULE `!a b n. ~(a = 0) ==> (n <= b DIV a <=> a * n <= b)`;;
1733
1734 (* ------------------------------------------------------------------------- *)
1735 (* This is particularly gratifying --- cf hideous manual proof in arith.ml   *)
1736 (* ------------------------------------------------------------------------- *)
1737
1738 (*** This doesn't now seem to work as well as it did; what changed?
1739
1740 time SOS_RULE
1741  `!a b c d. ~(b = 0) /\ b * c < (a + 1) * d ==> c DIV d <= a DIV b`;;
1742
1743  ***)
1744
1745 (* ------------------------------------------------------------------------- *)
1746 (* Key lemma for injectivity of Cantor-type pairing functions.               *)
1747 (* ------------------------------------------------------------------------- *)
1748
1749 time SOS_RULE
1750  `!x1 y1 x2 y2. ((x1 + y1) EXP 2 + x1 + 1 = (x2 + y2) EXP 2 + x2 + 1)
1751                 ==> (x1 + y1 = x2 + y2)`;;
1752
1753 time SOS_RULE
1754  `!x1 y1 x2 y2. ((x1 + y1) EXP 2 + x1 + 1 = (x2 + y2) EXP 2 + x2 + 1) /\
1755                 (x1 + y1 = x2 + y2)
1756                 ==> (x1 = x2) /\ (y1 = y2)`;;
1757
1758 time SOS_RULE
1759  `!x1 y1 x2 y2.
1760       (((x1 + y1) EXP 2 + 3 * x1 + y1) DIV 2 =
1761        ((x2 + y2) EXP 2 + 3 * x2 + y2) DIV 2)
1762        ==> (x1 + y1 = x2 + y2)`;;
1763
1764 time SOS_RULE
1765  `!x1 y1 x2 y2.
1766       (((x1 + y1) EXP 2 + 3 * x1 + y1) DIV 2 =
1767        ((x2 + y2) EXP 2 + 3 * x2 + y2) DIV 2) /\
1768       (x1 + y1 = x2 + y2)
1769       ==> (x1 = x2) /\ (y1 = y2)`;;
1770
1771 (* ------------------------------------------------------------------------- *)
1772 (* Reciprocal multiplication (actually just ARITH_RULE does these).          *)
1773 (* ------------------------------------------------------------------------- *)
1774
1775 time SOS_RULE `x <= 127 ==> ((86 * x) DIV 256 = x DIV 3)`;;
1776
1777 time SOS_RULE `x < 2 EXP 16 ==> ((104858 * x) DIV (2 EXP 20) = x DIV 10)`;;
1778
1779 (* ------------------------------------------------------------------------- *)
1780 (* This is more impressive since it's really nonlinear. See REMAINDER_DECODE *)
1781 (* ------------------------------------------------------------------------- *)
1782
1783 time SOS_RULE `0 < m /\ m < n  ==> ((m * ((n * x) DIV m + 1)) DIV n = x)`;;
1784
1785 (* ------------------------------------------------------------------------- *)
1786 (* Some conversion examples.                                                 *)
1787 (* ------------------------------------------------------------------------- *)
1788
1789 time SOS_CONV
1790  `&2 * x pow 4 + &2 * x pow 3 * y - x pow 2 * y pow 2 + &5 * y pow 4`;;
1791
1792 time SOS_CONV
1793  `x pow 4 - (&2 * y * z + &1) * x pow 2 +
1794   (y pow 2 * z pow 2 + &2 * y * z + &2)`;;
1795
1796 time SOS_CONV `&4 * x pow 4 +
1797           &4 * x pow 3 * y - &7 * x pow 2 * y pow 2 - &2 * x * y pow 3 +
1798           &10 * y pow 4`;;
1799
1800 time SOS_CONV `&4 * x pow 4 * y pow 6 + x pow 2 - x * y pow 2 + y pow 2`;;
1801
1802 time SOS_CONV
1803   `&4096 * (x pow 4 + x pow 2 + z pow 6 - &3 * x pow 2 * z pow 2) + &729`;;
1804
1805 time SOS_CONV
1806  `&120 * x pow 2 - &63 * x pow 4 + &10 * x pow 6 +
1807   &30 * x * y - &120 * y pow 2 + &120 * y pow 4 + &31`;;
1808
1809 time SOS_CONV
1810  `&9 * x pow 2 * y pow 4 + &9 * x pow 2 * z pow 4 + &36 * x pow 2 * y pow 3 +
1811   &36 * x pow 2 * y pow 2 - &48 * x * y * z pow 2 + &4 * y pow 4 +
1812   &4 * z pow 4 - &16 * y pow 3 + &16 * y pow 2`;;
1813
1814 time SOS_CONV
1815  `(x pow 2 + y pow 2 + z pow 2) *
1816   (x pow 4 * y pow 2 + x pow 2 * y pow 4 +
1817    z pow 6 - &3 * x pow 2 * y pow 2 *   z pow 2)`;;
1818
1819 time SOS_CONV
1820  `x pow 4 + y pow 4 + z pow 4 - &4 * x * y * z + x + y + z + &3`;;
1821
1822 (*** I think this will work, but normalization is slow
1823
1824 time SOS_CONV
1825  `&100 * (x pow 4 + y pow 4 + z pow 4 - &4 * x * y * z + x + y + z) + &212`;;
1826
1827  ***)
1828
1829 time SOS_CONV
1830  `&100 * ((&2 * x - &2) pow 2 + (x pow 3 - &8 * x - &2) pow 2) - &588`;;
1831
1832 time SOS_CONV
1833  `x pow 2 * (&120 - &63 * x pow 2 + &10 * x pow 4) + &30 * x * y +
1834     &30 * y pow 2 * (&4 * y pow 2 - &4) + &31`;;
1835
1836 (* ------------------------------------------------------------------------- *)
1837 (* Example of basic rule.                                                    *)
1838 (* ------------------------------------------------------------------------- *)
1839
1840 time PURE_SOS
1841   `!x. x pow 4 + y pow 4 + z pow 4 - &4 * x * y * z + x + y + z + &3
1842        >= &1 / &7`;;
1843
1844 time PURE_SOS
1845  `&0 <= &98 * x pow 12 +
1846         -- &980 * x pow 10 +
1847         &3038 * x pow 8 +
1848         -- &2968 * x pow 6 +
1849         &1022 * x pow 4 +
1850         -- &84 * x pow 2 +
1851         &2`;;
1852
1853 time PURE_SOS
1854  `!x. &0 <= &2 * x pow 14 +
1855             -- &84 * x pow 12 +
1856             &1022 * x pow 10 +
1857             -- &2968 * x pow 8 +
1858             &3038 * x pow 6 +
1859             -- &980 * x pow 4 +
1860             &98 * x pow 2`;;
1861
1862 (* ------------------------------------------------------------------------- *)
1863 (* From Zeng et al, JSC vol 37 (2004), p83-99.                               *)
1864 (* All of them work nicely with pure SOS_CONV, except (maybe) the one noted. *)
1865 (* ------------------------------------------------------------------------- *)
1866
1867 PURE_SOS
1868   `x pow 6 + y pow 6 + z pow 6 - &3 * x pow 2 * y pow 2 * z pow 2 >= &0`;;
1869
1870 PURE_SOS `x pow 4 + y pow 4 + z pow 4 + &1 - &4*x*y*z >= &0`;;
1871
1872 PURE_SOS `x pow 4 + &2*x pow 2*z + x pow 2 - &2*x*y*z + &2*y pow 2*z pow 2 +
1873 &2*y*z pow 2 + &2*z pow 2 - &2*x + &2* y*z + &1 >= &0`;;
1874
1875 (**** This is harder. Interestingly, this fails the pure SOS test, it seems.
1876       Yet only on rounding(!?) Poor Newton polytope optimization or something?
1877       But REAL_SOS does finally converge on the second run at level 12!
1878
1879 REAL_SOS
1880 `x pow 4*y pow 4 - &2*x pow 5*y pow 3*z pow 2 + x pow 6*y pow 2*z pow 4 + &2*x
1881 pow 2*y pow 3*z - &4* x pow 3*y pow 2*z pow 3 + &2*x pow 4*y*z pow 5 + z pow
1882 2*y pow 2 - &2*z pow 4*y*x + z pow 6*x pow 2 >= &0`;;
1883
1884  ****)
1885
1886 PURE_SOS
1887 `x pow 4 + &4*x pow 2*y pow 2 + &2*x*y*z pow 2 + &2*x*y*w pow 2 + y pow 4 + z
1888 pow 4 + w pow 4 + &2*z pow 2*w pow 2 + &2*x pow 2*w + &2*y pow 2*w + &2*x*y +
1889 &3*w pow 2 + &2*z pow 2 + &1 >= &0`;;
1890
1891 PURE_SOS
1892 `w pow 6 + &2*z pow 2*w pow 3 + x pow 4 + y pow 4 + z pow 4 + &2*x pow 2*w +
1893 &2*x pow 2*z + &3*x pow 2 + w pow 2 + &2*z*w + z pow 2 + &2*z + &2*w + &1 >=
1894 &0`;;
1895
1896 *****)