author  haftmann 
Tue, 13 Jul 2010 11:38:04 +0200  
changeset 37788  261c61fabc98 
parent 37764  3489daf839d5 
child 41491  a2ad5b824051 
permissions  rwrr 
37788  1 
(* Title: HOL/Matrix/FloatSparseMatrixBuilder.ML 
16784  2 
Author: Steven Obua 
3 
*) 

4 

22964  5 
signature FLOAT_SPARSE_MATIRX_BUILDER = 
16784  6 
sig 
22964  7 
include MATRIX_BUILDER 
16784  8 

22964  9 
structure cplex : CPLEX 
16784  10 

22964  11 
type float = Float.float 
12 
val approx_value : int > (float > float) > string > term * term 

13 
val approx_vector : int > (float > float) > vector > term * term 

14 
val approx_matrix : int > (float > float) > matrix > term * term 

16784  15 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24302
diff
changeset

16 
val mk_spvec_entry : int > float > term 
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24302
diff
changeset

17 
val mk_spvec_entry' : int > term > term 
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24302
diff
changeset

18 
val mk_spmat_entry : int > term > term 
22964  19 
val spvecT: typ 
20 
val spmatT: typ 

21 

22 
val v_elem_at : vector > int > string option 

23 
val m_elem_at : matrix > int > vector option 

24 
val v_only_elem : vector > int option 

25 
val v_fold : (int * string > 'a > 'a) > vector > 'a > 'a 

26 
val m_fold : (int * vector > 'a > 'a) > matrix > 'a > 'a 

20485  27 

22964  28 
val transpose_matrix : matrix > matrix 
16784  29 

22964  30 
val cut_vector : int > vector > vector 
31 
val cut_matrix : vector > int option > matrix > matrix 

22951  32 

23883  33 
val delete_matrix : int list > matrix > matrix 
34 
val cut_matrix' : int list > matrix > matrix 

35 
val delete_vector : int list > vector > vector 

36 
val cut_vector' : int list > vector > vector 

37 

38 
val indices_of_matrix : matrix > int list 

39 
val indices_of_vector : vector > int list 

40 

22964  41 
(* cplexProg c A b *) 
42 
val cplexProg : vector > matrix > vector > cplex.cplexProg * (string > int) 

43 
(* dual_cplexProg c A b *) 

44 
val dual_cplexProg : vector > matrix > vector > cplex.cplexProg * (string > int) 

45 
end; 

16784  46 

22964  47 
structure FloatSparseMatrixBuilder : FLOAT_SPARSE_MATIRX_BUILDER = 
16784  48 
struct 
49 

22964  50 
type float = Float.float 
31971
8c1b845ed105
renamed functor TableFun to Table, and GraphFun to Graph;
wenzelm
parents:
24630
diff
changeset

51 
structure Inttab = Table(type key = int val ord = rev_order o int_ord); 
16784  52 

53 
type vector = string Inttab.table 

54 
type matrix = vector Inttab.table 

20485  55 

22964  56 
val spvec_elemT = HOLogic.mk_prodT (HOLogic.natT, HOLogic.realT); 
57 
val spvecT = HOLogic.listT spvec_elemT; 

58 
val spmat_elemT = HOLogic.mk_prodT (HOLogic.natT, spvecT); 

59 
val spmatT = HOLogic.listT spmat_elemT; 

16784  60 

22964  61 
fun approx_value prec f = 
62 
FloatArith.approx_float prec (fn (x, y) => (f x, f y)); 

16784  63 

22951  64 
fun mk_spvec_entry i f = 
22964  65 
HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, FloatArith.mk_float f); 
16784  66 

24302  67 
fun mk_spvec_entry' i x = 
68 
HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, x); 

69 

22951  70 
fun mk_spmat_entry i e = 
22964  71 
HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, e); 
16784  72 

73 
fun approx_vector prec pprt vector = 

22964  74 
let 
75 
fun app (index, s) (lower, upper) = 

76 
let 

77 
val (flower, fupper) = approx_value prec pprt s 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24302
diff
changeset

78 
val index = HOLogic.mk_number HOLogic.natT index 
22964  79 
val elower = HOLogic.mk_prod (index, flower) 
80 
val eupper = HOLogic.mk_prod (index, fupper) 

81 
in (elower :: lower, eupper :: upper) end; 

82 
in 

23665
825bea0266db
adopted to new computing oracle and fixed bugs introduced by tuning
obua
parents:
23261
diff
changeset

83 
pairself (HOLogic.mk_list spvec_elemT) (Inttab.fold app vector ([], [])) 
22964  84 
end; 
16784  85 

