src/HOL/Matrix/cplex/FloatSparseMatrixBuilder.ML
author wenzelm
Thu Sep 15 17:16:56 2005 +0200 (2005-09-15)
changeset 17412 e26cb20ef0cc
parent 17261 193b84a70ca4
child 20485 3078fd2eec7b
permissions -rw-r--r--
TableFun/Symtab: curried lookup and update;
obua@16784
     1
(*  Title:      HOL/Matrix/cplex/FloatSparseMatrixBuilder.ML
obua@16784
     2
    ID:         $Id$
obua@16784
     3
    Author:     Steven Obua
obua@16784
     4
*)
obua@16784
     5
obua@16784
     6
structure FloatSparseMatrixBuilder :
obua@16784
     7
sig
obua@16784
     8
    include MATRIX_BUILDER
obua@16784
     9
obua@16784
    10
    structure cplex : CPLEX
obua@16784
    11
obua@16784
    12
    type float = IntInf.int*IntInf.int
obua@16784
    13
    type floatfunc = float -> float
obua@16784
    14
obua@16784
    15
obua@16784
    16
    val float2cterm : IntInf.int * IntInf.int -> cterm
obua@16784
    17
       
obua@16784
    18
    val approx_value : int -> floatfunc -> string -> cterm * cterm
obua@16784
    19
    val approx_vector : int -> floatfunc -> vector -> cterm * cterm
obua@16784
    20
    val approx_matrix : int -> floatfunc -> matrix -> cterm * cterm
obua@16784
    21
obua@16784
    22
    val mk_spvec_entry : int -> float -> term
obua@16784
    23
    val empty_spvec : term
obua@16784
    24
    val cons_spvec : term -> term -> term
obua@16784
    25
    val empty_spmat : term
obua@16784
    26
    val mk_spmat_entry : int -> term -> term
obua@16784
    27
    val cons_spmat : term -> term -> term
obua@16784
    28
    val sign_term : term -> cterm
obua@16784
    29
obua@16784
    30
    val v_elem_at : vector -> int -> string option
obua@16784
    31
    val m_elem_at : matrix -> int -> vector option
obua@16784
    32
    val v_only_elem : vector -> int option
obua@16784
    33
    val v_fold : ('a * (int * string) -> 'a) -> 'a -> vector -> 'a
obua@16784
    34
    val m_fold : ('a * (int * vector) -> 'a) -> 'a -> matrix -> 'a
obua@16784
    35
				   
obua@16784
    36
    val transpose_matrix : matrix -> matrix
obua@16784
    37
obua@16784
    38
    val cut_vector : int -> vector -> vector
obua@16784
    39
    val cut_matrix : vector -> (int option) -> matrix -> matrix
obua@16784
    40
    
obua@16784
    41
    (* cplexProg c A b *)
obua@16784
    42
    val cplexProg : vector -> matrix -> vector -> (cplex.cplexProg * (string -> int))
obua@16784
    43
    (* dual_cplexProg c A b *)
obua@16784
    44
    val dual_cplexProg : vector -> matrix -> vector -> (cplex.cplexProg * (string -> int))
obua@16784
    45
obua@16784
    46
    val real_spmatT : typ
obua@16784
    47
    val real_spvecT : typ
obua@16784
    48
end 
obua@16784
    49
= 
obua@16784
    50
struct
obua@16784
    51
obua@16784
    52
obua@16784
    53
structure Inttab = TableFun(type key = int val ord = (rev_order o int_ord));
obua@16784
    54
obua@16784
    55
type vector = string Inttab.table
obua@16784
    56
type matrix = vector Inttab.table
obua@16784
    57
type float = IntInf.int*IntInf.int
obua@16784
    58
type floatfunc = float -> float
obua@16784
    59
obua@16784
    60
val th = theory "FloatSparseMatrix"
obua@16784
    61
val sg = sign_of th
obua@16784
    62
	 
