15179
|
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
|
15196
|
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)
|
15179
|
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;
|