16784

1 
(* Title: HOL/Matrix/cplex/FloatSparseMatrixBuilder.ML


2 
ID: $Id$


3 
Author: Steven Obua


4 
*)


5 


6 
structure FloatSparseMatrixBuilder :


7 
sig


8 
include MATRIX_BUILDER


9 


10 
structure cplex : CPLEX


11 


12 
type float = IntInf.int*IntInf.int


13 
type floatfunc = float > float


14 


15 


16 
val float2cterm : IntInf.int * IntInf.int > cterm


17 


18 
val approx_value : int > floatfunc > string > cterm * cterm


19 
val approx_vector : int > floatfunc > vector > cterm * cterm


20 
val approx_matrix : int > floatfunc > matrix > cterm * cterm


21 


22 
val mk_spvec_entry : int > float > term


23 
val empty_spvec : term


24 
val cons_spvec : term > term > term


25 
val empty_spmat : term


26 
val mk_spmat_entry : int > term > term


27 
val cons_spmat : term > term > term


28 
val sign_term : term > cterm


29 


30 
val v_elem_at : vector > int > string option


31 
val m_elem_at : matrix > int > vector option


32 
val v_only_elem : vector > int option


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


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


35 


36 
val transpose_matrix : matrix > matrix


37 


38 
val cut_vector : int > vector > vector


39 
val cut_matrix : vector > (int option) > matrix > matrix


40 


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 


46 
val real_spmatT : typ


47 
val real_spvecT : typ


48 
end


49 
=


50 
struct


51 


52 


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


54 


55 
type vector = string Inttab.table


56 
type matrix = vector Inttab.table


57 
type float = IntInf.int*IntInf.int


58 
type floatfunc = float > float


59 


60 
val th = theory "FloatSparseMatrix"


61 
val sg = sign_of th


62 


63 
fun readtype s = Sign.intern_type sg s


64 
fun readterm s = Sign.intern_const sg s


65 


66 
val ty_list = readtype "list"


67 
val term_Nil = readterm "Nil"


68 
val term_Cons = readterm "Cons"


69 


70 
val spvec_elemT = HOLogic.mk_prodT (HOLogic.natT, HOLogic.realT)


71 
val spvecT = Type (ty_list, [spvec_elemT])


72 
val spmat_elemT = HOLogic.mk_prodT (HOLogic.natT, spvecT)


73 
val spmatT = Type (ty_list, [spmat_elemT])


74 


75 
val real_spmatT = spmatT


76 
val real_spvecT = spvecT


77 


78 
val empty_matrix_const = Const (term_Nil, spmatT)


79 
val empty_vector_const = Const (term_Nil, spvecT)


80 


81 
val Cons_spvec_const = Const (term_Cons, spvec_elemT > spvecT > spvecT)


82 
val Cons_spmat_const = Const (term_Cons, spmat_elemT > spmatT > spmatT)


83 


84 
val float_const = Const (readterm "float", HOLogic.mk_prodT (HOLogic.intT, HOLogic.intT) > HOLogic.realT)


85 


86 
val zero = IntInf.fromInt 0


87 
val minus_one = IntInf.fromInt ~1


88 
val two = IntInf.fromInt 2


89 


90 
fun mk_intinf ty n =


91 
let


92 
fun mk_bit n = if n = zero then HOLogic.false_const else HOLogic.true_const


93 


94 
fun bin_of n =


95 
if n = zero then HOLogic.pls_const


96 
else if n = minus_one then HOLogic.min_const


97 
else


98 
let


