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

22951

17 

16784

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

22951

40 

16784

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

22951

48 
end


49 
=

16784

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"

22951

61 

22578

62 
fun readtype s = Sign.intern_type th s


63 
fun readterm s = Sign.intern_const th s

22951

64 

16784

65 
val ty_list = readtype "list"


66 
val term_Nil = readterm "Nil"


67 
val term_Cons = readterm "Cons"

22951

68 

16784

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


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


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


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


73 


74 
val real_spmatT = spmatT


75 
val real_spvecT = spvecT


76 


77 
val empty_matrix_const = Const (term_Nil, spmatT)


78 
val empty_vector_const = Const (term_Nil, spvecT)


79 


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

22951

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

20485

82 

16873

83 
val float_const = Float.float_const


84 

16784

85 
val zero = IntInf.fromInt 0


86 
val minus_one = IntInf.fromInt ~1


87 
val two = IntInf.fromInt 2

20485

88 

16873

89 
val mk_intinf = Float.mk_intinf

16784

90 

22951

91 
val mk_float = Float.mk_float

16784

92 

22578

93 
fun float2cterm (a,b) = cterm_of th (mk_float (a,b))

16784

94 

16873

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

16784

96 

22951

97 
fun approx_value prec pprt value =

20485

98 
let


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


100 
in

22578

101 
(cterm_of th flower, cterm_of th fupper)

20485

102 
end

16784

103 

22578

104 
fun sign_term t = cterm_of th t

16784

105 


106 
val empty_spvec = empty_vector_const


107 


108 
val empty_spmat = empty_matrix_const


109 

22951

110 
fun mk_spvec_entry i f =

16784

111 
let

22951

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


113 
val term_f = mk_float f


114 
in


115 
HOLogic.mk_prod (term_i, term_f)

16784

116 
end


117 

22951

118 
fun mk_spmat_entry i e =

16784

119 
let

22951

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

16784

121 
in

22951

122 
HOLogic.mk_prod (term_i, e)

16784

123 
end


124 


125 
fun cons_spvec h t = Cons_spvec_const $ h $ t


126 

22951

127 
fun cons_spmat h t = Cons_spmat_const $ h $ t

16784

128 

22951

129 
fun approx_vector_term prec pprt vector =


130 
let


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


132 
let


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


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


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


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


137 
in


138 
(Cons_spvec_const $ elower $ vlower,


139 
Cons_spvec_const $ eupper $ vupper)


140 
end

16784

141 
in

22951

142 
Inttab.fold app vector (empty_vector_const, empty_vector_const)

16784

143 
end


144 


145 
fun approx_matrix_term prec pprt matrix =

21056

146 
let


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


148 
let

22951

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

21056

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


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


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


153 
in


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


155 
end


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


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

16784

158 


159 
fun approx_vector prec pprt vector =


160 
let

22951

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

16784

162 
in

22951

163 
(cterm_of th l, cterm_of th u)

16784

164 
end


165 

22951

166 
fun approx_matrix prec pprt matrix =

16784

167 
let

22951

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

16784

169 
in

22951

170 
(cterm_of th l, cterm_of th u)

16784

171 
end


172 


173 


174 
exception Nat_expected of int;


175 


176 
val zero_interval = approx_value_term 1 I "0"


177 

22951

178 
fun set_elem vector index str =


179 
if index < 0 then


180 
raise (Nat_expected index)


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


182 
vector


183 
else


184 
Inttab.update (index, str) vector

16784

185 

22951

186 
fun set_vector matrix index vector =

16784

187 
if index < 0 then

22951

188 
raise (Nat_expected index)

16784

189 
else if Inttab.is_empty vector then

22951

190 
matrix

16784

191 
else

22951

192 
Inttab.update (index, vector) matrix

16784

193 


194 
val empty_matrix = Inttab.empty


195 
val empty_vector = Inttab.empty


196 


197 
(* dual stuff *)


198 


199 
structure cplex = Cplex


200 

22951

201 
fun transpose_matrix matrix =

21056

202 
let


203 
fun upd j (i, s) =


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


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


206 
in Inttab.fold updm matrix empty_matrix end;

16784

207 


208 
exception No_name of string;


209 


210 
exception Superfluous_constr_right_hand_sides


211 

22951

212 
fun cplexProg c A b =

16784

213 
let

22951

214 
val ytable = ref Inttab.empty


215 
fun indexof s =


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


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


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

16784

219 

22951

220 
fun nameof i =


221 
let


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


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


224 
in


225 
s


226 
end


227 


228 
fun split_numstr s =


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


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


231 
else (true, s)


232 


233 
fun mk_term index s =


234 
let


235 
val (p, s) = split_numstr s


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


237 
in


238 
if p then prod else cplex.cplexNeg prod


239 
end

16784

240 

22951

241 
fun vec2sum vector =


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

16784

243 

22951

244 
fun mk_constr index vector c =


245 
let


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


247 
val (p, s) = split_numstr s


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


249 
in


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


251 
end


252 


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

16784

254 

22951

255 
val (list, b) = Inttab.fold


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


257 
A ([], b)


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

16784

259 

22951

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


261 
cplex.cplexVar y, cplex.cplexLeq,


262 
cplex.cplexInf)

16784

263 

22951

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

16784

265 

22951

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

16784

267 
in

22951

268 
(prog, indexof)

16784

269 
end


270 


271 

22951

272 
fun dual_cplexProg c A b =

16784

273 
let

22951

274 
fun indexof s =


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


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


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


278 


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

16784

280 

22951

281 
fun split_numstr s =


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


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


284 
else (true, s)


285 


286 
fun mk_term index s =


287 
let


288 
val (p, s) = split_numstr s


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


290 
in


291 
if p then prod else cplex.cplexNeg prod


292 
end

16784

293 

22951

294 
fun vec2sum vector =


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

16784

296 

22951

297 
fun mk_constr index vector c =


298 
let


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


300 
val (p, s) = split_numstr s


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


302 
in


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


304 
end

16784

305 

22951

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

16784

307 

22951

308 
val (list, c) = Inttab.fold


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


310 
(transpose_matrix A) ([], c)


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


312 


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

16784

314 
in

22951

315 
(prog, indexof)

16784

316 
end


317 

22951

318 
fun cut_vector size v =

21056

319 
let


320 
val count = ref 0;


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


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


323 
else I


324 
in


325 
Inttab.fold app v empty_vector


326 
end

16784

327 

22951

328 
fun cut_matrix vfilter vsize m =

21056

329 
let

22951

330 
fun app (i, v) =

21056

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


332 
else case vsize


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


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


335 
in Inttab.fold app m empty_matrix end


336 

17412

337 
fun v_elem_at v i = Inttab.lookup v i


338 
fun m_elem_at m i = Inttab.lookup m i

16784

339 

22951

340 
fun v_only_elem v =

16784

341 
case Inttab.min_key v of

22951

342 
NONE => NONE

16784

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

22951

344 
NONE => SOME vmin


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

16784

346 

21056

347 
fun v_fold f = Inttab.fold f;


348 
fun m_fold f = Inttab.fold f;

16784

349 


350 
end;