22964  86 
fun approx_matrix prec pprt vector = 
87 
let 

88 
fun app (index, v) (lower, upper) = 

89 
let 

90 
val (flower, fupper) = approx_vector prec pprt v 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24302
diff
changeset

91 
val index = HOLogic.mk_number HOLogic.natT index 
22964  92 
val elower = HOLogic.mk_prod (index, flower) 
93 
val eupper = HOLogic.mk_prod (index, fupper) 

94 
in (elower :: lower, eupper :: upper) end; 

95 
in 

23665
825bea0266db
adopted to new computing oracle and fixed bugs introduced by tuning
obua
parents:
23261
diff
changeset

96 
pairself (HOLogic.mk_list spmat_elemT) (Inttab.fold app vector ([], [])) 
22964  97 
end; 
16784  98 

99 
exception Nat_expected of int; 

100 

22964  101 
val zero_interval = approx_value 1 I "0" 
16784  102 

22951  103 
fun set_elem vector index str = 
104 
if index < 0 then 

105 
raise (Nat_expected index) 

22964  106 
else if (approx_value 1 I str) = zero_interval then 
22951  107 
vector 
108 
else 

109 
Inttab.update (index, str) vector 

16784  110 

22951  111 
fun set_vector matrix index vector = 
16784  112 
if index < 0 then 
22951  113 
raise (Nat_expected index) 
16784  114 
else if Inttab.is_empty vector then 
22951  115 
matrix 
16784  116 
else 
22951  117 
Inttab.update (index, vector) matrix 
16784  118 

119 
val empty_matrix = Inttab.empty 

120 
val empty_vector = Inttab.empty 

121 

122 
(* dual stuff *) 

123 

124 
structure cplex = Cplex 

125 

22951  126 
fun transpose_matrix matrix = 
21056  127 
let 
128 
fun upd j (i, s) = 

129 
Inttab.map_default (i, Inttab.empty) (Inttab.update (j, s)); 

130 
fun updm (j, v) = Inttab.fold (upd j) v; 

131 
in Inttab.fold updm matrix empty_matrix end; 

16784  132 

133 
exception No_name of string; 

134 

135 
exception Superfluous_constr_right_hand_sides 

136 

22951  137 
fun cplexProg c A b = 
16784  138 
let 
32740  139 
val ytable = Unsynchronized.ref Inttab.empty 
22951  140 
fun indexof s = 
141 
if String.size s = 0 then raise (No_name s) 

142 
else case Int.fromString (String.extract(s, 1, NONE)) of 

143 
SOME i => i  NONE => raise (No_name s) 

16784  144 

22951  145 
fun nameof i = 
146 
let 

147 
val s = "x"^(Int.toString i) 

32740  148 
val _ = Unsynchronized.change ytable (Inttab.update (i, s)) 
22951  149 
in 
150 
s 

151 
end 

152 

153 
fun split_numstr s = 

154 
if String.isPrefix "" s then (false,String.extract(s, 1, NONE)) 

155 
else if String.isPrefix "+" s then (true, String.extract(s, 1, NONE)) 

156 
else (true, s) 

157 

158 
fun mk_term index s = 

159 
let 

160 
val (p, s) = split_numstr s 

161 
val prod = cplex.cplexProd (cplex.cplexNum s, cplex.cplexVar (nameof index)) 

162 
in 

163 
if p then prod else cplex.cplexNeg prod 

164 
end 

16784  165 

22951  166 
fun vec2sum vector = 
167 
cplex.cplexSum (Inttab.fold (fn (index, s) => fn list => (mk_term index s) :: list) vector []) 

16784  168 

22951  169 
fun mk_constr index vector c = 
170 
let 

171 
val s = case Inttab.lookup c index of SOME s => s  NONE => "0" 

172 
val (p, s) = split_numstr s 

173 
val num = if p then cplex.cplexNum s else cplex.cplexNeg (cplex.cplexNum s) 

174 
in 

175 
(NONE, cplex.cplexConstr (cplex.cplexLeq, (vec2sum vector, num))) 

176 
end 

177 

178 
fun delete index c = Inttab.delete index c handle Inttab.UNDEF _ => c 

16784  179 

22951  180 
val (list, b) = Inttab.fold 
181 
(fn (index, v) => fn (list, c) => ((mk_constr index v c)::list, delete index c)) 

182 
A ([], b) 

183 
val _ = if Inttab.is_empty b then () else raise Superfluous_constr_right_hand_sides 

16784  184 

