15178

1 
signature FSPMLP =


2 
sig


3 
type linprog


4 


5 
val y : linprog > cterm


6 
val A : linprog > cterm * cterm


7 
val b : linprog > cterm


8 
val c : linprog > cterm * cterm


9 
val r : linprog > cterm


10 


11 
exception Load of string


12 


13 
val load : string > int > bool > linprog


14 
end


15 


16 
structure fspmlp : FSPMLP =


17 
struct


18 


19 
type linprog = cterm * (cterm * cterm) * cterm * (cterm * cterm) * cterm


20 


21 
fun y (c1, c2, c3, c4, c5) = c1


22 
fun A (c1, c2, c3, c4, c5) = c2


23 
fun b (c1, c2, c3, c4, c5) = c3


24 
fun c (c1, c2, c3, c4, c5) = c4


25 
fun r (c1, c2, c3, c4, c5) = c5


26 


27 
structure CplexFloatSparseMatrixConverter =


28 
MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = FloatSparseMatrixBuilder);


29 


30 
datatype bound_type = LOWER  UPPER


31 


32 
fun intbound_ord ((i1, b1),(i2,b2)) =


33 
if i1 < i2 then LESS


34 
else if i1 = i2 then


35 
(if b1 = b2 then EQUAL else if b1=LOWER then LESS else GREATER)


36 
else GREATER


37 


38 
structure Inttab = TableFun(type key = int val ord = (rev_order o int_ord));


39 


40 
structure VarGraph = TableFun(type key = int*bound_type val ord = intbound_ord);


41 
(* key > (float option) * (int > (float * (((float * float) * key) list)))) *)


42 
(* dest_key > (sure_bound * (row_index > (row_bound * (((coeff_lower * coeff_upper) * src_key) list)))) *)


43 


44 
exception Internal of string;


45 


46 
fun add_row_bound g dest_key row_index row_bound =


47 
let


48 
val x =


49 
case VarGraph.lookup (g, dest_key) of

15531

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


51 
 SOME (sure_bound, f) =>

15178

