author  wenzelm 
Tue, 18 Sep 2007 16:08:00 +0200  
changeset 24630  351a308ab58d 
parent 24302  3045683749af 
child 31971  8c1b845ed105 
permissions  rwrr 
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; 