author  wenzelm 
Sat, 07 Apr 2012 16:41:59 +0200  
changeset 47389  e8552cba702d 
parent 46988  9f492f5b0cec 
child 47455  26315a545e26 
permissions  rwrr 
37788  1 
(* Title: HOL/Matrix/fspmlp.ML 
16784  2 
Author: Steven Obua 
3 
*) 

4 

22951  5 
signature FSPMLP = 
16784  6 
sig 
7 
type linprog 

19404  8 
type vector = FloatSparseMatrixBuilder.vector 
9 
type matrix = FloatSparseMatrixBuilder.matrix 

16784  10 

22964  11 
val y : linprog > term 
12 
val A : linprog > term * term 

13 
val b : linprog > term 

14 
val c : linprog > term * term 

15 
val r12 : linprog > term * term 

16784  16 

17 
exception Load of string 

22951  18 

16784  19 
val load : string > int > bool > linprog 
20 
end 

21 

22964  22 
structure Fspmlp : FSPMLP = 
16784  23 
struct 
24 

19404  25 
type vector = FloatSparseMatrixBuilder.vector 
26 
type matrix = FloatSparseMatrixBuilder.matrix 

27 

24302  28 
type linprog = term * (term * term) * term * (term * term) * (term * term) 
16784  29 

46533  30 
fun y (c1, _, _, _, _) = c1 
31 
fun A (_, c2, _, _, _) = c2 

32 
fun b (_, _, c3, _, _) = c3 

33 
fun c (_, _, _, c4, _) = c4 

34 
fun r12 (_, _, _, _, c6) = c6 

16784  35 

22951  36 
structure CplexFloatSparseMatrixConverter = 
16784  37 
MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = FloatSparseMatrixBuilder); 
38 

39 
datatype bound_type = LOWER  UPPER 

40 

24135
e7fb7ef2a85e
added int type constraints to accomodate hacked SML/NJ;
wenzelm
parents:
23665
diff
changeset

41 
fun intbound_ord ((i1: int, b1),(i2,b2)) = 
16784  42 
if i1 < i2 then LESS 
22951  43 
else if i1 = i2 then 
44 
(if b1 = b2 then EQUAL else if b1=LOWER then LESS else GREATER) 

16784  45 
else GREATER 
46 

31971
8c1b845ed105
renamed functor TableFun to Table, and GraphFun to Graph;
wenzelm
parents:
24630
diff
changeset

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

31971
8c1b845ed105
renamed functor TableFun to Table, and GraphFun to Graph;
wenzelm
parents:
24630
diff
changeset

49 
structure VarGraph = Table(type key = int*bound_type val ord = intbound_ord); 
16784  50 
(* key > (float option) * (int > (float * (((float * float) * key) list)))) *) 
51 
(* dest_key > (sure_bound * (row_index > (row_bound * (((coeff_lower * coeff_upper) * src_key) list)))) *) 

52 

53 
exception Internal of string; 

54 

22951  55 
fun add_row_bound g dest_key row_index row_bound = 
56 
let 

57 
val x = 

58 
case VarGraph.lookup g dest_key of 

59 
NONE => (NONE, Inttab.update (row_index, (row_bound, [])) Inttab.empty) 

60 
 SOME (sure_bound, f) => 

61 
(sure_bound, 

62 
case Inttab.lookup f row_index of 

63 
NONE => Inttab.update (row_index, (row_bound, [])) f 

64 
 SOME _ => raise (Internal "add_row_bound")) 

16784  65 
in 
22951  66 
VarGraph.update (dest_key, x) g 
16784  67 
end 
68 

22951  69 
fun update_sure_bound g (key as (_, btype)) bound = 
70 
let 

71 
val x = 

72 
case VarGraph.lookup g key of 

73 
NONE => (SOME bound, Inttab.empty) 

74 
 SOME (NONE, f) => (SOME bound, f) 

75 
 SOME (SOME old_bound, f) => 

76 
(SOME ((case btype of 

22964  77 
UPPER => Float.min 
78 
 LOWER => Float.max) 

22951  79 
old_bound bound), f) 
80 
in 

81 
VarGraph.update (key, x) g 

82 
end 

83 

84 
fun get_sure_bound g key = 

85 
case VarGraph.lookup g key of 