99 
(*val (q,r) = IntInf.divMod (n, two): doesn't work in SML 10.0.7, but in newer versions!!!*)


100 
val q = IntInf.div (n, two)


101 
val r = IntInf.mod (n, two)


102 
in


103 
HOLogic.bit_const $ bin_of q $ mk_bit r


104 
end


105 
in


106 
HOLogic.number_of_const ty $ (bin_of n)


107 
end


108 


109 
fun mk_float (a,b) =


110 
float_const $ (HOLogic.mk_prod ((mk_intinf HOLogic.intT a), (mk_intinf HOLogic.intT b)))


111 


112 
fun float2cterm (a,b) = cterm_of sg (mk_float (a,b))


113 


114 
fun approx_value_term prec f value =


115 
let


116 
val (flower, fupper) = ExactFloatingPoint.approx_decstr_by_bin prec value


117 
val (flower, fupper) = (f flower, f fupper)


118 
in


119 
(mk_float flower, mk_float fupper)


120 
end


121 


122 
fun approx_value prec pprt value =


123 
let


124 
val (flower, fupper) = approx_value_term prec pprt value


125 
in


126 
(cterm_of sg flower, cterm_of sg fupper)


127 
end


128 


129 
fun sign_term t = cterm_of sg t


130 


131 
val empty_spvec = empty_vector_const


132 


133 
val empty_spmat = empty_matrix_const


134 


135 
fun mk_spvec_entry i f =


136 
let


137 
val term_i = mk_intinf HOLogic.natT (IntInf.fromInt i)


138 
val term_f = mk_float f


139 
in


140 
HOLogic.mk_prod (term_i, term_f)


141 
end


142 


143 
fun mk_spmat_entry i e =


144 
let


145 
val term_i = mk_intinf HOLogic.natT (IntInf.fromInt i)


146 
in


147 
HOLogic.mk_prod (term_i, e)


148 
end


149 


150 
fun cons_spvec h t = Cons_spvec_const $ h $ t


151 


152 
fun cons_spmat h t = Cons_spmat_const $ h $ t


153 


154 
fun approx_vector_term prec pprt vector =


155 
let


156 
fun app ((vlower, vupper), (index, s)) =


157 
let


158 
val (flower, fupper) = approx_value_term prec pprt s


159 
val index = mk_intinf HOLogic.natT (IntInf.fromInt index)


160 
val elower = HOLogic.mk_prod (index, flower)


161 
val eupper = HOLogic.mk_prod (index, fupper)


162 
in


163 
(Cons_spvec_const $ elower $ vlower,


164 
Cons_spvec_const $ eupper $ vupper)


165 
end


166 
in


167 
Inttab.foldl app ((empty_vector_const, empty_vector_const), vector)


168 
end


169 


170 
fun approx_matrix_term prec pprt matrix =


171 
let


172 
fun app ((mlower, mupper), (index, vector)) =


173 
let


174 
val (vlower, vupper) = approx_vector_term prec pprt vector


175 
val index = mk_intinf HOLogic.natT (IntInf.fromInt index)


176 
val elower = HOLogic.mk_prod (index, vlower)


177 
val eupper = HOLogic.mk_prod (index, vupper)


178 
in


179 
(Cons_spmat_const $ elower $ mlower,


180 
Cons_spmat_const $ eupper $ mupper)


181 
end


182 


183 
val (mlower, mupper) = Inttab.foldl app ((empty_matrix_const, empty_matrix_const), matrix)


184 
in


185 
Inttab.foldl app ((empty_matrix_const, empty_matrix_const), matrix)


186 
end


187 


188 
fun approx_vector prec pprt vector =


189 
let


190 
val (l, u) = approx_vector_term prec pprt vector


191 
in


192 
(cterm_of sg l, cterm_of sg u)


193 
end


194 


195 
fun approx_matrix prec pprt matrix =


196 
let


197 
val (l, u) = approx_matrix_term prec pprt matrix


198 
in


199 
(cterm_of sg l, cterm_of sg u)


200 
end


201 


202 


203 
exception Nat_expected of int;


204 


205 
val zero_interval = approx_value_term 1 I "0"


206 


207 
fun set_elem vector index str =


208 
if index < 0 then


209 
raise (Nat_expected index)


210 
else if (approx_value_term 1 I str) = zero_interval then


211 
vector


212 
else


213 
Inttab.update ((index, str), vector)


214 


215 
fun set_vector matrix index vector =


216 
if index < 0 then


217 
raise (Nat_expected index)


218 
else if Inttab.is_empty vector then


219 
matrix


220 
else


221 
Inttab.update ((index, vector), matrix)


222 


223 
val empty_matrix = Inttab.empty


224 
val empty_vector = Inttab.empty


225 


226 
(* dual stuff *)


227 


228 
structure cplex = Cplex


229 


230 
fun transpose_matrix matrix =


231 
let


232 
fun upd m j i x =


233 
case Inttab.lookup (m, j) of


234 
SOME v => Inttab.update ((j, Inttab.update ((i, x), v)), m)


235 
 NONE => Inttab.update ((j, Inttab.update ((i, x), Inttab.empty)), m)


236 


237 
fun updv j (m, (i, s)) = upd m i j s


238 


239 
fun updm (m, (j, v)) = Inttab.foldl (updv j) (m, v)


240 
in


241 
Inttab.foldl updm (empty_matrix, matrix)


242 
end


243 


244 
exception No_name of string;


245 


246 
exception Superfluous_constr_right_hand_sides


247 


248 
fun cplexProg c A b =


249 
let


250 
val ytable = ref Inttab.empty


251 
fun indexof s =


252 
if String.size s = 0 then raise (No_name s)


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


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


255 


256 
fun nameof i =


257 
let


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


259 
val _ = ytable := (Inttab.update ((i, s), !ytable))


260 
in


261 
s


262 
end


263 


264 
fun split_numstr s =


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


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


267 
else (true, s)


268 


269 
fun mk_term index s =


270 
let


271 
val (p, s) = split_numstr s


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


273 
in


274 
if p then prod else cplex.cplexNeg prod


275 
end


276 


277 
fun vec2sum vector =


278 
cplex.cplexSum (Inttab.foldl (fn (list, (index, s)) => (mk_term index s)::list) ([], vector))


279 


280 
fun mk_constr index vector c =


281 
let


282 
val s = case Inttab.lookup (c, index) of SOME s => s  NONE => "0"


283 
val (p, s) = split_numstr s


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


285 
in


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


287 
end


288 


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


290 


291 
val (list, b) = Inttab.foldl


292 
(fn ((list, c), (index, v)) => ((mk_constr index v c)::list, delete index c))


293 
(([], b), A)


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


295 


296 
fun mk_free y = cplex.cplexBounds (cplex.cplexNeg cplex.cplexInf, cplex.cplexLeq,


297 
cplex.cplexVar y, cplex.cplexLeq,


298 
cplex.cplexInf)


299 


300 
val yvars = Inttab.foldl (fn (l, (i, y)) => (mk_free y)::l) ([], !ytable)


301 


302 
val prog = cplex.cplexProg ("original", cplex.cplexMaximize (vec2sum c), list, yvars)


303 
in


304 
(prog, indexof)


305 
end


306 


307 


308 
fun dual_cplexProg c A b =


309 
let


310 
fun indexof s =


311 
if String.size s = 0 then raise (No_name s)


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


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


314 


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


316 


317 
fun split_numstr s =


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


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


320 
else (true, s)


321 


322 
fun mk_term index s =


323 
let


324 
val (p, s) = split_numstr s


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


326 
in


327 
if p then prod else cplex.cplexNeg prod


328 
end


329 


330 
fun vec2sum vector =


331 
cplex.cplexSum (Inttab.foldl (fn (list, (index, s)) => (mk_term index s)::list) ([], vector))


332 


333 
fun mk_constr index vector c =


334 
let


335 
val s = case Inttab.lookup (c, index) of SOME s => s  NONE => "0"


336 
val (p, s) = split_numstr s


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


338 
in


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


340 
end


341 


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


343 


344 
val (list, c) = Inttab.foldl


345 
(fn ((list, c), (index, v)) => ((mk_constr index v c)::list, delete index c))


346 
(([], c), transpose_matrix A)


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


348 


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


350 
in


351 
(prog, indexof)


352 
end


353 


354 
fun cut_vector size v =


355 
let


356 
val count = ref 0


357 
fun app (v, (i, s)) =


358 
if (!count < size) then


359 
(count := !count +1 ; Inttab.update ((i,s),v))


360 
else


361 
v


362 
in


363 
Inttab.foldl app (empty_vector, v)


364 
end


365 


366 
fun cut_matrix vfilter vsize m =


367 
let


368 
fun app (m, (i, v)) =


369 
if (Inttab.lookup (vfilter, i) = NONE) then


370 
m


371 
else


372 
case vsize of


373 
NONE => Inttab.update ((i,v), m)


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


375 
in


376 
Inttab.foldl app (empty_matrix, m)


377 
end


378 


379 
fun v_elem_at v i = Inttab.lookup (v,i)


380 
fun m_elem_at m i = Inttab.lookup (m,i)


381 


382 
fun v_only_elem v =


383 
case Inttab.min_key v of


384 
NONE => NONE


385 
 SOME vmin => (case Inttab.max_key v of


386 
NONE => SOME vmin


387 
 SOME vmax => if vmin = vmax then SOME vmin else NONE)


388 


389 
fun v_fold f a v = Inttab.foldl f (a,v)


390 


391 
fun m_fold f a m = Inttab.foldl f (a,m)


392 


393 
end;
