1 structure FloatSparseMatrixBuilder : |
|
2 sig |
|
3 include MATRIX_BUILDER |
|
4 |
|
5 structure cplex : CPLEX |
|
6 |
|
7 type float = IntInf.int*IntInf.int |
|
8 type floatfunc = float -> float |
|
9 |
|
10 |
|
11 val float2cterm : IntInf.int * IntInf.int -> cterm |
|
12 |
|
13 val approx_value : int -> floatfunc -> string -> cterm * cterm |
|
14 val approx_vector : int -> floatfunc -> vector -> cterm * cterm |
|
15 val approx_matrix : int -> floatfunc -> matrix -> cterm * cterm |
|
16 |
|
17 val mk_spvec_entry : int -> float -> term |
|
18 val empty_spvec : term |
|
19 val cons_spvec : term -> term -> term |
|
20 val empty_spmat : term |
|
21 val mk_spmat_entry : int -> term -> term |
|
22 val cons_spmat : term -> term -> term |
|
23 val sign_term : term -> cterm |
|
24 |
|
25 val v_elem_at : vector -> int -> string option |
|
26 val m_elem_at : matrix -> int -> vector option |
|
27 val v_only_elem : vector -> int option |
|
28 val v_fold : ('a * (int * string) -> 'a) -> 'a -> vector -> 'a |
|
29 val m_fold : ('a * (int * vector) -> 'a) -> 'a -> matrix -> 'a |
|
30 |
|
31 val transpose_matrix : matrix -> matrix |
|
32 |
|
33 val cut_vector : int -> vector -> vector |
|
34 val cut_matrix : vector -> (int option) -> matrix -> matrix |
|
35 |
|
36 (* cplexProg c A b *) |
|
37 val cplexProg : vector -> matrix -> vector -> (cplex.cplexProg * (string -> int)) |
|
38 (* dual_cplexProg c A b *) |
|
39 val dual_cplexProg : vector -> matrix -> vector -> (cplex.cplexProg * (string -> int)) |
|
40 |
|
41 val real_spmatT : typ |
|
42 val real_spvecT : typ |
|
43 end |
|
44 = |
|
45 struct |
|
46 |
|
47 |
|
48 structure Inttab = TableFun(type key = int val ord = (rev_order o int_ord)); |
|
49 |
|
50 type vector = string Inttab.table |
|
51 type matrix = vector Inttab.table |
|
52 type float = IntInf.int*IntInf.int |
|
53 type floatfunc = float -> float |
|
54 |
|
55 val th = theory "Float" |
|
56 val sg = sign_of th |
|
57 |
|
58 fun readtype s = Sign.intern_tycon sg s |
|
59 fun readterm s = Sign.intern_const sg s |
|
60 |
|
61 val ty_list = readtype "list" |
|
62 val term_Nil = readterm "Nil" |
|
63 val term_Cons = readterm "Cons" |
|
64 |
|
65 val spvec_elemT = HOLogic.mk_prodT (HOLogic.natT, HOLogic.realT) |
|
66 val spvecT = Type (ty_list, [spvec_elemT]) |
|
67 val spmat_elemT = HOLogic.mk_prodT (HOLogic.natT, spvecT) |
|
68 val spmatT = Type (ty_list, [spmat_elemT]) |
|
69 |
|
70 val real_spmatT = spmatT |
|
71 val real_spvecT = spvecT |
|
72 |
|
73 val empty_matrix_const = Const (term_Nil, spmatT) |
|
74 val empty_vector_const = Const (term_Nil, spvecT) |
|
75 |
|
76 val Cons_spvec_const = Const (term_Cons, spvec_elemT --> spvecT --> spvecT) |
|
77 val Cons_spmat_const = Const (term_Cons, spmat_elemT --> spmatT --> spmatT) |
|
78 |
|
79 val float_const = Const (readterm "float", HOLogic.mk_prodT (HOLogic.intT, HOLogic.intT) --> HOLogic.realT) |
|
80 |
|
81 val zero = IntInf.fromInt 0 |
|
82 val minus_one = IntInf.fromInt ~1 |
|
83 val two = IntInf.fromInt 2 |
|
84 |
|
85 fun mk_intinf ty n = |
|
86 let |
|
87 fun mk_bit n = if n = zero then HOLogic.false_const else HOLogic.true_const |
|
88 |
|
89 fun bin_of n = |
|
90 if n = zero then HOLogic.pls_const |
|
91 else if n = minus_one then HOLogic.min_const |
|
92 else |
|
93 let |
|
94 (*val (q,r) = IntInf.divMod (n, two): doesn't work in SML 10.0.7, but in newer versions!!!*) |
|
95 val q = IntInf.div (n, two) |
|
96 val r = IntInf.mod (n, two) |
|
97 in |
|
98 HOLogic.bit_const $ bin_of q $ mk_bit r |
|
99 end |
|
100 in |
|
101 HOLogic.number_of_const ty $ (bin_of n) |
|
102 end |
|
103 |
|
104 fun mk_float (a,b) = |
|
105 float_const $ (HOLogic.mk_prod ((mk_intinf HOLogic.intT a), (mk_intinf HOLogic.intT b))) |
|
106 |
|
107 fun float2cterm (a,b) = cterm_of sg (mk_float (a,b)) |
|
108 |
|
109 fun approx_value_term prec f value = |
|
110 let |
|
111 val (flower, fupper) = ExactFloatingPoint.approx_decstr_by_bin prec value |
|
112 val (flower, fupper) = (f flower, f fupper) |
|
113 in |
|
114 (mk_float flower, mk_float fupper) |
|
115 end |
|
116 |
|
117 fun approx_value prec pprt value = |
|
118 let |
|
119 val (flower, fupper) = approx_value_term prec pprt value |
|
120 in |
|
121 (cterm_of sg flower, cterm_of sg fupper) |
|
122 end |
|
123 |
|
124 fun sign_term t = cterm_of sg t |
|
125 |
|
126 val empty_spvec = empty_vector_const |
|
127 |
|
128 val empty_spmat = empty_matrix_const |
|
129 |
|
130 fun mk_spvec_entry i f = |
|
131 let |
|
132 val term_i = mk_intinf HOLogic.natT (IntInf.fromInt i) |
|
133 val term_f = mk_float f |
|
134 in |
|
135 HOLogic.mk_prod (term_i, term_f) |
|
136 end |
|
137 |
|
138 fun mk_spmat_entry i e = |
|
139 let |
|
140 val term_i = mk_intinf HOLogic.natT (IntInf.fromInt i) |
|
141 in |
|
142 HOLogic.mk_prod (term_i, e) |
|
143 end |
|
144 |
|
145 fun cons_spvec h t = Cons_spvec_const $ h $ t |
|
146 |
|
147 fun cons_spmat h t = Cons_spmat_const $ h $ t |
|
148 |
|
149 fun approx_vector_term prec pprt vector = |
|
150 let |
|
151 fun app ((vlower, vupper), (index, s)) = |
|
152 let |
|
153 val (flower, fupper) = approx_value_term prec pprt s |
|
154 val index = mk_intinf HOLogic.natT (IntInf.fromInt index) |
|
155 val elower = HOLogic.mk_prod (index, flower) |
|
156 val eupper = HOLogic.mk_prod (index, fupper) |
|
157 in |
|
158 (Cons_spvec_const $ elower $ vlower, |
|
159 Cons_spvec_const $ eupper $ vupper) |
|
160 end |
|
161 in |
|
162 Inttab.foldl app ((empty_vector_const, empty_vector_const), vector) |
|
163 end |
|
164 |
|
165 fun approx_matrix_term prec pprt matrix = |
|
166 let |
|
167 fun app ((mlower, mupper), (index, vector)) = |
|
168 let |
|
169 val (vlower, vupper) = approx_vector_term prec pprt vector |
|
170 val index = mk_intinf HOLogic.natT (IntInf.fromInt index) |
|
171 val elower = HOLogic.mk_prod (index, vlower) |
|
172 val eupper = HOLogic.mk_prod (index, vupper) |
|
173 in |
|
174 (Cons_spmat_const $ elower $ mlower, |
|
175 Cons_spmat_const $ eupper $ mupper) |
|
176 end |
|
177 |
|
178 val (mlower, mupper) = Inttab.foldl app ((empty_matrix_const, empty_matrix_const), matrix) |
|
179 in |
|
180 Inttab.foldl app ((empty_matrix_const, empty_matrix_const), matrix) |
|
181 end |
|
182 |
|
183 fun approx_vector prec pprt vector = |
|
184 let |
|
185 val (l, u) = approx_vector_term prec pprt vector |
|
186 in |
|
187 (cterm_of sg l, cterm_of sg u) |
|
188 end |
|
189 |
|
190 fun approx_matrix prec pprt matrix = |
|
191 let |
|
192 val (l, u) = approx_matrix_term prec pprt matrix |
|
193 in |
|
194 (cterm_of sg l, cterm_of sg u) |
|
195 end |
|
196 |
|
197 |
|
198 exception Nat_expected of int; |
|
199 |
|
200 val zero_interval = approx_value_term 1 I "0" |
|
201 |
|
202 fun set_elem vector index str = |
|
203 if index < 0 then |
|
204 raise (Nat_expected index) |
|
205 else if (approx_value_term 1 I str) = zero_interval then |
|
206 vector |
|
207 else |
|
208 Inttab.update ((index, str), vector) |
|
209 |
|
210 fun set_vector matrix index vector = |
|
211 if index < 0 then |
|
212 raise (Nat_expected index) |
|
213 else if Inttab.is_empty vector then |
|
214 matrix |
|
215 else |
|
216 Inttab.update ((index, vector), matrix) |
|
217 |
|
218 val empty_matrix = Inttab.empty |
|
219 val empty_vector = Inttab.empty |
|
220 |
|
221 (* dual stuff *) |
|
222 |
|
223 structure cplex = Cplex |
|
224 |
|
225 fun transpose_matrix matrix = |
|
226 let |
|
227 fun upd m j i x = |
|
228 case Inttab.lookup (m, j) of |
|
229 SOME v => Inttab.update ((j, Inttab.update ((i, x), v)), m) |
|
230 | NONE => Inttab.update ((j, Inttab.update ((i, x), Inttab.empty)), m) |
|
231 |
|
232 fun updv j (m, (i, s)) = upd m i j s |
|
233 |
|
234 fun updm (m, (j, v)) = Inttab.foldl (updv j) (m, v) |
|
235 in |
|
236 Inttab.foldl updm (empty_matrix, matrix) |
|
237 end |
|
238 |
|
239 exception No_name of string; |
|
240 |
|
241 exception Superfluous_constr_right_hand_sides |
|
242 |
|
243 fun cplexProg c A b = |
|
244 let |
|
245 val ytable = ref Inttab.empty |
|
246 fun indexof s = |
|
247 if String.size s = 0 then raise (No_name s) |
|
248 else case Int.fromString (String.extract(s, 1, NONE)) of |
|
249 SOME i => i | NONE => raise (No_name s) |
|
250 |
|
251 fun nameof i = |
|
252 let |
|
253 val s = "x"^(Int.toString i) |
|
254 val _ = ytable := (Inttab.update ((i, s), !ytable)) |
|
255 in |
|
256 s |
|
257 end |
|
258 |
|
259 fun split_numstr s = |
|
260 if String.isPrefix "-" s then (false,String.extract(s, 1, NONE)) |
|
261 else if String.isPrefix "+" s then (true, String.extract(s, 1, NONE)) |
|
262 else (true, s) |
|
263 |
|
264 fun mk_term index s = |
|
265 let |
|
266 val (p, s) = split_numstr s |
|
267 val prod = cplex.cplexProd (cplex.cplexNum s, cplex.cplexVar (nameof index)) |
|
268 in |
|
269 if p then prod else cplex.cplexNeg prod |
|
270 end |
|
271 |
|
272 fun vec2sum vector = |
|
273 cplex.cplexSum (Inttab.foldl (fn (list, (index, s)) => (mk_term index s)::list) ([], vector)) |
|
274 |
|
275 fun mk_constr index vector c = |
|
276 let |
|
277 val s = case Inttab.lookup (c, index) of SOME s => s | NONE => "0" |
|
278 val (p, s) = split_numstr s |
|
279 val num = if p then cplex.cplexNum s else cplex.cplexNeg (cplex.cplexNum s) |
|
280 in |
|
281 (NONE, cplex.cplexConstr (cplex.cplexLeq, (vec2sum vector, num))) |
|
282 end |
|
283 |
|
284 fun delete index c = Inttab.delete index c handle Inttab.UNDEF _ => c |
|
285 |
|
286 val (list, b) = Inttab.foldl |
|
287 (fn ((list, c), (index, v)) => ((mk_constr index v c)::list, delete index c)) |
|
288 (([], b), A) |
|
289 val _ = if Inttab.is_empty b then () else raise Superfluous_constr_right_hand_sides |
|
290 |
|
291 fun mk_free y = cplex.cplexBounds (cplex.cplexNeg cplex.cplexInf, cplex.cplexLeq, |
|
292 cplex.cplexVar y, cplex.cplexLeq, |
|
293 cplex.cplexInf) |
|
294 |
|
295 val yvars = Inttab.foldl (fn (l, (i, y)) => (mk_free y)::l) ([], !ytable) |
|
296 |
|
297 val prog = cplex.cplexProg ("original", cplex.cplexMaximize (vec2sum c), list, yvars) |
|
298 in |
|
299 (prog, indexof) |
|
300 end |
|
301 |
|
302 |
|
303 fun dual_cplexProg c A b = |
|
304 let |
|
305 fun indexof s = |
|
306 if String.size s = 0 then raise (No_name s) |
|
307 else case Int.fromString (String.extract(s, 1, NONE)) of |
|
308 SOME i => i | NONE => raise (No_name s) |
|
309 |
|
310 fun nameof i = "y"^(Int.toString i) |
|
311 |
|
312 fun split_numstr s = |
|
313 if String.isPrefix "-" s then (false,String.extract(s, 1, NONE)) |
|
314 else if String.isPrefix "+" s then (true, String.extract(s, 1, NONE)) |
|
315 else (true, s) |
|
316 |
|
317 fun mk_term index s = |
|
318 let |
|
319 val (p, s) = split_numstr s |
|
320 val prod = cplex.cplexProd (cplex.cplexNum s, cplex.cplexVar (nameof index)) |
|
321 in |
|
322 if p then prod else cplex.cplexNeg prod |
|
323 end |
|
324 |
|
325 fun vec2sum vector = |
|
326 cplex.cplexSum (Inttab.foldl (fn (list, (index, s)) => (mk_term index s)::list) ([], vector)) |
|
327 |
|
328 fun mk_constr index vector c = |
|
329 let |
|
330 val s = case Inttab.lookup (c, index) of SOME s => s | NONE => "0" |
|
331 val (p, s) = split_numstr s |
|
332 val num = if p then cplex.cplexNum s else cplex.cplexNeg (cplex.cplexNum s) |
|
333 in |
|
334 (NONE, cplex.cplexConstr (cplex.cplexEq, (vec2sum vector, num))) |
|
335 end |
|
336 |
|
337 fun delete index c = Inttab.delete index c handle Inttab.UNDEF _ => c |
|
338 |
|
339 val (list, c) = Inttab.foldl |
|
340 (fn ((list, c), (index, v)) => ((mk_constr index v c)::list, delete index c)) |
|
341 (([], c), transpose_matrix A) |
|
342 val _ = if Inttab.is_empty c then () else raise Superfluous_constr_right_hand_sides |
|
343 |
|
344 val prog = cplex.cplexProg ("dual", cplex.cplexMinimize (vec2sum b), list, []) |
|
345 in |
|
346 (prog, indexof) |
|
347 end |
|
348 |
|
349 fun cut_vector size v = |
|
350 let |
|
351 val count = ref 0 |
|
352 fun app (v, (i, s)) = |
|
353 if (!count < size) then |
|
354 (count := !count +1 ; Inttab.update ((i,s),v)) |
|
355 else |
|
356 v |
|
357 in |
|
358 Inttab.foldl app (empty_vector, v) |
|
359 end |
|
360 |
|
361 fun cut_matrix vfilter vsize m = |
|
362 let |
|
363 fun app (m, (i, v)) = |
|
364 if (Inttab.lookup (vfilter, i) = NONE) then |
|
365 m |
|
366 else |
|
367 case vsize of |
|
368 NONE => Inttab.update ((i,v), m) |
|
369 | SOME s => Inttab.update((i, cut_vector s v),m) |
|
370 in |
|
371 Inttab.foldl app (empty_matrix, m) |
|
372 end |
|
373 |
|
374 fun v_elem_at v i = Inttab.lookup (v,i) |
|
375 fun m_elem_at m i = Inttab.lookup (m,i) |
|
376 |
|
377 fun v_only_elem v = |
|
378 case Inttab.min_key v of |
|
379 NONE => NONE |
|
380 | SOME vmin => (case Inttab.max_key v of |
|
381 NONE => SOME vmin |
|
382 | SOME vmax => if vmin = vmax then SOME vmin else NONE) |
|
383 |
|
384 fun v_fold f a v = Inttab.foldl f (a,v) |
|
385 |
|
386 fun m_fold f a m = Inttab.foldl f (a,m) |
|
387 |
|
388 end; |
|