16784
|
1 |
(* Title: HOL/Matrix/cplex/fspmlp.ML
|
|
2 |
ID: $Id$
|
|
3 |
Author: Steven Obua
|
|
4 |
*)
|
|
5 |
|
22951
|
6 |
signature FSPMLP =
|
16784
|
7 |
sig
|
|
8 |
type linprog
|
19404
|
9 |
type vector = FloatSparseMatrixBuilder.vector
|
|
10 |
type matrix = FloatSparseMatrixBuilder.matrix
|
16784
|
11 |
|
22964
|
12 |
val y : linprog -> term
|
|
13 |
val A : linprog -> term * term
|
|
14 |
val b : linprog -> term
|
|
15 |
val c : linprog -> term * term
|
|
16 |
val r : linprog -> term
|
|
17 |
val r12 : linprog -> term * term
|
16784
|
18 |
|
|
19 |
exception Load of string
|
22951
|
20 |
|
16784
|
21 |
val load : string -> int -> bool -> linprog
|
|
22 |
end
|
|
23 |
|
22964
|
24 |
structure Fspmlp : FSPMLP =
|
16784
|
25 |
struct
|
|
26 |
|
19404
|
27 |
type vector = FloatSparseMatrixBuilder.vector
|
|
28 |
type matrix = FloatSparseMatrixBuilder.matrix
|
|
29 |
|
22964
|
30 |
type linprog = term * (term * term) * term * (term * term) * term * (term * term)
|
16784
|
31 |
|
|
32 |
fun y (c1, c2, c3, c4, c5, _) = c1
|
|
33 |
fun A (c1, c2, c3, c4, c5, _) = c2
|
|
34 |
fun b (c1, c2, c3, c4, c5, _) = c3
|
|
35 |
fun c (c1, c2, c3, c4, c5, _) = c4
|
|
36 |
fun r (c1, c2, c3, c4, c5, _) = c5
|
|
37 |
fun r12 (c1, c2, c3, c4, c5, c6) = c6
|
|
38 |
|
22951
|
39 |
structure CplexFloatSparseMatrixConverter =
|
16784
|
40 |
MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = FloatSparseMatrixBuilder);
|
|
41 |
|
|
42 |
datatype bound_type = LOWER | UPPER
|
|
43 |
|
22951
|
44 |
fun intbound_ord ((i1, b1),(i2,b2)) =
|
16784
|
45 |
if i1 < i2 then LESS
|
22951
|
46 |
else if i1 = i2 then
|
|
47 |
(if b1 = b2 then EQUAL else if b1=LOWER then LESS else GREATER)
|
16784
|
48 |
else GREATER
|
|
49 |
|
|
50 |
structure Inttab = TableFun(type key = int val ord = (rev_order o int_ord));
|
|
51 |
|
|
52 |
structure VarGraph = TableFun(type key = int*bound_type val ord = intbound_ord);
|
|
53 |
(* key -> (float option) * (int -> (float * (((float * float) * key) list)))) *)
|
|
54 |
(* dest_key -> (sure_bound * (row_index -> (row_bound * (((coeff_lower * coeff_upper) * src_key) list)))) *)
|
|
55 |
|
|
56 |
exception Internal of string;
|
|
57 |
|
22951
|
58 |
fun add_row_bound g dest_key row_index row_bound =
|
|
59 |
let
|
|
60 |
val x =
|
|
61 |
case VarGraph.lookup g dest_key of
|
|
62 |
NONE => (NONE, Inttab.update (row_index, (row_bound, [])) Inttab.empty)
|
|
63 |
| SOME (sure_bound, f) =>
|
|
64 |
(sure_bound,
|
|
65 |
case Inttab.lookup f row_index of
|
|
66 |
NONE => Inttab.update (row_index, (row_bound, [])) f
|
|
67 |
| SOME _ => raise (Internal "add_row_bound"))
|
16784
|
68 |
in
|
22951
|
69 |
VarGraph.update (dest_key, x) g
|
16784
|
70 |
end
|
|
71 |
|
22951
|
72 |
fun update_sure_bound g (key as (_, btype)) bound =
|
|
73 |
let
|
|
74 |
val x =
|
|
75 |
case VarGraph.lookup g key of
|
|
76 |
NONE => (SOME bound, Inttab.empty)
|
|
77 |
| SOME (NONE, f) => (SOME bound, f)
|
|
78 |
| SOME (SOME old_bound, f) =>
|
|
79 |
(SOME ((case btype of
|
22964
|
80 |
UPPER => Float.min
|
|
81 |
| LOWER => Float.max)
|
22951
|
82 |
old_bound bound), f)
|
|
83 |
in
|
|
84 |
VarGraph.update (key, x) g
|
|
85 |
end
|
|
86 |
|
|
87 |
fun get_sure_bound g key =
|
|
88 |
case VarGraph.lookup g key of
|
|
89 |
NONE => NONE
|
16784
|
90 |
| SOME (sure_bound, _) => sure_bound
|
|
91 |
|
22951
|
92 |
(*fun get_row_bound g key row_index =
|
17412
|
93 |
case VarGraph.lookup g key of
|
22951
|
94 |
NONE => NONE
|
16784
|
95 |
| SOME (sure_bound, f) =>
|
22951
|
96 |
(case Inttab.lookup f row_index of
|
|
97 |
NONE => NONE
|
|
98 |
| SOME (row_bound, _) => (sure_bound, row_bound))*)
|
|
99 |
|
|
100 |
fun add_edge g src_key dest_key row_index coeff =
|
17412
|
101 |
case VarGraph.lookup g dest_key of
|
22951
|
102 |
NONE => raise (Internal "add_edge: dest_key not found")
|
16784
|
103 |
| SOME (sure_bound, f) =>
|
22951
|
104 |
(case Inttab.lookup f row_index of
|
|
105 |
NONE => raise (Internal "add_edge: row_index not found")
|
|
106 |
| SOME (row_bound, sources) =>
|
|
107 |
VarGraph.update (dest_key, (sure_bound, Inttab.update (row_index, (row_bound, (coeff, src_key) :: sources)) f)) g)
|
16784
|
108 |
|
22951
|
109 |
fun split_graph g =
|
21056
|
110 |
let
|
|
111 |
fun split (key, (sure_bound, _)) (r1, r2) = case sure_bound
|
|
112 |
of NONE => (r1, r2)
|
|
113 |
| SOME bound => (case key
|
|
114 |
of (u, UPPER) => (r1, Inttab.update (u, bound) r2)
|
|
115 |
| (u, LOWER) => (Inttab.update (u, bound) r1, r2))
|
|
116 |
in VarGraph.fold split g (Inttab.empty, Inttab.empty) end
|
16784
|
117 |
|
21056
|
118 |
fun it2list t = Inttab.fold cons t [];
|
16784
|
119 |
|
|
120 |
(* If safe is true, termination is guaranteed, but the sure bounds may be not optimal (relative to the algorithm).
|
|
121 |
If safe is false, termination is not guaranteed, but on termination the sure bounds are optimal (relative to the algorithm) *)
|
22951
|
122 |
fun propagate_sure_bounds safe names g =
|
|
123 |
let
|
|
124 |
(* returns NONE if no new sure bound could be calculated, otherwise the new sure bound is returned *)
|
|
125 |
fun calc_sure_bound_from_sources g (key as (_, btype)) =
|
|
126 |
let
|
|
127 |
fun mult_upper x (lower, upper) =
|
23520
|
128 |
if Float.sign x = LESS then
|
22964
|
129 |
Float.mult x lower
|
22951
|
130 |
else
|
22964
|
131 |
Float.mult x upper
|
16784
|
132 |
|
22951
|
133 |
fun mult_lower x (lower, upper) =
|
23520
|
134 |
if Float.sign x = LESS then
|
22964
|
135 |
Float.mult x upper
|
22951
|
136 |
else
|
22964
|
137 |
Float.mult x lower
|
22951
|
138 |
|
|
139 |
val mult_btype = case btype of UPPER => mult_upper | LOWER => mult_lower
|
16784
|
140 |
|
22951
|
141 |
fun calc_sure_bound (row_index, (row_bound, sources)) sure_bound =
|
|
142 |
let
|
|
143 |
fun add_src_bound (coeff, src_key) sum =
|
|
144 |
case sum of
|
|
145 |
NONE => NONE
|
|
146 |
| SOME x =>
|
|
147 |
(case get_sure_bound g src_key of
|
|
148 |
NONE => NONE
|
22964
|
149 |
| SOME src_sure_bound => SOME (Float.add x (mult_btype src_sure_bound coeff)))
|
22951
|
150 |
in
|
|
151 |
case fold add_src_bound sources (SOME row_bound) of
|
|
152 |
NONE => sure_bound
|
|
153 |
| new_sure_bound as (SOME new_bound) =>
|
|
154 |
(case sure_bound of
|
|
155 |
NONE => new_sure_bound
|
|
156 |
| SOME old_bound =>
|
|
157 |
SOME (case btype of
|
22964
|
158 |
UPPER => Float.min old_bound new_bound
|
|
159 |
| LOWER => Float.max old_bound new_bound))
|
22951
|
160 |
end
|
|
161 |
in
|
|
162 |
case VarGraph.lookup g key of
|
|
163 |
NONE => NONE
|
|
164 |
| SOME (sure_bound, f) =>
|
|
165 |
let
|
|
166 |
val x = Inttab.fold calc_sure_bound f sure_bound
|
|
167 |
in
|
|
168 |
if x = sure_bound then NONE else x
|
|
169 |
end
|
|
170 |
end
|
16784
|
171 |
|
22951
|
172 |
fun propagate (key, _) (g, b) =
|
|
173 |
case calc_sure_bound_from_sources g key of
|
|
174 |
NONE => (g,b)
|
|
175 |
| SOME bound => (update_sure_bound g key bound,
|
|
176 |
if safe then
|
|
177 |
case get_sure_bound g key of
|
|
178 |
NONE => true
|
|
179 |
| _ => b
|
|
180 |
else
|
|
181 |
true)
|
16784
|
182 |
|
22951
|
183 |
val (g, b) = VarGraph.fold propagate g (g, false)
|
16784
|
184 |
in
|
22951
|
185 |
if b then propagate_sure_bounds safe names g else g
|
|
186 |
end
|
|
187 |
|
16784
|
188 |
exception Load of string;
|
|
189 |
|
22964
|
190 |
(*val empty_spvec = @ {term "Nil :: (nat * real) list"};
|
|
191 |
fun cons_spvec x xs = @ {term "Cons :: nat * real => nat * real => (nat * real) list"} $ x $ xs;
|
|
192 |
val empty_spmat = @ {term "Nil :: (nat * (nat * real) list) list"};
|
|
193 |
fun cons_spmat x xs = @ {term "Cons :: nat * (nat * real) list => (nat * (nat * real) list) list => (nat * (nat * real) list) list"} $ x $ xs;*)
|
|
194 |
val empty_spvec = Bound 0
|
|
195 |
fun cons_spvec x xs = Bound 0
|
|
196 |
val empty_spmat = Bound 0
|
|
197 |
fun cons_spmat x xs = Bound 0
|
|
198 |
|
22951
|
199 |
fun calcr safe_propagation xlen names prec A b =
|
16784
|
200 |
let
|
22951
|
201 |
val empty = Inttab.empty
|
|
202 |
|
|
203 |
fun instab t i x = Inttab.update (i, x) t
|
|
204 |
|
|
205 |
fun isnegstr x = String.isPrefix "-" x
|
|
206 |
fun negstr x = if isnegstr x then String.extract (x, 1, NONE) else "-"^x
|
16784
|
207 |
|
22951
|
208 |
fun test_1 (lower, upper) =
|
|
209 |
if lower = upper then
|
23297
|
210 |
(if Float.eq (lower, (~1, 0)) then ~1
|
|
211 |
else if Float.eq (lower, (1, 0)) then 1
|
22951
|
212 |
else 0)
|
|
213 |
else 0
|
16784
|
214 |
|
22951
|
215 |
fun calcr (row_index, a) g =
|
|
216 |
let
|
|
217 |
val b = FloatSparseMatrixBuilder.v_elem_at b row_index
|
22964
|
218 |
val (_, b2) = FloatArith.approx_decstr_by_bin prec (case b of NONE => "0" | SOME b => b)
|
22951
|
219 |
val approx_a = FloatSparseMatrixBuilder.v_fold (fn (i, s) => fn l =>
|
22964
|
220 |
(i, FloatArith.approx_decstr_by_bin prec s)::l) a []
|
16784
|
221 |
|
22951
|
222 |
fun fold_dest_nodes (dest_index, dest_value) g =
|
|
223 |
let
|
|
224 |
val dest_test = test_1 dest_value
|
|
225 |
in
|
|
226 |
if dest_test = 0 then
|
|
227 |
g
|
|
228 |
else let
|
|
229 |
val (dest_key as (_, dest_btype), row_bound) =
|
|
230 |
if dest_test = ~1 then
|
22964
|
231 |
((dest_index, LOWER), Float.neg b2)
|
22951
|
232 |
else
|
|
233 |
((dest_index, UPPER), b2)
|
16784
|
234 |
|
22951
|
235 |
fun fold_src_nodes (src_index, src_value as (src_lower, src_upper)) g =
|
|
236 |
if src_index = dest_index then g
|
|
237 |
else
|
|
238 |
let
|
|
239 |
val coeff = case dest_btype of
|
22964
|
240 |
UPPER => (Float.neg src_upper, Float.neg src_lower)
|
22951
|
241 |
| LOWER => src_value
|
|
242 |
in
|
23520
|
243 |
if Float.sign src_lower = LESS then
|
22951
|
244 |
add_edge g (src_index, UPPER) dest_key row_index coeff
|
|
245 |
else
|
|
246 |
add_edge g (src_index, LOWER) dest_key row_index coeff
|
|
247 |
end
|
|
248 |
in
|
|
249 |
fold fold_src_nodes approx_a (add_row_bound g dest_key row_index row_bound)
|
|
250 |
end
|
|
251 |
end
|
|
252 |
in
|
|
253 |
case approx_a of
|
|
254 |
[] => g
|
|
255 |
| [(u, a)] =>
|
|
256 |
let
|
|
257 |
val atest = test_1 a
|
|
258 |
in
|
|
259 |
if atest = ~1 then
|
22964
|
260 |
update_sure_bound g (u, LOWER) (Float.neg b2)
|
22951
|
261 |
else if atest = 1 then
|
|
262 |
update_sure_bound g (u, UPPER) b2
|
|
263 |
else
|
|
264 |
g
|
|
265 |
end
|
|
266 |
| _ => fold fold_dest_nodes approx_a g
|
|
267 |
end
|
16784
|
268 |
|
22951
|
269 |
val g = FloatSparseMatrixBuilder.m_fold calcr A VarGraph.empty
|
|
270 |
val g = propagate_sure_bounds safe_propagation names g
|
|
271 |
|
|
272 |
val (r1, r2) = split_graph g
|
|
273 |
|
|
274 |
fun add_row_entry m index value =
|
|
275 |
let
|
22964
|
276 |
val vec = cons_spvec (FloatSparseMatrixBuilder.mk_spvec_entry 0 value) empty_spvec
|
22951
|
277 |
in
|
22964
|
278 |
cons_spmat (FloatSparseMatrixBuilder.mk_spmat_entry index vec) m
|
22951
|
279 |
end
|
16784
|
280 |
|
22951
|
281 |
fun abs_estimate i r1 r2 =
|
|
282 |
if i = 0 then
|
22964
|
283 |
let val e = empty_spmat in (e, (e, e)) end
|
22951
|
284 |
else
|
|
285 |
let
|
|
286 |
val index = xlen-i
|
|
287 |
val (r, (r12_1, r12_2)) = abs_estimate (i-1) r1 r2
|
|
288 |
val b1 = case Inttab.lookup r1 index of NONE => raise (Load ("x-value not bounded from below: "^(names index))) | SOME x => x
|
|
289 |
val b2 = case Inttab.lookup r2 index of NONE => raise (Load ("x-value not bounded from above: "^(names index))) | SOME x => x
|
22964
|
290 |
val abs_max = Float.max (Float.neg (Float.negative_part b1)) (Float.positive_part b2)
|
23261
|
291 |
val i' = Integer.int index
|
22951
|
292 |
in
|
22964
|
293 |
(add_row_entry r i' abs_max, (add_row_entry r12_1 i' b1, add_row_entry r12_2 i' b2))
|
22951
|
294 |
end
|
|
295 |
val (r, (r1, r2)) = abs_estimate xlen r1 r2
|
16784
|
296 |
in
|
22964
|
297 |
(r, (r1, r2))
|
16784
|
298 |
end
|
22951
|
299 |
|
16784
|
300 |
fun load filename prec safe_propagation =
|
|
301 |
let
|
22951
|
302 |
val prog = Cplex.load_cplexFile filename
|
|
303 |
val prog = Cplex.elim_nonfree_bounds prog
|
|
304 |
val prog = Cplex.relax_strict_ineqs prog
|
|
305 |
val (maximize, c, A, b, (xlen, names, _)) = CplexFloatSparseMatrixConverter.convert_prog prog
|
|
306 |
val (r, (r1, r2)) = calcr safe_propagation xlen names prec A b
|
|
307 |
val _ = if maximize then () else raise Load "sorry, cannot handle minimization problems"
|
|
308 |
val (dualprog, indexof) = FloatSparseMatrixBuilder.dual_cplexProg c A b
|
|
309 |
val results = Cplex.solve dualprog
|
|
310 |
val (optimal,v) = CplexFloatSparseMatrixConverter.convert_results results indexof
|
|
311 |
val A = FloatSparseMatrixBuilder.cut_matrix v NONE A
|
|
312 |
fun id x = x
|
|
313 |
val v = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 v
|
|
314 |
val b = FloatSparseMatrixBuilder.transpose_matrix (FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 b)
|
|
315 |
val c = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 c
|
22964
|
316 |
val (y1, _) = FloatSparseMatrixBuilder.approx_matrix prec Float.positive_part v
|
22951
|
317 |
val A = FloatSparseMatrixBuilder.approx_matrix prec id A
|
|
318 |
val (_,b2) = FloatSparseMatrixBuilder.approx_matrix prec id b
|
|
319 |
val c = FloatSparseMatrixBuilder.approx_matrix prec id c
|
16784
|
320 |
in
|
22951
|
321 |
(y1, A, b2, c, r, (r1, r2))
|
|
322 |
end handle CplexFloatSparseMatrixConverter.Converter s => (raise (Load ("Converter: "^s)))
|
16784
|
323 |
|
|
324 |
end
|