# HG changeset patch # User nipkow # Date 1248851209 -7200 # Node ID b642e4813e538fe3adbcbfc4dd593a481307d5b7 # Parent 99711ef9d582249ab73e9e93374fe8c0b8c3b05d# Parent d50f0cb67578330da12d0f460cc148262461e724 Added remote-SOS changes by Philipp Meyer diff -r 99711ef9d582 -r b642e4813e53 CONTRIBUTORS --- 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 ----------------------------- diff -r 99711ef9d582 -r b642e4813e53 lib/scripts/neos/NeosCSDPClient.py --- /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 \n") + sys.exit(1) + +neos=xmlrpclib.Server("http://%s:%d" % (Variables.NEOS_HOST, Variables.NEOS_PORT)) + +xmlfile = open(sys.argv[1],"r") +xml_pre = "\nsdp\ncsdp\nSPARSE_SDPA\n\n\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) + + diff -r 99711ef9d582 -r b642e4813e53 lib/scripts/neos/config.py --- /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 diff -r 99711ef9d582 -r b642e4813e53 src/HOL/Library/Sum_Of_Squares.thy --- 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 \ a < 0" by sos @@ -115,3 +118,4 @@ *) end + diff -r 99711ef9d582 -r b642e4813e53 src/HOL/Library/sos_wrapper.ML --- /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 diff -r 99711ef9d582 -r b642e4813e53 src/HOL/Library/sum_of_squares.ML --- 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;