16784
|
1 |
(* Title: HOL/Matrix/cplex/FloatSparseMatrixBuilder.ML
|
|
2 |
ID: $Id$
|
|
3 |
Author: Steven Obua
|
|
4 |
*)
|
|
5 |
|
22964
|
6 |
signature FLOAT_SPARSE_MATIRX_BUILDER =
|
16784
|
7 |
sig
|
22964
|
8 |
include MATRIX_BUILDER
|
16784
|
9 |
|
22964
|
10 |
structure cplex : CPLEX
|
16784
|
11 |
|
22964
|
12 |
type float = Float.float
|
|
13 |
val approx_value : int -> (float -> float) -> string -> term * term
|
|
14 |
val approx_vector : int -> (float -> float) -> vector -> term * term
|
|
15 |
val approx_matrix : int -> (float -> float) -> matrix -> term * term
|
16784
|
16 |
|
22964
|
17 |
val mk_spvec_entry : Intt.int -> float -> term
|
|
18 |
val mk_spmat_entry : Intt.int -> term -> term
|
|
19 |
val spvecT: typ
|
|
20 |
val spmatT: typ
|
|
21 |
|
|
22 |
val v_elem_at : vector -> int -> string option
|
|
23 |
val m_elem_at : matrix -> int -> vector option
|
|
24 |
val v_only_elem : vector -> int option
|
|
25 |
val v_fold : (int * string -> 'a -> 'a) -> vector -> 'a -> 'a
|
|
26 |
val m_fold : (int * vector -> 'a -> 'a) -> matrix -> 'a -> 'a
|
20485
|
27 |
|
22964
|
28 |
val transpose_matrix : matrix -> matrix
|
16784
|
29 |
|
22964
|
30 |
val cut_vector : int -> vector -> vector
|
|
31 |
val cut_matrix : vector -> int option -> matrix -> matrix
|
22951
|
32 |
|
22964
|
33 |
(* cplexProg c A b *)
|
|
34 |
val cplexProg : vector -> matrix -> vector -> cplex.cplexProg * (string -> int)
|
|
35 |
(* dual_cplexProg c A b *)
|
|
36 |
val dual_cplexProg : vector -> matrix -> vector -> cplex.cplexProg * (string -> int)
|
|
37 |
end;
|
16784
|
38 |
|
22964
|
39 |
structure FloatSparseMatrixBuilder : FLOAT_SPARSE_MATIRX_BUILDER =
|
16784
|
40 |
struct
|
|
41 |
|
22964
|
42 |
type float = Float.float
|
|
43 |
structure Inttab = TableFun(type key = int val ord = rev_order o int_ord);
|
16784
|
44 |
|
|
45 |
type vector = string Inttab.table
|
|
46 |
type matrix = vector Inttab.table
|
20485
|
47 |
|
22964
|
48 |
val spvec_elemT = HOLogic.mk_prodT (HOLogic.natT, HOLogic.realT);
|
|
49 |
val spvecT = HOLogic.listT spvec_elemT;
|
|
50 |
val spmat_elemT = HOLogic.mk_prodT (HOLogic.natT, spvecT);
|
|
51 |
val spmatT = HOLogic.listT spmat_elemT;
|
16784
|
52 |
|
22964
|
53 |
fun approx_value prec f =
|
|
54 |
FloatArith.approx_float prec (fn (x, y) => (f x, f y));
|
16784
|
55 |
|
22951
|
56 |
fun mk_spvec_entry i f =
|
22964
|
57 |
HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, FloatArith.mk_float f);
|
16784
|
58 |
|
22951
|
59 |
fun mk_spmat_entry i e =
|
22964
|
60 |
HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, e);
|
16784
|
61 |
|
|
62 |
fun approx_vector prec pprt vector =
|
22964
|
63 |
let
|
|
64 |
fun app (index, s) (lower, upper) =
|
|
65 |
let
|
|
66 |
val (flower, fupper) = approx_value prec pprt s
|
|
67 |
val index = HOLogic.mk_number HOLogic.natT (Intt.int index)
|
|
68 |
val elower = HOLogic.mk_prod (index, flower)
|
|
69 |
val eupper = HOLogic.mk_prod (index, fupper)
|
|
70 |
in (elower :: lower, eupper :: upper) end;
|
|
71 |
in
|
|
72 |
pairself (HOLogic.mk_list spvecT) (Inttab.fold app vector ([], []))
|
|
73 |
end;
|
16784
|
74 |
|
22964
|
75 |
fun approx_matrix prec pprt vector =
|
|
76 |
let
|
|
77 |
fun app (index, v) (lower, upper) =
|
|
78 |
let
|
|
79 |
val (flower, fupper) = approx_vector prec pprt v
|
|
80 |
val index = HOLogic.mk_number HOLogic.natT (Intt.int index)
|
|
81 |
val elower = HOLogic.mk_prod (index, flower)
|
|
82 |
val eupper = HOLogic.mk_prod (index, fupper)
|
|
83 |
in (elower :: lower, eupper :: upper) end;
|
|
84 |
in
|
|
85 |
pairself (HOLogic.mk_list spmatT) (Inttab.fold app vector ([], []))
|
|
86 |
end;
|
16784
|
87 |
|
|
88 |
exception Nat_expected of int;
|
|
89 |
|
22964
|
90 |
val zero_interval = approx_value 1 I "0"
|
16784
|
91 |
|
22951
|
92 |
fun set_elem vector index str =
|
|
93 |
if index < 0 then
|
|
94 |
raise (Nat_expected index)
|
22964
|
95 |
else if (approx_value 1 I str) = zero_interval then
|
22951
|
96 |
vector
|
|
97 |
else
|
|
98 |
Inttab.update (index, str) vector
|
16784
|
99 |
|
22951
|
100 |
fun set_vector matrix index vector =
|
16784
|
101 |
if index < 0 then
|
22951
|
102 |
raise (Nat_expected index)
|
16784
|
103 |
else if Inttab.is_empty vector then
|
22951
|
104 |
matrix
|
16784
|
105 |
else
|
22951
|
106 |
Inttab.update (index, vector) matrix
|
16784
|
107 |
|
|
108 |
val empty_matrix = Inttab.empty
|
|
109 |
val empty_vector = Inttab.empty
|
|
110 |
|
|
111 |
(* dual stuff *)
|
|
112 |
|
|
113 |
structure cplex = Cplex
|
|
114 |
|
22951
|
115 |
fun transpose_matrix matrix =
|
21056
|
116 |
let
|
|
117 |
fun upd j (i, s) =
|
|
118 |
Inttab.map_default (i, Inttab.empty) (Inttab.update (j, s));
|
|
119 |
fun updm (j, v) = Inttab.fold (upd j) v;
|
|
120 |
in Inttab.fold updm matrix empty_matrix end;
|
16784
|
121 |
|
|
122 |
exception No_name of string;
|
|
123 |
|
|
124 |
exception Superfluous_constr_right_hand_sides
|
|
125 |
|
22951
|
126 |
fun cplexProg c A b =
|
16784
|
127 |
let
|
22951
|
128 |
val ytable = ref Inttab.empty
|
|
129 |
fun indexof s =
|
|
130 |
if String.size s = 0 then raise (No_name s)
|
|
131 |
else case Int.fromString (String.extract(s, 1, NONE)) of
|
|
132 |
SOME i => i | NONE => raise (No_name s)
|
16784
|
133 |
|
22951
|
134 |
fun nameof i =
|
|
135 |
let
|
|
136 |
val s = "x"^(Int.toString i)
|
|
137 |
val _ = change ytable (Inttab.update (i, s))
|
|
138 |
in
|
|
139 |
s
|
|
140 |
end
|
|
141 |
|
|
142 |
fun split_numstr s =
|
|
143 |
if String.isPrefix "-" s then (false,String.extract(s, 1, NONE))
|
|
144 |
else if String.isPrefix "+" s then (true, String.extract(s, 1, NONE))
|
|
145 |
else (true, s)
|
|
146 |
|
|
147 |
fun mk_term index s =
|
|
148 |
let
|
|
149 |
val (p, s) = split_numstr s
|
|
150 |
val prod = cplex.cplexProd (cplex.cplexNum s, cplex.cplexVar (nameof index))
|
|
151 |
in
|
|
152 |
if p then prod else cplex.cplexNeg prod
|
|
153 |
end
|
16784
|
154 |
|
22951
|
155 |
fun vec2sum vector =
|
|
156 |
cplex.cplexSum (Inttab.fold (fn (index, s) => fn list => (mk_term index s) :: list) vector [])
|
16784
|
157 |
|
22951
|
158 |
fun mk_constr index vector c =
|
|
159 |
let
|
|
160 |
val s = case Inttab.lookup c index of SOME s => s | NONE => "0"
|
|
161 |
val (p, s) = split_numstr s
|
|
162 |
val num = if p then cplex.cplexNum s else cplex.cplexNeg (cplex.cplexNum s)
|
|
163 |
in
|
|
164 |
(NONE, cplex.cplexConstr (cplex.cplexLeq, (vec2sum vector, num)))
|
|
165 |
end
|
|
166 |
|
|
167 |
fun delete index c = Inttab.delete index c handle Inttab.UNDEF _ => c
|
16784
|
168 |
|
22951
|
169 |
val (list, b) = Inttab.fold
|
|
170 |
(fn (index, v) => fn (list, c) => ((mk_constr index v c)::list, delete index c))
|
|
171 |
A ([], b)
|
|
172 |
val _ = if Inttab.is_empty b then () else raise Superfluous_constr_right_hand_sides
|
16784
|
173 |
|
22951
|
174 |
fun mk_free y = cplex.cplexBounds (cplex.cplexNeg cplex.cplexInf, cplex.cplexLeq,
|
|
175 |
cplex.cplexVar y, cplex.cplexLeq,
|
|
176 |
cplex.cplexInf)
|
16784
|
177 |
|
22951
|
178 |
val yvars = Inttab.fold (fn (i, y) => fn l => (mk_free y)::l) (!ytable) []
|
16784
|
179 |
|
22951
|
180 |
val prog = cplex.cplexProg ("original", cplex.cplexMaximize (vec2sum c), list, yvars)
|
16784
|
181 |
in
|
22951
|
182 |
(prog, indexof)
|
16784
|
183 |
end
|
|
184 |
|
|
185 |
|
22951
|
186 |
fun dual_cplexProg c A b =
|
16784
|
187 |
let
|
22951
|
188 |
fun indexof s =
|
|
189 |
if String.size s = 0 then raise (No_name s)
|
|
190 |
else case Int.fromString (String.extract(s, 1, NONE)) of
|
|
191 |
SOME i => i | NONE => raise (No_name s)
|
|
192 |
|
|
193 |
fun nameof i = "y"^(Int.toString i)
|
16784
|
194 |
|
22951
|
195 |
fun split_numstr s =
|
|
196 |
if String.isPrefix "-" s then (false,String.extract(s, 1, NONE))
|
|
197 |
else if String.isPrefix "+" s then (true, String.extract(s, 1, NONE))
|
|
198 |
else (true, s)
|
|
199 |
|
|
200 |
fun mk_term index s =
|
|
201 |
let
|
|
202 |
val (p, s) = split_numstr s
|
|
203 |
val prod = cplex.cplexProd (cplex.cplexNum s, cplex.cplexVar (nameof index))
|
|
204 |
in
|
|
205 |
if p then prod else cplex.cplexNeg prod
|
|
206 |
end
|
16784
|
207 |
|
22951
|
208 |
fun vec2sum vector =
|
|
209 |
cplex.cplexSum (Inttab.fold (fn (index, s) => fn list => (mk_term index s)::list) vector [])
|
16784
|
210 |
|
22951
|
211 |
fun mk_constr index vector c =
|
|
212 |
let
|
|
213 |
val s = case Inttab.lookup c index of SOME s => s | NONE => "0"
|
|
214 |
val (p, s) = split_numstr s
|
|
215 |
val num = if p then cplex.cplexNum s else cplex.cplexNeg (cplex.cplexNum s)
|
|
216 |
in
|
|
217 |
(NONE, cplex.cplexConstr (cplex.cplexEq, (vec2sum vector, num)))
|
|
218 |
end
|
16784
|
219 |
|
22951
|
220 |
fun delete index c = Inttab.delete index c handle Inttab.UNDEF _ => c
|
16784
|
221 |
|
22951
|
222 |
val (list, c) = Inttab.fold
|
|
223 |
(fn (index, v) => fn (list, c) => ((mk_constr index v c)::list, delete index c))
|
|
224 |
(transpose_matrix A) ([], c)
|
|
225 |
val _ = if Inttab.is_empty c then () else raise Superfluous_constr_right_hand_sides
|
|
226 |
|
|
227 |
val prog = cplex.cplexProg ("dual", cplex.cplexMinimize (vec2sum b), list, [])
|
16784
|
228 |
in
|
22951
|
229 |
(prog, indexof)
|
16784
|
230 |
end
|
|
231 |
|
22951
|
232 |
fun cut_vector size v =
|
21056
|
233 |
let
|
|
234 |
val count = ref 0;
|
|
235 |
fun app (i, s) = if (!count < size) then
|
|
236 |
(count := !count +1 ; Inttab.update (i, s))
|
|
237 |
else I
|
|
238 |
in
|
|
239 |
Inttab.fold app v empty_vector
|
|
240 |
end
|
16784
|
241 |
|
22951
|
242 |
fun cut_matrix vfilter vsize m =
|
21056
|
243 |
let
|
22951
|
244 |
fun app (i, v) =
|
21056
|
245 |
if is_none (Inttab.lookup vfilter i) then I
|
|
246 |
else case vsize
|
|
247 |
of NONE => Inttab.update (i, v)
|
|
248 |
| SOME s => Inttab.update (i, cut_vector s v)
|
|
249 |
in Inttab.fold app m empty_matrix end
|
|
250 |
|
17412
|
251 |
fun v_elem_at v i = Inttab.lookup v i
|
|
252 |
fun m_elem_at m i = Inttab.lookup m i
|
16784
|
253 |
|
22951
|
254 |
fun v_only_elem v =
|
16784
|
255 |
case Inttab.min_key v of
|
22951
|
256 |
NONE => NONE
|
16784
|
257 |
| SOME vmin => (case Inttab.max_key v of
|
22951
|
258 |
NONE => SOME vmin
|
|
259 |
| SOME vmax => if vmin = vmax then SOME vmin else NONE)
|
16784
|
260 |
|
21056
|
261 |
fun v_fold f = Inttab.fold f;
|
|
262 |
fun m_fold f = Inttab.fold f;
|
16784
|
263 |
|
|
264 |
end;
|