86 
NONE => NONE 

16784  87 
 SOME (sure_bound, _) => sure_bound 
88 

22951  89 
(*fun get_row_bound g key row_index = 
17412  90 
case VarGraph.lookup g key of 
22951  91 
NONE => NONE 
16784  92 
 SOME (sure_bound, f) => 
22951  93 
(case Inttab.lookup f row_index of 
94 
NONE => NONE 

95 
 SOME (row_bound, _) => (sure_bound, row_bound))*) 

96 

97 
fun add_edge g src_key dest_key row_index coeff = 

17412  98 
case VarGraph.lookup g dest_key of 
22951  99 
NONE => raise (Internal "add_edge: dest_key not found") 
16784  100 
 SOME (sure_bound, f) => 
22951  101 
(case Inttab.lookup f row_index of 
102 
NONE => raise (Internal "add_edge: row_index not found") 

103 
 SOME (row_bound, sources) => 

104 
VarGraph.update (dest_key, (sure_bound, Inttab.update (row_index, (row_bound, (coeff, src_key) :: sources)) f)) g) 

16784  105 

22951  106 
fun split_graph g = 
21056  107 
let 
108 
fun split (key, (sure_bound, _)) (r1, r2) = case sure_bound 

109 
of NONE => (r1, r2) 

110 
 SOME bound => (case key 

111 
of (u, UPPER) => (r1, Inttab.update (u, bound) r2) 

112 
 (u, LOWER) => (Inttab.update (u, bound) r1, r2)) 

113 
in VarGraph.fold split g (Inttab.empty, Inttab.empty) end 

16784  114 

115 
(* If safe is true, termination is guaranteed, but the sure bounds may be not optimal (relative to the algorithm). 

116 
If safe is false, termination is not guaranteed, but on termination the sure bounds are optimal (relative to the algorithm) *) 

22951  117 
fun propagate_sure_bounds safe names g = 
118 
let 

119 
(* returns NONE if no new sure bound could be calculated, otherwise the new sure bound is returned *) 

120 
fun calc_sure_bound_from_sources g (key as (_, btype)) = 

121 
let 

122 
fun mult_upper x (lower, upper) = 

23520  123 
if Float.sign x = LESS then 
22964  124 
Float.mult x lower 
22951  125 
else 
22964  126 
Float.mult x upper 
16784  127 

22951  128 
fun mult_lower x (lower, upper) = 
23520  129 
if Float.sign x = LESS then 
22964  130 
Float.mult x upper 
22951  131 
else 
22964  132 
Float.mult x lower 
22951  133 

134 
val mult_btype = case btype of UPPER => mult_upper  LOWER => mult_lower 

16784  135 

46533  136 
fun calc_sure_bound (_, (row_bound, sources)) sure_bound = 
22951  137 
let 
138 
fun add_src_bound (coeff, src_key) sum = 

139 
case sum of 

140 
NONE => NONE 

141 
 SOME x => 

142 
(case get_sure_bound g src_key of 

143 
NONE => NONE 

22964  144 
 SOME src_sure_bound => SOME (Float.add x (mult_btype src_sure_bound coeff))) 
22951  145 
in 
146 
case fold add_src_bound sources (SOME row_bound) of 

147 
NONE => sure_bound 

148 
 new_sure_bound as (SOME new_bound) => 

149 
(case sure_bound of 

150 
NONE => new_sure_bound 

151 
 SOME old_bound => 

152 
SOME (case btype of 

22964  153 
UPPER => Float.min old_bound new_bound 
154 
 LOWER => Float.max old_bound new_bound)) 

22951  155 
end 
156 
in 

157 
case VarGraph.lookup g key of 

158 
NONE => NONE 

159 
 SOME (sure_bound, f) => 

160 
let 

161 
val x = Inttab.fold calc_sure_bound f sure_bound 

162 
in 

163 
if x = sure_bound then NONE else x 

164 
end 

165 
end 

16784  166 

22951  167 
fun propagate (key, _) (g, b) = 
168 
case calc_sure_bound_from_sources g key of 

169 
NONE => (g,b) 

170 
 SOME bound => (update_sure_bound g key bound, 

171 
if safe then 

172 
case get_sure_bound g key of 

173 
NONE => true 

174 
 _ => b 

175 
else 

176 
true) 