22951  185 
fun mk_free y = cplex.cplexBounds (cplex.cplexNeg cplex.cplexInf, cplex.cplexLeq, 
186 
cplex.cplexVar y, cplex.cplexLeq, 

187 
cplex.cplexInf) 

16784  188 

22951  189 
val yvars = Inttab.fold (fn (i, y) => fn l => (mk_free y)::l) (!ytable) [] 
16784  190 

22951  191 
val prog = cplex.cplexProg ("original", cplex.cplexMaximize (vec2sum c), list, yvars) 
16784  192 
in 
22951  193 
(prog, indexof) 
16784  194 
end 
195 

196 

22951  197 
fun dual_cplexProg c A b = 
16784  198 
let 
22951  199 
fun indexof s = 
200 
if String.size s = 0 then raise (No_name s) 

201 
else case Int.fromString (String.extract(s, 1, NONE)) of 

202 
SOME i => i  NONE => raise (No_name s) 

203 

204 
fun nameof i = "y"^(Int.toString i) 

16784  205 

22951  206 
fun split_numstr s = 
207 
if String.isPrefix "" s then (false,String.extract(s, 1, NONE)) 

208 
else if String.isPrefix "+" s then (true, String.extract(s, 1, NONE)) 

209 
else (true, s) 

210 

211 
fun mk_term index s = 

212 
let 

213 
val (p, s) = split_numstr s 

214 
val prod = cplex.cplexProd (cplex.cplexNum s, cplex.cplexVar (nameof index)) 

215 
in 

216 
if p then prod else cplex.cplexNeg prod 

217 
end 

16784  218 

22951  219 
fun vec2sum vector = 
220 
cplex.cplexSum (Inttab.fold (fn (index, s) => fn list => (mk_term index s)::list) vector []) 

16784  221 

22951  222 
fun mk_constr index vector c = 
223 
let 

224 
val s = case Inttab.lookup c index of SOME s => s  NONE => "0" 

225 
val (p, s) = split_numstr s 

226 
val num = if p then cplex.cplexNum s else cplex.cplexNeg (cplex.cplexNum s) 

227 
in 

228 
(NONE, cplex.cplexConstr (cplex.cplexEq, (vec2sum vector, num))) 

229 
end 

16784  230 

22951  231 
fun delete index c = Inttab.delete index c handle Inttab.UNDEF _ => c 
16784  232 

22951  233 
val (list, c) = Inttab.fold 
234 
(fn (index, v) => fn (list, c) => ((mk_constr index v c)::list, delete index c)) 

235 
(transpose_matrix A) ([], c) 

236 
val _ = if Inttab.is_empty c then () else raise Superfluous_constr_right_hand_sides 

237 

238 
val prog = cplex.cplexProg ("dual", cplex.cplexMinimize (vec2sum b), list, []) 

16784  239 
in 
22951  240 
(prog, indexof) 
16784  241 
end 
242 

22951  243 
fun cut_vector size v = 
21056  244 
let 
32740  245 
val count = Unsynchronized.ref 0; 
21056  246 
fun app (i, s) = if (!count < size) then 
247 
(count := !count +1 ; Inttab.update (i, s)) 

248 
else I 

249 
in 

250 
Inttab.fold app v empty_vector 

251 
end 

16784  252 

22951  253 
fun cut_matrix vfilter vsize m = 
21056  254 
let 
22951  255 
fun app (i, v) = 
21056  256 
if is_none (Inttab.lookup vfilter i) then I 
257 
else case vsize 

258 
of NONE => Inttab.update (i, v) 

259 
 SOME s => Inttab.update (i, cut_vector s v) 

260 
in Inttab.fold app m empty_matrix end 

261 

17412  262 
fun v_elem_at v i = Inttab.lookup v i 
263 
fun m_elem_at m i = Inttab.lookup m i 

16784  264 

22951  265 
fun v_only_elem v = 
16784  266 
case Inttab.min_key v of 
22951  267 
NONE => NONE 
16784  268 
 SOME vmin => (case Inttab.max_key v of 
22951  269 
NONE => SOME vmin 
270 
 SOME vmax => if vmin = vmax then SOME vmin else NONE) 

16784  271 

21056  272 
fun v_fold f = Inttab.fold f; 
273 
fun m_fold f = Inttab.fold f; 

16784  274 

23883  275 
fun indices_of_vector v = Inttab.keys v 
276 
fun indices_of_matrix m = Inttab.keys m 

277 
fun delete_vector indices v = fold Inttab.delete indices v 

278 
fun delete_matrix indices m = fold Inttab.delete indices m 

279 
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 

280 
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 

281 

282 

283 

16784  284 
end; 