1.1 --- a/etc/components Wed Aug 05 17:10:10 2009 +0200
1.2 +++ b/etc/components Thu Aug 06 19:51:59 2009 +0200
1.3 @@ -12,4 +12,5 @@
1.4 src/Sequents
1.5 #misc components
1.6 src/HOL/Tools/ATP_Manager
1.7 +src/HOL/Library/Sum_Of_Squares
1.8
2.1 --- a/etc/settings Wed Aug 05 17:10:10 2009 +0200
2.2 +++ b/etc/settings Thu Aug 06 19:51:59 2009 +0200
2.3 @@ -222,7 +222,6 @@
2.4 #JEDIT_JAVA_OPTIONS="-server -Xms128m -Xmx512m"
2.5 JEDIT_OPTIONS="-reuseview -noserver -nobackground"
2.6
2.7 -
2.8 ###
2.9 ### External reasoning tools
2.10 ###
2.11 @@ -274,6 +273,9 @@
2.12 # Jerusat 1.3 (SAT Solver, cf. Isabelle/src/HOL/Tools/sat_solver.ML)
2.13 #JERUSAT_HOME=/usr/local/bin
2.14
2.15 +# CSDP (SDP Solver, cf. Isabelle/src/HOL/Library/Sum_of_Squares/sos_wrapper.ML)
2.16 +#CSDP_EXE=csdp
2.17 +
2.18 # For configuring HOL/Matrix/cplex
2.19 # LP_SOLVER is the default solver. It can be changed during runtime via Cplex.set_solver.
2.20 # First option: use the commercial cplex solver
3.1 --- a/lib/scripts/neos/NeosCSDPClient.py Wed Aug 05 17:10:10 2009 +0200
3.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
3.3 @@ -1,56 +0,0 @@
3.4 -#!/usr/bin/env python
3.5 -import sys
3.6 -import xmlrpclib
3.7 -import time
3.8 -import re
3.9 -
3.10 -from config import Variables
3.11 -
3.12 -if len(sys.argv) < 3 or len(sys.argv) > 3:
3.13 - sys.stderr.write("Usage: NeosCSDPClient <input_filename> <output_filename>\n")
3.14 - sys.exit(1)
3.15 -
3.16 -neos=xmlrpclib.Server("http://%s:%d" % (Variables.NEOS_HOST, Variables.NEOS_PORT))
3.17 -
3.18 -xmlfile = open(sys.argv[1],"r")
3.19 -xml_pre = "<document>\n<category>sdp</category>\n<solver>csdp</solver>\n<inputMethod>SPARSE_SDPA</inputMethod>\n<dat><![CDATA["
3.20 -xml_post = "]]></dat>\n</document>\n"
3.21 -xml = xml_pre
3.22 -buffer = 1
3.23 -while buffer:
3.24 - buffer = xmlfile.read()
3.25 - xml += buffer
3.26 -xmlfile.close()
3.27 -xml += xml_post
3.28 -
3.29 -(jobNumber,password) = neos.submitJob(xml)
3.30 -
3.31 -if jobNumber == 0:
3.32 - sys.stdout.write("error submitting job: %s" % password)
3.33 - sys.exit(-1)
3.34 -else:
3.35 - sys.stdout.write("jobNumber = %d\tpassword = %s\n" % (jobNumber,password))
3.36 -
3.37 -offset=0
3.38 -status="Waiting"
3.39 -while status == "Running" or status=="Waiting":
3.40 - time.sleep(1)
3.41 - (msg,offset) = neos.getIntermediateResults(jobNumber,password,offset)
3.42 - sys.stdout.write(msg.data)
3.43 - status = neos.getJobStatus(jobNumber, password)
3.44 -
3.45 -msg = neos.getFinalResults(jobNumber, password).data
3.46 -result = msg.split("Solution:")
3.47 -
3.48 -sys.stdout.write(result[0])
3.49 -if len(result) > 1:
3.50 - plain_msg = result[1].strip()
3.51 - if plain_msg != "":
3.52 - output = open(sys.argv[2],"w")
3.53 - output.write(plain_msg)
3.54 - output.close()
3.55 - sys.exit(0)
3.56 -
3.57 -sys.exit(2)
3.58 -
3.59 -
4.1 --- a/lib/scripts/neos/config.py Wed Aug 05 17:10:10 2009 +0200
4.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
4.3 @@ -1,3 +0,0 @@
4.4 -class Variables:
4.5 - NEOS_HOST="neos.mcs.anl.gov"
4.6 - NEOS_PORT=3332
5.1 --- a/src/HOL/IsaMakefile Wed Aug 05 17:10:10 2009 +0200
5.2 +++ b/src/HOL/IsaMakefile Thu Aug 06 19:51:59 2009 +0200
5.3 @@ -316,49 +316,45 @@
5.4
5.5 HOL-Library: HOL $(LOG)/HOL-Library.gz
5.6
5.7 -
5.8 $(LOG)/HOL-Library.gz: $(OUT)/HOL Library/SetsAndFunctions.thy \
5.9 - Library/Abstract_Rat.thy \
5.10 - Library/BigO.thy Library/ContNotDenum.thy Library/Efficient_Nat.thy \
5.11 - Library/Euclidean_Space.thy Library/Sum_Of_Squares.thy Library/positivstellensatz.ML \
5.12 - Library/Fset.thy Library/Convex_Euclidean_Space.thy \
5.13 - Library/sum_of_squares.ML Library/Glbs.thy Library/normarith.ML \
5.14 - Library/Executable_Set.thy Library/Infinite_Set.thy \
5.15 - Library/FuncSet.thy Library/Permutations.thy Library/Determinants.thy\
5.16 - Library/Bit.thy Library/Topology_Euclidean_Space.thy \
5.17 - Library/Finite_Cartesian_Product.thy \
5.18 - Library/FrechetDeriv.thy Library/Fraction_Field.thy\
5.19 - Library/Fundamental_Theorem_Algebra.thy \
5.20 - Library/Inner_Product.thy Library/Kleene_Algebra.thy Library/Lattice_Syntax.thy \
5.21 - Library/Legacy_GCD.thy \
5.22 - Library/Library.thy Library/List_Prefix.thy Library/List_Set.thy Library/State_Monad.thy \
5.23 - Library/Nat_Int_Bij.thy Library/Multiset.thy Library/Permutation.thy \
5.24 - Library/Primes.thy Library/Pocklington.thy Library/Quotient.thy \
5.25 - Library/Quicksort.thy Library/Nat_Infinity.thy Library/Word.thy \
5.26 - Library/README.html Library/Continuity.thy Library/Order_Relation.thy \
5.27 - Library/Nested_Environment.thy Library/Ramsey.thy Library/Zorn.thy \
5.28 - Library/Library/ROOT.ML Library/Library/document/root.tex \
5.29 - Library/Library/document/root.bib Library/While_Combinator.thy \
5.30 - Library/Product_ord.thy Library/Char_nat.thy Library/Char_ord.thy \
5.31 - Library/Option_ord.thy Library/Sublist_Order.thy \
5.32 - Library/List_lexord.thy Library/Commutative_Ring.thy \
5.33 - Library/comm_ring.ML Library/Coinductive_List.thy \
5.34 - Library/AssocList.thy Library/Formal_Power_Series.thy \
5.35 - Library/Binomial.thy Library/Eval_Witness.thy \
5.36 - Library/Code_Char.thy \
5.37 + Library/Abstract_Rat.thy Library/BigO.thy Library/ContNotDenum.thy \
5.38 + Library/Efficient_Nat.thy Library/Euclidean_Space.thy \
5.39 + Library/Sum_Of_Squares.thy Library/Sum_Of_Squares/sos_wrapper.ML \
5.40 + Library/Sum_Of_Squares/sum_of_squares.ML Library/Fset.thy \
5.41 + Library/Convex_Euclidean_Space.thy Library/Glbs.thy \
5.42 + Library/normarith.ML Library/Executable_Set.thy \
5.43 + Library/Infinite_Set.thy Library/FuncSet.thy \
5.44 + Library/Permutations.thy Library/Determinants.thy Library/Bit.thy \
5.45 + Library/Topology_Euclidean_Space.thy \
5.46 + Library/Finite_Cartesian_Product.thy Library/FrechetDeriv.thy \
5.47 + Library/Fraction_Field.thy Library/Fundamental_Theorem_Algebra.thy \
5.48 + Library/Inner_Product.thy Library/Kleene_Algebra.thy \
5.49 + Library/Lattice_Syntax.thy Library/Legacy_GCD.thy \
5.50 + Library/Library.thy Library/List_Prefix.thy Library/List_Set.thy \
5.51 + Library/State_Monad.thy Library/Nat_Int_Bij.thy Library/Multiset.thy \
5.52 + Library/Permutation.thy Library/Primes.thy Library/Pocklington.thy \
5.53 + Library/Quotient.thy Library/Quicksort.thy Library/Nat_Infinity.thy \
5.54 + Library/Word.thy Library/README.html Library/Continuity.thy \
5.55 + Library/Order_Relation.thy Library/Nested_Environment.thy \
5.56 + Library/Ramsey.thy Library/Zorn.thy Library/Library/ROOT.ML \
5.57 + Library/Library/document/root.tex Library/Library/document/root.bib \
5.58 + Library/While_Combinator.thy Library/Product_ord.thy \
5.59 + Library/Char_nat.thy Library/Char_ord.thy Library/Option_ord.thy \
5.60 + Library/Sublist_Order.thy Library/List_lexord.thy \
5.61 + Library/Commutative_Ring.thy Library/comm_ring.ML \
5.62 + Library/Coinductive_List.thy Library/AssocList.thy \
5.63 + Library/Formal_Power_Series.thy Library/Binomial.thy \
5.64 + Library/Eval_Witness.thy Library/Code_Char.thy \
5.65 Library/Code_Char_chr.thy Library/Code_Integer.thy \
5.66 Library/Mapping.thy Library/Numeral_Type.thy Library/Reflection.thy \
5.67 Library/Boolean_Algebra.thy Library/Countable.thy \
5.68 Library/Diagonalize.thy Library/RBT.thy Library/Univ_Poly.thy \
5.69 - Library/Poly_Deriv.thy \
5.70 - Library/Polynomial.thy \
5.71 - Library/Preorder.thy \
5.72 - Library/Product_plus.thy \
5.73 - Library/Product_Vector.thy \
5.74 - Library/Tree.thy \
5.75 - Library/Enum.thy Library/Float.thy $(SRC)/Tools/float.ML $(SRC)/HOL/Tools/float_arith.ML \
5.76 - Library/reify_data.ML Library/reflection.ML \
5.77 - Library/LaTeXsugar.thy Library/OptionalSugar.thy
5.78 + Library/Poly_Deriv.thy Library/Polynomial.thy Library/Preorder.thy \
5.79 + Library/Product_plus.thy Library/Product_Vector.thy Library/Tree.thy \
5.80 + Library/Enum.thy Library/Float.thy $(SRC)/Tools/float.ML \
5.81 + $(SRC)/HOL/Tools/float_arith.ML Library/positivstellensatz.ML \
5.82 + Library/reify_data.ML Library/reflection.ML Library/LaTeXsugar.thy \
5.83 + Library/OptionalSugar.thy
5.84 @cd Library; $(ISABELLE_TOOL) usedir $(OUT)/HOL Library
5.85
5.86
6.1 --- a/src/HOL/Library/Sum_Of_Squares.thy Wed Aug 05 17:10:10 2009 +0200
6.2 +++ b/src/HOL/Library/Sum_Of_Squares.thy Thu Aug 06 19:51:59 2009 +0200
6.3 @@ -2,9 +2,9 @@
6.4 Author: Amine Chaieb, University of Cambridge
6.5
6.6 In order to use the method sos, call it with (sos remote_csdp) to use the remote solver
6.7 -or install CSDP (https://projects.coin-or.org/Csdp/), put the executable csdp on your path,
6.8 -and call it with (sos csdp). By default, sos calls remote_csdp. This can take of the order
6.9 -of a minute for one sos call, because sos calls CSDP repeatedly.
6.10 +or install CSDP (https://projects.coin-or.org/Csdp/), set the Isabelle environment
6.11 +variable CSDP_EXE and call it with (sos csdp). By default, sos calls remote_csdp.
6.12 +This can take of the order of a minute for one sos call, because sos calls CSDP repeatedly.
6.13 If you install CSDP locally, sos calls typically takes only a few seconds.
6.14
6.15 *)
6.16 @@ -13,11 +13,19 @@
6.17 multiplication and ordering using semidefinite programming*}
6.18
6.19 theory Sum_Of_Squares
6.20 - imports Complex_Main (* "~~/src/HOL/Decision_Procs/Dense_Linear_Order" *)
6.21 - uses "positivstellensatz.ML" "sum_of_squares.ML" "sos_wrapper.ML"
6.22 - begin
6.23 +imports Complex_Main (* "~~/src/HOL/Decision_Procs/Dense_Linear_Order" *)
6.24 +uses
6.25 + ("positivstellensatz.ML")
6.26 + ("Sum_Of_Squares/sum_of_squares.ML")
6.27 + ("Sum_Of_Squares/sos_wrapper.ML")
6.28 +begin
6.29
6.30 (* setup sos tactic *)
6.31 +
6.32 +use "positivstellensatz.ML"
6.33 +use "Sum_Of_Squares/sum_of_squares.ML"
6.34 +use "Sum_Of_Squares/sos_wrapper.ML"
6.35 +
6.36 setup SosWrapper.setup
6.37
6.38 text{* Tests -- commented since they work only when csdp is installed or take too long with remote csdps *}
7.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
7.2 +++ b/src/HOL/Library/Sum_Of_Squares/etc/settings Thu Aug 06 19:51:59 2009 +0200
7.3 @@ -0,0 +1,1 @@
7.4 +ISABELLE_SUM_OF_SQUARES="$COMPONENT"
8.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
8.2 +++ b/src/HOL/Library/Sum_Of_Squares/neos_csdp_client Thu Aug 06 19:51:59 2009 +0200
8.3 @@ -0,0 +1,65 @@
8.4 +#!/usr/bin/env python
8.5 +import sys
8.6 +import xmlrpclib
8.7 +import time
8.8 +import re
8.9 +
8.10 +# Neos server config
8.11 +NEOS_HOST="neos.mcs.anl.gov"
8.12 +NEOS_PORT=3332
8.13 +
8.14 +if len(sys.argv) < 3 or len(sys.argv) > 3:
8.15 + sys.stderr.write("Usage: NeosCSDPClient <input_filename> <output_filename>\n")
8.16 + sys.exit(1)
8.17 +
8.18 +neos=xmlrpclib.Server("http://%s:%d" % (NEOS_HOST, NEOS_PORT))
8.19 +
8.20 +inputfile = open(sys.argv[1],"r")
8.21 +xml_pre = "<document>\n<category>sdp</category>\n<solver>csdp</solver>\n<inputMethod>SPARSE_SDPA</inputMethod>\n<dat><![CDATA["
8.22 +xml_post = "]]></dat>\n</document>\n"
8.23 +xml = xml_pre
8.24 +buffer = 1
8.25 +while buffer:
8.26 + buffer = inputfile.read()
8.27 + xml += buffer
8.28 +inputfile.close()
8.29 +xml += xml_post
8.30 +
8.31 +(jobNumber,password) = neos.submitJob(xml)
8.32 +
8.33 +if jobNumber == 0:
8.34 + sys.stdout.write("error submitting job: %s" % password)
8.35 + sys.exit(20)
8.36 +else:
8.37 + sys.stdout.write("jobNumber = %d\tpassword = %s\n" % (jobNumber,password))
8.38 +
8.39 +offset=0
8.40 +messages = ""
8.41 +status="Waiting"
8.42 +while status == "Running" or status=="Waiting":
8.43 + time.sleep(1)
8.44 + (msg,offset) = neos.getIntermediateResults(jobNumber,password,offset)
8.45 + messages += msg.data
8.46 + sys.stdout.write(msg.data)
8.47 + status = neos.getJobStatus(jobNumber, password)
8.48 +
8.49 +msg = neos.getFinalResults(jobNumber, password).data
8.50 +sys.stdout.write("---------- Begin CSDP Output -------------\n");
8.51 +sys.stdout.write(msg)
8.52 +
8.53 +# extract solution
8.54 +result = msg.split("Solution:")
8.55 +if len(result) > 1:
8.56 + output = open(sys.argv[2],"w")
8.57 + output.write(result[1].strip())
8.58 + output.close()
8.59 +
8.60 +# extract return code
8.61 +p = re.compile(r"^Error: Command exited with non-zero status (\d+)$", re.MULTILINE)
8.62 +m = p.search(messages)
8.63 +if m:
8.64 + sys.exit(int(m.group(1)))
8.65 +else:
8.66 + sys.exit(0)
8.67 +
8.68 +
9.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
9.2 +++ b/src/HOL/Library/Sum_Of_Squares/sos_wrapper.ML Thu Aug 06 19:51:59 2009 +0200
9.3 @@ -0,0 +1,140 @@
9.4 +(* Title: sos_wrapper.ML
9.5 + Author: Philipp Meyer, TU Muenchen
9.6 +
9.7 +Added functionality for sums of squares, e.g. calling a remote prover
9.8 +*)
9.9 +
9.10 +signature SOS_WRAPPER =
9.11 +sig
9.12 +
9.13 + datatype prover_result = Success | Failure | Error
9.14 +
9.15 + val setup: theory -> theory
9.16 + val destdir: string ref
9.17 +end
9.18 +
9.19 +structure SosWrapper : SOS_WRAPPER=
9.20 +struct
9.21 +
9.22 +datatype prover_result = Success | Failure | Error
9.23 +fun str_of_result Success = "Success"
9.24 + | str_of_result Failure = "Failure"
9.25 + | str_of_result Error = "Error"
9.26 +
9.27 +(*** output control ***)
9.28 +
9.29 +fun debug s = if ! Sos.debugging then Output.writeln s else ()
9.30 +val write = Output.priority
9.31 +
9.32 +(*** calling provers ***)
9.33 +
9.34 +val destdir = ref ""
9.35 +
9.36 +fun filename dir name =
9.37 + let
9.38 + val probfile = Path.basic (name ^ serial_string ())
9.39 + val dir_path = Path.explode dir
9.40 + in
9.41 + if dir = "" then
9.42 + File.tmp_path probfile
9.43 + else
9.44 + if File.exists dir_path then
9.45 + Path.append dir_path probfile
9.46 + else
9.47 + error ("No such directory: " ^ dir)
9.48 + end
9.49 +
9.50 +fun run_solver name cmd find_failure input =
9.51 + let
9.52 + val _ = write ("Calling solver: " ^ name)
9.53 +
9.54 + (* create input file *)
9.55 + val dir = ! destdir
9.56 + val input_file = filename dir "sos_in"
9.57 + val _ = File.write input_file input
9.58 +
9.59 + (* call solver *)
9.60 + val output_file = filename dir "sos_out"
9.61 + val (output, rv) = system_out (
9.62 + if File.exists cmd then space_implode " "
9.63 + [File.shell_path cmd, File.platform_path input_file, File.platform_path output_file]
9.64 + else error ("Bad executable: " ^ File.shell_path cmd))
9.65 +
9.66 + (* read and analysize output *)
9.67 + val (res, res_msg) = find_failure rv
9.68 + val result = if File.exists output_file then File.read output_file else ""
9.69 +
9.70 + (* remove temporary files *)
9.71 + val _ = if dir = "" then
9.72 + (File.rm input_file ; if File.exists output_file then File.rm output_file else ())
9.73 + else ()
9.74 +
9.75 + val _ = debug ("Solver output:\n" ^ output)
9.76 +
9.77 + val _ = write (str_of_result res ^ ": " ^ res_msg)
9.78 + in
9.79 + case res of
9.80 + Success => result
9.81 + | Failure => raise Sos.Failure res_msg
9.82 + | Error => error ("Prover failed: " ^ res_msg)
9.83 + end
9.84 +
9.85 +(*** various provers ***)
9.86 +
9.87 +(* local csdp client *)
9.88 +
9.89 +fun find_csdp_failure rv =
9.90 + case rv of
9.91 + 0 => (Success, "SDP solved")
9.92 + | 1 => (Failure, "SDP is primal infeasible")
9.93 + | 2 => (Failure, "SDP is dual infeasible")
9.94 + | 3 => (Success, "SDP solved with reduced accuracy")
9.95 + | 4 => (Failure, "Maximum iterations reached")
9.96 + | 5 => (Failure, "Stuck at edge of primal feasibility")
9.97 + | 6 => (Failure, "Stuck at edge of dual infeasibility")
9.98 + | 7 => (Failure, "Lack of progress")
9.99 + | 8 => (Failure, "X, Z, or O was singular")
9.100 + | 9 => (Failure, "Detected NaN or Inf values")
9.101 + | _ => (Error, "return code is " ^ string_of_int rv)
9.102 +
9.103 +val csdp = ("$CSDP_EXE", find_csdp_failure)
9.104 +
9.105 +(* remote neos server *)
9.106 +
9.107 +fun find_neos_failure rv =
9.108 + case rv of
9.109 + 20 => (Error, "error submitting job")
9.110 + | 21 => (Error, "no solution")
9.111 + | _ => find_csdp_failure rv
9.112 +
9.113 +val neos_csdp = ("$ISABELLE_SUM_OF_SQUARES/neos_csdp_client", find_neos_failure)
9.114 +
9.115 +(* save provers in table *)
9.116 +
9.117 +val provers =
9.118 + Symtab.make [("remote_csdp", neos_csdp),("csdp", csdp)]
9.119 +
9.120 +fun get_prover name =
9.121 + case Symtab.lookup provers name of
9.122 + SOME prover => prover
9.123 + | NONE => error ("unknown prover: " ^ name)
9.124 +
9.125 +fun call_solver name =
9.126 + let
9.127 + val (cmd, find_failure) = get_prover name
9.128 + in
9.129 + run_solver name (Path.explode cmd) find_failure
9.130 + end
9.131 +
9.132 +(* setup tactic *)
9.133 +
9.134 +val def_solver = "remote_csdp"
9.135 +
9.136 +fun sos_solver name = (SIMPLE_METHOD' o (Sos.sos_tac (call_solver name)))
9.137 +
9.138 +val sos_method = Scan.optional (Scan.lift OuterParse.xname) def_solver >> sos_solver
9.139 +
9.140 +val setup = Method.setup @{binding sos} sos_method
9.141 + "Prove universal problems over the reals using sums of squares"
9.142 +
9.143 +end
10.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
10.2 +++ b/src/HOL/Library/Sum_Of_Squares/sum_of_squares.ML Thu Aug 06 19:51:59 2009 +0200
10.3 @@ -0,0 +1,1512 @@
10.4 +(* Title: sum_of_squares.ML
10.5 + Authors: Amine Chaieb, University of Cambridge
10.6 + Philipp Meyer, TU Muenchen
10.7 +
10.8 +A tactic for proving nonlinear inequalities
10.9 +*)
10.10 +
10.11 +signature SOS =
10.12 +sig
10.13 +
10.14 + val sos_tac : (string -> string) -> Proof.context -> int -> Tactical.tactic
10.15 +
10.16 + val debugging : bool ref;
10.17 +
10.18 + exception Failure of string;
10.19 +end
10.20 +
10.21 +structure Sos : SOS =
10.22 +struct
10.23 +
10.24 +val rat_0 = Rat.zero;
10.25 +val rat_1 = Rat.one;
10.26 +val rat_2 = Rat.two;
10.27 +val rat_10 = Rat.rat_of_int 10;
10.28 +val rat_1_2 = rat_1 // rat_2;
10.29 +val max = curry IntInf.max;
10.30 +val min = curry IntInf.min;
10.31 +
10.32 +val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int;
10.33 +val numerator_rat = Rat.quotient_of_rat #> fst #> Rat.rat_of_int;
10.34 +fun int_of_rat a =
10.35 + case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int";
10.36 +fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y));
10.37 +
10.38 +fun rat_pow r i =
10.39 + let fun pow r i =
10.40 + if i = 0 then rat_1 else
10.41 + let val d = pow r (i div 2)
10.42 + in d */ d */ (if i mod 2 = 0 then rat_1 else r)
10.43 + end
10.44 + in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end;
10.45 +
10.46 +fun round_rat r =
10.47 + let val (a,b) = Rat.quotient_of_rat (Rat.abs r)
10.48 + val d = a div b
10.49 + val s = if r </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
10.50 + val x2 = 2 * (a - (b * d))
10.51 + in s (if x2 >= b then d + 1 else d) end
10.52 +
10.53 +val abs_rat = Rat.abs;
10.54 +val pow2 = rat_pow rat_2;
10.55 +val pow10 = rat_pow rat_10;
10.56 +
10.57 +val debugging = ref false;
10.58 +
10.59 +exception Sanity;
10.60 +
10.61 +exception Unsolvable;
10.62 +
10.63 +exception Failure of string;
10.64 +
10.65 +(* Turn a rational into a decimal string with d sig digits. *)
10.66 +
10.67 +local
10.68 +fun normalize y =
10.69 + if abs_rat y </ (rat_1 // rat_10) then normalize (rat_10 */ y) - 1
10.70 + else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1
10.71 + else 0
10.72 + in
10.73 +fun decimalize d x =
10.74 + if x =/ rat_0 then "0.0" else
10.75 + let
10.76 + val y = Rat.abs x
10.77 + val e = normalize y
10.78 + val z = pow10(~ e) */ y +/ rat_1
10.79 + val k = int_of_rat (round_rat(pow10 d */ z))
10.80 + in (if x </ rat_0 then "-0." else "0.") ^
10.81 + implode(tl(explode(string_of_int k))) ^
10.82 + (if e = 0 then "" else "e"^string_of_int e)
10.83 + end
10.84 +end;
10.85 +
10.86 +(* Iterations over numbers, and lists indexed by numbers. *)
10.87 +
10.88 +fun itern k l f a =
10.89 + case l of
10.90 + [] => a
10.91 + | h::t => itern (k + 1) t f (f h k a);
10.92 +
10.93 +fun iter (m,n) f a =
10.94 + if n < m then a
10.95 + else iter (m+1,n) f (f m a);
10.96 +
10.97 +(* The main types. *)
10.98 +
10.99 +fun strict_ord ord (x,y) = case ord (x,y) of LESS => LESS | _ => GREATER
10.100 +
10.101 +structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord);
10.102 +
10.103 +type vector = int* Rat.rat Intfunc.T;
10.104 +
10.105 +type matrix = (int*int)*(Rat.rat Intpairfunc.T);
10.106 +
10.107 +type monomial = int Ctermfunc.T;
10.108 +
10.109 +val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t))
10.110 + fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2)
10.111 +structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord)
10.112 +
10.113 +type poly = Rat.rat Monomialfunc.T;
10.114 +
10.115 + fun iszero (k,r) = r =/ rat_0;
10.116 +
10.117 +fun fold_rev2 f l1 l2 b =
10.118 + case (l1,l2) of
10.119 + ([],[]) => b
10.120 + | (h1::t1,h2::t2) => f h1 h2 (fold_rev2 f t1 t2 b)
10.121 + | _ => error "fold_rev2";
10.122 +
10.123 +(* Vectors. Conventionally indexed 1..n. *)
10.124 +
10.125 +fun vector_0 n = (n,Intfunc.undefined):vector;
10.126 +
10.127 +fun dim (v:vector) = fst v;
10.128 +
10.129 +fun vector_const c n =
10.130 + if c =/ rat_0 then vector_0 n
10.131 + else (n,fold_rev (fn k => Intfunc.update (k,c)) (1 upto n) Intfunc.undefined) :vector;
10.132 +
10.133 +val vector_1 = vector_const rat_1;
10.134 +
10.135 +fun vector_cmul c (v:vector) =
10.136 + let val n = dim v
10.137 + in if c =/ rat_0 then vector_0 n
10.138 + else (n,Intfunc.mapf (fn x => c */ x) (snd v))
10.139 + end;
10.140 +
10.141 +fun vector_neg (v:vector) = (fst v,Intfunc.mapf Rat.neg (snd v)) :vector;
10.142 +
10.143 +fun vector_add (v1:vector) (v2:vector) =
10.144 + let val m = dim v1
10.145 + val n = dim v2
10.146 + in if m <> n then error "vector_add: incompatible dimensions"
10.147 + else (n,Intfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd v1) (snd v2)) :vector
10.148 + end;
10.149 +
10.150 +fun vector_sub v1 v2 = vector_add v1 (vector_neg v2);
10.151 +
10.152 +fun vector_dot (v1:vector) (v2:vector) =
10.153 + let val m = dim v1
10.154 + val n = dim v2
10.155 + in if m <> n then error "vector_dot: incompatible dimensions"
10.156 + else Intfunc.fold (fn (i,x) => fn a => x +/ a)
10.157 + (Intfunc.combine (curry op */) (fn x => x =/ rat_0) (snd v1) (snd v2)) rat_0
10.158 + end;
10.159 +
10.160 +fun vector_of_list l =
10.161 + let val n = length l
10.162 + in (n,fold_rev2 (curry Intfunc.update) (1 upto n) l Intfunc.undefined) :vector
10.163 + end;
10.164 +
10.165 +(* Matrices; again rows and columns indexed from 1. *)
10.166 +
10.167 +fun matrix_0 (m,n) = ((m,n),Intpairfunc.undefined):matrix;
10.168 +
10.169 +fun dimensions (m:matrix) = fst m;
10.170 +
10.171 +fun matrix_const c (mn as (m,n)) =
10.172 + if m <> n then error "matrix_const: needs to be square"
10.173 + else if c =/ rat_0 then matrix_0 mn
10.174 + else (mn,fold_rev (fn k => Intpairfunc.update ((k,k), c)) (1 upto n) Intpairfunc.undefined) :matrix;;
10.175 +
10.176 +val matrix_1 = matrix_const rat_1;
10.177 +
10.178 +fun matrix_cmul c (m:matrix) =
10.179 + let val (i,j) = dimensions m
10.180 + in if c =/ rat_0 then matrix_0 (i,j)
10.181 + else ((i,j),Intpairfunc.mapf (fn x => c */ x) (snd m))
10.182 + end;
10.183 +
10.184 +fun matrix_neg (m:matrix) =
10.185 + (dimensions m, Intpairfunc.mapf Rat.neg (snd m)) :matrix;
10.186 +
10.187 +fun matrix_add (m1:matrix) (m2:matrix) =
10.188 + let val d1 = dimensions m1
10.189 + val d2 = dimensions m2
10.190 + in if d1 <> d2
10.191 + then error "matrix_add: incompatible dimensions"
10.192 + else (d1,Intpairfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd m1) (snd m2)) :matrix
10.193 + end;;
10.194 +
10.195 +fun matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2);
10.196 +
10.197 +fun row k (m:matrix) =
10.198 + let val (i,j) = dimensions m
10.199 + in (j,
10.200 + Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then Intfunc.update (j,c) a else a) (snd m) Intfunc.undefined ) : vector
10.201 + end;
10.202 +
10.203 +fun column k (m:matrix) =
10.204 + let val (i,j) = dimensions m
10.205 + in (i,
10.206 + Intpairfunc.fold (fn ((i,j), c) => fn a => if j = k then Intfunc.update (i,c) a else a) (snd m) Intfunc.undefined)
10.207 + : vector
10.208 + end;
10.209 +
10.210 +fun transp (m:matrix) =
10.211 + let val (i,j) = dimensions m
10.212 + in
10.213 + ((j,i),Intpairfunc.fold (fn ((i,j), c) => fn a => Intpairfunc.update ((j,i), c) a) (snd m) Intpairfunc.undefined) :matrix
10.214 + end;
10.215 +
10.216 +fun diagonal (v:vector) =
10.217 + let val n = dim v
10.218 + in ((n,n),Intfunc.fold (fn (i, c) => fn a => Intpairfunc.update ((i,i), c) a) (snd v) Intpairfunc.undefined) : matrix
10.219 + end;
10.220 +
10.221 +fun matrix_of_list l =
10.222 + let val m = length l
10.223 + in if m = 0 then matrix_0 (0,0) else
10.224 + let val n = length (hd l)
10.225 + in ((m,n),itern 1 l (fn v => fn i => itern 1 v (fn c => fn j => Intpairfunc.update ((i,j), c))) Intpairfunc.undefined)
10.226 + end
10.227 + end;
10.228 +
10.229 +(* Monomials. *)
10.230 +
10.231 +fun monomial_eval assig (m:monomial) =
10.232 + Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (Ctermfunc.apply assig x) k)
10.233 + m rat_1;
10.234 +val monomial_1 = (Ctermfunc.undefined:monomial);
10.235 +
10.236 +fun monomial_var x = Ctermfunc.onefunc (x, 1) :monomial;
10.237 +
10.238 +val (monomial_mul:monomial->monomial->monomial) =
10.239 + Ctermfunc.combine (curry op +) (K false);
10.240 +
10.241 +fun monomial_pow (m:monomial) k =
10.242 + if k = 0 then monomial_1
10.243 + else Ctermfunc.mapf (fn x => k * x) m;
10.244 +
10.245 +fun monomial_divides (m1:monomial) (m2:monomial) =
10.246 + Ctermfunc.fold (fn (x, k) => fn a => Ctermfunc.tryapplyd m2 x 0 >= k andalso a) m1 true;;
10.247 +
10.248 +fun monomial_div (m1:monomial) (m2:monomial) =
10.249 + let val m = Ctermfunc.combine (curry op +)
10.250 + (fn x => x = 0) m1 (Ctermfunc.mapf (fn x => ~ x) m2)
10.251 + in if Ctermfunc.fold (fn (x, k) => fn a => k >= 0 andalso a) m true then m
10.252 + else error "monomial_div: non-divisible"
10.253 + end;
10.254 +
10.255 +fun monomial_degree x (m:monomial) =
10.256 + Ctermfunc.tryapplyd m x 0;;
10.257 +
10.258 +fun monomial_lcm (m1:monomial) (m2:monomial) =
10.259 + fold_rev (fn x => Ctermfunc.update (x, max (monomial_degree x m1) (monomial_degree x m2)))
10.260 + (gen_union (is_equal o cterm_ord) (Ctermfunc.dom m1, Ctermfunc.dom m2)) (Ctermfunc.undefined :monomial);
10.261 +
10.262 +fun monomial_multidegree (m:monomial) =
10.263 + Ctermfunc.fold (fn (x, k) => fn a => k + a) m 0;;
10.264 +
10.265 +fun monomial_variables m = Ctermfunc.dom m;;
10.266 +
10.267 +(* Polynomials. *)
10.268 +
10.269 +fun eval assig (p:poly) =
10.270 + Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0;
10.271 +
10.272 +val poly_0 = (Monomialfunc.undefined:poly);
10.273 +
10.274 +fun poly_isconst (p:poly) =
10.275 + Monomialfunc.fold (fn (m, c) => fn a => Ctermfunc.is_undefined m andalso a) p true;
10.276 +
10.277 +fun poly_var x = Monomialfunc.onefunc (monomial_var x,rat_1) :poly;
10.278 +
10.279 +fun poly_const c =
10.280 + if c =/ rat_0 then poly_0 else Monomialfunc.onefunc(monomial_1, c);
10.281 +
10.282 +fun poly_cmul c (p:poly) =
10.283 + if c =/ rat_0 then poly_0
10.284 + else Monomialfunc.mapf (fn x => c */ x) p;
10.285 +
10.286 +fun poly_neg (p:poly) = (Monomialfunc.mapf Rat.neg p :poly);;
10.287 +
10.288 +fun poly_add (p1:poly) (p2:poly) =
10.289 + (Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2 :poly);
10.290 +
10.291 +fun poly_sub p1 p2 = poly_add p1 (poly_neg p2);
10.292 +
10.293 +fun poly_cmmul (c,m) (p:poly) =
10.294 + if c =/ rat_0 then poly_0
10.295 + else if Ctermfunc.is_undefined m
10.296 + then Monomialfunc.mapf (fn d => c */ d) p
10.297 + else Monomialfunc.fold (fn (m', d) => fn a => (Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0;
10.298 +
10.299 +fun poly_mul (p1:poly) (p2:poly) =
10.300 + Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0;
10.301 +
10.302 +fun poly_div (p1:poly) (p2:poly) =
10.303 + if not(poly_isconst p2)
10.304 + then error "poly_div: non-constant" else
10.305 + let val c = eval Ctermfunc.undefined p2
10.306 + in if c =/ rat_0 then error "poly_div: division by zero"
10.307 + else poly_cmul (Rat.inv c) p1
10.308 + end;
10.309 +
10.310 +fun poly_square p = poly_mul p p;
10.311 +
10.312 +fun poly_pow p k =
10.313 + if k = 0 then poly_const rat_1
10.314 + else if k = 1 then p
10.315 + else let val q = poly_square(poly_pow p (k div 2)) in
10.316 + if k mod 2 = 1 then poly_mul p q else q end;
10.317 +
10.318 +fun poly_exp p1 p2 =
10.319 + if not(poly_isconst p2)
10.320 + then error "poly_exp: not a constant"
10.321 + else poly_pow p1 (int_of_rat (eval Ctermfunc.undefined p2));
10.322 +
10.323 +fun degree x (p:poly) =
10.324 + Monomialfunc.fold (fn (m,c) => fn a => max (monomial_degree x m) a) p 0;
10.325 +
10.326 +fun multidegree (p:poly) =
10.327 + Monomialfunc.fold (fn (m, c) => fn a => max (monomial_multidegree m) a) p 0;
10.328 +
10.329 +fun poly_variables (p:poly) =
10.330 + sort cterm_ord (Monomialfunc.fold_rev (fn (m, c) => curry (gen_union (is_equal o cterm_ord)) (monomial_variables m)) p []);;
10.331 +
10.332 +(* Order monomials for human presentation. *)
10.333 +
10.334 +fun cterm_ord (t,t') = TermOrd.fast_term_ord (term_of t, term_of t');
10.335 +
10.336 +val humanorder_varpow = prod_ord cterm_ord (rev_order o int_ord);
10.337 +
10.338 +local
10.339 + fun ord (l1,l2) = case (l1,l2) of
10.340 + (_,[]) => LESS
10.341 + | ([],_) => GREATER
10.342 + | (h1::t1,h2::t2) =>
10.343 + (case humanorder_varpow (h1, h2) of
10.344 + LESS => LESS
10.345 + | EQUAL => ord (t1,t2)
10.346 + | GREATER => GREATER)
10.347 +in fun humanorder_monomial m1 m2 =
10.348 + ord (sort humanorder_varpow (Ctermfunc.graph m1),
10.349 + sort humanorder_varpow (Ctermfunc.graph m2))
10.350 +end;
10.351 +
10.352 +fun fold1 f l = case l of
10.353 + [] => error "fold1"
10.354 + | [x] => x
10.355 + | (h::t) => f h (fold1 f t);
10.356 +
10.357 +(* Conversions to strings. *)
10.358 +
10.359 +fun string_of_vector min_size max_size (v:vector) =
10.360 + let val n_raw = dim v
10.361 + in if n_raw = 0 then "[]" else
10.362 + let
10.363 + val n = max min_size (min n_raw max_size)
10.364 + val xs = map (Rat.string_of_rat o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
10.365 + in "[" ^ fold1 (fn s => fn t => s ^ ", " ^ t) xs ^
10.366 + (if n_raw > max_size then ", ...]" else "]")
10.367 + end
10.368 + end;
10.369 +
10.370 +fun string_of_matrix max_size (m:matrix) =
10.371 + let
10.372 + val (i_raw,j_raw) = dimensions m
10.373 + val i = min max_size i_raw
10.374 + val j = min max_size j_raw
10.375 + val rstr = map (fn k => string_of_vector j j (row k m)) (1 upto i)
10.376 + in "["^ fold1 (fn s => fn t => s^";\n "^t) rstr ^
10.377 + (if j > max_size then "\n ...]" else "]")
10.378 + end;
10.379 +
10.380 +fun string_of_term t =
10.381 + case t of
10.382 + a$b => "("^(string_of_term a)^" "^(string_of_term b)^")"
10.383 + | Abs x =>
10.384 + let val (xn, b) = Term.dest_abs x
10.385 + in "(\\"^xn^"."^(string_of_term b)^")"
10.386 + end
10.387 + | Const(s,_) => s
10.388 + | Free (s,_) => s
10.389 + | Var((s,_),_) => s
10.390 + | _ => error "string_of_term";
10.391 +
10.392 +val string_of_cterm = string_of_term o term_of;
10.393 +
10.394 +fun string_of_varpow x k =
10.395 + if k = 1 then string_of_cterm x
10.396 + else string_of_cterm x^"^"^string_of_int k;
10.397 +
10.398 +fun string_of_monomial m =
10.399 + if Ctermfunc.is_undefined m then "1" else
10.400 + let val vps = fold_rev (fn (x,k) => fn a => string_of_varpow x k :: a)
10.401 + (sort humanorder_varpow (Ctermfunc.graph m)) []
10.402 + in fold1 (fn s => fn t => s^"*"^t) vps
10.403 + end;
10.404 +
10.405 +fun string_of_cmonomial (c,m) =
10.406 + if Ctermfunc.is_undefined m then Rat.string_of_rat c
10.407 + else if c =/ rat_1 then string_of_monomial m
10.408 + else Rat.string_of_rat c ^ "*" ^ string_of_monomial m;;
10.409 +
10.410 +fun string_of_poly (p:poly) =
10.411 + if Monomialfunc.is_undefined p then "<<0>>" else
10.412 + let
10.413 + val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1 m2) (Monomialfunc.graph p)
10.414 + val s = fold (fn (m,c) => fn a =>
10.415 + if c </ rat_0 then a ^ " - " ^ string_of_cmonomial(Rat.neg c,m)
10.416 + else a ^ " + " ^ string_of_cmonomial(c,m))
10.417 + cms ""
10.418 + val s1 = String.substring (s, 0, 3)
10.419 + val s2 = String.substring (s, 3, String.size s - 3)
10.420 + in "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>"
10.421 + end;
10.422 +
10.423 +(* Conversion from HOL term. *)
10.424 +
10.425 +local
10.426 + val neg_tm = @{cterm "uminus :: real => _"}
10.427 + val add_tm = @{cterm "op + :: real => _"}
10.428 + val sub_tm = @{cterm "op - :: real => _"}
10.429 + val mul_tm = @{cterm "op * :: real => _"}
10.430 + val inv_tm = @{cterm "inverse :: real => _"}
10.431 + val div_tm = @{cterm "op / :: real => _"}
10.432 + val pow_tm = @{cterm "op ^ :: real => _"}
10.433 + val zero_tm = @{cterm "0:: real"}
10.434 + val is_numeral = can (HOLogic.dest_number o term_of)
10.435 + fun is_comb t = case t of _$_ => true | _ => false
10.436 + fun poly_of_term tm =
10.437 + if tm aconvc zero_tm then poly_0
10.438 + else if RealArith.is_ratconst tm
10.439 + then poly_const(RealArith.dest_ratconst tm)
10.440 + else
10.441 + (let val (lop,r) = Thm.dest_comb tm
10.442 + in if lop aconvc neg_tm then poly_neg(poly_of_term r)
10.443 + else if lop aconvc inv_tm then
10.444 + let val p = poly_of_term r
10.445 + in if poly_isconst p
10.446 + then poly_const(Rat.inv (eval Ctermfunc.undefined p))
10.447 + else error "poly_of_term: inverse of non-constant polyomial"
10.448 + end
10.449 + else (let val (opr,l) = Thm.dest_comb lop
10.450 + in
10.451 + if opr aconvc pow_tm andalso is_numeral r
10.452 + then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r)
10.453 + else if opr aconvc add_tm
10.454 + then poly_add (poly_of_term l) (poly_of_term r)
10.455 + else if opr aconvc sub_tm
10.456 + then poly_sub (poly_of_term l) (poly_of_term r)
10.457 + else if opr aconvc mul_tm
10.458 + then poly_mul (poly_of_term l) (poly_of_term r)
10.459 + else if opr aconvc div_tm
10.460 + then let
10.461 + val p = poly_of_term l
10.462 + val q = poly_of_term r
10.463 + in if poly_isconst q then poly_cmul (Rat.inv (eval Ctermfunc.undefined q)) p
10.464 + else error "poly_of_term: division by non-constant polynomial"
10.465 + end
10.466 + else poly_var tm
10.467 +
10.468 + end
10.469 + handle CTERM ("dest_comb",_) => poly_var tm)
10.470 + end
10.471 + handle CTERM ("dest_comb",_) => poly_var tm)
10.472 +in
10.473 +val poly_of_term = fn tm =>
10.474 + if type_of (term_of tm) = @{typ real} then poly_of_term tm
10.475 + else error "poly_of_term: term does not have real type"
10.476 +end;
10.477 +
10.478 +(* String of vector (just a list of space-separated numbers). *)
10.479 +
10.480 +fun sdpa_of_vector (v:vector) =
10.481 + let
10.482 + val n = dim v
10.483 + val strs = map (decimalize 20 o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
10.484 + in fold1 (fn x => fn y => x ^ " " ^ y) strs ^ "\n"
10.485 + end;
10.486 +
10.487 +fun increasing f ord (x,y) = ord (f x, f y);
10.488 +fun triple_int_ord ((a,b,c),(a',b',c')) =
10.489 + prod_ord int_ord (prod_ord int_ord int_ord)
10.490 + ((a,(b,c)),(a',(b',c')));
10.491 +structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord);
10.492 +
10.493 +(* String for block diagonal matrix numbered k. *)
10.494 +
10.495 +fun sdpa_of_blockdiagonal k m =
10.496 + let
10.497 + val pfx = string_of_int k ^" "
10.498 + val ents =
10.499 + Inttriplefunc.fold (fn ((b,i,j), c) => fn a => if i > j then a else ((b,i,j),c)::a) m []
10.500 + val entss = sort (increasing fst triple_int_ord ) ents
10.501 + in fold_rev (fn ((b,i,j),c) => fn a =>
10.502 + pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
10.503 + " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
10.504 + end;
10.505 +
10.506 +(* String for a matrix numbered k, in SDPA sparse format. *)
10.507 +
10.508 +fun sdpa_of_matrix k (m:matrix) =
10.509 + let
10.510 + val pfx = string_of_int k ^ " 1 "
10.511 + val ms = Intpairfunc.fold (fn ((i,j), c) => fn a => if i > j then a else ((i,j),c)::a) (snd m) []
10.512 + val mss = sort (increasing fst (prod_ord int_ord int_ord)) ms
10.513 + in fold_rev (fn ((i,j),c) => fn a =>
10.514 + pfx ^ string_of_int i ^ " " ^ string_of_int j ^
10.515 + " " ^ decimalize 20 c ^ "\n" ^ a) mss ""
10.516 + end;;
10.517 +
10.518 +(* ------------------------------------------------------------------------- *)
10.519 +(* String in SDPA sparse format for standard SDP problem: *)
10.520 +(* *)
10.521 +(* X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD *)
10.522 +(* Minimize obj_1 * v_1 + ... obj_m * v_m *)
10.523 +(* ------------------------------------------------------------------------- *)
10.524 +
10.525 +fun sdpa_of_problem obj mats =
10.526 + let
10.527 + val m = length mats - 1
10.528 + val (n,_) = dimensions (hd mats)
10.529 + in
10.530 + string_of_int m ^ "\n" ^
10.531 + "1\n" ^
10.532 + string_of_int n ^ "\n" ^
10.533 + sdpa_of_vector obj ^
10.534 + fold_rev2 (fn k => fn m => fn a => sdpa_of_matrix (k - 1) m ^ a) (1 upto length mats) mats ""
10.535 + end;
10.536 +
10.537 +fun index_char str chr pos =
10.538 + if pos >= String.size str then ~1
10.539 + else if String.sub(str,pos) = chr then pos
10.540 + else index_char str chr (pos + 1);
10.541 +fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b);
10.542 +fun rat_of_string s =
10.543 + let val n = index_char s #"/" 0 in
10.544 + if n = ~1 then s |> IntInf.fromString |> valOf |> Rat.rat_of_int
10.545 + else
10.546 + let val SOME numer = IntInf.fromString(String.substring(s,0,n))
10.547 + val SOME den = IntInf.fromString (String.substring(s,n+1,String.size s - n - 1))
10.548 + in rat_of_quotient(numer, den)
10.549 + end
10.550 + end;
10.551 +
10.552 +fun isspace x = x = " " ;
10.553 +fun isnum x = x mem_string ["0","1","2","3","4","5","6","7","8","9"]
10.554 +
10.555 +(* More parser basics. *)
10.556 +
10.557 +local
10.558 + open Scan
10.559 +in
10.560 + val word = this_string
10.561 + fun token s =
10.562 + repeat ($$ " ") |-- word s --| repeat ($$ " ")
10.563 + val numeral = one isnum
10.564 + val decimalint = bulk numeral >> (rat_of_string o implode)
10.565 + val decimalfrac = bulk numeral
10.566 + >> (fn s => rat_of_string(implode s) // pow10 (length s))
10.567 + val decimalsig =
10.568 + decimalint -- option (Scan.$$ "." |-- decimalfrac)
10.569 + >> (fn (h,NONE) => h | (h,SOME x) => h +/ x)
10.570 + fun signed prs =
10.571 + $$ "-" |-- prs >> Rat.neg
10.572 + || $$ "+" |-- prs
10.573 + || prs;
10.574 +
10.575 +fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs
10.576 +
10.577 + val exponent = ($$ "e" || $$ "E") |-- signed decimalint;
10.578 +
10.579 + val decimal = signed decimalsig -- (emptyin rat_0|| exponent)
10.580 + >> (fn (h, x) => h */ pow10 (int_of_rat x));
10.581 +end;
10.582 +
10.583 + fun mkparser p s =
10.584 + let val (x,rst) = p (explode s)
10.585 + in if null rst then x
10.586 + else error "mkparser: unparsed input"
10.587 + end;;
10.588 +
10.589 +(* Parse back csdp output. *)
10.590 +
10.591 + fun ignore inp = ((),[])
10.592 + fun csdpoutput inp = ((decimal -- Scan.bulk (Scan.$$ " " |-- Scan.option decimal) >> (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp
10.593 + val parse_csdpoutput = mkparser csdpoutput
10.594 +
10.595 +(* Run prover on a problem in linear form. *)
10.596 +
10.597 +fun run_problem prover obj mats =
10.598 + parse_csdpoutput (prover (sdpa_of_problem obj mats))
10.599 +
10.600 +(* Try some apparently sensible scaling first. Note that this is purely to *)
10.601 +(* get a cleaner translation to floating-point, and doesn't affect any of *)
10.602 +(* the results, in principle. In practice it seems a lot better when there *)
10.603 +(* are extreme numbers in the original problem. *)
10.604 +
10.605 + (* Version for (int*int) keys *)
10.606 +local
10.607 + fun max_rat x y = if x </ y then y else x
10.608 + fun common_denominator fld amat acc =
10.609 + fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
10.610 + fun maximal_element fld amat acc =
10.611 + fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
10.612 +fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
10.613 + in Real.fromLargeInt a / Real.fromLargeInt b end;
10.614 +in
10.615 +
10.616 +fun pi_scale_then solver (obj:vector) mats =
10.617 + let
10.618 + val cd1 = fold_rev (common_denominator Intpairfunc.fold) mats (rat_1)
10.619 + val cd2 = common_denominator Intfunc.fold (snd obj) (rat_1)
10.620 + val mats' = map (Intpairfunc.mapf (fn x => cd1 */ x)) mats
10.621 + val obj' = vector_cmul cd2 obj
10.622 + val max1 = fold_rev (maximal_element Intpairfunc.fold) mats' (rat_0)
10.623 + val max2 = maximal_element Intfunc.fold (snd obj') (rat_0)
10.624 + val scal1 = pow2 (20 - trunc(Math.ln (float_of_rat max1) / Math.ln 2.0))
10.625 + val scal2 = pow2 (20 - trunc(Math.ln (float_of_rat max2) / Math.ln 2.0))
10.626 + val mats'' = map (Intpairfunc.mapf (fn x => x */ scal1)) mats'
10.627 + val obj'' = vector_cmul scal2 obj'
10.628 + in solver obj'' mats''
10.629 + end
10.630 +end;
10.631 +
10.632 +(* Try some apparently sensible scaling first. Note that this is purely to *)
10.633 +(* get a cleaner translation to floating-point, and doesn't affect any of *)
10.634 +(* the results, in principle. In practice it seems a lot better when there *)
10.635 +(* are extreme numbers in the original problem. *)
10.636 +
10.637 + (* Version for (int*int*int) keys *)
10.638 +local
10.639 + fun max_rat x y = if x </ y then y else x
10.640 + fun common_denominator fld amat acc =
10.641 + fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
10.642 + fun maximal_element fld amat acc =
10.643 + fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
10.644 +fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
10.645 + in Real.fromLargeInt a / Real.fromLargeInt b end;
10.646 +fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0)
10.647 +in
10.648 +
10.649 +fun tri_scale_then solver (obj:vector) mats =
10.650 + let
10.651 + val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1)
10.652 + val cd2 = common_denominator Intfunc.fold (snd obj) (rat_1)
10.653 + val mats' = map (Inttriplefunc.mapf (fn x => cd1 */ x)) mats
10.654 + val obj' = vector_cmul cd2 obj
10.655 + val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0)
10.656 + val max2 = maximal_element Intfunc.fold (snd obj') (rat_0)
10.657 + val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0))
10.658 + val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0))
10.659 + val mats'' = map (Inttriplefunc.mapf (fn x => x */ scal1)) mats'
10.660 + val obj'' = vector_cmul scal2 obj'
10.661 + in solver obj'' mats''
10.662 + end
10.663 +end;
10.664 +
10.665 +(* Round a vector to "nice" rationals. *)
10.666 +
10.667 +fun nice_rational n x = round_rat (n */ x) // n;;
10.668 +fun nice_vector n ((d,v) : vector) =
10.669 + (d, Intfunc.fold (fn (i,c) => fn a =>
10.670 + let val y = nice_rational n c
10.671 + in if c =/ rat_0 then a
10.672 + else Intfunc.update (i,y) a end) v Intfunc.undefined):vector
10.673 +
10.674 +fun dest_ord f x = is_equal (f x);
10.675 +
10.676 +(* Stuff for "equations" ((int*int*int)->num functions). *)
10.677 +
10.678 +fun tri_equation_cmul c eq =
10.679 + if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
10.680 +
10.681 +fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
10.682 +
10.683 +fun tri_equation_eval assig eq =
10.684 + let fun value v = Inttriplefunc.apply assig v
10.685 + in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
10.686 + end;
10.687 +
10.688 +(* Eliminate among linear equations: return unconstrained variables and *)
10.689 +(* assignments for the others in terms of them. We give one pseudo-variable *)
10.690 +(* "one" that's used for a constant term. *)
10.691 +
10.692 +local
10.693 + fun extract_first p l = case l of (* FIXME : use find_first instead *)
10.694 + [] => error "extract_first"
10.695 + | h::t => if p h then (h,t) else
10.696 + let val (k,s) = extract_first p t in (k,h::s) end
10.697 +fun eliminate vars dun eqs = case vars of
10.698 + [] => if forall Inttriplefunc.is_undefined eqs then dun
10.699 + else raise Unsolvable
10.700 + | v::vs =>
10.701 + ((let
10.702 + val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs
10.703 + val a = Inttriplefunc.apply eq v
10.704 + val eq' = tri_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
10.705 + fun elim e =
10.706 + let val b = Inttriplefunc.tryapplyd e v rat_0
10.707 + in if b =/ rat_0 then e else
10.708 + tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
10.709 + end
10.710 + in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
10.711 + end)
10.712 + handle Failure _ => eliminate vs dun eqs)
10.713 +in
10.714 +fun tri_eliminate_equations one vars eqs =
10.715 + let
10.716 + val assig = eliminate vars Inttriplefunc.undefined eqs
10.717 + val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
10.718 + in (distinct (dest_ord triple_int_ord) vs, assig)
10.719 + end
10.720 +end;
10.721 +
10.722 +(* Eliminate all variables, in an essentially arbitrary order. *)
10.723 +
10.724 +fun tri_eliminate_all_equations one =
10.725 + let
10.726 + fun choose_variable eq =
10.727 + let val (v,_) = Inttriplefunc.choose eq
10.728 + in if is_equal (triple_int_ord(v,one)) then
10.729 + let val eq' = Inttriplefunc.undefine v eq
10.730 + in if Inttriplefunc.is_undefined eq' then error "choose_variable"
10.731 + else fst (Inttriplefunc.choose eq')
10.732 + end
10.733 + else v
10.734 + end
10.735 + fun eliminate dun eqs = case eqs of
10.736 + [] => dun
10.737 + | eq::oeqs =>
10.738 + if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
10.739 + let val v = choose_variable eq
10.740 + val a = Inttriplefunc.apply eq v
10.741 + val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a)
10.742 + (Inttriplefunc.undefine v eq)
10.743 + fun elim e =
10.744 + let val b = Inttriplefunc.tryapplyd e v rat_0
10.745 + in if b =/ rat_0 then e
10.746 + else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
10.747 + end
10.748 + in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun))
10.749 + (map elim oeqs)
10.750 + end
10.751 +in fn eqs =>
10.752 + let
10.753 + val assig = eliminate Inttriplefunc.undefined eqs
10.754 + val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
10.755 + in (distinct (dest_ord triple_int_ord) vs,assig)
10.756 + end
10.757 +end;
10.758 +
10.759 +(* Solve equations by assigning arbitrary numbers. *)
10.760 +
10.761 +fun tri_solve_equations one eqs =
10.762 + let
10.763 + val (vars,assigs) = tri_eliminate_all_equations one eqs
10.764 + val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars
10.765 + (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
10.766 + val ass =
10.767 + Inttriplefunc.combine (curry op +/) (K false)
10.768 + (Inttriplefunc.mapf (tri_equation_eval vfn) assigs) vfn
10.769 + in if forall (fn e => tri_equation_eval ass e =/ rat_0) eqs
10.770 + then Inttriplefunc.undefine one ass else raise Sanity
10.771 + end;
10.772 +
10.773 +(* Multiply equation-parametrized poly by regular poly and add accumulator. *)
10.774 +
10.775 +fun tri_epoly_pmul p q acc =
10.776 + Monomialfunc.fold (fn (m1, c) => fn a =>
10.777 + Monomialfunc.fold (fn (m2,e) => fn b =>
10.778 + let val m = monomial_mul m1 m2
10.779 + val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined
10.780 + in Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b
10.781 + end) q a) p acc ;
10.782 +
10.783 +(* Usual operations on equation-parametrized poly. *)
10.784 +
10.785 +fun tri_epoly_cmul c l =
10.786 + if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (tri_equation_cmul c) l;;
10.787 +
10.788 +val tri_epoly_neg = tri_epoly_cmul (Rat.rat_of_int ~1);
10.789 +
10.790 +val tri_epoly_add = Inttriplefunc.combine tri_equation_add Inttriplefunc.is_undefined;
10.791 +
10.792 +fun tri_epoly_sub p q = tri_epoly_add p (tri_epoly_neg q);;
10.793 +
10.794 +(* Stuff for "equations" ((int*int)->num functions). *)
10.795 +
10.796 +fun pi_equation_cmul c eq =
10.797 + if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
10.798 +
10.799 +fun pi_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
10.800 +
10.801 +fun pi_equation_eval assig eq =
10.802 + let fun value v = Inttriplefunc.apply assig v
10.803 + in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
10.804 + end;
10.805 +
10.806 +(* Eliminate among linear equations: return unconstrained variables and *)
10.807 +(* assignments for the others in terms of them. We give one pseudo-variable *)
10.808 +(* "one" that's used for a constant term. *)
10.809 +
10.810 +local
10.811 +fun extract_first p l = case l of
10.812 + [] => error "extract_first"
10.813 + | h::t => if p h then (h,t) else
10.814 + let val (k,s) = extract_first p t in (k,h::s) end
10.815 +fun eliminate vars dun eqs = case vars of
10.816 + [] => if forall Inttriplefunc.is_undefined eqs then dun
10.817 + else raise Unsolvable
10.818 + | v::vs =>
10.819 + let
10.820 + val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs
10.821 + val a = Inttriplefunc.apply eq v
10.822 + val eq' = pi_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
10.823 + fun elim e =
10.824 + let val b = Inttriplefunc.tryapplyd e v rat_0
10.825 + in if b =/ rat_0 then e else
10.826 + pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
10.827 + end
10.828 + in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
10.829 + end
10.830 + handle Failure _ => eliminate vs dun eqs
10.831 +in
10.832 +fun pi_eliminate_equations one vars eqs =
10.833 + let
10.834 + val assig = eliminate vars Inttriplefunc.undefined eqs
10.835 + val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
10.836 + in (distinct (dest_ord triple_int_ord) vs, assig)
10.837 + end
10.838 +end;
10.839 +
10.840 +(* Eliminate all variables, in an essentially arbitrary order. *)
10.841 +
10.842 +fun pi_eliminate_all_equations one =
10.843 + let
10.844 + fun choose_variable eq =
10.845 + let val (v,_) = Inttriplefunc.choose eq
10.846 + in if is_equal (triple_int_ord(v,one)) then
10.847 + let val eq' = Inttriplefunc.undefine v eq
10.848 + in if Inttriplefunc.is_undefined eq' then error "choose_variable"
10.849 + else fst (Inttriplefunc.choose eq')
10.850 + end
10.851 + else v
10.852 + end
10.853 + fun eliminate dun eqs = case eqs of
10.854 + [] => dun
10.855 + | eq::oeqs =>
10.856 + if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
10.857 + let val v = choose_variable eq
10.858 + val a = Inttriplefunc.apply eq v
10.859 + val eq' = pi_equation_cmul ((Rat.rat_of_int ~1) // a)
10.860 + (Inttriplefunc.undefine v eq)
10.861 + fun elim e =
10.862 + let val b = Inttriplefunc.tryapplyd e v rat_0
10.863 + in if b =/ rat_0 then e
10.864 + else pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
10.865 + end
10.866 + in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun))
10.867 + (map elim oeqs)
10.868 + end
10.869 +in fn eqs =>
10.870 + let
10.871 + val assig = eliminate Inttriplefunc.undefined eqs
10.872 + val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
10.873 + in (distinct (dest_ord triple_int_ord) vs,assig)
10.874 + end
10.875 +end;
10.876 +
10.877 +(* Solve equations by assigning arbitrary numbers. *)
10.878 +
10.879 +fun pi_solve_equations one eqs =
10.880 + let
10.881 + val (vars,assigs) = pi_eliminate_all_equations one eqs
10.882 + val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars
10.883 + (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
10.884 + val ass =
10.885 + Inttriplefunc.combine (curry op +/) (K false)
10.886 + (Inttriplefunc.mapf (pi_equation_eval vfn) assigs) vfn
10.887 + in if forall (fn e => pi_equation_eval ass e =/ rat_0) eqs
10.888 + then Inttriplefunc.undefine one ass else raise Sanity
10.889 + end;
10.890 +
10.891 +(* Multiply equation-parametrized poly by regular poly and add accumulator. *)
10.892 +
10.893 +fun pi_epoly_pmul p q acc =
10.894 + Monomialfunc.fold (fn (m1, c) => fn a =>
10.895 + Monomialfunc.fold (fn (m2,e) => fn b =>
10.896 + let val m = monomial_mul m1 m2
10.897 + val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined
10.898 + in Monomialfunc.update (m,pi_equation_add (pi_equation_cmul c e) es) b
10.899 + end) q a) p acc ;
10.900 +
10.901 +(* Usual operations on equation-parametrized poly. *)
10.902 +
10.903 +fun pi_epoly_cmul c l =
10.904 + if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (pi_equation_cmul c) l;;
10.905 +
10.906 +val pi_epoly_neg = pi_epoly_cmul (Rat.rat_of_int ~1);
10.907 +
10.908 +val pi_epoly_add = Inttriplefunc.combine pi_equation_add Inttriplefunc.is_undefined;
10.909 +
10.910 +fun pi_epoly_sub p q = pi_epoly_add p (pi_epoly_neg q);;
10.911 +
10.912 +fun allpairs f l1 l2 = fold_rev (fn x => (curry (op @)) (map (f x) l2)) l1 [];
10.913 +
10.914 +(* Hence produce the "relevant" monomials: those whose squares lie in the *)
10.915 +(* Newton polytope of the monomials in the input. (This is enough according *)
10.916 +(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal, *)
10.917 +(* vol 45, pp. 363--374, 1978. *)
10.918 +(* *)
10.919 +(* These are ordered in sort of decreasing degree. In particular the *)
10.920 +(* constant monomial is last; this gives an order in diagonalization of the *)
10.921 +(* quadratic form that will tend to display constants. *)
10.922 +
10.923 +(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form. *)
10.924 +
10.925 +local
10.926 +fun diagonalize n i m =
10.927 + if Intpairfunc.is_undefined (snd m) then []
10.928 + else
10.929 + let val a11 = Intpairfunc.tryapplyd (snd m) (i,i) rat_0
10.930 + in if a11 </ rat_0 then raise Failure "diagonalize: not PSD"
10.931 + else if a11 =/ rat_0 then
10.932 + if Intfunc.is_undefined (snd (row i m)) then diagonalize n (i + 1) m
10.933 + else raise Failure "diagonalize: not PSD ___ "
10.934 + else
10.935 + let
10.936 + val v = row i m
10.937 + val v' = (fst v, Intfunc.fold (fn (i, c) => fn a =>
10.938 + let val y = c // a11
10.939 + in if y = rat_0 then a else Intfunc.update (i,y) a
10.940 + end) (snd v) Intfunc.undefined)
10.941 + fun upt0 x y a = if y = rat_0 then a else Intpairfunc.update (x,y) a
10.942 + val m' =
10.943 + ((n,n),
10.944 + iter (i+1,n) (fn j =>
10.945 + iter (i+1,n) (fn k =>
10.946 + (upt0 (j,k) (Intpairfunc.tryapplyd (snd m) (j,k) rat_0 -/ Intfunc.tryapplyd (snd v) j rat_0 */ Intfunc.tryapplyd (snd v') k rat_0))))
10.947 + Intpairfunc.undefined)
10.948 + in (a11,v')::diagonalize n (i + 1) m'
10.949 + end
10.950 + end
10.951 +in
10.952 +fun diag m =
10.953 + let
10.954 + val nn = dimensions m
10.955 + val n = fst nn
10.956 + in if snd nn <> n then error "diagonalize: non-square matrix"
10.957 + else diagonalize n 1 m
10.958 + end
10.959 +end;
10.960 +
10.961 +fun gcd_rat a b = Rat.rat_of_int (Integer.gcd (int_of_rat a) (int_of_rat b));
10.962 +
10.963 +(* Adjust a diagonalization to collect rationals at the start. *)
10.964 + (* FIXME : Potentially polymorphic keys, but here only: integers!! *)
10.965 +local
10.966 + fun upd0 x y a = if y =/ rat_0 then a else Intfunc.update(x,y) a;
10.967 + fun mapa f (d,v) =
10.968 + (d, Intfunc.fold (fn (i,c) => fn a => upd0 i (f c) a) v Intfunc.undefined)
10.969 + fun adj (c,l) =
10.970 + let val a =
10.971 + Intfunc.fold (fn (i,c) => fn a => lcm_rat a (denominator_rat c))
10.972 + (snd l) rat_1 //
10.973 + Intfunc.fold (fn (i,c) => fn a => gcd_rat a (numerator_rat c))
10.974 + (snd l) rat_0
10.975 + in ((c // (a */ a)),mapa (fn x => a */ x) l)
10.976 + end
10.977 +in
10.978 +fun deration d = if null d then (rat_0,d) else
10.979 + let val d' = map adj d
10.980 + val a = fold (lcm_rat o denominator_rat o fst) d' rat_1 //
10.981 + fold (gcd_rat o numerator_rat o fst) d' rat_0
10.982 + in ((rat_1 // a),map (fn (c,l) => (a */ c,l)) d')
10.983 + end
10.984 +end;
10.985 +
10.986 +(* Enumeration of monomials with given multidegree bound. *)
10.987 +
10.988 +fun enumerate_monomials d vars =
10.989 + if d < 0 then []
10.990 + else if d = 0 then [Ctermfunc.undefined]
10.991 + else if null vars then [monomial_1] else
10.992 + let val alts =
10.993 + map (fn k => let val oths = enumerate_monomials (d - k) (tl vars)
10.994 + in map (fn ks => if k = 0 then ks else Ctermfunc.update (hd vars, k) ks) oths end) (0 upto d)
10.995 + in fold1 (curry op @) alts
10.996 + end;
10.997 +
10.998 +(* Enumerate products of distinct input polys with degree <= d. *)
10.999 +(* We ignore any constant input polynomials. *)
10.1000 +(* Give the output polynomial and a record of how it was derived. *)
10.1001 +
10.1002 +local
10.1003 + open RealArith
10.1004 +in
10.1005 +fun enumerate_products d pols =
10.1006 +if d = 0 then [(poly_const rat_1,Rational_lt rat_1)]
10.1007 +else if d < 0 then [] else
10.1008 +case pols of
10.1009 + [] => [(poly_const rat_1,Rational_lt rat_1)]
10.1010 + | (p,b)::ps =>
10.1011 + let val e = multidegree p
10.1012 + in if e = 0 then enumerate_products d ps else
10.1013 + enumerate_products d ps @
10.1014 + map (fn (q,c) => (poly_mul p q,Product(b,c)))
10.1015 + (enumerate_products (d - e) ps)
10.1016 + end
10.1017 +end;
10.1018 +
10.1019 +(* Convert regular polynomial. Note that we treat (0,0,0) as -1. *)
10.1020 +
10.1021 +fun epoly_of_poly p =
10.1022 + Monomialfunc.fold (fn (m,c) => fn a => Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p Monomialfunc.undefined;
10.1023 +
10.1024 +(* String for block diagonal matrix numbered k. *)
10.1025 +
10.1026 +fun sdpa_of_blockdiagonal k m =
10.1027 + let
10.1028 + val pfx = string_of_int k ^" "
10.1029 + val ents =
10.1030 + Inttriplefunc.fold
10.1031 + (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a)
10.1032 + m []
10.1033 + val entss = sort (increasing fst triple_int_ord) ents
10.1034 + in fold_rev (fn ((b,i,j),c) => fn a =>
10.1035 + pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
10.1036 + " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
10.1037 + end;
10.1038 +
10.1039 +(* SDPA for problem using block diagonal (i.e. multiple SDPs) *)
10.1040 +
10.1041 +fun sdpa_of_blockproblem nblocks blocksizes obj mats =
10.1042 + let val m = length mats - 1
10.1043 + in
10.1044 + string_of_int m ^ "\n" ^
10.1045 + string_of_int nblocks ^ "\n" ^
10.1046 + (fold1 (fn s => fn t => s^" "^t) (map string_of_int blocksizes)) ^
10.1047 + "\n" ^
10.1048 + sdpa_of_vector obj ^
10.1049 + fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a)
10.1050 + (1 upto length mats) mats ""
10.1051 + end;
10.1052 +
10.1053 +(* Run prover on a problem in block diagonal form. *)
10.1054 +
10.1055 +fun run_blockproblem prover nblocks blocksizes obj mats=
10.1056 + parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats))
10.1057 +
10.1058 +(* 3D versions of matrix operations to consider blocks separately. *)
10.1059 +
10.1060 +val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0);
10.1061 +fun bmatrix_cmul c bm =
10.1062 + if c =/ rat_0 then Inttriplefunc.undefined
10.1063 + else Inttriplefunc.mapf (fn x => c */ x) bm;
10.1064 +
10.1065 +val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1);
10.1066 +fun bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);;
10.1067 +
10.1068 +(* Smash a block matrix into components. *)
10.1069 +
10.1070 +fun blocks blocksizes bm =
10.1071 + map (fn (bs,b0) =>
10.1072 + let val m = Inttriplefunc.fold
10.1073 + (fn ((b,i,j),c) => fn a => if b = b0 then Intpairfunc.update ((i,j),c) a else a) bm Intpairfunc.undefined
10.1074 + val d = Intpairfunc.fold (fn ((i,j),c) => fn a => max a (max i j)) m 0
10.1075 + in (((bs,bs),m):matrix) end)
10.1076 + (blocksizes ~~ (1 upto length blocksizes));;
10.1077 +
10.1078 +(* FIXME : Get rid of this !!!*)
10.1079 +local
10.1080 + fun tryfind_with msg f [] = raise Failure msg
10.1081 + | tryfind_with msg f (x::xs) = (f x handle Failure s => tryfind_with s f xs);
10.1082 +in
10.1083 + fun tryfind f = tryfind_with "tryfind" f
10.1084 +end
10.1085 +
10.1086 +(*
10.1087 +fun tryfind f [] = error "tryfind"
10.1088 + | tryfind f (x::xs) = (f x handle ERROR _ => tryfind f xs);
10.1089 +*)
10.1090 +
10.1091 +(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
10.1092 +
10.1093 +
10.1094 +local
10.1095 + open RealArith
10.1096 +in
10.1097 +fun real_positivnullstellensatz_general prover linf d eqs leqs pol =
10.1098 +let
10.1099 + val vars = fold_rev (curry (gen_union (op aconvc)) o poly_variables)
10.1100 + (pol::eqs @ map fst leqs) []
10.1101 + val monoid = if linf then
10.1102 + (poly_const rat_1,Rational_lt rat_1)::
10.1103 + (filter (fn (p,c) => multidegree p <= d) leqs)
10.1104 + else enumerate_products d leqs
10.1105 + val nblocks = length monoid
10.1106 + fun mk_idmultiplier k p =
10.1107 + let
10.1108 + val e = d - multidegree p
10.1109 + val mons = enumerate_monomials e vars
10.1110 + val nons = mons ~~ (1 upto length mons)
10.1111 + in (mons,
10.1112 + fold_rev (fn (m,n) => Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons Monomialfunc.undefined)
10.1113 + end
10.1114 +
10.1115 + fun mk_sqmultiplier k (p,c) =
10.1116 + let
10.1117 + val e = (d - multidegree p) div 2
10.1118 + val mons = enumerate_monomials e vars
10.1119 + val nons = mons ~~ (1 upto length mons)
10.1120 + in (mons,
10.1121 + fold_rev (fn (m1,n1) =>
10.1122 + fold_rev (fn (m2,n2) => fn a =>
10.1123 + let val m = monomial_mul m1 m2
10.1124 + in if n1 > n2 then a else
10.1125 + let val c = if n1 = n2 then rat_1 else rat_2
10.1126 + val e = Monomialfunc.tryapplyd a m Inttriplefunc.undefined
10.1127 + in Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a
10.1128 + end
10.1129 + end) nons)
10.1130 + nons Monomialfunc.undefined)
10.1131 + end
10.1132 +
10.1133 + val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid)
10.1134 + val (idmonlist,ids) = split_list(map2 mk_idmultiplier (1 upto length eqs) eqs)
10.1135 + val blocksizes = map length sqmonlist
10.1136 + val bigsum =
10.1137 + fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids
10.1138 + (fold_rev2 (fn (p,c) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs
10.1139 + (epoly_of_poly(poly_neg pol)))
10.1140 + val eqns = Monomialfunc.fold (fn (m,e) => fn a => e::a) bigsum []
10.1141 + val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns
10.1142 + val qvars = (0,0,0)::pvs
10.1143 + val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig
10.1144 + fun mk_matrix v =
10.1145 + Inttriplefunc.fold (fn ((b,i,j), ass) => fn m =>
10.1146 + if b < 0 then m else
10.1147 + let val c = Inttriplefunc.tryapplyd ass v rat_0
10.1148 + in if c = rat_0 then m else
10.1149 + Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m)
10.1150 + end)
10.1151 + allassig Inttriplefunc.undefined
10.1152 + val diagents = Inttriplefunc.fold
10.1153 + (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a)
10.1154 + allassig Inttriplefunc.undefined
10.1155 +
10.1156 + val mats = map mk_matrix qvars
10.1157 + val obj = (length pvs,
10.1158 + itern 1 pvs (fn v => fn i => Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0))
10.1159 + Intfunc.undefined)
10.1160 + val raw_vec = if null pvs then vector_0 0
10.1161 + else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats
10.1162 + fun int_element (d,v) i = Intfunc.tryapplyd v i rat_0
10.1163 + fun cterm_element (d,v) i = Ctermfunc.tryapplyd v i rat_0
10.1164 +
10.1165 + fun find_rounding d =
10.1166 + let
10.1167 + val _ = if !debugging
10.1168 + then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n")
10.1169 + else ()
10.1170 + val vec = nice_vector d raw_vec
10.1171 + val blockmat = iter (1,dim vec)
10.1172 + (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a)
10.1173 + (bmatrix_neg (nth mats 0))
10.1174 + val allmats = blocks blocksizes blockmat
10.1175 + in (vec,map diag allmats)
10.1176 + end
10.1177 + val (vec,ratdias) =
10.1178 + if null pvs then find_rounding rat_1
10.1179 + else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @
10.1180 + map pow2 (5 upto 66))
10.1181 + val newassigs =
10.1182 + fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k))
10.1183 + (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1))
10.1184 + val finalassigs =
10.1185 + Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs
10.1186 + fun poly_of_epoly p =
10.1187 + Monomialfunc.fold (fn (v,e) => fn a => Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a)
10.1188 + p Monomialfunc.undefined
10.1189 + fun mk_sos mons =
10.1190 + let fun mk_sq (c,m) =
10.1191 + (c,fold_rev (fn k=> fn a => Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a)
10.1192 + (1 upto length mons) Monomialfunc.undefined)
10.1193 + in map mk_sq
10.1194 + end
10.1195 + val sqs = map2 mk_sos sqmonlist ratdias
10.1196 + val cfs = map poly_of_epoly ids
10.1197 + val msq = filter (fn (a,b) => not (null b)) (map2 pair monoid sqs)
10.1198 + fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0
10.1199 + val sanity =
10.1200 + fold_rev (fn ((p,c),s) => poly_add (poly_mul p (eval_sq s))) msq
10.1201 + (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs
10.1202 + (poly_neg pol))
10.1203 +
10.1204 +in if not(Monomialfunc.is_undefined sanity) then raise Sanity else
10.1205 + (cfs,map (fn (a,b) => (snd a,b)) msq)
10.1206 + end
10.1207 +
10.1208 +
10.1209 +end;
10.1210 +
10.1211 +(* Iterative deepening. *)
10.1212 +
10.1213 +fun deepen f n =
10.1214 + (writeln ("Searching with depth limit " ^ string_of_int n) ; (f n handle Failure s => (writeln ("failed with message: " ^ s) ; deepen f (n+1))))
10.1215 +
10.1216 +(* The ordering so we can create canonical HOL polynomials. *)
10.1217 +
10.1218 +fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon);
10.1219 +
10.1220 +fun monomial_order (m1,m2) =
10.1221 + if Ctermfunc.is_undefined m2 then LESS
10.1222 + else if Ctermfunc.is_undefined m1 then GREATER
10.1223 + else
10.1224 + let val mon1 = dest_monomial m1
10.1225 + val mon2 = dest_monomial m2
10.1226 + val deg1 = fold (curry op + o snd) mon1 0
10.1227 + val deg2 = fold (curry op + o snd) mon2 0
10.1228 + in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS
10.1229 + else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2)
10.1230 + end;
10.1231 +
10.1232 +fun dest_poly p =
10.1233 + map (fn (m,c) => (c,dest_monomial m))
10.1234 + (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p));
10.1235 +
10.1236 +(* Map back polynomials and their composites to HOL. *)
10.1237 +
10.1238 +local
10.1239 + open Thm Numeral RealArith
10.1240 +in
10.1241 +
10.1242 +fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x)
10.1243 + (mk_cnumber @{ctyp nat} k)
10.1244 +
10.1245 +fun cterm_of_monomial m =
10.1246 + if Ctermfunc.is_undefined m then @{cterm "1::real"}
10.1247 + else
10.1248 + let
10.1249 + val m' = dest_monomial m
10.1250 + val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' []
10.1251 + in fold1 (fn s => fn t => capply (capply @{cterm "op * :: real => _"} s) t) vps
10.1252 + end
10.1253 +
10.1254 +fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c
10.1255 + else if c = Rat.one then cterm_of_monomial m
10.1256 + else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
10.1257 +
10.1258 +fun cterm_of_poly p =
10.1259 + if Monomialfunc.is_undefined p then @{cterm "0::real"}
10.1260 + else
10.1261 + let
10.1262 + val cms = map cterm_of_cmonomial
10.1263 + (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
10.1264 + in fold1 (fn t1 => fn t2 => capply(capply @{cterm "op + :: real => _"} t1) t2) cms
10.1265 + end;
10.1266 +
10.1267 +fun cterm_of_sqterm (c,p) = Product(Rational_lt c,Square(cterm_of_poly p));
10.1268 +
10.1269 +fun cterm_of_sos (pr,sqs) = if null sqs then pr
10.1270 + else Product(pr,fold1 (fn a => fn b => Sum(a,b)) (map cterm_of_sqterm sqs));
10.1271 +
10.1272 +end
10.1273 +
10.1274 +(* Interface to HOL. *)
10.1275 +local
10.1276 + open Thm Conv RealArith
10.1277 + val concl = dest_arg o cprop_of
10.1278 + fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
10.1279 +in
10.1280 + (* FIXME: Replace tryfind by get_first !! *)
10.1281 +fun real_nonlinear_prover prover ctxt =
10.1282 + let
10.1283 + val {add,mul,neg,pow,sub,main} = Normalizer.semiring_normalizers_ord_wrapper ctxt
10.1284 + (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"}))
10.1285 + simple_cterm_ord
10.1286 + val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
10.1287 + real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
10.1288 + fun mainf translator (eqs,les,lts) =
10.1289 + let
10.1290 + val eq0 = map (poly_of_term o dest_arg1 o concl) eqs
10.1291 + val le0 = map (poly_of_term o dest_arg o concl) les
10.1292 + val lt0 = map (poly_of_term o dest_arg o concl) lts
10.1293 + val eqp0 = map (fn (t,i) => (t,Axiom_eq i)) (eq0 ~~ (0 upto (length eq0 - 1)))
10.1294 + val lep0 = map (fn (t,i) => (t,Axiom_le i)) (le0 ~~ (0 upto (length le0 - 1)))
10.1295 + val ltp0 = map (fn (t,i) => (t,Axiom_lt i)) (lt0 ~~ (0 upto (length lt0 - 1)))
10.1296 + val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0
10.1297 + val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0
10.1298 + val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0
10.1299 + fun trivial_axiom (p,ax) =
10.1300 + case ax of
10.1301 + Axiom_eq n => if eval Ctermfunc.undefined p <>/ Rat.zero then nth eqs n
10.1302 + else raise Failure "trivial_axiom: Not a trivial axiom"
10.1303 + | Axiom_le n => if eval Ctermfunc.undefined p </ Rat.zero then nth les n
10.1304 + else raise Failure "trivial_axiom: Not a trivial axiom"
10.1305 + | Axiom_lt n => if eval Ctermfunc.undefined p <=/ Rat.zero then nth lts n
10.1306 + else raise Failure "trivial_axiom: Not a trivial axiom"
10.1307 + | _ => error "trivial_axiom: Not a trivial axiom"
10.1308 + in
10.1309 + ((let val th = tryfind trivial_axiom (keq @ klep @ kltp)
10.1310 + in fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv field_comp_conv) th end)
10.1311 + handle Failure _ => (
10.1312 + let
10.1313 + val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
10.1314 + val leq = lep @ ltp
10.1315 + fun tryall d =
10.1316 + let val e = multidegree pol
10.1317 + val k = if e = 0 then 0 else d div e
10.1318 + val eq' = map fst eq
10.1319 + in tryfind (fn i => (d,i,real_positivnullstellensatz_general prover false d eq' leq
10.1320 + (poly_neg(poly_pow pol i))))
10.1321 + (0 upto k)
10.1322 + end
10.1323 + val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
10.1324 + val proofs_ideal =
10.1325 + map2 (fn q => fn (p,ax) => Eqmul(cterm_of_poly q,ax)) cert_ideal eq
10.1326 + val proofs_cone = map cterm_of_sos cert_cone
10.1327 + val proof_ne = if null ltp then Rational_lt Rat.one else
10.1328 + let val p = fold1 (fn s => fn t => Product(s,t)) (map snd ltp)
10.1329 + in funpow i (fn q => Product(p,q)) (Rational_lt Rat.one)
10.1330 + end
10.1331 + val proof = fold1 (fn s => fn t => Sum(s,t))
10.1332 + (proof_ne :: proofs_ideal @ proofs_cone)
10.1333 + in writeln "Translating proof certificate to HOL";
10.1334 + translator (eqs,les,lts) proof
10.1335 + end))
10.1336 + end
10.1337 + in mainf end
10.1338 +end
10.1339 +
10.1340 +fun C f x y = f y x;
10.1341 + (* FIXME : This is very bad!!!*)
10.1342 +fun subst_conv eqs t =
10.1343 + let
10.1344 + val t' = fold (Thm.cabs o Thm.lhs_of) eqs t
10.1345 + in Conv.fconv_rule (Thm.beta_conversion true) (fold (C combination) eqs (reflexive t'))
10.1346 + end
10.1347 +
10.1348 +(* A wrapper that tries to substitute away variables first. *)
10.1349 +
10.1350 +local
10.1351 + open Thm Conv RealArith
10.1352 + fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
10.1353 + val concl = dest_arg o cprop_of
10.1354 + val shuffle1 =
10.1355 + fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps) })
10.1356 + val shuffle2 =
10.1357 + fconv_rule (rewr_conv @{lemma "(x + a == y) == (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps)})
10.1358 + fun substitutable_monomial fvs tm = case term_of tm of
10.1359 + Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm)
10.1360 + else raise Failure "substitutable_monomial"
10.1361 + | @{term "op * :: real => _"}$c$(t as Free _ ) =>
10.1362 + if is_ratconst (dest_arg1 tm) andalso not (member (op aconvc) fvs (dest_arg tm))
10.1363 + then (dest_ratconst (dest_arg1 tm),dest_arg tm) else raise Failure "substitutable_monomial"
10.1364 + | @{term "op + :: real => _"}$s$t =>
10.1365 + (substitutable_monomial (add_cterm_frees (dest_arg tm) fvs) (dest_arg1 tm)
10.1366 + handle Failure _ => substitutable_monomial (add_cterm_frees (dest_arg1 tm) fvs) (dest_arg tm))
10.1367 + | _ => raise Failure "substitutable_monomial"
10.1368 +
10.1369 + fun isolate_variable v th =
10.1370 + let val w = dest_arg1 (cprop_of th)
10.1371 + in if v aconvc w then th
10.1372 + else case term_of w of
10.1373 + @{term "op + :: real => _"}$s$t =>
10.1374 + if dest_arg1 w aconvc v then shuffle2 th
10.1375 + else isolate_variable v (shuffle1 th)
10.1376 + | _ => error "isolate variable : This should not happen?"
10.1377 + end
10.1378 +in
10.1379 +
10.1380 +fun real_nonlinear_subst_prover prover ctxt =
10.1381 + let
10.1382 + val {add,mul,neg,pow,sub,main} = Normalizer.semiring_normalizers_ord_wrapper ctxt
10.1383 + (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"}))
10.1384 + simple_cterm_ord
10.1385 +
10.1386 + val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
10.1387 + real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
10.1388 +
10.1389 + fun make_substitution th =
10.1390 + let
10.1391 + val (c,v) = substitutable_monomial [] (dest_arg1(concl th))
10.1392 + val th1 = Drule.arg_cong_rule (capply @{cterm "op * :: real => _"} (cterm_of_rat (Rat.inv c))) (mk_meta_eq th)
10.1393 + val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1
10.1394 + in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2)
10.1395 + end
10.1396 + fun oprconv cv ct =
10.1397 + let val g = Thm.dest_fun2 ct
10.1398 + in if g aconvc @{cterm "op <= :: real => _"}
10.1399 + orelse g aconvc @{cterm "op < :: real => _"}
10.1400 + then arg_conv cv ct else arg1_conv cv ct
10.1401 + end
10.1402 + fun mainf translator =
10.1403 + let
10.1404 + fun substfirst(eqs,les,lts) =
10.1405 + ((let
10.1406 + val eth = tryfind make_substitution eqs
10.1407 + val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv)))
10.1408 + in substfirst
10.1409 + (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t
10.1410 + aconvc @{cterm "0::real"}) (map modify eqs),
10.1411 + map modify les,map modify lts)
10.1412 + end)
10.1413 + handle Failure _ => real_nonlinear_prover prover ctxt translator (rev eqs, rev les, rev lts))
10.1414 + in substfirst
10.1415 + end
10.1416 +
10.1417 +
10.1418 + in mainf
10.1419 + end
10.1420 +
10.1421 +(* Overall function. *)
10.1422 +
10.1423 +fun real_sos prover ctxt t = gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) t;
10.1424 +end;
10.1425 +
10.1426 +(* A tactic *)
10.1427 +fun strip_all ct =
10.1428 + case term_of ct of
10.1429 + Const("all",_) $ Abs (xn,xT,p) =>
10.1430 + let val (a,(v,t')) = (apsnd (Thm.dest_abs (SOME xn)) o Thm.dest_comb) ct
10.1431 + in apfst (cons v) (strip_all t')
10.1432 + end
10.1433 +| _ => ([],ct)
10.1434 +
10.1435 +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})
10.1436 +
10.1437 +val known_sos_constants =
10.1438 + [@{term "op ==>"}, @{term "Trueprop"},
10.1439 + @{term "op -->"}, @{term "op &"}, @{term "op |"},
10.1440 + @{term "Not"}, @{term "op = :: bool => _"},
10.1441 + @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"},
10.1442 + @{term "op = :: real => _"}, @{term "op < :: real => _"},
10.1443 + @{term "op <= :: real => _"},
10.1444 + @{term "op + :: real => _"}, @{term "op - :: real => _"},
10.1445 + @{term "op * :: real => _"}, @{term "uminus :: real => _"},
10.1446 + @{term "op / :: real => _"}, @{term "inverse :: real => _"},
10.1447 + @{term "op ^ :: real => _"}, @{term "abs :: real => _"},
10.1448 + @{term "min :: real => _"}, @{term "max :: real => _"},
10.1449 + @{term "0::real"}, @{term "1::real"}, @{term "number_of :: int => real"},
10.1450 + @{term "number_of :: int => nat"},
10.1451 + @{term "Int.Bit0"}, @{term "Int.Bit1"},
10.1452 + @{term "Int.Pls"}, @{term "Int.Min"}];
10.1453 +
10.1454 +fun check_sos kcts ct =
10.1455 + let
10.1456 + val t = term_of ct
10.1457 + val _ = if not (null (Term.add_tfrees t [])
10.1458 + andalso null (Term.add_tvars t []))
10.1459 + then error "SOS: not sos. Additional type varables" else ()
10.1460 + val fs = Term.add_frees t []
10.1461 + val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs
10.1462 + then error "SOS: not sos. Variables with type not real" else ()
10.1463 + val vs = Term.add_vars t []
10.1464 + val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs
10.1465 + then error "SOS: not sos. Variables with type not real" else ()
10.1466 + val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
10.1467 + val _ = if null ukcs then ()
10.1468 + else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
10.1469 +in () end
10.1470 +
10.1471 +fun core_sos_tac prover ctxt = CSUBGOAL (fn (ct, i) =>
10.1472 + let val _ = check_sos known_sos_constants ct
10.1473 + val (avs, p) = strip_all ct
10.1474 + val th = standard (fold_rev forall_intr avs (real_sos prover ctxt (Thm.dest_arg p)))
10.1475 + in rtac th i end);
10.1476 +
10.1477 +fun default_SOME f NONE v = SOME v
10.1478 + | default_SOME f (SOME v) _ = SOME v;
10.1479 +
10.1480 +fun lift_SOME f NONE a = f a
10.1481 + | lift_SOME f (SOME a) _ = SOME a;
10.1482 +
10.1483 +
10.1484 +local
10.1485 + val is_numeral = can (HOLogic.dest_number o term_of)
10.1486 +in
10.1487 +fun get_denom b ct = case term_of ct of
10.1488 + @{term "op / :: real => _"} $ _ $ _ =>
10.1489 + if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct)
10.1490 + else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct)) (Thm.dest_arg ct, b)
10.1491 + | @{term "op < :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
10.1492 + | @{term "op <= :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
10.1493 + | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct)
10.1494 + | _ => NONE
10.1495 +end;
10.1496 +
10.1497 +fun elim_one_denom_tac ctxt =
10.1498 +CSUBGOAL (fn (P,i) =>
10.1499 + case get_denom false P of
10.1500 + NONE => no_tac
10.1501 + | SOME (d,ord) =>
10.1502 + let
10.1503 + val ss = simpset_of ctxt addsimps @{thms field_simps}
10.1504 + addsimps [@{thm nonzero_power_divide}, @{thm power_divide}]
10.1505 + val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)]
10.1506 + (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto}
10.1507 + else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast})
10.1508 + in (rtac th i THEN Simplifier.asm_full_simp_tac ss i) end);
10.1509 +
10.1510 +fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
10.1511 +
10.1512 +fun sos_tac prover ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac prover ctxt
10.1513 +
10.1514 +
10.1515 +end;
11.1 --- a/src/HOL/Library/sos_wrapper.ML Wed Aug 05 17:10:10 2009 +0200
11.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
11.3 @@ -1,158 +0,0 @@
11.4 -(* Title: sos_wrapper.ML
11.5 - Author: Philipp Meyer, TU Muenchen
11.6 -
11.7 -Added functionality for sums of squares, e.g. calling a remote prover
11.8 -*)
11.9 -
11.10 -signature SOS_WRAPPER =
11.11 -sig
11.12 -
11.13 - datatype prover_result = Success | PartialSuccess | Failure | Error
11.14 - type prover = string * (int -> string -> prover_result * string)
11.15 -
11.16 - val setup: theory -> theory
11.17 -end
11.18 -
11.19 -structure SosWrapper : SOS_WRAPPER=
11.20 -struct
11.21 -
11.22 -datatype prover_result = Success | PartialSuccess | Failure | Error
11.23 -type prover =
11.24 - string * (* command name *)
11.25 - (int -> string ->prover_result * string) (* function to find failure from return value and output *)
11.26 -
11.27 -
11.28 -(*** output control ***)
11.29 -
11.30 -fun debug s = Output.debug (fn () => s)
11.31 -val answer = Output.priority
11.32 -val write = Output.writeln
11.33 -
11.34 -(*** calling provers ***)
11.35 -
11.36 -val destdir = ref ""
11.37 -
11.38 -fun filename dir name =
11.39 - let
11.40 - val probfile = Path.basic (name ^ serial_string ())
11.41 - in
11.42 - if dir = "" then
11.43 - File.tmp_path probfile
11.44 - else
11.45 - if File.exists (Path.explode dir) then
11.46 - Path.append (Path.explode dir) probfile
11.47 - else
11.48 - error ("No such directory: " ^ dir)
11.49 - end
11.50 -
11.51 -fun is_success Success = true
11.52 - | is_success PartialSuccess = true
11.53 - | is_success _ = false
11.54 -fun str_of_status Success = "Success"
11.55 - | str_of_status PartialSuccess = "Partial Success"
11.56 - | str_of_status Failure= "Failure"
11.57 - | str_of_status Error= "Error"
11.58 -
11.59 -fun run_solver name (cmd, find_failure) input =
11.60 - let
11.61 - val _ = answer ("Calling solver: " ^ name)
11.62 -
11.63 - (* create input file *)
11.64 - val dir = ! destdir
11.65 - val input_file = filename dir "sos_in"
11.66 - val _ = File.write input_file input
11.67 -
11.68 - val _ = debug "Solver input:"
11.69 - val _ = debug input
11.70 -
11.71 - (* call solver *)
11.72 - val output_file = filename dir "sos_out"
11.73 - val (output, rv) = system_out (cmd ^ " " ^ (Path.implode input_file) ^
11.74 - " " ^ (Path.implode output_file))
11.75 -
11.76 - (* read and analysize output *)
11.77 - val (res, res_msg) = find_failure rv output
11.78 - val result = if is_success res then File.read output_file else ""
11.79 -
11.80 - (* remove temporary files *)
11.81 - val _ = if dir = "" then (File.rm input_file ; if File.exists output_file then File.rm output_file else ()) else ()
11.82 -
11.83 - val _ = debug "Solver output:"
11.84 - val _ = debug output
11.85 - val _ = debug "Solver result:"
11.86 - val _ = debug result
11.87 -
11.88 - val _ = answer (str_of_status res ^ ": " ^ res_msg)
11.89 -
11.90 - in
11.91 - if is_success res then
11.92 - result
11.93 - else
11.94 - error ("Prover failed: " ^ res_msg)
11.95 - end
11.96 -
11.97 -(*** various provers ***)
11.98 -
11.99 -(* local csdp client *)
11.100 -
11.101 -fun find_csdp_run_failure rv _ =
11.102 - case rv of
11.103 - 0 => (Success, "SDP solved")
11.104 - | 1 => (Failure, "SDP is primal infeasible")
11.105 - | 2 => (Failure, "SDP is dual infeasible")
11.106 - | 3 => (PartialSuccess, "SDP solved with reduced accuracy")
11.107 - | _ => (Failure, "return code is " ^ string_of_int rv)
11.108 -
11.109 -val csdp = ("csdp", find_csdp_run_failure)
11.110 -
11.111 -(* remote neos server *)
11.112 -
11.113 -fun find_neos_failure rv output =
11.114 - if rv = 2 then (Failure, "no solution") else
11.115 - if rv <> 0 then (Error, "return code is " ^ string_of_int rv) else
11.116 - let
11.117 - fun find_success str =
11.118 - if String.isPrefix "Success: " str then
11.119 - SOME (Success, unprefix "Success: " str)
11.120 - else if String.isPrefix "Partial Success: " str then
11.121 - SOME (PartialSuccess, unprefix "Partial Success: " str)
11.122 - else if String.isPrefix "Failure: " str then
11.123 - SOME (Failure, unprefix "Failure: " str)
11.124 - else
11.125 - NONE
11.126 - val exit_line = get_first find_success (split_lines output)
11.127 - in
11.128 - case exit_line of
11.129 - SOME (status, msg) =>
11.130 - if String.isPrefix "SDP solved" msg then
11.131 - (status, msg)
11.132 - else (Failure, msg)
11.133 - | NONE => (Failure, "no success")
11.134 - end
11.135 -
11.136 -val neos_csdp = ("$ISABELLE_HOME/lib/scripts/neos/NeosCSDPClient.py", find_neos_failure)
11.137 -
11.138 -(* save provers in table *)
11.139 -
11.140 -val provers =
11.141 - Symtab.make [("remote_csdp", neos_csdp),("csdp", csdp)]
11.142 -
11.143 -fun get_prover name =
11.144 - case Symtab.lookup provers name of
11.145 - SOME prover => prover
11.146 - | NONE => error ("unknown prover: " ^ name)
11.147 -
11.148 -fun call_solver name =
11.149 - run_solver name (get_prover name)
11.150 -
11.151 -(* setup tactic *)
11.152 -
11.153 -val def_solver = "remote_csdp"
11.154 -
11.155 -fun sos_solver name = (SIMPLE_METHOD' o (Sos.sos_tac (call_solver name)))
11.156 -
11.157 -val sos_method = Scan.optional (Scan.lift OuterParse.xname) def_solver >> sos_solver
11.158 -
11.159 -val setup = Method.setup @{binding sos} sos_method "Prove universal problems over the reals using sums of squares"
11.160 -
11.161 -end
12.1 --- a/src/HOL/Library/sum_of_squares.ML Wed Aug 05 17:10:10 2009 +0200
12.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
12.3 @@ -1,1754 +0,0 @@
12.4 -(* Title: sum_of_squares.ML
12.5 - Authors: Amine Chaieb, University of Cambridge
12.6 - Philipp Meyer, TU Muenchen
12.7 -
12.8 -A tactic for proving nonlinear inequalities
12.9 -*)
12.10 -
12.11 -signature SOS =
12.12 -sig
12.13 -
12.14 - val sos_tac : (string -> string) -> Proof.context -> int -> Tactical.tactic
12.15 -
12.16 -end
12.17 -
12.18 -structure Sos : SOS =
12.19 -struct
12.20 -
12.21 -
12.22 -val rat_0 = Rat.zero;
12.23 -val rat_1 = Rat.one;
12.24 -val rat_2 = Rat.two;
12.25 -val rat_10 = Rat.rat_of_int 10;
12.26 -val rat_1_2 = rat_1 // rat_2;
12.27 -val max = curry IntInf.max;
12.28 -val min = curry IntInf.min;
12.29 -
12.30 -val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int;
12.31 -val numerator_rat = Rat.quotient_of_rat #> fst #> Rat.rat_of_int;
12.32 -fun int_of_rat a =
12.33 - case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int";
12.34 -fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y));
12.35 -
12.36 -fun rat_pow r i =
12.37 - let fun pow r i =
12.38 - if i = 0 then rat_1 else
12.39 - let val d = pow r (i div 2)
12.40 - in d */ d */ (if i mod 2 = 0 then rat_1 else r)
12.41 - end
12.42 - in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end;
12.43 -
12.44 -fun round_rat r =
12.45 - let val (a,b) = Rat.quotient_of_rat (Rat.abs r)
12.46 - val d = a div b
12.47 - val s = if r </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
12.48 - val x2 = 2 * (a - (b * d))
12.49 - in s (if x2 >= b then d + 1 else d) end
12.50 -
12.51 -val abs_rat = Rat.abs;
12.52 -val pow2 = rat_pow rat_2;
12.53 -val pow10 = rat_pow rat_10;
12.54 -
12.55 -val debugging = ref false;
12.56 -
12.57 -exception Sanity;
12.58 -
12.59 -exception Unsolvable;
12.60 -
12.61 -(* Turn a rational into a decimal string with d sig digits. *)
12.62 -
12.63 -local
12.64 -fun normalize y =
12.65 - if abs_rat y </ (rat_1 // rat_10) then normalize (rat_10 */ y) - 1
12.66 - else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1
12.67 - else 0
12.68 - in
12.69 -fun decimalize d x =
12.70 - if x =/ rat_0 then "0.0" else
12.71 - let
12.72 - val y = Rat.abs x
12.73 - val e = normalize y
12.74 - val z = pow10(~ e) */ y +/ rat_1
12.75 - val k = int_of_rat (round_rat(pow10 d */ z))
12.76 - in (if x </ rat_0 then "-0." else "0.") ^
12.77 - implode(tl(explode(string_of_int k))) ^
12.78 - (if e = 0 then "" else "e"^string_of_int e)
12.79 - end
12.80 -end;
12.81 -
12.82 -(* Iterations over numbers, and lists indexed by numbers. *)
12.83 -
12.84 -fun itern k l f a =
12.85 - case l of
12.86 - [] => a
12.87 - | h::t => itern (k + 1) t f (f h k a);
12.88 -
12.89 -fun iter (m,n) f a =
12.90 - if n < m then a
12.91 - else iter (m+1,n) f (f m a);
12.92 -
12.93 -(* The main types. *)
12.94 -
12.95 -fun strict_ord ord (x,y) = case ord (x,y) of LESS => LESS | _ => GREATER
12.96 -
12.97 -structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord);
12.98 -
12.99 -type vector = int* Rat.rat Intfunc.T;
12.100 -
12.101 -type matrix = (int*int)*(Rat.rat Intpairfunc.T);
12.102 -
12.103 -type monomial = int Ctermfunc.T;
12.104 -
12.105 -val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t))
12.106 - fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2)
12.107 -structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord)
12.108 -
12.109 -type poly = Rat.rat Monomialfunc.T;
12.110 -
12.111 - fun iszero (k,r) = r =/ rat_0;
12.112 -
12.113 -fun fold_rev2 f l1 l2 b =
12.114 - case (l1,l2) of
12.115 - ([],[]) => b
12.116 - | (h1::t1,h2::t2) => f h1 h2 (fold_rev2 f t1 t2 b)
12.117 - | _ => error "fold_rev2";
12.118 -
12.119 -(* Vectors. Conventionally indexed 1..n. *)
12.120 -
12.121 -fun vector_0 n = (n,Intfunc.undefined):vector;
12.122 -
12.123 -fun dim (v:vector) = fst v;
12.124 -
12.125 -fun vector_const c n =
12.126 - if c =/ rat_0 then vector_0 n
12.127 - else (n,fold_rev (fn k => Intfunc.update (k,c)) (1 upto n) Intfunc.undefined) :vector;
12.128 -
12.129 -val vector_1 = vector_const rat_1;
12.130 -
12.131 -fun vector_cmul c (v:vector) =
12.132 - let val n = dim v
12.133 - in if c =/ rat_0 then vector_0 n
12.134 - else (n,Intfunc.mapf (fn x => c */ x) (snd v))
12.135 - end;
12.136 -
12.137 -fun vector_neg (v:vector) = (fst v,Intfunc.mapf Rat.neg (snd v)) :vector;
12.138 -
12.139 -fun vector_add (v1:vector) (v2:vector) =
12.140 - let val m = dim v1
12.141 - val n = dim v2
12.142 - in if m <> n then error "vector_add: incompatible dimensions"
12.143 - else (n,Intfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd v1) (snd v2)) :vector
12.144 - end;
12.145 -
12.146 -fun vector_sub v1 v2 = vector_add v1 (vector_neg v2);
12.147 -
12.148 -fun vector_dot (v1:vector) (v2:vector) =
12.149 - let val m = dim v1
12.150 - val n = dim v2
12.151 - in if m <> n then error "vector_dot: incompatible dimensions"
12.152 - else Intfunc.fold (fn (i,x) => fn a => x +/ a)
12.153 - (Intfunc.combine (curry op */) (fn x => x =/ rat_0) (snd v1) (snd v2)) rat_0
12.154 - end;
12.155 -
12.156 -fun vector_of_list l =
12.157 - let val n = length l
12.158 - in (n,fold_rev2 (curry Intfunc.update) (1 upto n) l Intfunc.undefined) :vector
12.159 - end;
12.160 -
12.161 -(* Matrices; again rows and columns indexed from 1. *)
12.162 -
12.163 -fun matrix_0 (m,n) = ((m,n),Intpairfunc.undefined):matrix;
12.164 -
12.165 -fun dimensions (m:matrix) = fst m;
12.166 -
12.167 -fun matrix_const c (mn as (m,n)) =
12.168 - if m <> n then error "matrix_const: needs to be square"
12.169 - else if c =/ rat_0 then matrix_0 mn
12.170 - else (mn,fold_rev (fn k => Intpairfunc.update ((k,k), c)) (1 upto n) Intpairfunc.undefined) :matrix;;
12.171 -
12.172 -val matrix_1 = matrix_const rat_1;
12.173 -
12.174 -fun matrix_cmul c (m:matrix) =
12.175 - let val (i,j) = dimensions m
12.176 - in if c =/ rat_0 then matrix_0 (i,j)
12.177 - else ((i,j),Intpairfunc.mapf (fn x => c */ x) (snd m))
12.178 - end;
12.179 -
12.180 -fun matrix_neg (m:matrix) =
12.181 - (dimensions m, Intpairfunc.mapf Rat.neg (snd m)) :matrix;
12.182 -
12.183 -fun matrix_add (m1:matrix) (m2:matrix) =
12.184 - let val d1 = dimensions m1
12.185 - val d2 = dimensions m2
12.186 - in if d1 <> d2
12.187 - then error "matrix_add: incompatible dimensions"
12.188 - else (d1,Intpairfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd m1) (snd m2)) :matrix
12.189 - end;;
12.190 -
12.191 -fun matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2);
12.192 -
12.193 -fun row k (m:matrix) =
12.194 - let val (i,j) = dimensions m
12.195 - in (j,
12.196 - Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then Intfunc.update (j,c) a else a) (snd m) Intfunc.undefined ) : vector
12.197 - end;
12.198 -
12.199 -fun column k (m:matrix) =
12.200 - let val (i,j) = dimensions m
12.201 - in (i,
12.202 - Intpairfunc.fold (fn ((i,j), c) => fn a => if j = k then Intfunc.update (i,c) a else a) (snd m) Intfunc.undefined)
12.203 - : vector
12.204 - end;
12.205 -
12.206 -fun transp (m:matrix) =
12.207 - let val (i,j) = dimensions m
12.208 - in
12.209 - ((j,i),Intpairfunc.fold (fn ((i,j), c) => fn a => Intpairfunc.update ((j,i), c) a) (snd m) Intpairfunc.undefined) :matrix
12.210 - end;
12.211 -
12.212 -fun diagonal (v:vector) =
12.213 - let val n = dim v
12.214 - in ((n,n),Intfunc.fold (fn (i, c) => fn a => Intpairfunc.update ((i,i), c) a) (snd v) Intpairfunc.undefined) : matrix
12.215 - end;
12.216 -
12.217 -fun matrix_of_list l =
12.218 - let val m = length l
12.219 - in if m = 0 then matrix_0 (0,0) else
12.220 - let val n = length (hd l)
12.221 - in ((m,n),itern 1 l (fn v => fn i => itern 1 v (fn c => fn j => Intpairfunc.update ((i,j), c))) Intpairfunc.undefined)
12.222 - end
12.223 - end;
12.224 -
12.225 -(* Monomials. *)
12.226 -
12.227 -fun monomial_eval assig (m:monomial) =
12.228 - Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (Ctermfunc.apply assig x) k)
12.229 - m rat_1;
12.230 -val monomial_1 = (Ctermfunc.undefined:monomial);
12.231 -
12.232 -fun monomial_var x = Ctermfunc.onefunc (x, 1) :monomial;
12.233 -
12.234 -val (monomial_mul:monomial->monomial->monomial) =
12.235 - Ctermfunc.combine (curry op +) (K false);
12.236 -
12.237 -fun monomial_pow (m:monomial) k =
12.238 - if k = 0 then monomial_1
12.239 - else Ctermfunc.mapf (fn x => k * x) m;
12.240 -
12.241 -fun monomial_divides (m1:monomial) (m2:monomial) =
12.242 - Ctermfunc.fold (fn (x, k) => fn a => Ctermfunc.tryapplyd m2 x 0 >= k andalso a) m1 true;;
12.243 -
12.244 -fun monomial_div (m1:monomial) (m2:monomial) =
12.245 - let val m = Ctermfunc.combine (curry op +)
12.246 - (fn x => x = 0) m1 (Ctermfunc.mapf (fn x => ~ x) m2)
12.247 - in if Ctermfunc.fold (fn (x, k) => fn a => k >= 0 andalso a) m true then m
12.248 - else error "monomial_div: non-divisible"
12.249 - end;
12.250 -
12.251 -fun monomial_degree x (m:monomial) =
12.252 - Ctermfunc.tryapplyd m x 0;;
12.253 -
12.254 -fun monomial_lcm (m1:monomial) (m2:monomial) =
12.255 - fold_rev (fn x => Ctermfunc.update (x, max (monomial_degree x m1) (monomial_degree x m2)))
12.256 - (gen_union (is_equal o cterm_ord) (Ctermfunc.dom m1, Ctermfunc.dom m2)) (Ctermfunc.undefined :monomial);
12.257 -
12.258 -fun monomial_multidegree (m:monomial) =
12.259 - Ctermfunc.fold (fn (x, k) => fn a => k + a) m 0;;
12.260 -
12.261 -fun monomial_variables m = Ctermfunc.dom m;;
12.262 -
12.263 -(* Polynomials. *)
12.264 -
12.265 -fun eval assig (p:poly) =
12.266 - Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0;
12.267 -
12.268 -val poly_0 = (Monomialfunc.undefined:poly);
12.269 -
12.270 -fun poly_isconst (p:poly) =
12.271 - Monomialfunc.fold (fn (m, c) => fn a => Ctermfunc.is_undefined m andalso a) p true;
12.272 -
12.273 -fun poly_var x = Monomialfunc.onefunc (monomial_var x,rat_1) :poly;
12.274 -
12.275 -fun poly_const c =
12.276 - if c =/ rat_0 then poly_0 else Monomialfunc.onefunc(monomial_1, c);
12.277 -
12.278 -fun poly_cmul c (p:poly) =
12.279 - if c =/ rat_0 then poly_0
12.280 - else Monomialfunc.mapf (fn x => c */ x) p;
12.281 -
12.282 -fun poly_neg (p:poly) = (Monomialfunc.mapf Rat.neg p :poly);;
12.283 -
12.284 -fun poly_add (p1:poly) (p2:poly) =
12.285 - (Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2 :poly);
12.286 -
12.287 -fun poly_sub p1 p2 = poly_add p1 (poly_neg p2);
12.288 -
12.289 -fun poly_cmmul (c,m) (p:poly) =
12.290 - if c =/ rat_0 then poly_0
12.291 - else if Ctermfunc.is_undefined m
12.292 - then Monomialfunc.mapf (fn d => c */ d) p
12.293 - else Monomialfunc.fold (fn (m', d) => fn a => (Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0;
12.294 -
12.295 -fun poly_mul (p1:poly) (p2:poly) =
12.296 - Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0;
12.297 -
12.298 -fun poly_div (p1:poly) (p2:poly) =
12.299 - if not(poly_isconst p2)
12.300 - then error "poly_div: non-constant" else
12.301 - let val c = eval Ctermfunc.undefined p2
12.302 - in if c =/ rat_0 then error "poly_div: division by zero"
12.303 - else poly_cmul (Rat.inv c) p1
12.304 - end;
12.305 -
12.306 -fun poly_square p = poly_mul p p;
12.307 -
12.308 -fun poly_pow p k =
12.309 - if k = 0 then poly_const rat_1
12.310 - else if k = 1 then p
12.311 - else let val q = poly_square(poly_pow p (k div 2)) in
12.312 - if k mod 2 = 1 then poly_mul p q else q end;
12.313 -
12.314 -fun poly_exp p1 p2 =
12.315 - if not(poly_isconst p2)
12.316 - then error "poly_exp: not a constant"
12.317 - else poly_pow p1 (int_of_rat (eval Ctermfunc.undefined p2));
12.318 -
12.319 -fun degree x (p:poly) =
12.320 - Monomialfunc.fold (fn (m,c) => fn a => max (monomial_degree x m) a) p 0;
12.321 -
12.322 -fun multidegree (p:poly) =
12.323 - Monomialfunc.fold (fn (m, c) => fn a => max (monomial_multidegree m) a) p 0;
12.324 -
12.325 -fun poly_variables (p:poly) =
12.326 - sort cterm_ord (Monomialfunc.fold_rev (fn (m, c) => curry (gen_union (is_equal o cterm_ord)) (monomial_variables m)) p []);;
12.327 -
12.328 -(* Order monomials for human presentation. *)
12.329 -
12.330 -fun cterm_ord (t,t') = TermOrd.fast_term_ord (term_of t, term_of t');
12.331 -
12.332 -val humanorder_varpow = prod_ord cterm_ord (rev_order o int_ord);
12.333 -
12.334 -local
12.335 - fun ord (l1,l2) = case (l1,l2) of
12.336 - (_,[]) => LESS
12.337 - | ([],_) => GREATER
12.338 - | (h1::t1,h2::t2) =>
12.339 - (case humanorder_varpow (h1, h2) of
12.340 - LESS => LESS
12.341 - | EQUAL => ord (t1,t2)
12.342 - | GREATER => GREATER)
12.343 -in fun humanorder_monomial m1 m2 =
12.344 - ord (sort humanorder_varpow (Ctermfunc.graph m1),
12.345 - sort humanorder_varpow (Ctermfunc.graph m2))
12.346 -end;
12.347 -
12.348 -fun fold1 f l = case l of
12.349 - [] => error "fold1"
12.350 - | [x] => x
12.351 - | (h::t) => f h (fold1 f t);
12.352 -
12.353 -(* Conversions to strings. *)
12.354 -
12.355 -fun string_of_vector min_size max_size (v:vector) =
12.356 - let val n_raw = dim v
12.357 - in if n_raw = 0 then "[]" else
12.358 - let
12.359 - val n = max min_size (min n_raw max_size)
12.360 - val xs = map (Rat.string_of_rat o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
12.361 - in "[" ^ fold1 (fn s => fn t => s ^ ", " ^ t) xs ^
12.362 - (if n_raw > max_size then ", ...]" else "]")
12.363 - end
12.364 - end;
12.365 -
12.366 -fun string_of_matrix max_size (m:matrix) =
12.367 - let
12.368 - val (i_raw,j_raw) = dimensions m
12.369 - val i = min max_size i_raw
12.370 - val j = min max_size j_raw
12.371 - val rstr = map (fn k => string_of_vector j j (row k m)) (1 upto i)
12.372 - in "["^ fold1 (fn s => fn t => s^";\n "^t) rstr ^
12.373 - (if j > max_size then "\n ...]" else "]")
12.374 - end;
12.375 -
12.376 -fun string_of_term t =
12.377 - case t of
12.378 - a$b => "("^(string_of_term a)^" "^(string_of_term b)^")"
12.379 - | Abs x =>
12.380 - let val (xn, b) = Term.dest_abs x
12.381 - in "(\\"^xn^"."^(string_of_term b)^")"
12.382 - end
12.383 - | Const(s,_) => s
12.384 - | Free (s,_) => s
12.385 - | Var((s,_),_) => s
12.386 - | _ => error "string_of_term";
12.387 -
12.388 -val string_of_cterm = string_of_term o term_of;
12.389 -
12.390 -fun string_of_varpow x k =
12.391 - if k = 1 then string_of_cterm x
12.392 - else string_of_cterm x^"^"^string_of_int k;
12.393 -
12.394 -fun string_of_monomial m =
12.395 - if Ctermfunc.is_undefined m then "1" else
12.396 - let val vps = fold_rev (fn (x,k) => fn a => string_of_varpow x k :: a)
12.397 - (sort humanorder_varpow (Ctermfunc.graph m)) []
12.398 - in fold1 (fn s => fn t => s^"*"^t) vps
12.399 - end;
12.400 -
12.401 -fun string_of_cmonomial (c,m) =
12.402 - if Ctermfunc.is_undefined m then Rat.string_of_rat c
12.403 - else if c =/ rat_1 then string_of_monomial m
12.404 - else Rat.string_of_rat c ^ "*" ^ string_of_monomial m;;
12.405 -
12.406 -fun string_of_poly (p:poly) =
12.407 - if Monomialfunc.is_undefined p then "<<0>>" else
12.408 - let
12.409 - val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1 m2) (Monomialfunc.graph p)
12.410 - val s = fold (fn (m,c) => fn a =>
12.411 - if c </ rat_0 then a ^ " - " ^ string_of_cmonomial(Rat.neg c,m)
12.412 - else a ^ " + " ^ string_of_cmonomial(c,m))
12.413 - cms ""
12.414 - val s1 = String.substring (s, 0, 3)
12.415 - val s2 = String.substring (s, 3, String.size s - 3)
12.416 - in "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>"
12.417 - end;
12.418 -
12.419 -(* Conversion from HOL term. *)
12.420 -
12.421 -local
12.422 - val neg_tm = @{cterm "uminus :: real => _"}
12.423 - val add_tm = @{cterm "op + :: real => _"}
12.424 - val sub_tm = @{cterm "op - :: real => _"}
12.425 - val mul_tm = @{cterm "op * :: real => _"}
12.426 - val inv_tm = @{cterm "inverse :: real => _"}
12.427 - val div_tm = @{cterm "op / :: real => _"}
12.428 - val pow_tm = @{cterm "op ^ :: real => _"}
12.429 - val zero_tm = @{cterm "0:: real"}
12.430 - val is_numeral = can (HOLogic.dest_number o term_of)
12.431 - fun is_comb t = case t of _$_ => true | _ => false
12.432 - fun poly_of_term tm =
12.433 - if tm aconvc zero_tm then poly_0
12.434 - else if RealArith.is_ratconst tm
12.435 - then poly_const(RealArith.dest_ratconst tm)
12.436 - else
12.437 - (let val (lop,r) = Thm.dest_comb tm
12.438 - in if lop aconvc neg_tm then poly_neg(poly_of_term r)
12.439 - else if lop aconvc inv_tm then
12.440 - let val p = poly_of_term r
12.441 - in if poly_isconst p
12.442 - then poly_const(Rat.inv (eval Ctermfunc.undefined p))
12.443 - else error "poly_of_term: inverse of non-constant polyomial"
12.444 - end
12.445 - else (let val (opr,l) = Thm.dest_comb lop
12.446 - in
12.447 - if opr aconvc pow_tm andalso is_numeral r
12.448 - then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r)
12.449 - else if opr aconvc add_tm
12.450 - then poly_add (poly_of_term l) (poly_of_term r)
12.451 - else if opr aconvc sub_tm
12.452 - then poly_sub (poly_of_term l) (poly_of_term r)
12.453 - else if opr aconvc mul_tm
12.454 - then poly_mul (poly_of_term l) (poly_of_term r)
12.455 - else if opr aconvc div_tm
12.456 - then let
12.457 - val p = poly_of_term l
12.458 - val q = poly_of_term r
12.459 - in if poly_isconst q then poly_cmul (Rat.inv (eval Ctermfunc.undefined q)) p
12.460 - else error "poly_of_term: division by non-constant polynomial"
12.461 - end
12.462 - else poly_var tm
12.463 -
12.464 - end
12.465 - handle CTERM ("dest_comb",_) => poly_var tm)
12.466 - end
12.467 - handle CTERM ("dest_comb",_) => poly_var tm)
12.468 -in
12.469 -val poly_of_term = fn tm =>
12.470 - if type_of (term_of tm) = @{typ real} then poly_of_term tm
12.471 - else error "poly_of_term: term does not have real type"
12.472 -end;
12.473 -
12.474 -(* String of vector (just a list of space-separated numbers). *)
12.475 -
12.476 -fun sdpa_of_vector (v:vector) =
12.477 - let
12.478 - val n = dim v
12.479 - val strs = map (decimalize 20 o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
12.480 - in fold1 (fn x => fn y => x ^ " " ^ y) strs ^ "\n"
12.481 - end;
12.482 -
12.483 -fun increasing f ord (x,y) = ord (f x, f y);
12.484 -fun triple_int_ord ((a,b,c),(a',b',c')) =
12.485 - prod_ord int_ord (prod_ord int_ord int_ord)
12.486 - ((a,(b,c)),(a',(b',c')));
12.487 -structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord);
12.488 -
12.489 -(* String for block diagonal matrix numbered k. *)
12.490 -
12.491 -fun sdpa_of_blockdiagonal k m =
12.492 - let
12.493 - val pfx = string_of_int k ^" "
12.494 - val ents =
12.495 - Inttriplefunc.fold (fn ((b,i,j), c) => fn a => if i > j then a else ((b,i,j),c)::a) m []
12.496 - val entss = sort (increasing fst triple_int_ord ) ents
12.497 - in fold_rev (fn ((b,i,j),c) => fn a =>
12.498 - pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
12.499 - " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
12.500 - end;
12.501 -
12.502 -(* String for a matrix numbered k, in SDPA sparse format. *)
12.503 -
12.504 -fun sdpa_of_matrix k (m:matrix) =
12.505 - let
12.506 - val pfx = string_of_int k ^ " 1 "
12.507 - val ms = Intpairfunc.fold (fn ((i,j), c) => fn a => if i > j then a else ((i,j),c)::a) (snd m) []
12.508 - val mss = sort (increasing fst (prod_ord int_ord int_ord)) ms
12.509 - in fold_rev (fn ((i,j),c) => fn a =>
12.510 - pfx ^ string_of_int i ^ " " ^ string_of_int j ^
12.511 - " " ^ decimalize 20 c ^ "\n" ^ a) mss ""
12.512 - end;;
12.513 -
12.514 -(* ------------------------------------------------------------------------- *)
12.515 -(* String in SDPA sparse format for standard SDP problem: *)
12.516 -(* *)
12.517 -(* X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD *)
12.518 -(* Minimize obj_1 * v_1 + ... obj_m * v_m *)
12.519 -(* ------------------------------------------------------------------------- *)
12.520 -
12.521 -fun sdpa_of_problem obj mats =
12.522 - let
12.523 - val m = length mats - 1
12.524 - val (n,_) = dimensions (hd mats)
12.525 - in
12.526 - string_of_int m ^ "\n" ^
12.527 - "1\n" ^
12.528 - string_of_int n ^ "\n" ^
12.529 - sdpa_of_vector obj ^
12.530 - fold_rev2 (fn k => fn m => fn a => sdpa_of_matrix (k - 1) m ^ a) (1 upto length mats) mats ""
12.531 - end;
12.532 -
12.533 -fun index_char str chr pos =
12.534 - if pos >= String.size str then ~1
12.535 - else if String.sub(str,pos) = chr then pos
12.536 - else index_char str chr (pos + 1);
12.537 -fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b);
12.538 -fun rat_of_string s =
12.539 - let val n = index_char s #"/" 0 in
12.540 - if n = ~1 then s |> IntInf.fromString |> valOf |> Rat.rat_of_int
12.541 - else
12.542 - let val SOME numer = IntInf.fromString(String.substring(s,0,n))
12.543 - val SOME den = IntInf.fromString (String.substring(s,n+1,String.size s - n - 1))
12.544 - in rat_of_quotient(numer, den)
12.545 - end
12.546 - end;
12.547 -
12.548 -fun isspace x = x = " " ;
12.549 -fun isnum x = x mem_string ["0","1","2","3","4","5","6","7","8","9"]
12.550 -
12.551 -(* More parser basics. *)
12.552 -
12.553 -local
12.554 - open Scan
12.555 -in
12.556 - val word = this_string
12.557 - fun token s =
12.558 - repeat ($$ " ") |-- word s --| repeat ($$ " ")
12.559 - val numeral = one isnum
12.560 - val decimalint = bulk numeral >> (rat_of_string o implode)
12.561 - val decimalfrac = bulk numeral
12.562 - >> (fn s => rat_of_string(implode s) // pow10 (length s))
12.563 - val decimalsig =
12.564 - decimalint -- option (Scan.$$ "." |-- decimalfrac)
12.565 - >> (fn (h,NONE) => h | (h,SOME x) => h +/ x)
12.566 - fun signed prs =
12.567 - $$ "-" |-- prs >> Rat.neg
12.568 - || $$ "+" |-- prs
12.569 - || prs;
12.570 -
12.571 -fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs
12.572 -
12.573 - val exponent = ($$ "e" || $$ "E") |-- signed decimalint;
12.574 -
12.575 - val decimal = signed decimalsig -- (emptyin rat_0|| exponent)
12.576 - >> (fn (h, x) => h */ pow10 (int_of_rat x));
12.577 -end;
12.578 -
12.579 - fun mkparser p s =
12.580 - let val (x,rst) = p (explode s)
12.581 - in if null rst then x
12.582 - else error "mkparser: unparsed input"
12.583 - end;;
12.584 -val parse_decimal = mkparser decimal;
12.585 -
12.586 -fun fix err prs =
12.587 - prs || (fn x=> error err);
12.588 -
12.589 -fun listof prs sep err =
12.590 - prs -- Scan.bulk (sep |-- fix err prs) >> uncurry cons;
12.591 -
12.592 -(* Parse back a vector. *)
12.593 -
12.594 - val vector =
12.595 - token "{" |-- listof decimal (token ",") "decimal" --| token "}"
12.596 - >> vector_of_list
12.597 - val parse_vector = mkparser vector
12.598 - fun skipupto dscr prs inp =
12.599 - (dscr |-- prs
12.600 - || Scan.one (K true) |-- skipupto dscr prs) inp
12.601 - fun ignore inp = ((),[])
12.602 - fun sdpaoutput inp = skipupto (word "xVec" -- token "=")
12.603 - (vector --| ignore) inp
12.604 - fun csdpoutput inp = ((decimal -- Scan.bulk (Scan.$$ " " |-- Scan.option decimal) >> (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp
12.605 - val parse_sdpaoutput = mkparser sdpaoutput
12.606 - val parse_csdpoutput = mkparser csdpoutput
12.607 -
12.608 -(* Run prover on a problem in linear form. *)
12.609 -
12.610 -fun run_problem prover obj mats =
12.611 - parse_csdpoutput (prover (sdpa_of_problem obj mats))
12.612 -
12.613 -(*
12.614 -UNUSED
12.615 -
12.616 -(* Also parse the SDPA output to test success (CSDP yields a return code). *)
12.617 -
12.618 -local
12.619 - val prs =
12.620 - skipupto (word "phase.value" -- token "=")
12.621 - (Scan.option (Scan.$$ "p") -- Scan.option (Scan.$$ "d")
12.622 - -- (word "OPT" || word "FEAS"))
12.623 -in
12.624 - fun sdpa_run_succeeded s =
12.625 - (prs (explode s); true) handle _ => false
12.626 -end;
12.627 -
12.628 -(* The default parameters. Unfortunately this goes to a fixed file. *)
12.629 -
12.630 -val sdpa_default_parameters =
12.631 -"100 unsigned int maxIteration; \n1.0E-7 double 0.0 < epsilonStar;\n1.0E2 double 0.0 < lambdaStar;\n2.0 double 1.0 < omegaStar;\n-1.0E5 double lowerBound;\n1.0E5 double upperBound;\n0.1 double 0.0 <= betaStar < 1.0;\n0.2 double 0.0 <= betaBar < 1.0, betaStar <= betaBar;\n0.9 double 0.0 < gammaStar < 1.0;\n1.0E-7 double 0.0 < epsilonDash;\n";;
12.632 -
12.633 -(* These were suggested by Makoto Yamashita for problems where we are *)
12.634 -(* right at the edge of the semidefinite cone, as sometimes happens. *)
12.635 -
12.636 -val sdpa_alt_parameters =
12.637 -"1000 unsigned int maxIteration;\n1.0E-7 double 0.0 < epsilonStar;\n1.0E4 double 0.0 < lambdaStar;\n2.0 double 1.0 < omegaStar;\n-1.0E5 double lowerBound;\n1.0E5 double upperBound;\n0.1 double 0.0 <= betaStar < 1.0;\n0.2 double 0.0 <= betaBar < 1.0, betaStar <= betaBar;\n0.9 double 0.0 < gammaStar < 1.0;\n1.0E-7 double 0.0 < epsilonDash;\n";;
12.638 -
12.639 -val sdpa_params = sdpa_alt_parameters;;
12.640 -
12.641 -(* CSDP parameters; so far I'm sticking with the defaults. *)
12.642 -
12.643 -val csdp_default_parameters =
12.644 -"axtol=1.0e-8\natytol=1.0e-8\nobjtol=1.0e-8\npinftol=1.0e8\ndinftol=1.0e8\nmaxiter=100\nminstepfrac=0.9\nmaxstepfrac=0.97\nminstepp=1.0e-8\nminstepd=1.0e-8\nusexzgap=1\ntweakgap=0\naffine=0\nprintlevel=1\n";;
12.645 -
12.646 -val csdp_params = csdp_default_parameters;;
12.647 -
12.648 -fun tmp_file pre suf =
12.649 - let val i = string_of_int (round (random()))
12.650 - val name = Path.append (Path.variable "ISABELLE_TMP") (Path.explode (pre ^ i ^ suf))
12.651 - in
12.652 - if File.exists name then tmp_file pre suf
12.653 - else name
12.654 - end;
12.655 -
12.656 -(* Now call SDPA on a problem and parse back the output. *)
12.657 -
12.658 -fun run_sdpa dbg obj mats =
12.659 - let
12.660 - val input_file = tmp_file "sos" ".dat-s"
12.661 - val output_file = tmp_file "sos" ".out"
12.662 - val params_file = tmp_file "param" ".sdpa"
12.663 - val current_dir = File.pwd()
12.664 - val _ = File.write input_file
12.665 - (sdpa_of_problem "" obj mats)
12.666 - val _ = File.write params_file sdpa_params
12.667 - val _ = File.cd (Path.variable "ISABELLE_TMP")
12.668 - val _ = File.system_command ("sdpa "^ (Path.implode input_file) ^ " " ^
12.669 - (Path.implode output_file) ^
12.670 - (if dbg then "" else "> /dev/null"))
12.671 - val opr = File.read output_file
12.672 - in if not(sdpa_run_succeeded opr) then error "sdpa: call failed"
12.673 - else
12.674 - let val res = parse_sdpaoutput opr
12.675 - in ((if dbg then ()
12.676 - else (File.rm input_file; File.rm output_file ; File.cd current_dir));
12.677 - res)
12.678 - end
12.679 - end;
12.680 -
12.681 -fun sdpa obj mats = run_sdpa (!debugging) obj mats;
12.682 -
12.683 -(* The same thing with CSDP. *)
12.684 -
12.685 -fun run_csdp dbg obj mats =
12.686 - let
12.687 - val input_file = tmp_file "sos" ".dat-s"
12.688 - val output_file = tmp_file "sos" ".out"
12.689 - val params_file = tmp_file "param" ".csdp"
12.690 - val current_dir = File.pwd()
12.691 - val _ = File.write input_file (sdpa_of_problem "" obj mats)
12.692 - val _ = File.write params_file csdp_params
12.693 - val _ = File.cd (Path.variable "ISABELLE_TMP")
12.694 - val rv = system ("csdp "^(Path.implode input_file) ^ " "
12.695 - ^ (Path.implode output_file) ^
12.696 - (if dbg then "" else "> /dev/null"))
12.697 - val opr = File.read output_file
12.698 - val res = parse_csdpoutput opr
12.699 - in
12.700 - ((if dbg then ()
12.701 - else (File.rm input_file; File.rm output_file ; File.cd current_dir));
12.702 - (rv,res))
12.703 - end;
12.704 -
12.705 -fun csdp obj mats =
12.706 - let
12.707 - val (rv,res) = run_csdp (!debugging) obj mats
12.708 - in
12.709 - ((if rv = 1 orelse rv = 2 then error "csdp: Problem is infeasible"
12.710 - else if rv = 3 then writeln "csdp warning: Reduced accuracy"
12.711 - else if rv <> 0 then error ("csdp: error "^string_of_int rv)
12.712 - else ());
12.713 - res)
12.714 - end;
12.715 -
12.716 -*)
12.717 -
12.718 -(* Try some apparently sensible scaling first. Note that this is purely to *)
12.719 -(* get a cleaner translation to floating-point, and doesn't affect any of *)
12.720 -(* the results, in principle. In practice it seems a lot better when there *)
12.721 -(* are extreme numbers in the original problem. *)
12.722 -
12.723 - (* Version for (int*int) keys *)
12.724 -local
12.725 - fun max_rat x y = if x </ y then y else x
12.726 - fun common_denominator fld amat acc =
12.727 - fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
12.728 - fun maximal_element fld amat acc =
12.729 - fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
12.730 -fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
12.731 - in Real.fromLargeInt a / Real.fromLargeInt b end;
12.732 -in
12.733 -
12.734 -fun pi_scale_then solver (obj:vector) mats =
12.735 - let
12.736 - val cd1 = fold_rev (common_denominator Intpairfunc.fold) mats (rat_1)
12.737 - val cd2 = common_denominator Intfunc.fold (snd obj) (rat_1)
12.738 - val mats' = map (Intpairfunc.mapf (fn x => cd1 */ x)) mats
12.739 - val obj' = vector_cmul cd2 obj
12.740 - val max1 = fold_rev (maximal_element Intpairfunc.fold) mats' (rat_0)
12.741 - val max2 = maximal_element Intfunc.fold (snd obj') (rat_0)
12.742 - val scal1 = pow2 (20 - trunc(Math.ln (float_of_rat max1) / Math.ln 2.0))
12.743 - val scal2 = pow2 (20 - trunc(Math.ln (float_of_rat max2) / Math.ln 2.0))
12.744 - val mats'' = map (Intpairfunc.mapf (fn x => x */ scal1)) mats'
12.745 - val obj'' = vector_cmul scal2 obj'
12.746 - in solver obj'' mats''
12.747 - end
12.748 -end;
12.749 -
12.750 -(* Try some apparently sensible scaling first. Note that this is purely to *)
12.751 -(* get a cleaner translation to floating-point, and doesn't affect any of *)
12.752 -(* the results, in principle. In practice it seems a lot better when there *)
12.753 -(* are extreme numbers in the original problem. *)
12.754 -
12.755 - (* Version for (int*int*int) keys *)
12.756 -local
12.757 - fun max_rat x y = if x </ y then y else x
12.758 - fun common_denominator fld amat acc =
12.759 - fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
12.760 - fun maximal_element fld amat acc =
12.761 - fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
12.762 -fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
12.763 - in Real.fromLargeInt a / Real.fromLargeInt b end;
12.764 -fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0)
12.765 -in
12.766 -
12.767 -fun tri_scale_then solver (obj:vector) mats =
12.768 - let
12.769 - val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1)
12.770 - val cd2 = common_denominator Intfunc.fold (snd obj) (rat_1)
12.771 - val mats' = map (Inttriplefunc.mapf (fn x => cd1 */ x)) mats
12.772 - val obj' = vector_cmul cd2 obj
12.773 - val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0)
12.774 - val max2 = maximal_element Intfunc.fold (snd obj') (rat_0)
12.775 - val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0))
12.776 - val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0))
12.777 - val mats'' = map (Inttriplefunc.mapf (fn x => x */ scal1)) mats'
12.778 - val obj'' = vector_cmul scal2 obj'
12.779 - in solver obj'' mats''
12.780 - end
12.781 -end;
12.782 -
12.783 -(* Round a vector to "nice" rationals. *)
12.784 -
12.785 -fun nice_rational n x = round_rat (n */ x) // n;;
12.786 -fun nice_vector n ((d,v) : vector) =
12.787 - (d, Intfunc.fold (fn (i,c) => fn a =>
12.788 - let val y = nice_rational n c
12.789 - in if c =/ rat_0 then a
12.790 - else Intfunc.update (i,y) a end) v Intfunc.undefined):vector
12.791 -
12.792 -(*
12.793 -UNUSED
12.794 -
12.795 -(* Reduce linear program to SDP (diagonal matrices) and test with CSDP. This *)
12.796 -(* one tests A [-1;x1;..;xn] >= 0 (i.e. left column is negated constants). *)
12.797 -
12.798 -fun linear_program_basic a =
12.799 - let
12.800 - val (m,n) = dimensions a
12.801 - val mats = map (fn j => diagonal (column j a)) (1 upto n)
12.802 - val obj = vector_const rat_1 m
12.803 - val (rv,res) = run_csdp false obj mats
12.804 - in if rv = 1 orelse rv = 2 then false
12.805 - else if rv = 0 then true
12.806 - else error "linear_program: An error occurred in the SDP solver"
12.807 - end;
12.808 -
12.809 -(* Alternative interface testing A x >= b for matrix A, vector b. *)
12.810 -
12.811 -fun linear_program a b =
12.812 - let val (m,n) = dimensions a
12.813 - in if dim b <> m then error "linear_program: incompatible dimensions"
12.814 - else
12.815 - let
12.816 - val mats = diagonal b :: map (fn j => diagonal (column j a)) (1 upto n)
12.817 - val obj = vector_const rat_1 m
12.818 - val (rv,res) = run_csdp false obj mats
12.819 - in if rv = 1 orelse rv = 2 then false
12.820 - else if rv = 0 then true
12.821 - else error "linear_program: An error occurred in the SDP solver"
12.822 - end
12.823 - end;
12.824 -
12.825 -(* Test whether a point is in the convex hull of others. Rather than use *)
12.826 -(* computational geometry, express as linear inequalities and call CSDP. *)
12.827 -(* This is a bit lazy of me, but it's easy and not such a bottleneck so far. *)
12.828 -
12.829 -fun in_convex_hull pts pt =
12.830 - let
12.831 - val pts1 = (1::pt) :: map (fn x => 1::x) pts
12.832 - val pts2 = map (fn p => map (fn x => ~x) p @ p) pts1
12.833 - val n = length pts + 1
12.834 - val v = 2 * (length pt + 1)
12.835 - val m = v + n - 1
12.836 - val mat = ((m,n),
12.837 - itern 1 pts2 (fn pts => fn j => itern 1 pts
12.838 - (fn x => fn i => Intpairfunc.update ((i,j), Rat.rat_of_int x)))
12.839 - (iter (1,n) (fn i => Intpairfunc.update((v + i,i+1), rat_1))
12.840 - Intpairfunc.undefined))
12.841 - in linear_program_basic mat
12.842 - end;
12.843 -
12.844 -(* Filter down a set of points to a minimal set with the same convex hull. *)
12.845 -
12.846 -local
12.847 - fun augment1 (m::ms) = if in_convex_hull ms m then ms else ms@[m]
12.848 - fun augment m ms = funpow 3 augment1 (m::ms)
12.849 -in
12.850 -fun minimal_convex_hull mons =
12.851 - let val mons' = fold_rev augment (tl mons) [hd mons]
12.852 - in funpow (length mons') augment1 mons'
12.853 - end
12.854 -end;
12.855 -
12.856 -*)
12.857 -
12.858 -fun dest_ord f x = is_equal (f x);
12.859 -
12.860 -
12.861 -
12.862 -(* Stuff for "equations" ((int*int*int)->num functions). *)
12.863 -
12.864 -fun tri_equation_cmul c eq =
12.865 - if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
12.866 -
12.867 -fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
12.868 -
12.869 -fun tri_equation_eval assig eq =
12.870 - let fun value v = Inttriplefunc.apply assig v
12.871 - in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
12.872 - end;
12.873 -
12.874 -(* Eliminate among linear equations: return unconstrained variables and *)
12.875 -(* assignments for the others in terms of them. We give one pseudo-variable *)
12.876 -(* "one" that's used for a constant term. *)
12.877 -
12.878 -local
12.879 - fun extract_first p l = case l of (* FIXME : use find_first instead *)
12.880 - [] => error "extract_first"
12.881 - | h::t => if p h then (h,t) else
12.882 - let val (k,s) = extract_first p t in (k,h::s) end
12.883 -fun eliminate vars dun eqs = case vars of
12.884 - [] => if forall Inttriplefunc.is_undefined eqs then dun
12.885 - else raise Unsolvable
12.886 - | v::vs =>
12.887 - ((let
12.888 - val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs
12.889 - val a = Inttriplefunc.apply eq v
12.890 - val eq' = tri_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
12.891 - fun elim e =
12.892 - let val b = Inttriplefunc.tryapplyd e v rat_0
12.893 - in if b =/ rat_0 then e else
12.894 - tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
12.895 - end
12.896 - in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
12.897 - end)
12.898 - handle ERROR _ => eliminate vs dun eqs)
12.899 -in
12.900 -fun tri_eliminate_equations one vars eqs =
12.901 - let
12.902 - val assig = eliminate vars Inttriplefunc.undefined eqs
12.903 - val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
12.904 - in (distinct (dest_ord triple_int_ord) vs, assig)
12.905 - end
12.906 -end;
12.907 -
12.908 -(* Eliminate all variables, in an essentially arbitrary order. *)
12.909 -
12.910 -fun tri_eliminate_all_equations one =
12.911 - let
12.912 - fun choose_variable eq =
12.913 - let val (v,_) = Inttriplefunc.choose eq
12.914 - in if is_equal (triple_int_ord(v,one)) then
12.915 - let val eq' = Inttriplefunc.undefine v eq
12.916 - in if Inttriplefunc.is_undefined eq' then error "choose_variable"
12.917 - else fst (Inttriplefunc.choose eq')
12.918 - end
12.919 - else v
12.920 - end
12.921 - fun eliminate dun eqs = case eqs of
12.922 - [] => dun
12.923 - | eq::oeqs =>
12.924 - if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
12.925 - let val v = choose_variable eq
12.926 - val a = Inttriplefunc.apply eq v
12.927 - val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a)
12.928 - (Inttriplefunc.undefine v eq)
12.929 - fun elim e =
12.930 - let val b = Inttriplefunc.tryapplyd e v rat_0
12.931 - in if b =/ rat_0 then e
12.932 - else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
12.933 - end
12.934 - in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun))
12.935 - (map elim oeqs)
12.936 - end
12.937 -in fn eqs =>
12.938 - let
12.939 - val assig = eliminate Inttriplefunc.undefined eqs
12.940 - val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
12.941 - in (distinct (dest_ord triple_int_ord) vs,assig)
12.942 - end
12.943 -end;
12.944 -
12.945 -(* Solve equations by assigning arbitrary numbers. *)
12.946 -
12.947 -fun tri_solve_equations one eqs =
12.948 - let
12.949 - val (vars,assigs) = tri_eliminate_all_equations one eqs
12.950 - val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars
12.951 - (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
12.952 - val ass =
12.953 - Inttriplefunc.combine (curry op +/) (K false)
12.954 - (Inttriplefunc.mapf (tri_equation_eval vfn) assigs) vfn
12.955 - in if forall (fn e => tri_equation_eval ass e =/ rat_0) eqs
12.956 - then Inttriplefunc.undefine one ass else raise Sanity
12.957 - end;
12.958 -
12.959 -(* Multiply equation-parametrized poly by regular poly and add accumulator. *)
12.960 -
12.961 -fun tri_epoly_pmul p q acc =
12.962 - Monomialfunc.fold (fn (m1, c) => fn a =>
12.963 - Monomialfunc.fold (fn (m2,e) => fn b =>
12.964 - let val m = monomial_mul m1 m2
12.965 - val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined
12.966 - in Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b
12.967 - end) q a) p acc ;
12.968 -
12.969 -(* Usual operations on equation-parametrized poly. *)
12.970 -
12.971 -fun tri_epoly_cmul c l =
12.972 - if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (tri_equation_cmul c) l;;
12.973 -
12.974 -val tri_epoly_neg = tri_epoly_cmul (Rat.rat_of_int ~1);
12.975 -
12.976 -val tri_epoly_add = Inttriplefunc.combine tri_equation_add Inttriplefunc.is_undefined;
12.977 -
12.978 -fun tri_epoly_sub p q = tri_epoly_add p (tri_epoly_neg q);;
12.979 -
12.980 -(* Stuff for "equations" ((int*int)->num functions). *)
12.981 -
12.982 -fun pi_equation_cmul c eq =
12.983 - if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
12.984 -
12.985 -fun pi_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
12.986 -
12.987 -fun pi_equation_eval assig eq =
12.988 - let fun value v = Inttriplefunc.apply assig v
12.989 - in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
12.990 - end;
12.991 -
12.992 -(* Eliminate among linear equations: return unconstrained variables and *)
12.993 -(* assignments for the others in terms of them. We give one pseudo-variable *)
12.994 -(* "one" that's used for a constant term. *)
12.995 -
12.996 -local
12.997 -fun extract_first p l = case l of
12.998 - [] => error "extract_first"
12.999 - | h::t => if p h then (h,t) else
12.1000 - let val (k,s) = extract_first p t in (k,h::s) end
12.1001 -fun eliminate vars dun eqs = case vars of
12.1002 - [] => if forall Inttriplefunc.is_undefined eqs then dun
12.1003 - else raise Unsolvable
12.1004 - | v::vs =>
12.1005 - let
12.1006 - val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs
12.1007 - val a = Inttriplefunc.apply eq v
12.1008 - val eq' = pi_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
12.1009 - fun elim e =
12.1010 - let val b = Inttriplefunc.tryapplyd e v rat_0
12.1011 - in if b =/ rat_0 then e else
12.1012 - pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
12.1013 - end
12.1014 - in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
12.1015 - end
12.1016 - handle ERROR _ => eliminate vs dun eqs
12.1017 -in
12.1018 -fun pi_eliminate_equations one vars eqs =
12.1019 - let
12.1020 - val assig = eliminate vars Inttriplefunc.undefined eqs
12.1021 - val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
12.1022 - in (distinct (dest_ord triple_int_ord) vs, assig)
12.1023 - end
12.1024 -end;
12.1025 -
12.1026 -(* Eliminate all variables, in an essentially arbitrary order. *)
12.1027 -
12.1028 -fun pi_eliminate_all_equations one =
12.1029 - let
12.1030 - fun choose_variable eq =
12.1031 - let val (v,_) = Inttriplefunc.choose eq
12.1032 - in if is_equal (triple_int_ord(v,one)) then
12.1033 - let val eq' = Inttriplefunc.undefine v eq
12.1034 - in if Inttriplefunc.is_undefined eq' then error "choose_variable"
12.1035 - else fst (Inttriplefunc.choose eq')
12.1036 - end
12.1037 - else v
12.1038 - end
12.1039 - fun eliminate dun eqs = case eqs of
12.1040 - [] => dun
12.1041 - | eq::oeqs =>
12.1042 - if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
12.1043 - let val v = choose_variable eq
12.1044 - val a = Inttriplefunc.apply eq v
12.1045 - val eq' = pi_equation_cmul ((Rat.rat_of_int ~1) // a)
12.1046 - (Inttriplefunc.undefine v eq)
12.1047 - fun elim e =
12.1048 - let val b = Inttriplefunc.tryapplyd e v rat_0
12.1049 - in if b =/ rat_0 then e
12.1050 - else pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
12.1051 - end
12.1052 - in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun))
12.1053 - (map elim oeqs)
12.1054 - end
12.1055 -in fn eqs =>
12.1056 - let
12.1057 - val assig = eliminate Inttriplefunc.undefined eqs
12.1058 - val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
12.1059 - in (distinct (dest_ord triple_int_ord) vs,assig)
12.1060 - end
12.1061 -end;
12.1062 -
12.1063 -(* Solve equations by assigning arbitrary numbers. *)
12.1064 -
12.1065 -fun pi_solve_equations one eqs =
12.1066 - let
12.1067 - val (vars,assigs) = pi_eliminate_all_equations one eqs
12.1068 - val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars
12.1069 - (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
12.1070 - val ass =
12.1071 - Inttriplefunc.combine (curry op +/) (K false)
12.1072 - (Inttriplefunc.mapf (pi_equation_eval vfn) assigs) vfn
12.1073 - in if forall (fn e => pi_equation_eval ass e =/ rat_0) eqs
12.1074 - then Inttriplefunc.undefine one ass else raise Sanity
12.1075 - end;
12.1076 -
12.1077 -(* Multiply equation-parametrized poly by regular poly and add accumulator. *)
12.1078 -
12.1079 -fun pi_epoly_pmul p q acc =
12.1080 - Monomialfunc.fold (fn (m1, c) => fn a =>
12.1081 - Monomialfunc.fold (fn (m2,e) => fn b =>
12.1082 - let val m = monomial_mul m1 m2
12.1083 - val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined
12.1084 - in Monomialfunc.update (m,pi_equation_add (pi_equation_cmul c e) es) b
12.1085 - end) q a) p acc ;
12.1086 -
12.1087 -(* Usual operations on equation-parametrized poly. *)
12.1088 -
12.1089 -fun pi_epoly_cmul c l =
12.1090 - if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (pi_equation_cmul c) l;;
12.1091 -
12.1092 -val pi_epoly_neg = pi_epoly_cmul (Rat.rat_of_int ~1);
12.1093 -
12.1094 -val pi_epoly_add = Inttriplefunc.combine pi_equation_add Inttriplefunc.is_undefined;
12.1095 -
12.1096 -fun pi_epoly_sub p q = pi_epoly_add p (pi_epoly_neg q);;
12.1097 -
12.1098 -fun allpairs f l1 l2 = fold_rev (fn x => (curry (op @)) (map (f x) l2)) l1 [];
12.1099 -
12.1100 -(* Hence produce the "relevant" monomials: those whose squares lie in the *)
12.1101 -(* Newton polytope of the monomials in the input. (This is enough according *)
12.1102 -(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal, *)
12.1103 -(* vol 45, pp. 363--374, 1978. *)
12.1104 -(* *)
12.1105 -(* These are ordered in sort of decreasing degree. In particular the *)
12.1106 -(* constant monomial is last; this gives an order in diagonalization of the *)
12.1107 -(* quadratic form that will tend to display constants. *)
12.1108 -
12.1109 -(*
12.1110 -UNUSED
12.1111 -
12.1112 -fun newton_polytope pol =
12.1113 - let
12.1114 - val vars = poly_variables pol
12.1115 - val mons = map (fn m => map (fn x => monomial_degree x m) vars)
12.1116 - (Monomialfunc.dom pol)
12.1117 - val ds = map (fn x => (degree x pol + 1) div 2) vars
12.1118 - val all = fold_rev (fn n => allpairs cons (0 upto n)) ds [[]]
12.1119 - val mons' = minimal_convex_hull mons
12.1120 - val all' =
12.1121 - filter (fn m => in_convex_hull mons' (map (fn x => 2 * x) m)) all
12.1122 - in map (fn m => fold_rev2 (fn v => fn i => fn a => if i = 0 then a else Ctermfunc.update (v,i) a)
12.1123 - vars m monomial_1) (rev all')
12.1124 - end;
12.1125 -
12.1126 -*)
12.1127 -
12.1128 -(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form. *)
12.1129 -
12.1130 -local
12.1131 -fun diagonalize n i m =
12.1132 - if Intpairfunc.is_undefined (snd m) then []
12.1133 - else
12.1134 - let val a11 = Intpairfunc.tryapplyd (snd m) (i,i) rat_0
12.1135 - in if a11 </ rat_0 then error "diagonalize: not PSD"
12.1136 - else if a11 =/ rat_0 then
12.1137 - if Intfunc.is_undefined (snd (row i m)) then diagonalize n (i + 1) m
12.1138 - else error "diagonalize: not PSD ___ "
12.1139 - else
12.1140 - let
12.1141 - val v = row i m
12.1142 - val v' = (fst v, Intfunc.fold (fn (i, c) => fn a =>
12.1143 - let val y = c // a11
12.1144 - in if y = rat_0 then a else Intfunc.update (i,y) a
12.1145 - end) (snd v) Intfunc.undefined)
12.1146 - fun upt0 x y a = if y = rat_0 then a else Intpairfunc.update (x,y) a
12.1147 - val m' =
12.1148 - ((n,n),
12.1149 - iter (i+1,n) (fn j =>
12.1150 - iter (i+1,n) (fn k =>
12.1151 - (upt0 (j,k) (Intpairfunc.tryapplyd (snd m) (j,k) rat_0 -/ Intfunc.tryapplyd (snd v) j rat_0 */ Intfunc.tryapplyd (snd v') k rat_0))))
12.1152 - Intpairfunc.undefined)
12.1153 - in (a11,v')::diagonalize n (i + 1) m'
12.1154 - end
12.1155 - end
12.1156 -in
12.1157 -fun diag m =
12.1158 - let
12.1159 - val nn = dimensions m
12.1160 - val n = fst nn
12.1161 - in if snd nn <> n then error "diagonalize: non-square matrix"
12.1162 - else diagonalize n 1 m
12.1163 - end
12.1164 -end;
12.1165 -
12.1166 -fun gcd_rat a b = Rat.rat_of_int (Integer.gcd (int_of_rat a) (int_of_rat b));
12.1167 -
12.1168 -(* Adjust a diagonalization to collect rationals at the start. *)
12.1169 - (* FIXME : Potentially polymorphic keys, but here only: integers!! *)
12.1170 -local
12.1171 - fun upd0 x y a = if y =/ rat_0 then a else Intfunc.update(x,y) a;
12.1172 - fun mapa f (d,v) =
12.1173 - (d, Intfunc.fold (fn (i,c) => fn a => upd0 i (f c) a) v Intfunc.undefined)
12.1174 - fun adj (c,l) =
12.1175 - let val a =
12.1176 - Intfunc.fold (fn (i,c) => fn a => lcm_rat a (denominator_rat c))
12.1177 - (snd l) rat_1 //
12.1178 - Intfunc.fold (fn (i,c) => fn a => gcd_rat a (numerator_rat c))
12.1179 - (snd l) rat_0
12.1180 - in ((c // (a */ a)),mapa (fn x => a */ x) l)
12.1181 - end
12.1182 -in
12.1183 -fun deration d = if null d then (rat_0,d) else
12.1184 - let val d' = map adj d
12.1185 - val a = fold (lcm_rat o denominator_rat o fst) d' rat_1 //
12.1186 - fold (gcd_rat o numerator_rat o fst) d' rat_0
12.1187 - in ((rat_1 // a),map (fn (c,l) => (a */ c,l)) d')
12.1188 - end
12.1189 -end;
12.1190 -
12.1191 -(* Enumeration of monomials with given multidegree bound. *)
12.1192 -
12.1193 -fun enumerate_monomials d vars =
12.1194 - if d < 0 then []
12.1195 - else if d = 0 then [Ctermfunc.undefined]
12.1196 - else if null vars then [monomial_1] else
12.1197 - let val alts =
12.1198 - map (fn k => let val oths = enumerate_monomials (d - k) (tl vars)
12.1199 - in map (fn ks => if k = 0 then ks else Ctermfunc.update (hd vars, k) ks) oths end) (0 upto d)
12.1200 - in fold1 (curry op @) alts
12.1201 - end;
12.1202 -
12.1203 -(* Enumerate products of distinct input polys with degree <= d. *)
12.1204 -(* We ignore any constant input polynomials. *)
12.1205 -(* Give the output polynomial and a record of how it was derived. *)
12.1206 -
12.1207 -local
12.1208 - open RealArith
12.1209 -in
12.1210 -fun enumerate_products d pols =
12.1211 -if d = 0 then [(poly_const rat_1,Rational_lt rat_1)]
12.1212 -else if d < 0 then [] else
12.1213 -case pols of
12.1214 - [] => [(poly_const rat_1,Rational_lt rat_1)]
12.1215 - | (p,b)::ps =>
12.1216 - let val e = multidegree p
12.1217 - in if e = 0 then enumerate_products d ps else
12.1218 - enumerate_products d ps @
12.1219 - map (fn (q,c) => (poly_mul p q,Product(b,c)))
12.1220 - (enumerate_products (d - e) ps)
12.1221 - end
12.1222 -end;
12.1223 -
12.1224 -(* Convert regular polynomial. Note that we treat (0,0,0) as -1. *)
12.1225 -
12.1226 -fun epoly_of_poly p =
12.1227 - Monomialfunc.fold (fn (m,c) => fn a => Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p Monomialfunc.undefined;
12.1228 -
12.1229 -(* String for block diagonal matrix numbered k. *)
12.1230 -
12.1231 -fun sdpa_of_blockdiagonal k m =
12.1232 - let
12.1233 - val pfx = string_of_int k ^" "
12.1234 - val ents =
12.1235 - Inttriplefunc.fold
12.1236 - (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a)
12.1237 - m []
12.1238 - val entss = sort (increasing fst triple_int_ord) ents
12.1239 - in fold_rev (fn ((b,i,j),c) => fn a =>
12.1240 - pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
12.1241 - " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
12.1242 - end;
12.1243 -
12.1244 -(* SDPA for problem using block diagonal (i.e. multiple SDPs) *)
12.1245 -
12.1246 -fun sdpa_of_blockproblem nblocks blocksizes obj mats =
12.1247 - let val m = length mats - 1
12.1248 - in
12.1249 - string_of_int m ^ "\n" ^
12.1250 - string_of_int nblocks ^ "\n" ^
12.1251 - (fold1 (fn s => fn t => s^" "^t) (map string_of_int blocksizes)) ^
12.1252 - "\n" ^
12.1253 - sdpa_of_vector obj ^
12.1254 - fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a)
12.1255 - (1 upto length mats) mats ""
12.1256 - end;
12.1257 -
12.1258 -(* Run prover on a problem in block diagonal form. *)
12.1259 -
12.1260 -fun run_blockproblem prover nblocks blocksizes obj mats=
12.1261 - parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats))
12.1262 -
12.1263 -(*
12.1264 -UNUSED
12.1265 -
12.1266 -(* Hence run CSDP on a problem in block diagonal form. *)
12.1267 -
12.1268 -fun run_csdp dbg nblocks blocksizes obj mats =
12.1269 - let
12.1270 - val input_file = tmp_file "sos" ".dat-s"
12.1271 - val output_file = tmp_file "sos" ".out"
12.1272 - val params_file = tmp_file "param" ".csdp"
12.1273 - val _ = File.write input_file
12.1274 - (sdpa_of_blockproblem "" nblocks blocksizes obj mats)
12.1275 - val _ = File.write params_file csdp_params
12.1276 - val current_dir = File.pwd()
12.1277 - val _ = File.cd (Path.variable "ISABELLE_TMP")
12.1278 - val rv = system ("csdp "^(Path.implode input_file) ^ " "
12.1279 - ^ (Path.implode output_file) ^
12.1280 - (if dbg then "" else "> /dev/null"))
12.1281 - val opr = File.read output_file
12.1282 - val res = parse_csdpoutput opr
12.1283 - in
12.1284 - ((if dbg then ()
12.1285 - else (File.rm input_file ; File.rm output_file ; File.cd current_dir));
12.1286 - (rv,res))
12.1287 - end;
12.1288 -
12.1289 -fun csdp nblocks blocksizes obj mats =
12.1290 - let
12.1291 - val (rv,res) = run_csdp (!debugging) nblocks blocksizes obj mats
12.1292 - in ((if rv = 1 orelse rv = 2 then error "csdp: Problem is infeasible"
12.1293 - else if rv = 3 then writeln "csdp warning: Reduced accuracy"
12.1294 - else if rv <> 0 then error ("csdp: error "^string_of_int rv)
12.1295 - else ());
12.1296 - res)
12.1297 - end;
12.1298 -*)
12.1299 -
12.1300 -(* 3D versions of matrix operations to consider blocks separately. *)
12.1301 -
12.1302 -val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0);
12.1303 -fun bmatrix_cmul c bm =
12.1304 - if c =/ rat_0 then Inttriplefunc.undefined
12.1305 - else Inttriplefunc.mapf (fn x => c */ x) bm;
12.1306 -
12.1307 -val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1);
12.1308 -fun bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);;
12.1309 -
12.1310 -(* Smash a block matrix into components. *)
12.1311 -
12.1312 -fun blocks blocksizes bm =
12.1313 - map (fn (bs,b0) =>
12.1314 - let val m = Inttriplefunc.fold
12.1315 - (fn ((b,i,j),c) => fn a => if b = b0 then Intpairfunc.update ((i,j),c) a else a) bm Intpairfunc.undefined
12.1316 - val d = Intpairfunc.fold (fn ((i,j),c) => fn a => max a (max i j)) m 0
12.1317 - in (((bs,bs),m):matrix) end)
12.1318 - (blocksizes ~~ (1 upto length blocksizes));;
12.1319 -
12.1320 -(* FIXME : Get rid of this !!!*)
12.1321 -local
12.1322 - fun tryfind_with msg f [] = error msg
12.1323 - | tryfind_with msg f (x::xs) = (f x handle ERROR s => tryfind_with s f xs);
12.1324 -in
12.1325 - fun tryfind f = tryfind_with "tryfind" f
12.1326 -end
12.1327 -
12.1328 -(*
12.1329 -fun tryfind f [] = error "tryfind"
12.1330 - | tryfind f (x::xs) = (f x handle ERROR _ => tryfind f xs);
12.1331 -*)
12.1332 -
12.1333 -(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
12.1334 -
12.1335 -
12.1336 -local
12.1337 - open RealArith
12.1338 -in
12.1339 -fun real_positivnullstellensatz_general prover linf d eqs leqs pol =
12.1340 -let
12.1341 - val vars = fold_rev (curry (gen_union (op aconvc)) o poly_variables)
12.1342 - (pol::eqs @ map fst leqs) []
12.1343 - val monoid = if linf then
12.1344 - (poly_const rat_1,Rational_lt rat_1)::
12.1345 - (filter (fn (p,c) => multidegree p <= d) leqs)
12.1346 - else enumerate_products d leqs
12.1347 - val nblocks = length monoid
12.1348 - fun mk_idmultiplier k p =
12.1349 - let
12.1350 - val e = d - multidegree p
12.1351 - val mons = enumerate_monomials e vars
12.1352 - val nons = mons ~~ (1 upto length mons)
12.1353 - in (mons,
12.1354 - fold_rev (fn (m,n) => Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons Monomialfunc.undefined)
12.1355 - end
12.1356 -
12.1357 - fun mk_sqmultiplier k (p,c) =
12.1358 - let
12.1359 - val e = (d - multidegree p) div 2
12.1360 - val mons = enumerate_monomials e vars
12.1361 - val nons = mons ~~ (1 upto length mons)
12.1362 - in (mons,
12.1363 - fold_rev (fn (m1,n1) =>
12.1364 - fold_rev (fn (m2,n2) => fn a =>
12.1365 - let val m = monomial_mul m1 m2
12.1366 - in if n1 > n2 then a else
12.1367 - let val c = if n1 = n2 then rat_1 else rat_2
12.1368 - val e = Monomialfunc.tryapplyd a m Inttriplefunc.undefined
12.1369 - in Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a
12.1370 - end
12.1371 - end) nons)
12.1372 - nons Monomialfunc.undefined)
12.1373 - end
12.1374 -
12.1375 - val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid)
12.1376 - val (idmonlist,ids) = split_list(map2 mk_idmultiplier (1 upto length eqs) eqs)
12.1377 - val blocksizes = map length sqmonlist
12.1378 - val bigsum =
12.1379 - fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids
12.1380 - (fold_rev2 (fn (p,c) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs
12.1381 - (epoly_of_poly(poly_neg pol)))
12.1382 - val eqns = Monomialfunc.fold (fn (m,e) => fn a => e::a) bigsum []
12.1383 - val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns
12.1384 - val qvars = (0,0,0)::pvs
12.1385 - val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig
12.1386 - fun mk_matrix v =
12.1387 - Inttriplefunc.fold (fn ((b,i,j), ass) => fn m =>
12.1388 - if b < 0 then m else
12.1389 - let val c = Inttriplefunc.tryapplyd ass v rat_0
12.1390 - in if c = rat_0 then m else
12.1391 - Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m)
12.1392 - end)
12.1393 - allassig Inttriplefunc.undefined
12.1394 - val diagents = Inttriplefunc.fold
12.1395 - (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a)
12.1396 - allassig Inttriplefunc.undefined
12.1397 -
12.1398 - val mats = map mk_matrix qvars
12.1399 - val obj = (length pvs,
12.1400 - itern 1 pvs (fn v => fn i => Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0))
12.1401 - Intfunc.undefined)
12.1402 - val raw_vec = if null pvs then vector_0 0
12.1403 - else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats
12.1404 - fun int_element (d,v) i = Intfunc.tryapplyd v i rat_0
12.1405 - fun cterm_element (d,v) i = Ctermfunc.tryapplyd v i rat_0
12.1406 -
12.1407 - fun find_rounding d =
12.1408 - let
12.1409 - val _ = if !debugging
12.1410 - then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n")
12.1411 - else ()
12.1412 - val vec = nice_vector d raw_vec
12.1413 - val blockmat = iter (1,dim vec)
12.1414 - (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a)
12.1415 - (bmatrix_neg (nth mats 0))
12.1416 - val allmats = blocks blocksizes blockmat
12.1417 - in (vec,map diag allmats)
12.1418 - end
12.1419 - val (vec,ratdias) =
12.1420 - if null pvs then find_rounding rat_1
12.1421 - else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @
12.1422 - map pow2 (5 upto 66))
12.1423 - val newassigs =
12.1424 - fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k))
12.1425 - (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1))
12.1426 - val finalassigs =
12.1427 - Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs
12.1428 - fun poly_of_epoly p =
12.1429 - Monomialfunc.fold (fn (v,e) => fn a => Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a)
12.1430 - p Monomialfunc.undefined
12.1431 - fun mk_sos mons =
12.1432 - let fun mk_sq (c,m) =
12.1433 - (c,fold_rev (fn k=> fn a => Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a)
12.1434 - (1 upto length mons) Monomialfunc.undefined)
12.1435 - in map mk_sq
12.1436 - end
12.1437 - val sqs = map2 mk_sos sqmonlist ratdias
12.1438 - val cfs = map poly_of_epoly ids
12.1439 - val msq = filter (fn (a,b) => not (null b)) (map2 pair monoid sqs)
12.1440 - fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0
12.1441 - val sanity =
12.1442 - fold_rev (fn ((p,c),s) => poly_add (poly_mul p (eval_sq s))) msq
12.1443 - (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs
12.1444 - (poly_neg pol))
12.1445 -
12.1446 -in if not(Monomialfunc.is_undefined sanity) then raise Sanity else
12.1447 - (cfs,map (fn (a,b) => (snd a,b)) msq)
12.1448 - end
12.1449 -
12.1450 -
12.1451 -end;
12.1452 -
12.1453 -(* Iterative deepening. *)
12.1454 -
12.1455 -fun deepen f n =
12.1456 - (writeln ("Searching with depth limit " ^ string_of_int n) ; (f n handle ERROR s => (writeln ("failed with message: " ^ s) ; deepen f (n+1))))
12.1457 -
12.1458 -(* The ordering so we can create canonical HOL polynomials. *)
12.1459 -
12.1460 -fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon);
12.1461 -
12.1462 -fun monomial_order (m1,m2) =
12.1463 - if Ctermfunc.is_undefined m2 then LESS
12.1464 - else if Ctermfunc.is_undefined m1 then GREATER
12.1465 - else
12.1466 - let val mon1 = dest_monomial m1
12.1467 - val mon2 = dest_monomial m2
12.1468 - val deg1 = fold (curry op + o snd) mon1 0
12.1469 - val deg2 = fold (curry op + o snd) mon2 0
12.1470 - in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS
12.1471 - else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2)
12.1472 - end;
12.1473 -
12.1474 -fun dest_poly p =
12.1475 - map (fn (m,c) => (c,dest_monomial m))
12.1476 - (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p));
12.1477 -
12.1478 -(* Map back polynomials and their composites to HOL. *)
12.1479 -
12.1480 -local
12.1481 - open Thm Numeral RealArith
12.1482 -in
12.1483 -
12.1484 -fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x)
12.1485 - (mk_cnumber @{ctyp nat} k)
12.1486 -
12.1487 -fun cterm_of_monomial m =
12.1488 - if Ctermfunc.is_undefined m then @{cterm "1::real"}
12.1489 - else
12.1490 - let
12.1491 - val m' = dest_monomial m
12.1492 - val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' []
12.1493 - in fold1 (fn s => fn t => capply (capply @{cterm "op * :: real => _"} s) t) vps
12.1494 - end
12.1495 -
12.1496 -fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c
12.1497 - else if c = Rat.one then cterm_of_monomial m
12.1498 - else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
12.1499 -
12.1500 -fun cterm_of_poly p =
12.1501 - if Monomialfunc.is_undefined p then @{cterm "0::real"}
12.1502 - else
12.1503 - let
12.1504 - val cms = map cterm_of_cmonomial
12.1505 - (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
12.1506 - in fold1 (fn t1 => fn t2 => capply(capply @{cterm "op + :: real => _"} t1) t2) cms
12.1507 - end;
12.1508 -
12.1509 -fun cterm_of_sqterm (c,p) = Product(Rational_lt c,Square(cterm_of_poly p));
12.1510 -
12.1511 -fun cterm_of_sos (pr,sqs) = if null sqs then pr
12.1512 - else Product(pr,fold1 (fn a => fn b => Sum(a,b)) (map cterm_of_sqterm sqs));
12.1513 -
12.1514 -end
12.1515 -
12.1516 -(* Interface to HOL. *)
12.1517 -local
12.1518 - open Thm Conv RealArith
12.1519 - val concl = dest_arg o cprop_of
12.1520 - fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
12.1521 -in
12.1522 - (* FIXME: Replace tryfind by get_first !! *)
12.1523 -fun real_nonlinear_prover prover ctxt =
12.1524 - let
12.1525 - val {add,mul,neg,pow,sub,main} = Normalizer.semiring_normalizers_ord_wrapper ctxt
12.1526 - (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"}))
12.1527 - simple_cterm_ord
12.1528 - val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
12.1529 - real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
12.1530 - fun mainf translator (eqs,les,lts) =
12.1531 - let
12.1532 - val eq0 = map (poly_of_term o dest_arg1 o concl) eqs
12.1533 - val le0 = map (poly_of_term o dest_arg o concl) les
12.1534 - val lt0 = map (poly_of_term o dest_arg o concl) lts
12.1535 - val eqp0 = map (fn (t,i) => (t,Axiom_eq i)) (eq0 ~~ (0 upto (length eq0 - 1)))
12.1536 - val lep0 = map (fn (t,i) => (t,Axiom_le i)) (le0 ~~ (0 upto (length le0 - 1)))
12.1537 - val ltp0 = map (fn (t,i) => (t,Axiom_lt i)) (lt0 ~~ (0 upto (length lt0 - 1)))
12.1538 - val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0
12.1539 - val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0
12.1540 - val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0
12.1541 - fun trivial_axiom (p,ax) =
12.1542 - case ax of
12.1543 - Axiom_eq n => if eval Ctermfunc.undefined p <>/ Rat.zero then nth eqs n
12.1544 - else error "trivial_axiom: Not a trivial axiom"
12.1545 - | Axiom_le n => if eval Ctermfunc.undefined p </ Rat.zero then nth les n
12.1546 - else error "trivial_axiom: Not a trivial axiom"
12.1547 - | Axiom_lt n => if eval Ctermfunc.undefined p <=/ Rat.zero then nth lts n
12.1548 - else error "trivial_axiom: Not a trivial axiom"
12.1549 - | _ => error "trivial_axiom: Not a trivial axiom"
12.1550 - in
12.1551 - ((let val th = tryfind trivial_axiom (keq @ klep @ kltp)
12.1552 - in fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv field_comp_conv) th end)
12.1553 - handle ERROR _ => (
12.1554 - let
12.1555 - val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
12.1556 - val leq = lep @ ltp
12.1557 - fun tryall d =
12.1558 - let val e = multidegree pol
12.1559 - val k = if e = 0 then 0 else d div e
12.1560 - val eq' = map fst eq
12.1561 - in tryfind (fn i => (d,i,real_positivnullstellensatz_general prover false d eq' leq
12.1562 - (poly_neg(poly_pow pol i))))
12.1563 - (0 upto k)
12.1564 - end
12.1565 - val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
12.1566 - val proofs_ideal =
12.1567 - map2 (fn q => fn (p,ax) => Eqmul(cterm_of_poly q,ax)) cert_ideal eq
12.1568 - val proofs_cone = map cterm_of_sos cert_cone
12.1569 - val proof_ne = if null ltp then Rational_lt Rat.one else
12.1570 - let val p = fold1 (fn s => fn t => Product(s,t)) (map snd ltp)
12.1571 - in funpow i (fn q => Product(p,q)) (Rational_lt Rat.one)
12.1572 - end
12.1573 - val proof = fold1 (fn s => fn t => Sum(s,t))
12.1574 - (proof_ne :: proofs_ideal @ proofs_cone)
12.1575 - in writeln "Translating proof certificate to HOL";
12.1576 - translator (eqs,les,lts) proof
12.1577 - end))
12.1578 - end
12.1579 - in mainf end
12.1580 -end
12.1581 -
12.1582 -fun C f x y = f y x;
12.1583 - (* FIXME : This is very bad!!!*)
12.1584 -fun subst_conv eqs t =
12.1585 - let
12.1586 - val t' = fold (Thm.cabs o Thm.lhs_of) eqs t
12.1587 - in Conv.fconv_rule (Thm.beta_conversion true) (fold (C combination) eqs (reflexive t'))
12.1588 - end
12.1589 -
12.1590 -(* A wrapper that tries to substitute away variables first. *)
12.1591 -
12.1592 -local
12.1593 - open Thm Conv RealArith
12.1594 - fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
12.1595 - val concl = dest_arg o cprop_of
12.1596 - val shuffle1 =
12.1597 - fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps) })
12.1598 - val shuffle2 =
12.1599 - fconv_rule (rewr_conv @{lemma "(x + a == y) == (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps)})
12.1600 - fun substitutable_monomial fvs tm = case term_of tm of
12.1601 - Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm)
12.1602 - else error "substitutable_monomial"
12.1603 - | @{term "op * :: real => _"}$c$(t as Free _ ) =>
12.1604 - if is_ratconst (dest_arg1 tm) andalso not (member (op aconvc) fvs (dest_arg tm))
12.1605 - then (dest_ratconst (dest_arg1 tm),dest_arg tm) else error "substitutable_monomial"
12.1606 - | @{term "op + :: real => _"}$s$t =>
12.1607 - (substitutable_monomial (add_cterm_frees (dest_arg tm) fvs) (dest_arg1 tm)
12.1608 - handle ERROR _ => substitutable_monomial (add_cterm_frees (dest_arg1 tm) fvs) (dest_arg tm))
12.1609 - | _ => error "substitutable_monomial"
12.1610 -
12.1611 - fun isolate_variable v th =
12.1612 - let val w = dest_arg1 (cprop_of th)
12.1613 - in if v aconvc w then th
12.1614 - else case term_of w of
12.1615 - @{term "op + :: real => _"}$s$t =>
12.1616 - if dest_arg1 w aconvc v then shuffle2 th
12.1617 - else isolate_variable v (shuffle1 th)
12.1618 - | _ => error "isolate variable : This should not happen?"
12.1619 - end
12.1620 -in
12.1621 -
12.1622 -fun real_nonlinear_subst_prover prover ctxt =
12.1623 - let
12.1624 - val {add,mul,neg,pow,sub,main} = Normalizer.semiring_normalizers_ord_wrapper ctxt
12.1625 - (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"}))
12.1626 - simple_cterm_ord
12.1627 -
12.1628 - val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
12.1629 - real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
12.1630 -
12.1631 - fun make_substitution th =
12.1632 - let
12.1633 - val (c,v) = substitutable_monomial [] (dest_arg1(concl th))
12.1634 - val th1 = Drule.arg_cong_rule (capply @{cterm "op * :: real => _"} (cterm_of_rat (Rat.inv c))) (mk_meta_eq th)
12.1635 - val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1
12.1636 - in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2)
12.1637 - end
12.1638 - fun oprconv cv ct =
12.1639 - let val g = Thm.dest_fun2 ct
12.1640 - in if g aconvc @{cterm "op <= :: real => _"}
12.1641 - orelse g aconvc @{cterm "op < :: real => _"}
12.1642 - then arg_conv cv ct else arg1_conv cv ct
12.1643 - end
12.1644 - fun mainf translator =
12.1645 - let
12.1646 - fun substfirst(eqs,les,lts) =
12.1647 - ((let
12.1648 - val eth = tryfind make_substitution eqs
12.1649 - val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv)))
12.1650 - in substfirst
12.1651 - (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t
12.1652 - aconvc @{cterm "0::real"}) (map modify eqs),
12.1653 - map modify les,map modify lts)
12.1654 - end)
12.1655 - handle ERROR _ => real_nonlinear_prover prover ctxt translator (rev eqs, rev les, rev lts))
12.1656 - in substfirst
12.1657 - end
12.1658 -
12.1659 -
12.1660 - in mainf
12.1661 - end
12.1662 -
12.1663 -(* Overall function. *)
12.1664 -
12.1665 -fun real_sos prover ctxt t = gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) t;
12.1666 -end;
12.1667 -
12.1668 -(* A tactic *)
12.1669 -fun strip_all ct =
12.1670 - case term_of ct of
12.1671 - Const("all",_) $ Abs (xn,xT,p) =>
12.1672 - let val (a,(v,t')) = (apsnd (Thm.dest_abs (SOME xn)) o Thm.dest_comb) ct
12.1673 - in apfst (cons v) (strip_all t')
12.1674 - end
12.1675 -| _ => ([],ct)
12.1676 -
12.1677 -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})
12.1678 -
12.1679 -val known_sos_constants =
12.1680 - [@{term "op ==>"}, @{term "Trueprop"},
12.1681 - @{term "op -->"}, @{term "op &"}, @{term "op |"},
12.1682 - @{term "Not"}, @{term "op = :: bool => _"},
12.1683 - @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"},
12.1684 - @{term "op = :: real => _"}, @{term "op < :: real => _"},
12.1685 - @{term "op <= :: real => _"},
12.1686 - @{term "op + :: real => _"}, @{term "op - :: real => _"},
12.1687 - @{term "op * :: real => _"}, @{term "uminus :: real => _"},
12.1688 - @{term "op / :: real => _"}, @{term "inverse :: real => _"},
12.1689 - @{term "op ^ :: real => _"}, @{term "abs :: real => _"},
12.1690 - @{term "min :: real => _"}, @{term "max :: real => _"},
12.1691 - @{term "0::real"}, @{term "1::real"}, @{term "number_of :: int => real"},
12.1692 - @{term "number_of :: int => nat"},
12.1693 - @{term "Int.Bit0"}, @{term "Int.Bit1"},
12.1694 - @{term "Int.Pls"}, @{term "Int.Min"}];
12.1695 -
12.1696 -fun check_sos kcts ct =
12.1697 - let
12.1698 - val t = term_of ct
12.1699 - val _ = if not (null (Term.add_tfrees t [])
12.1700 - andalso null (Term.add_tvars t []))
12.1701 - then error "SOS: not sos. Additional type varables" else ()
12.1702 - val fs = Term.add_frees t []
12.1703 - val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs
12.1704 - then error "SOS: not sos. Variables with type not real" else ()
12.1705 - val vs = Term.add_vars t []
12.1706 - val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs
12.1707 - then error "SOS: not sos. Variables with type not real" else ()
12.1708 - val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
12.1709 - val _ = if null ukcs then ()
12.1710 - else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
12.1711 -in () end
12.1712 -
12.1713 -fun core_sos_tac prover ctxt = CSUBGOAL (fn (ct, i) =>
12.1714 - let val _ = check_sos known_sos_constants ct
12.1715 - val (avs, p) = strip_all ct
12.1716 - val th = standard (fold_rev forall_intr avs (real_sos prover ctxt (Thm.dest_arg p)))
12.1717 - in rtac th i end);
12.1718 -
12.1719 -fun default_SOME f NONE v = SOME v
12.1720 - | default_SOME f (SOME v) _ = SOME v;
12.1721 -
12.1722 -fun lift_SOME f NONE a = f a
12.1723 - | lift_SOME f (SOME a) _ = SOME a;
12.1724 -
12.1725 -
12.1726 -local
12.1727 - val is_numeral = can (HOLogic.dest_number o term_of)
12.1728 -in
12.1729 -fun get_denom b ct = case term_of ct of
12.1730 - @{term "op / :: real => _"} $ _ $ _ =>
12.1731 - if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct)
12.1732 - else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct)) (Thm.dest_arg ct, b)
12.1733 - | @{term "op < :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
12.1734 - | @{term "op <= :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
12.1735 - | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct)
12.1736 - | _ => NONE
12.1737 -end;
12.1738 -
12.1739 -fun elim_one_denom_tac ctxt =
12.1740 -CSUBGOAL (fn (P,i) =>
12.1741 - case get_denom false P of
12.1742 - NONE => no_tac
12.1743 - | SOME (d,ord) =>
12.1744 - let
12.1745 - val ss = simpset_of ctxt addsimps @{thms field_simps}
12.1746 - addsimps [@{thm nonzero_power_divide}, @{thm power_divide}]
12.1747 - val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)]
12.1748 - (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto}
12.1749 - else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast})
12.1750 - in (rtac th i THEN Simplifier.asm_full_simp_tac ss i) end);
12.1751 -
12.1752 -fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
12.1753 -
12.1754 -fun sos_tac prover ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac prover ctxt
12.1755 -
12.1756 -
12.1757 -end;