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

21056

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


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

20485

35 

16784

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)

20485

83 

16873

84 
val float_const = Float.float_const


85 

16784

86 
val zero = IntInf.fromInt 0


87 
val minus_one = IntInf.fromInt ~1


88 
val two = IntInf.fromInt 2

20485

89 

16873

90 
val mk_intinf = Float.mk_intinf

16784

91 

16873

92 
val mk_float = Float.mk_float

16784

93 


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


95 

16873

96 
fun approx_value_term prec f = Float.approx_float prec (fn (x, y) => (f x, f y))

16784

97 


98 
fun approx_value prec pprt value =

20485

99 
let


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


101 
in


102 
(cterm_of sg flower, cterm_of sg fupper)


103 
end

16784

104 


105 
fun sign_term t = cterm_of sg t


106 


107 
val empty_spvec = empty_vector_const


108 


109 
val empty_spmat = empty_matrix_const


110 


111 
fun mk_spvec_entry i f =


112 
let


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


114 
val term_f = mk_float f


115 
in


116 
HOLogic.mk_prod (term_i, term_f)


117 
end


118 


119 
fun mk_spmat_entry i e =


120 
let


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


122 
in


123 
HOLogic.mk_prod (term_i, e)


124 
end


125 


126 
fun cons_spvec h t = Cons_spvec_const $ h $ t


127 


128 
fun cons_spmat h t = Cons_spmat_const $ h $ t


129 


130 
fun approx_vector_term prec pprt vector =


131 
let

21056

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

16784

133 
let


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


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


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


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


138 
in


139 
(Cons_spvec_const $ elower $ vlower,


140 
Cons_spvec_const $ eupper $ vupper)


141 
end


142 
in

21056

143 
Inttab.fold app vector (empty_vector_const, empty_vector_const)

16784

144 
end


145 


146 
fun approx_matrix_term prec pprt matrix =

21056

147 
let


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


149 
let


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


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


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


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


154 
in


155 
(Cons_spmat_const $ elower $ mlower, Cons_spmat_const $ eupper $ mupper)


156 
end


157 
val (mlower, mupper) = Inttab.fold app matrix (empty_matrix_const, empty_matrix_const)


158 
in Inttab.fold app matrix (empty_matrix_const, empty_matrix_const) end;

16784

159 


160 
fun approx_vector prec pprt vector =


161 
let


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


163 
in


164 
(cterm_of sg l, cterm_of sg u)


165 
end


166 


167 
fun approx_matrix prec pprt matrix =


168 
let


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


170 
in


171 
(cterm_of sg l, cterm_of sg u)


172 
end


173 


174 


175 
exception Nat_expected of int;


176 


177 
val zero_interval = approx_value_term 1 I "0"


178 


179 
fun set_elem vector index str =


180 
if index < 0 then


181 
raise (Nat_expected index)


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


183 
vector


184 
else

17412

185 
Inttab.update (index, str) vector

16784

186 


187 
fun set_vector matrix index vector =


188 
if index < 0 then


189 
raise (Nat_expected index)


190 
else if Inttab.is_empty vector then


191 
matrix


192 
else

17412

193 
Inttab.update (index, vector) matrix

16784

194 


195 
val empty_matrix = Inttab.empty


196 
val empty_vector = Inttab.empty


197 


198 
(* dual stuff *)


199 


200 
structure cplex = Cplex


201 


202 
fun transpose_matrix matrix =

21056

203 
let


204 
fun upd j (i, s) =


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


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


207 
in Inttab.fold updm matrix empty_matrix end;

16784

208 


209 
exception No_name of string;


210 


211 
exception Superfluous_constr_right_hand_sides


212 


213 
fun cplexProg c A b =


214 
let


215 
val ytable = ref Inttab.empty


216 
fun indexof s =


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


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


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


220 


221 
fun nameof i =


222 
let


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

17412

224 
val _ = change ytable (Inttab.update (i, s))

16784

225 
in


226 
s


227 
end


228 


229 
fun split_numstr s =


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


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


232 
else (true, s)


233 


234 
fun mk_term index s =


235 
let


236 
val (p, s) = split_numstr s


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


238 
in


239 
if p then prod else cplex.cplexNeg prod


240 
end


241 


242 
fun vec2sum vector =

21056

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

16784

244 


245 
fun mk_constr index vector c =


246 
let

17412

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

16784

248 
val (p, s) = split_numstr s


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


250 
in


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


252 
end


253 


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


255 

21056

256 
val (list, b) = Inttab.fold


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


258 
A ([], b)

16784

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


260 


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


262 
cplex.cplexVar y, cplex.cplexLeq,


263 
cplex.cplexInf)


264 

21056

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

16784

266 


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


268 
in


269 
(prog, indexof)


270 
end


271 


272 


273 
fun dual_cplexProg c A b =


274 
let


275 
fun indexof s =


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


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


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


279 


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


281 


282 
fun split_numstr s =


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


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


285 
else (true, s)


286 


287 
fun mk_term index s =


288 
let


289 
val (p, s) = split_numstr s


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


291 
in


292 
if p then prod else cplex.cplexNeg prod


293 
end


294 


295 
fun vec2sum vector =

21056

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

16784

297 


298 
fun mk_constr index vector c =


299 
let

17412

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

16784

301 
val (p, s) = split_numstr s


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


303 
in


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


305 
end


306 


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


308 

21056

309 
val (list, c) = Inttab.fold


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


311 
(transpose_matrix A) ([], c)

16784

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


313 


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


315 
in


316 
(prog, indexof)


317 
end


318 


319 
fun cut_vector size v =

21056

320 
let


321 
val count = ref 0;


322 
fun app (i, s) = if (!count < size) then


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


324 
else I


325 
in


326 
Inttab.fold app v empty_vector


327 
end

16784

328 


329 
fun cut_matrix vfilter vsize m =

21056

330 
let


331 
fun app (i, v) =


332 
if is_none (Inttab.lookup vfilter i) then I


333 
else case vsize


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


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


336 
in Inttab.fold app m empty_matrix end


337 

17412

338 
fun v_elem_at v i = Inttab.lookup v i


339 
fun m_elem_at m i = Inttab.lookup m i

16784

340 


341 
fun v_only_elem v =


342 
case Inttab.min_key v of


343 
NONE => NONE


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


345 
NONE => SOME vmin


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


347 

21056

348 
fun v_fold f = Inttab.fold f;


349 
fun m_fold f = Inttab.fold f;

16784

350 


351 
end;
