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