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