52 
(sure_bound,


53 
case Inttab.lookup (f, row_index) of

15531

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


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

15178

56 
in


57 
VarGraph.update ((dest_key, x), g)


58 
end


59 


60 
fun update_sure_bound g (key as (_, btype)) bound =


61 
let


62 
val x =


63 
case VarGraph.lookup (g, key) of

15531

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


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


66 
 SOME (SOME old_bound, f) =>


67 
(SOME ((case btype of

15178

68 
UPPER => FloatArith.min


69 
 LOWER => FloatArith.max)


70 
old_bound bound), f)


71 
in


72 
VarGraph.update ((key, x), g)


73 
end


74 


75 
fun get_sure_bound g key =


76 
case VarGraph.lookup (g, key) of

15531

77 
NONE => NONE


78 
 SOME (sure_bound, _) => sure_bound

15178

79 


80 
(*fun get_row_bound g key row_index =


81 
case VarGraph.lookup (g, key) of

15531

82 
NONE => NONE


83 
 SOME (sure_bound, f) =>

15178

84 
(case Inttab.lookup (f, row_index) of

15531

85 
NONE => NONE


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

15178

87 


88 
fun add_edge g src_key dest_key row_index coeff =


89 
case VarGraph.lookup (g, dest_key) of

15531

90 
NONE => raise (Internal "add_edge: dest_key not found")


91 
 SOME (sure_bound, f) =>

15178

92 
(case Inttab.lookup (f, row_index) of

15531

93 
NONE => raise (Internal "add_edge: row_index not found")


94 
 SOME (row_bound, sources) =>

15178

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


96 


97 
fun split_graph g =


98 
let


99 
fun split ((r1, r2), (key, (sure_bound, _))) =


100 
case sure_bound of

15531

101 
NONE => (r1, r2)


102 
 SOME bound =>

15178

103 
(case key of


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


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


106 
in


107 
VarGraph.foldl split ((Inttab.empty, Inttab.empty), g)


108 
end


109 


110 
fun it2list t =


111 
let


112 
fun tolist (l, a) = a::l


113 
in


114 
Inttab.foldl tolist ([], t)


115 
end


116 


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


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


119 
fun propagate_sure_bounds safe names g =


120 
let

15531

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

15178

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


123 
let


124 
fun mult_upper x (lower, upper) =


125 
if FloatArith.is_negative x then


126 
FloatArith.mul x lower


127 
else


128 
FloatArith.mul x upper


129 


130 
fun mult_lower x (lower, upper) =


131 
if FloatArith.is_negative x then


132 
FloatArith.mul x upper


133 
else


134 
FloatArith.mul x lower


135 


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


137 


138 
fun calc_sure_bound (sure_bound, (row_index, (row_bound, sources))) =


139 
let


140 
fun add_src_bound (sum, (coeff, src_key)) =


141 
case sum of

15531

142 
NONE => NONE


143 
 SOME x =>

15178

144 
(case get_sure_bound g src_key of

15531

145 
NONE => NONE


146 
 SOME src_sure_bound => SOME (FloatArith.add x (mult_btype src_sure_bound coeff)))

15178

147 
in

15570

148 
case Library.foldl add_src_bound (SOME row_bound, sources) of

15531

149 
NONE => sure_bound


150 
 new_sure_bound as (SOME new_bound) =>

15178

151 
(case sure_bound of

15531

152 
NONE => new_sure_bound


153 
 SOME old_bound =>


154 
SOME (case btype of

15178

155 
UPPER => FloatArith.min old_bound new_bound


156 
 LOWER => FloatArith.max old_bound new_bound))


157 
end


158 
in


159 
case VarGraph.lookup (g, key) of

15531

160 
NONE => NONE


161 
 SOME (sure_bound, f) =>

15178

162 
let


163 
val x = Inttab.foldl calc_sure_bound (sure_bound, f)


164 
in

15531

165 
if x = sure_bound then NONE else x

15178

166 
end


167 
end


168 


169 
fun propagate ((g, b), (key, _)) =


170 
case calc_sure_bound_from_sources g key of

15531

171 
NONE => (g,b)


172 
 SOME bound => (update_sure_bound g key bound,

15178

173 
if safe then


174 
case get_sure_bound g key of

15531

175 
NONE => true

15178

176 
 _ => b


177 
else


178 
true)


179 


180 
val (g, b) = VarGraph.foldl propagate ((g, false), g)


181 
in


182 
if b then propagate_sure_bounds safe names g else g


183 
end


184 


185 
exception Load of string;


186 


187 
fun calcr safe_propagation xlen names prec A b =


188 
let


189 
val empty = Inttab.empty


190 


191 
fun instab t i x = Inttab.update ((i,x), t)


192 


193 
fun isnegstr x = String.isPrefix "" x


194 
fun negstr x = if isnegstr x then String.extract (x, 1, NONE) else ""^x


195 


196 
fun test_1 (lower, upper) =


197 
if lower = upper then


198 
(if FloatArith.is_equal lower (IntInf.fromInt ~1, FloatArith.izero) then ~1


199 
else if FloatArith.is_equal lower (IntInf.fromInt 1, FloatArith.izero) then 1


200 
else 0)


201 
else 0


202 


203 
fun calcr (g, (row_index, a)) =


204 
let


205 
val b = FloatSparseMatrixBuilder.v_elem_at b row_index

15531

206 
val (_, b2) = ExactFloatingPoint.approx_decstr_by_bin prec (case b of NONE => "0"  SOME b => b)

15178

207 
val approx_a = FloatSparseMatrixBuilder.v_fold (fn (l, (i,s)) =>


208 
(i, ExactFloatingPoint.approx_decstr_by_bin prec s)::l) [] a


209 


210 
fun fold_dest_nodes (g, (dest_index, dest_value)) =


211 
let


212 
val dest_test = test_1 dest_value


213 
in


214 
if dest_test = 0 then


215 
g


216 
else let


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


218 
if dest_test = ~1 then


219 
((dest_index, LOWER), FloatArith.neg b2)


220 
else


221 
((dest_index, UPPER), b2)


222 


223 
fun fold_src_nodes (g, (src_index, src_value as (src_lower, src_upper))) =


224 
if src_index = dest_index then g


225 
else


226 
let


227 
val coeff = case dest_btype of


228 
UPPER => (FloatArith.neg src_upper, FloatArith.neg src_lower)


229 
 LOWER => src_value


230 
in


231 
if FloatArith.is_negative src_lower then


232 
add_edge g (src_index, UPPER) dest_key row_index coeff


233 
else


234 
add_edge g (src_index, LOWER) dest_key row_index coeff


235 
end


236 
in

15570

237 
Library.foldl fold_src_nodes ((add_row_bound g dest_key row_index row_bound), approx_a)

15178

238 
end


239 
end


240 
in


241 
case approx_a of


242 
[] => g


243 
 [(u, a)] =>


244 
let


245 
val atest = test_1 a


246 
in


247 
if atest = ~1 then


248 
update_sure_bound g (u, LOWER) (FloatArith.neg b2)


249 
else if atest = 1 then


250 
update_sure_bound g (u, UPPER) b2


251 
else


252 
g


253 
end

15570

254 
 _ => Library.foldl fold_dest_nodes (g, approx_a)

15178

255 
end


256 


257 
val g = FloatSparseMatrixBuilder.m_fold calcr VarGraph.empty A


258 
val g = propagate_sure_bounds safe_propagation names g


259 


260 
val (r1, r2) = split_graph g


261 


262 
fun abs_estimate i r1 r2 =


263 
if i = 0 then FloatSparseMatrixBuilder.empty_spmat


264 
else


265 
let


266 
val index = xleni


267 
val r = abs_estimate (i1) r1 r2

15531

268 
val b1 = case Inttab.lookup (r1, index) of NONE => raise (Load ("xvalue not bounded from below: "^(names index)))  SOME x => x


269 
val b2 = case Inttab.lookup (r2, index) of NONE => raise (Load ("xvalue not bounded from above: "^(names index)))  SOME x => x

15178

270 
val abs_max = FloatArith.max (FloatArith.neg (FloatArith.negative_part b1)) (FloatArith.positive_part b2)


271 
val vec = FloatSparseMatrixBuilder.cons_spvec (FloatSparseMatrixBuilder.mk_spvec_entry 0 abs_max) FloatSparseMatrixBuilder.empty_spvec


272 
in


273 
FloatSparseMatrixBuilder.cons_spmat (FloatSparseMatrixBuilder.mk_spmat_entry index vec) r


274 
end


275 
in


276 
FloatSparseMatrixBuilder.sign_term (abs_estimate xlen r1 r2)


277 
end


278 


279 
fun load filename prec safe_propagation =


280 
let


281 
val prog = Cplex.load_cplexFile filename


282 
val prog = Cplex.elim_nonfree_bounds prog


283 
val prog = Cplex.relax_strict_ineqs prog


284 
val (maximize, c, A, b, (xlen, names, _)) = CplexFloatSparseMatrixConverter.convert_prog prog


285 
val r = calcr safe_propagation xlen names prec A b


286 
val _ = if maximize then () else raise Load "sorry, cannot handle minimization problems"


287 
val (dualprog, indexof) = FloatSparseMatrixBuilder.dual_cplexProg c A b


288 
val results = Cplex.solve dualprog


289 
val (optimal,v) = CplexFloatSparseMatrixConverter.convert_results results indexof

15531

290 
val A = FloatSparseMatrixBuilder.cut_matrix v NONE A

15178

291 
fun id x = x


292 
val v = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 v


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


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


295 
val (y1, _) = FloatSparseMatrixBuilder.approx_matrix prec FloatArith.positive_part v


296 
val A = FloatSparseMatrixBuilder.approx_matrix prec id A


297 
val (_,b2) = FloatSparseMatrixBuilder.approx_matrix prec id b


298 
val c = FloatSparseMatrixBuilder.approx_matrix prec id c


299 
in


300 
(y1, A, b2, c, r)


301 
end handle CplexFloatSparseMatrixConverter.Converter s => (raise (Load ("Converter: "^s)))


302 


303 
end 