obua@16784
    63
fun readtype s = Sign.intern_type sg s
obua@16784
    64
fun readterm s = Sign.intern_const sg s
obua@16784
    65
		 
obua@16784
    66
val ty_list = readtype "list"
obua@16784
    67
val term_Nil = readterm "Nil"
obua@16784
    68
val term_Cons = readterm "Cons"
obua@16784
    69
		
obua@16784
    70
val spvec_elemT = HOLogic.mk_prodT (HOLogic.natT, HOLogic.realT)
obua@16784
    71
val spvecT = Type (ty_list, [spvec_elemT])
obua@16784
    72
val spmat_elemT = HOLogic.mk_prodT (HOLogic.natT, spvecT)
obua@16784
    73
val spmatT = Type (ty_list, [spmat_elemT])
obua@16784
    74
obua@16784
    75
val real_spmatT = spmatT
obua@16784
    76
val real_spvecT = spvecT
obua@16784
    77
obua@16784
    78
val empty_matrix_const = Const (term_Nil, spmatT)
obua@16784
    79
val empty_vector_const = Const (term_Nil, spvecT)
obua@16784
    80
obua@16784
    81
val Cons_spvec_const = Const (term_Cons, spvec_elemT --> spvecT --> spvecT)
obua@16784
    82
val Cons_spmat_const = Const (term_Cons, spmat_elemT --> spmatT --> spmatT) 
obua@16784
    83
			 
obua@16873
    84
val float_const = Float.float_const
obua@16873
    85
obua@16873
    86
(*val float_const = Const (readterm "float", HOLogic.mk_prodT (HOLogic.intT, HOLogic.intT) --> HOLogic.realT)*)
obua@16784
    87
		  
obua@16784
    88
val zero = IntInf.fromInt 0
obua@16784
    89
val minus_one = IntInf.fromInt ~1
obua@16784
    90
val two = IntInf.fromInt 2
obua@16784
    91
	  
obua@16873
    92
val mk_intinf = Float.mk_intinf
obua@16873
    93
