--- a/CONTRIBUTORS Tue Jul 28 20:26:39 2009 +0200
+++ b/CONTRIBUTORS Wed Jul 29 09:06:49 2009 +0200
@@ -7,6 +7,9 @@
Contributions to this Isabelle version
--------------------------------------
+* July 2009: Philipp Meyer, TUM
+ HOL/Library/Sum_of_Squares: functionality to call a remote csdp prover
+
* July 2009: Florian Haftmann, TUM
New quickcheck implementation using new code generator
@@ -19,6 +22,9 @@
* June 2009: Florian Haftmann, TUM
HOL/Library/Tree: searchtrees implementing mappings, ready to use for code generation
+* March 2009: Philipp Meyer, TUM
+ minimalization algorithm for results from sledgehammer call
+
Contributions to Isabelle2009
-----------------------------
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/lib/scripts/neos/NeosCSDPClient.py Wed Jul 29 09:06:49 2009 +0200
@@ -0,0 +1,56 @@
+#!/usr/bin/env python
+import sys
+import xmlrpclib
+import time
+import re
+
+from config import Variables
+
+if len(sys.argv) < 3 or len(sys.argv) > 3:
+ sys.stderr.write("Usage: NeosCSDPClient <input_filename> <output_filename>\n")
+ sys.exit(1)
+
+neos=xmlrpclib.Server("http://%s:%d" % (Variables.NEOS_HOST, Variables.NEOS_PORT))
+
+xmlfile = open(sys.argv[1],"r")
+xml_pre = "<document>\n<category>sdp</category>\n<solver>csdp</solver>\n<inputMethod>SPARSE_SDPA</inputMethod>\n<dat><![CDATA["
+xml_post = "]]></dat>\n</document>\n"
+xml = xml_pre
+buffer = 1
+while buffer:
+ buffer = xmlfile.read()
+ xml += buffer
+xmlfile.close()
+xml += xml_post
+
+(jobNumber,password) = neos.submitJob(xml)
+
+if jobNumber == 0:
+ sys.stdout.write("error submitting job: %s" % password)
+ sys.exit(-1)
+else:
+ sys.stdout.write("jobNumber = %d\tpassword = %s\n" % (jobNumber,password))
+
+offset=0
+status="Waiting"
+while status == "Running" or status=="Waiting":
+ time.sleep(1)
+ (msg,offset) = neos.getIntermediateResults(jobNumber,password,offset)
+ sys.stdout.write(msg.data)
+ status = neos.getJobStatus(jobNumber, password)
+
+msg = neos.getFinalResults(jobNumber, password).data
+result = msg.split("Solution:")
+
+sys.stdout.write(result[0])
+if len(result) > 1:
+ plain_msg = result[1].strip()
+ if plain_msg != "":
+ output = open(sys.argv[2],"w")
+ output.write(plain_msg)
+ output.close()
+ sys.exit(0)
+
+sys.exit(2)
+
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/lib/scripts/neos/config.py Wed Jul 29 09:06:49 2009 +0200
@@ -0,0 +1,3 @@
+class Variables:
+ NEOS_HOST="neos.mcs.anl.gov"
+ NEOS_PORT=3332
--- a/src/HOL/Library/Sum_Of_Squares.thy Tue Jul 28 20:26:39 2009 +0200
+++ b/src/HOL/Library/Sum_Of_Squares.thy Wed Jul 29 09:06:49 2009 +0200
@@ -6,19 +6,22 @@
multiplication and ordering using semidefinite programming*}
theory Sum_Of_Squares
imports Complex_Main (* "~~/src/HOL/Decision_Procs/Dense_Linear_Order" *)
- uses "positivstellensatz.ML" "sum_of_squares.ML"
+ uses "positivstellensatz.ML" "sum_of_squares.ML" "sos_wrapper.ML"
begin
(* Note:
-In order to use the method sos, install CSDP (https://projects.coin-or.org/Csdp/) and put the executable csdp on your path.
+In order to use the method sos, call it with (sos remote_csdp) to use the remote solver
+or install CSDP (https://projects.coin-or.org/Csdp/), put the executable csdp on your path,
+and call it with (sos csdp)
*)
-method_setup sos = {* Scan.succeed (SIMPLE_METHOD' o Sos.sos_tac) *}
- "Prove universal problems over the reals using sums of squares"
+(* setup sos tactic *)
+setup SosWrapper.setup
-text{* Tests -- commented since they work only when csdp is installed -- see above *}
+text{* Tests -- commented since they work only when csdp is installed or take too long with remote csdps *}
+
(*
lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x \<Longrightarrow> a < 0" by sos
@@ -115,3 +118,4 @@
*)
end
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Library/sos_wrapper.ML Wed Jul 29 09:06:49 2009 +0200
@@ -0,0 +1,158 @@
+(* Title: sos_wrapper.ML
+ Author: Philipp Meyer, TU Muenchen
+
+Added functionality for sums of squares, e.g. calling a remote prover
+*)
+
+signature SOS_WRAPPER =
+sig
+
+ datatype prover_result = Success | PartialSuccess | Failure | Error
+ type prover = string * (int -> string -> prover_result * string)
+
+ val setup: theory -> theory
+end
+
+structure SosWrapper : SOS_WRAPPER=
+struct
+
+datatype prover_result = Success | PartialSuccess | Failure | Error
+type prover =
+ string * (* command name *)
+ (int -> string ->prover_result * string) (* function to find failure from return value and output *)
+
+
+(*** output control ***)
+
+fun debug s = Output.debug (fn () => s)
+val answer = Output.priority
+val write = Output.writeln
+
+(*** calling provers ***)
+
+val destdir = ref ""
+
+fun filename dir name =
+ let
+ val probfile = Path.basic (name ^ serial_string ())
+ in
+ if dir = "" then
+ File.tmp_path probfile
+ else
+ if File.exists (Path.explode dir) then
+ Path.append (Path.explode dir) probfile
+ else
+ error ("No such directory: " ^ dir)
+ end
+
+fun is_success Success = true
+ | is_success PartialSuccess = true
+ | is_success _ = false
+fun str_of_status Success = "Success"
+ | str_of_status PartialSuccess = "Partial Success"
+ | str_of_status Failure= "Failure"
+ | str_of_status Error= "Error"
+
+fun run_solver name (cmd, find_failure) input =
+ let
+ val _ = answer ("Calling solver: " ^ name)
+
+ (* create input file *)
+ val dir = ! destdir
+ val input_file = filename dir "sos_in"
+ val _ = File.write input_file input
+
+ val _ = debug "Solver input:"
+ val _ = debug input
+
+ (* call solver *)
+ val output_file = filename dir "sos_out"
+ val (output, rv) = system_out (cmd ^ " " ^ (Path.implode input_file) ^
+ " " ^ (Path.implode output_file))
+
+ (* read and analysize output *)
+ val (res, res_msg) = find_failure rv output
+ val result = if is_success res then File.read output_file else ""
+
+ (* remove temporary files *)
+ val _ = if dir = "" then (File.rm input_file ; if File.exists output_file then File.rm output_file else ()) else ()
+
+ val _ = debug "Solver output:"
+ val _ = debug output
+ val _ = debug "Solver result:"
+ val _ = debug result
+
+ val _ = answer (str_of_status res ^ ": " ^ res_msg)
+
+ in
+ if is_success res then
+ result
+ else
+ error ("Prover failed: " ^ res_msg)
+ end
+
+(*** various provers ***)
+
+(* local csdp client *)
+
+fun find_csdp_run_failure rv _ =
+ case rv of
+ 0 => (Success, "SDP solved")
+ | 1 => (Failure, "SDP is primal infeasible")
+ | 2 => (Failure, "SDP is dual infeasible")
+ | 3 => (PartialSuccess, "SDP solved with reduced accuracy")
+ | _ => (Failure, "return code is " ^ string_of_int rv)
+
+val csdp = ("csdp", find_csdp_run_failure)
+
+(* remote neos server *)
+
+fun find_neos_failure rv output =
+ if rv = 2 then (Failure, "no solution") else
+ if rv <> 0 then (Error, "return code is " ^ string_of_int rv) else
+ let
+ fun find_success str =
+ if String.isPrefix "Success: " str then
+ SOME (Success, unprefix "Success: " str)
+ else if String.isPrefix "Partial Success: " str then
+ SOME (PartialSuccess, unprefix "Partial Success: " str)
+ else if String.isPrefix "Failure: " str then
+ SOME (Failure, unprefix "Failure: " str)
+ else
+ NONE
+ val exit_line = get_first find_success (split_lines output)
+ in
+ case exit_line of
+ SOME (status, msg) =>
+ if String.isPrefix "SDP solved" msg then
+ (status, msg)
+ else (Failure, msg)
+ | NONE => (Failure, "no success")
+ end
+
+val neos_csdp = ("$ISABELLE_HOME/lib/scripts/neos/NeosCSDPClient.py", find_neos_failure)
+
+(* save provers in table *)
+
+val provers =
+ Symtab.make [("remote_csdp", neos_csdp),("csdp", csdp)]
+
+fun get_prover name =
+ case Symtab.lookup provers name of
+ SOME prover => prover
+ | NONE => error ("unknown prover: " ^ name)
+
+fun call_solver name =
+ run_solver name (get_prover name)
+
+(* setup tactic *)
+
+val def_solver = "remote_csdp"
+
+fun sos_solver name = (SIMPLE_METHOD' o (Sos.sos_tac (call_solver name)))
+
+val sos_method = Scan.optional (Scan.lift OuterParse.xname) def_solver >> sos_solver
+
+val setup = Method.setup @{binding sos} sos_method "Prove universal problems over the reals using sums of squares"
+
+end
--- a/src/HOL/Library/sum_of_squares.ML Tue Jul 28 20:26:39 2009 +0200
+++ b/src/HOL/Library/sum_of_squares.ML Wed Jul 29 09:06:49 2009 +0200
@@ -1,6 +1,21 @@
-structure Sos =
+(* Title: sum_of_squares.ML
+ Authors: Amine Chaieb, University of Cambridge
+ Philipp Meyer, TU Muenchen
+
+A tactic for proving nonlinear inequalities
+*)
+
+signature SOS =
+sig
+
+ val sos_tac : (string -> string) -> Proof.context -> int -> Tactical.tactic
+
+end
+
+structure Sos : SOS =
struct
+
val rat_0 = Rat.zero;
val rat_1 = Rat.one;
val rat_2 = Rat.two;
@@ -500,11 +515,11 @@
(* Minimize obj_1 * v_1 + ... obj_m * v_m *)
(* ------------------------------------------------------------------------- *)
-fun sdpa_of_problem comment obj mats =
+fun sdpa_of_problem obj mats =
let
val m = length mats - 1
val (n,_) = dimensions (hd mats)
- in "\"" ^ comment ^ "\"\n" ^
+ in
string_of_int m ^ "\n" ^
"1\n" ^
string_of_int n ^ "\n" ^
@@ -587,6 +602,14 @@
val parse_sdpaoutput = mkparser sdpaoutput
val parse_csdpoutput = mkparser csdpoutput
+(* Run prover on a problem in linear form. *)
+
+fun run_problem prover obj mats =
+ parse_csdpoutput (prover (sdpa_of_problem obj mats))
+
+(*
+UNUSED
+
(* Also parse the SDPA output to test success (CSDP yields a return code). *)
local
@@ -687,6 +710,8 @@
res)
end;
+*)
+
(* Try some apparently sensible scaling first. Note that this is purely to *)
(* get a cleaner translation to floating-point, and doesn't affect any of *)
(* the results, in principle. In practice it seems a lot better when there *)
@@ -761,6 +786,8 @@
in if c =/ rat_0 then a
else Intfunc.update (i,y) a end) v Intfunc.undefined):vector
+(*
+UNUSED
(* Reduce linear program to SDP (diagonal matrices) and test with CSDP. This *)
(* one tests A [-1;x1;..;xn] >= 0 (i.e. left column is negated constants). *)
@@ -823,8 +850,12 @@
end
end;
+*)
+
fun dest_ord f x = is_equal (f x);
+
+
(* Stuff for "equations" ((int*int*int)->num functions). *)
fun tri_equation_cmul c eq =
@@ -1072,6 +1103,9 @@
(* constant monomial is last; this gives an order in diagonalization of the *)
(* quadratic form that will tend to display constants. *)
+(*
+UNUSED
+
fun newton_polytope pol =
let
val vars = poly_variables pol
@@ -1086,6 +1120,8 @@
vars m monomial_1) (rev all')
end;
+*)
+
(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form. *)
local
@@ -1204,9 +1240,9 @@
(* SDPA for problem using block diagonal (i.e. multiple SDPs) *)
-fun sdpa_of_blockproblem comment nblocks blocksizes obj mats =
+fun sdpa_of_blockproblem nblocks blocksizes obj mats =
let val m = length mats - 1
- in "\"" ^ comment ^ "\"\n" ^
+ in
string_of_int m ^ "\n" ^
string_of_int nblocks ^ "\n" ^
(fold1 (fn s => fn t => s^" "^t) (map string_of_int blocksizes)) ^
@@ -1216,6 +1252,14 @@
(1 upto length mats) mats ""
end;
+(* Run prover on a problem in block diagonal form. *)
+
+fun run_blockproblem prover nblocks blocksizes obj mats=
+ parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats))
+
+(*
+UNUSED
+
(* Hence run CSDP on a problem in block diagonal form. *)
fun run_csdp dbg nblocks blocksizes obj mats =
@@ -1248,6 +1292,7 @@
else ());
res)
end;
+*)
(* 3D versions of matrix operations to consider blocks separately. *)
@@ -1270,17 +1315,25 @@
(blocksizes ~~ (1 upto length blocksizes));;
(* FIXME : Get rid of this !!!*)
+local
+ fun tryfind_with msg f [] = error msg
+ | tryfind_with msg f (x::xs) = (f x handle ERROR s => tryfind_with s f xs);
+in
+ fun tryfind f = tryfind_with "tryfind" f
+end
+
+(*
fun tryfind f [] = error "tryfind"
| tryfind f (x::xs) = (f x handle ERROR _ => tryfind f xs);
-
+*)
(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
-
+
local
open RealArith
in
-fun real_positivnullstellensatz_general linf d eqs leqs pol =
+fun real_positivnullstellensatz_general prover linf d eqs leqs pol =
let
val vars = fold_rev (curry (gen_union (op aconvc)) o poly_variables)
(pol::eqs @ map fst leqs) []
@@ -1344,7 +1397,7 @@
itern 1 pvs (fn v => fn i => Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0))
Intfunc.undefined)
val raw_vec = if null pvs then vector_0 0
- else tri_scale_then (csdp nblocks blocksizes) obj mats
+ else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats
fun int_element (d,v) i = Intfunc.tryapplyd v i rat_0
fun cterm_element (d,v) i = Ctermfunc.tryapplyd v i rat_0
@@ -1464,7 +1517,7 @@
fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
in
(* FIXME: Replace tryfind by get_first !! *)
-fun real_nonlinear_prover ctxt =
+fun real_nonlinear_prover prover ctxt =
let
val {add,mul,neg,pow,sub,main} = Normalizer.semiring_normalizers_ord_wrapper ctxt
(valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"}))
@@ -1502,7 +1555,7 @@
let val e = multidegree pol
val k = if e = 0 then 0 else d div e
val eq' = map fst eq
- in tryfind (fn i => (d,i,real_positivnullstellensatz_general false d eq' leq
+ in tryfind (fn i => (d,i,real_positivnullstellensatz_general prover false d eq' leq
(poly_neg(poly_pow pol i))))
(0 upto k)
end
@@ -1563,7 +1616,7 @@
end
in
-fun real_nonlinear_subst_prover ctxt =
+fun real_nonlinear_subst_prover prover ctxt =
let
val {add,mul,neg,pow,sub,main} = Normalizer.semiring_normalizers_ord_wrapper ctxt
(valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"}))
@@ -1596,7 +1649,7 @@
aconvc @{cterm "0::real"}) (map modify eqs),
map modify les,map modify lts)
end)
- handle ERROR _ => real_nonlinear_prover ctxt translator (rev eqs, rev les, rev lts))
+ handle ERROR _ => real_nonlinear_prover prover ctxt translator (rev eqs, rev les, rev lts))
in substfirst
end
@@ -1606,7 +1659,7 @@
(* Overall function. *)
-fun real_sos ctxt t = gen_prover_real_arith ctxt (real_nonlinear_subst_prover ctxt) t;
+fun real_sos prover ctxt t = gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) t;
end;
(* A tactic *)
@@ -1618,7 +1671,7 @@
end
| _ => ([],ct)
-fun core_sos_conv ctxt t = Drule.arg_cong_rule @{cterm Trueprop} (real_sos ctxt (Thm.dest_arg t) RS @{thm Eq_TrueI})
+fun core_sos_conv prover ctxt t = Drule.arg_cong_rule @{cterm Trueprop} (real_sos prover ctxt (Thm.dest_arg t) RS @{thm Eq_TrueI})
val known_sos_constants =
[@{term "op ==>"}, @{term "Trueprop"},
@@ -1654,10 +1707,10 @@
else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
in () end
-fun core_sos_tac ctxt = CSUBGOAL (fn (ct, i) =>
+fun core_sos_tac prover ctxt = CSUBGOAL (fn (ct, i) =>
let val _ = check_sos known_sos_constants ct
val (avs, p) = strip_all ct
- val th = standard (fold_rev forall_intr avs (real_sos ctxt (Thm.dest_arg p)))
+ val th = standard (fold_rev forall_intr avs (real_sos prover ctxt (Thm.dest_arg p)))
in rtac th i end);
fun default_SOME f NONE v = SOME v
@@ -1695,7 +1748,7 @@
fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
-fun sos_tac ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac ctxt
+fun sos_tac prover ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac prover ctxt
end;