16784  177 

22951  178 
val (g, b) = VarGraph.fold propagate g (g, false) 
16784  179 
in 
22951  180 
if b then propagate_sure_bounds safe names g else g 
181 
end 

182 

16784  183 
exception Load of string; 
184 

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

185 
val empty_spvec = @{term "Nil :: real spvec"}; 
825bea0266db
adopted to new computing oracle and fixed bugs introduced by tuning
obua
parents:
23520
diff
changeset

186 
fun cons_spvec x xs = @{term "Cons :: nat * real => real spvec => real spvec"} $ x $ xs; 
825bea0266db
adopted to new computing oracle and fixed bugs introduced by tuning
obua
parents:
23520
diff
changeset

187 
val empty_spmat = @{term "Nil :: real spmat"}; 
825bea0266db
adopted to new computing oracle and fixed bugs introduced by tuning
obua
parents:
23520
diff
changeset

188 
fun cons_spmat x xs = @{term "Cons :: nat * real spvec => real spmat => real spmat"} $ x $ xs; 
22964  189 

22951  190 
fun calcr safe_propagation xlen names prec A b = 
16784  191 
let 
22951  192 
fun test_1 (lower, upper) = 
193 
if lower = upper then 

23297  194 
(if Float.eq (lower, (~1, 0)) then ~1 
195 
else if Float.eq (lower, (1, 0)) then 1 

22951  196 
else 0) 
197 
else 0 

16784  198 

22951  199 
fun calcr (row_index, a) g = 
200 
let 

201 
val b = FloatSparseMatrixBuilder.v_elem_at b row_index 

22964  202 
val (_, b2) = FloatArith.approx_decstr_by_bin prec (case b of NONE => "0"  SOME b => b) 
22951  203 
val approx_a = FloatSparseMatrixBuilder.v_fold (fn (i, s) => fn l => 
22964  204 
(i, FloatArith.approx_decstr_by_bin prec s)::l) a [] 
16784  205 

22951  206 
fun fold_dest_nodes (dest_index, dest_value) g = 
207 
let 

208 
val dest_test = test_1 dest_value 

209 
in 

210 
if dest_test = 0 then 

211 
g 

212 
else let 

213 
val (dest_key as (_, dest_btype), row_bound) = 

214 
if dest_test = ~1 then 

22964  215 
((dest_index, LOWER), Float.neg b2) 
22951  216 
else 
217 
((dest_index, UPPER), b2) 

16784  218 

22951  219 
fun fold_src_nodes (src_index, src_value as (src_lower, src_upper)) g = 
220 
if src_index = dest_index then g 

221 
else 

222 
let 

223 
val coeff = case dest_btype of 

22964  224 
UPPER => (Float.neg src_upper, Float.neg src_lower) 
22951  225 
 LOWER => src_value 
226 
in 

23520  227 
if Float.sign src_lower = LESS then 
22951  228 
add_edge g (src_index, UPPER) dest_key row_index coeff 
229 
else 

230 
add_edge g (src_index, LOWER) dest_key row_index coeff 

231 
end 

232 
in 

233 
fold fold_src_nodes approx_a (add_row_bound g dest_key row_index row_bound) 

234 
end 

235 
end 

236 
in 

237 
case approx_a of 

238 
[] => g 

239 
 [(u, a)] => 

240 
let 

241 
val atest = test_1 a 

242 
in 

243 
if atest = ~1 then 

22964  244 
update_sure_bound g (u, LOWER) (Float.neg b2) 
22951  245 
else if atest = 1 then 
246 
update_sure_bound g (u, UPPER) b2 

247 
else 

248 
g 

249 
end 

250 
 _ => fold fold_dest_nodes approx_a g 

251 
end 

16784  252 

22951  253 
val g = FloatSparseMatrixBuilder.m_fold calcr A VarGraph.empty 
23665
825bea0266db
adopted to new computing oracle and fixed bugs introduced by tuning
obua
parents:
23520
diff
changeset

254 

22951  255 
val g = propagate_sure_bounds safe_propagation names g 
256 

32960
69916a850301
eliminated hard tabulators, guessing at each author's individual tabwidth;
wenzelm
parents:
31971
diff
changeset