(*    let
obua@16784
    94
	fun mk_bit n = if n = zero then HOLogic.false_const else HOLogic.true_const
obua@16784
    95
								 
obua@16784
    96
	fun bin_of n = 
obua@16784
    97
	    if n = zero then HOLogic.pls_const
obua@16784
    98
	    else if n = minus_one then HOLogic.min_const
obua@16784
    99
	    else 
obua@16784
   100
		let 
obua@16784
   101
		    (*val (q,r) = IntInf.divMod (n, two): doesn't work in SML 10.0.7, but in newer versions!!!*)
obua@16784
   102
	            val q = IntInf.div (n, two)
obua@16784
   103
		    val r = IntInf.mod (n, two)
obua@16784
   104
		in
obua@16784
   105
		    HOLogic.bit_const $ bin_of q $ mk_bit r
obua@16784
   106
		end
obua@16784
   107
    in 
obua@16784
   108
	HOLogic.number_of_const ty $ (bin_of n)
obua@16873
   109
    end*)
obua@16784
   110
obua@16873
   111
val mk_float  = Float.mk_float 
obua@16873
   112
(*    float_const $ (HOLogic.mk_prod ((mk_intinf HOLogic.intT a), (mk_intinf HOLogic.intT b)))*)
obua@16784
   113
obua@16784
   114
fun float2cterm (a,b) = cterm_of sg (mk_float (a,b))
obua@16784
   115
obua@16873
   116
fun approx_value_term prec f = Float.approx_float prec (fn (x, y) => (f x, f y))
obua@16873
   117
(*    let
obua@16784
   118
	val (flower, fupper) = ExactFloatingPoint.approx_decstr_by_bin prec value
obua@16784
   119
	val (flower, fupper) = (f flower, f fupper)
obua@16784
   120
    in
obua@16784
   121
	(mk_float flower, mk_float fupper)
obua@16873
   122
    end*)
obua@16784
   123
obua@16784
   124
fun approx_value prec pprt value = 
obua@16784
   125
    let
obua@16784
   126
	val (flower, fupper) = approx_value_term prec pprt value			       
obua@16784
   127
    in
obua@16784
   128
	(cterm_of sg flower, cterm_of sg fupper)
obua@16784
   129
    end
obua@16784
   130
obua@16784
   131
fun sign_term t = cterm_of sg t
obua@16784
   132
obua@16784
   133
val empty_spvec = empty_vector_const
obua@16784
   134
obua@16784
   135
val empty_spmat = empty_matrix_const
obua@16784
   136
obua@16784
   137
fun mk_spvec_entry i f = 
obua@16784
   138
    let
obua@16784
   139
	val term_i = mk_intinf HOLogic.natT (IntInf.fromInt i)
obua@16784
   140
	val term_f = mk_float f
obua@16784
   141
    in 
obua@16784
   142
	HOLogic.mk_prod (term_i, term_f)
obua@16784
   143
    end
obua@16784
   144
obua@16784
   145
fun mk_spmat_entry i e = 
obua@16784
   146
    let
obua@16784
   147
	val term_i = mk_intinf HOLogic.natT (IntInf.fromInt i)
obua@16784
   148
    in
obua@16784
   149
	HOLogic.mk_prod (term_i, e)
obua@16784
   150
    end
obua@16784
   151
obua@16784
   152
fun cons_spvec h t = Cons_spvec_const $ h $ t
obua@16784
   153
obua@16784
   154
fun cons_spmat h t = Cons_spmat_const $ h $ t 
obua@16784
   155
obua@16784
   156
fun approx_vector_term prec pprt vector = 
obua@16784
   157
    let 	 
obua@16784
   158
	fun app ((vlower, vupper), (index, s)) = 
obua@16784
   159
	    let
obua@16784
   160
		val (flower, fupper) = approx_value_term prec pprt s
obua@16784
   161
		val index = mk_intinf HOLogic.natT (IntInf.fromInt index)
obua@16784
   162
		val elower = HOLogic.mk_prod (index, flower)
obua@16784
   163
		val eupper = HOLogic.mk_prod (index, fupper)
obua@16784
   164
	    in
obua@16784
   165
		(Cons_spvec_const $ elower $ vlower,
obua@16784
   166
		 Cons_spvec_const $ eupper $ vupper)	
obua@16784
   167
	    end
obua@16784
   168
    in
obua@16784
   169
	Inttab.foldl app ((empty_vector_const, empty_vector_const), vector)	
obua@16784
   170
    end
obua@16784
   171
obua@16784
   172
fun approx_matrix_term prec pprt matrix =
obua@16784
   173
    let 	 
obua@16784
   174
	fun app ((mlower, mupper), (index, vector)) = 
obua@16784
   175
	    let
obua@16784
   176
		val (vlower, vupper) = approx_vector_term prec pprt vector  
obua@16784
   177
		val index = mk_intinf HOLogic.natT (IntInf.fromInt index)
obua@16784
   178
		val elower = HOLogic.mk_prod (index, vlower)
obua@16784
   179
		val eupper = HOLogic.mk_prod (index, vupper)
obua@16784
   180
	    in
obua@16784
   181
		(Cons_spmat_const $ elower $ mlower,
obua@16784
   182
		 Cons_spmat_const $ eupper $ mupper)		
obua@16784
   183
	    end
obua@16784
   184
		
obua@16784
   185
	val (mlower, mupper) = Inttab.foldl app ((empty_matrix_const, empty_matrix_const), matrix)
obua@16784
   186
    in
obua@16784
   187
	Inttab.foldl app ((empty_matrix_const, empty_matrix_const), matrix)	
obua@16784
   188
    end
obua@16784
   189
obua@16784
   190
fun approx_vector prec pprt vector =
obua@16784
   191
    let
obua@16784
   192
	val (l, u) = approx_vector_term prec pprt vector
obua@16784
   193
    in
obua@16784
   194
	(cterm_of sg l, cterm_of sg u)
obua@16784
   195
    end
obua@16784
   196
obua@16784
   197
fun approx_matrix prec pprt matrix = 
obua@16784
   198
    let
obua@16784
   199
	val (l, u) = approx_matrix_term prec pprt matrix
obua@16784
   200
    in
obua@16784
   201
	(cterm_of sg l, cterm_of sg u)
obua@16784
   202
    end
obua@16784
   203
obua@16784
   204
obua@16784
   205
exception Nat_expected of int;
obua@16784
   206
obua@16784
   207
val zero_interval = approx_value_term 1 I "0"
obua@16784
   208
obua@16784
   209
fun set_elem vector index str = 
obua@16784
   210
    if index < 0 then 
obua@16784
   211
	raise (Nat_expected index)
obua@16784
   212
    else if (approx_value_term 1 I str) = zero_interval then 
obua@16784
   213
	vector
obua@16784
   214
    else  
wenzelm@17412
   215
	Inttab.update (index, str) vector
obua@16784
   216
obua@16784
   217
fun set_vector matrix index vector = 
obua@16784
   218
    if index < 0 then
obua@16784
   219
	raise (Nat_expected index)
obua@16784
   220
    else if Inttab.is_empty vector then
obua@16784
   221
	matrix
obua@16784
   222
    else
wenzelm@17412
   223
	Inttab.update (index, vector) matrix
obua@16784
   224
obua@16784
   225
val empty_matrix = Inttab.empty
obua@16784
   226
val empty_vector = Inttab.empty
obua@16784
   227
obua@16784
   228
(* dual stuff *)
obua@16784
   229
obua@16784
   230
structure cplex = Cplex
obua@16784
   231
obua@16784
   232
fun transpose_matrix matrix = 
obua@16784
   233
    let
obua@16784
   234
	fun upd m j i x =
wenzelm@17412
   235
	    case Inttab.lookup m j of
wenzelm@17412
   236
		SOME v => Inttab.update (j, Inttab.update (i, x) v) m
wenzelm@17412
   237
	      | NONE => Inttab.update (j, Inttab.update (i, x) Inttab.empty) m
obua@16784
   238
obua@16784
   239
	fun updv j (m, (i, s)) = upd m i j s
obua@16784
   240
obua@16784
   241
	fun updm (m, (j, v)) = Inttab.foldl (updv j) (m, v)
obua@16784
   242
    in
obua@16784
   243
	Inttab.foldl updm (empty_matrix, matrix)
obua@16784
   244
    end
obua@16784
   245
obua@16784
   246
exception No_name of string;
obua@16784
   247
obua@16784
   248
exception Superfluous_constr_right_hand_sides
obua@16784
   249
obua@16784
   250
fun cplexProg c A b = 
obua@16784
   251
    let
obua@16784
   252
	val ytable = ref Inttab.empty
obua@16784
   253
	fun indexof s = 
obua@16784
   254
	    if String.size s = 0 then raise (No_name s)
obua@16784
   255
	    else case Int.fromString (String.extract(s, 1, NONE)) of 
obua@16784
   256
		     SOME i => i | NONE => raise (No_name s)
obua@16784
   257
				     
obua@16784
   258
	fun nameof i = 
obua@16784
   259
	    let 
obua@16784
   260
		val s = "x"^(Int.toString i)
wenzelm@17412
   261
		val _ = change ytable (Inttab.update (i, s))
obua@16784
   262
	    in
obua@16784
   263
		s
obua@16784
   264
	    end
obua@16784
   265
	    
obua@16784
   266
	fun split_numstr s = 
obua@16784
   267
	    if String.isPrefix "-" s then (false,String.extract(s, 1, NONE))
obua@16784
   268
	    else if String.isPrefix "+" s then (true, String.extract(s, 1, NONE))
obua@16784
   269
	    else (true, s)
obua@16784
   270
obua@16784
   271
	fun mk_term index s =
obua@16784
   272
	    let 
obua@16784
   273
		val (p, s) = split_numstr s
obua@16784
   274
		val prod = cplex.cplexProd (cplex.cplexNum s, cplex.cplexVar (nameof index))
obua@16784
   275
	    in
obua@16784
   276
		if p then prod else cplex.cplexNeg prod
obua@16784
   277
	    end
obua@16784
   278
obua@16784
   279
	fun vec2sum vector =
obua@16784
   280
	    cplex.cplexSum (Inttab.foldl (fn (list, (index, s)) => (mk_term index s)::list) ([], vector))
obua@16784
   281
		       		       
obua@16784
   282
	fun mk_constr index vector c = 
obua@16784
   283
	    let 
wenzelm@17412
   284
		val s = case Inttab.lookup c index of SOME s => s | NONE => "0"
obua@16784
   285
		val (p, s) = split_numstr s
obua@16784
   286
		val num = if p then cplex.cplexNum s else cplex.cplexNeg (cplex.cplexNum s)
obua@16784
   287
	    in
obua@16784
   288
		(NONE, cplex.cplexConstr (cplex.cplexLeq, (vec2sum vector, num)))
obua@16784
   289
	    end
obua@16784
   290
obua@16784
   291
	fun delete index c = Inttab.delete index c handle Inttab.UNDEF _ => c
obua@16784
   292
obua@16784
   293
	val (list, b) = Inttab.foldl 
obua@16784
   294
			    (fn ((list, c), (index, v)) => ((mk_constr index v c)::list, delete index c))
obua@16784
   295
			    (([], b), A)
obua@16784
   296
	val _ = if Inttab.is_empty b then () else raise Superfluous_constr_right_hand_sides
obua@16784
   297
obua@16784
   298
	fun mk_free y = cplex.cplexBounds (cplex.cplexNeg cplex.cplexInf, cplex.cplexLeq, 
obua@16784
   299
					   cplex.cplexVar y, cplex.cplexLeq,
obua@16784
   300
					   cplex.cplexInf)
obua@16784
   301
obua@16784
   302
	val yvars = Inttab.foldl (fn (l, (i, y)) => (mk_free y)::l) ([], !ytable)
obua@16784
   303
obua@16784
   304
	val prog = cplex.cplexProg ("original", cplex.cplexMaximize (vec2sum c), list, yvars)
obua@16784
   305
    in
obua@16784
   306
	(prog, indexof)
obua@16784
   307
    end
obua@16784
   308
obua@16784
   309
obua@16784
   310
fun dual_cplexProg c A b = 
obua@16784
   311
    let
obua@16784
   312
	fun indexof s = 
obua@16784
   313
	    if String.size s = 0 then raise (No_name s)
obua@16784
   314
	    else case Int.fromString (String.extract(s, 1, NONE)) of 
obua@16784
   315
		     SOME i => i | NONE => raise (No_name s)
obua@16784
   316
				     
obua@16784
   317
	fun nameof i = "y"^(Int.toString i)
obua@16784
   318
	    
obua@16784
   319
	fun split_numstr s = 
obua@16784
   320
	    if String.isPrefix "-" s then (false,String.extract(s, 1, NONE))
obua@16784
   321
	    else if String.isPrefix "+" s then (true, String.extract(s, 1, NONE))
obua@16784
   322
	    else (true, s)
obua@16784
   323
obua@16784
   324
	fun mk_term index s =
obua@16784
   325
	    let 
obua@16784
   326
		val (p, s) = split_numstr s
obua@16784
   327
		val prod = cplex.cplexProd (cplex.cplexNum s, cplex.cplexVar (nameof index))
obua@16784
   328
	    in
obua@16784
   329
		if p then prod else cplex.cplexNeg prod
obua@16784
   330
	    end
obua@16784
   331
obua@16784
   332
	fun vec2sum vector =
obua@16784
   333
	    cplex.cplexSum (Inttab.foldl (fn (list, (index, s)) => (mk_term index s)::list) ([], vector))
obua@16784
   334
		       		       
obua@16784
   335
	fun mk_constr index vector c = 
obua@16784
   336
	    let 
wenzelm@17412
   337
		val s = case Inttab.lookup c index of SOME s => s | NONE => "0"
obua@16784
   338
		val (p, s) = split_numstr s
obua@16784
   339
		val num = if p then cplex.cplexNum s else cplex.cplexNeg (cplex.cplexNum s)
obua@16784
   340
	    in
obua@16784
   341
		(NONE, cplex.cplexConstr (cplex.cplexEq, (vec2sum vector, num)))
obua@16784
   342
	    end
obua@16784
   343
obua@16784
   344
	fun delete index c = Inttab.delete index c handle Inttab.UNDEF _ => c
obua@16784
   345
obua@16784
   346
	val (list, c) = Inttab.foldl 
obua@16784
   347
			    (fn ((list, c), (index, v)) => ((mk_constr index v c)::list, delete index c))
obua@16784
   348
			    (([], c), transpose_matrix A)
obua@16784
   349
	val _ = if Inttab.is_empty c then () else raise Superfluous_constr_right_hand_sides
obua@16784
   350
obua@16784
   351
	val prog = cplex.cplexProg ("dual", cplex.cplexMinimize (vec2sum b), list, [])
obua@16784
   352
    in
obua@16784
   353
	(prog, indexof)
obua@16784
   354
    end
obua@16784
   355
obua@16784
   356
fun cut_vector size v = 
obua@16784
   357
    let
obua@16784
   358
	val count = ref 0 
obua@16784
   359
	fun app (v, (i, s)) = 
obua@16784
   360
	    if (!count < size) then
wenzelm@17412
   361
		(count := !count +1 ; Inttab.update (i,s) v)
obua@16784
   362
	    else
obua@16784
   363
		v
obua@16784
   364
    in
obua@16784
   365
	Inttab.foldl app (empty_vector, v)
obua@16784
   366
    end
obua@16784
   367
obua@16784
   368
fun cut_matrix vfilter vsize m = 
obua@16784
   369
    let 
obua@16784
   370
	fun app (m, (i, v)) = 
wenzelm@17412
   371
	    if Inttab.lookup vfilter i = NONE then 
obua@16784
   372
		m 
obua@16784
   373
	    else 
obua@16784
   374
		case vsize of
wenzelm@17412
   375
		    NONE => Inttab.update (i,v) m
wenzelm@17412
   376
		  | SOME s => Inttab.update (i, cut_vector s v) m
obua@16784
   377
    in
obua@16784
   378
	Inttab.foldl app (empty_matrix, m)
obua@16784
   379
    end
obua@16784
   380
		 
wenzelm@17412
   381
fun v_elem_at v i = Inttab.lookup v i
wenzelm@17412
   382
fun m_elem_at m i = Inttab.lookup m i
obua@16784
   383
obua@16784
   384
fun v_only_elem v = 
obua@16784
   385
    case Inttab.min_key v of
obua@16784
   386
	NONE => NONE
obua@16784
   387
      | SOME vmin => (case Inttab.max_key v of
obua@16784
   388
			  NONE => SOME vmin
obua@16784
   389
			| SOME vmax => if vmin = vmax then SOME vmin else NONE)
obua@16784
   390
obua@16784
   391
fun v_fold f a v = Inttab.foldl f (a,v) 
obua@16784
   392
obua@16784
   393
fun m_fold f a m = Inttab.foldl f (a,m)
obua@16784
   394
obua@16784
   395
end;