257 
val (r1, r2) = split_graph g 
22951  258 

24302  259 
fun add_row_entry m index f vname value = 
22951  260 
let 
32960
69916a850301
eliminated hard tabulators, guessing at each author's individual tabwidth;
wenzelm
parents:
31971
diff
changeset

261 
val v = (case value of 
69916a850301
eliminated hard tabulators, guessing at each author's individual tabwidth;
wenzelm
parents:
31971
diff
changeset

262 
SOME value => FloatSparseMatrixBuilder.mk_spvec_entry 0 value 
69916a850301
eliminated hard tabulators, guessing at each author's individual tabwidth;
wenzelm
parents:
31971
diff
changeset

263 
 NONE => FloatSparseMatrixBuilder.mk_spvec_entry' 0 (f $ (Var ((vname,0), HOLogic.realT)))) 
24302  264 
val vec = cons_spvec v empty_spvec 
22951  265 
in 
22964  266 
cons_spmat (FloatSparseMatrixBuilder.mk_spmat_entry index vec) m 
22951  267 
end 
16784  268 

22951  269 
fun abs_estimate i r1 r2 = 
270 
if i = 0 then 

24302  271 
let val e = empty_spmat in (e, e) end 
22951  272 
else 
273 
let 

274 
val index = xleni 

24302  275 
val (r12_1, r12_2) = abs_estimate (i1) r1 r2 
276 
val b1 = Inttab.lookup r1 index 

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

277 
val b2 = Inttab.lookup r2 index 
22951  278 
in 
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24302
diff
changeset

279 
(add_row_entry r12_1 index @{term "lbound :: real => real"} ((names index)^"l") b1, 
32960
69916a850301
eliminated hard tabulators, guessing at each author's individual tabwidth;
wenzelm
parents:
31971
diff
changeset

280 
add_row_entry r12_2 index @{term "ubound :: real => real"} ((names index)^"u") b2) 
22951  281 
end 
23665
825bea0266db
adopted to new computing oracle and fixed bugs introduced by tuning
obua
parents:
23520
diff
changeset

282 

24302  283 
val (r1, r2) = abs_estimate xlen r1 r2 
23665
825bea0266db
adopted to new computing oracle and fixed bugs introduced by tuning
obua
parents:
23520
diff
changeset

284 

16784  285 
in 
24302  286 
(r1, r2) 
16784  287 
end 
22951  288 

16784  289 
fun load filename prec safe_propagation = 
290 
let 

22951  291 
val prog = Cplex.load_cplexFile filename 
292 
val prog = Cplex.elim_nonfree_bounds prog 

293 
val prog = Cplex.relax_strict_ineqs prog 

32960
69916a850301
eliminated hard tabulators, guessing at each author's individual tabwidth;
wenzelm
parents:
31971
diff
changeset

294 
val (maximize, c, A, b, (xlen, names, _)) = CplexFloatSparseMatrixConverter.convert_prog prog 
24302  295 
val (r1, r2) = calcr safe_propagation xlen names prec A b 
22951  296 
val _ = if maximize then () else raise Load "sorry, cannot handle minimization problems" 
297 
val (dualprog, indexof) = FloatSparseMatrixBuilder.dual_cplexProg c A b 

298 
val results = Cplex.solve dualprog 

46533  299 
val (_, v) = CplexFloatSparseMatrixConverter.convert_results results indexof 
24302  300 
(*val A = FloatSparseMatrixBuilder.cut_matrix v NONE A*) 
22951  301 
fun id x = x 
302 
val v = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 v 

303 
val b = FloatSparseMatrixBuilder.transpose_matrix (FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 b) 

304 
val c = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 c 

22964  305 
val (y1, _) = FloatSparseMatrixBuilder.approx_matrix prec Float.positive_part v 
22951  306 
val A = FloatSparseMatrixBuilder.approx_matrix prec id A 
307 
val (_,b2) = FloatSparseMatrixBuilder.approx_matrix prec id b 

308 
val c = FloatSparseMatrixBuilder.approx_matrix prec id c 

16784  309 
in 
24302  310 
(y1, A, b2, c, (r1, r2)) 
22951  311 
end handle CplexFloatSparseMatrixConverter.Converter s => (raise (Load ("Converter: "^s))) 
16784  312 

313 
end 