src/HOL/Matrix/cplex/FloatSparseMatrixBuilder.ML
author obua
Mon Jul 09 17:39:55 2007 +0200 (2007-07-09)
changeset 23665 825bea0266db
parent 23261 85f27f79232f
child 23883 7d5aa704454e
permissions -rw-r--r--
adopted to new computing oracle and fixed bugs introduced by tuning
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
haftmann@22964
     6
signature FLOAT_SPARSE_MATIRX_BUILDER =
obua@16784
     7
sig
haftmann@22964
     8
  include MATRIX_BUILDER
obua@16784
     9
haftmann@22964
    10
  structure cplex : CPLEX
obua@16784
    11
haftmann@22964
    12
  type float = Float.float
haftmann@22964
    13
  val approx_value : int -> (float -> float) -> string -> term * term
haftmann@22964
    14
  val approx_vector : int -> (float -> float) -> vector -> term * term
haftmann@22964
    15
  val approx_matrix : int -> (float -> float) -> matrix -> term * term
obua@16784
    16
haftmann@23261
    17
  val mk_spvec_entry : integer -> float -> term
haftmann@23261
    18
  val mk_spmat_entry : integer -> term -> term
haftmann@22964
    19
  val spvecT: typ
haftmann@22964
    20
  val spmatT: typ
haftmann@22964
    21
  
haftmann@22964
    22
  val v_elem_at : vector -> int -> string option
haftmann@22964
    23
  val m_elem_at : matrix -> int -> vector option
haftmann@22964
    24
  val v_only_elem : vector -> int option
haftmann@22964
    25
  val v_fold : (int * string -> 'a -> 'a) -> vector -> 'a -> 'a
haftmann@22964
    26
  val m_fold : (int * vector -> 'a -> 'a) -> matrix -> 'a -> 'a
haftmann@20485
    27
haftmann@22964
    28
  val transpose_matrix : matrix -> matrix
obua@16784
    29
haftmann@22964
    30
  val cut_vector : int -> vector -> vector
haftmann@22964
    31
  val cut_matrix : vector -> int option -> matrix -> matrix
haftmann@22951
    32
haftmann@22964
    33
  (* cplexProg c A b *)
haftmann@22964
    34
  val cplexProg : vector -> matrix -> vector -> cplex.cplexProg * (string -> int)
haftmann@22964
    35
  (* dual_cplexProg c A b *)
haftmann@22964
    36
  val dual_cplexProg : vector -> matrix -> vector -> cplex.cplexProg * (string -> int)
haftmann@22964
    37
end;
obua@16784
    38
haftmann@22964
    39
structure FloatSparseMatrixBuilder : FLOAT_SPARSE_MATIRX_BUILDER =
obua@16784
    40
struct
obua@16784
    41
haftmann@22964
    42
type float = Float.float
haftmann@22964
    43
structure Inttab = TableFun(type key = int val ord = rev_order o int_ord);
obua@16784
    44
obua@16784
    45
type vector = string Inttab.table
obua@16784
    46
type matrix = vector Inttab.table
haftmann@20485
    47
haftmann@22964
    48
val spvec_elemT = HOLogic.mk_prodT (HOLogic.natT, HOLogic.realT);
haftmann@22964
    49
val spvecT = HOLogic.listT spvec_elemT;
haftmann@22964
    50
val spmat_elemT = HOLogic.mk_prodT (HOLogic.natT, spvecT);
haftmann@22964
    51
val spmatT = HOLogic.listT spmat_elemT;
obua@16784
    52
haftmann@22964
    53
fun approx_value prec f =
haftmann@22964
    54
  FloatArith.approx_float prec (fn (x, y) => (f x, f y));
obua@16784
    55
haftmann@22951
    56
fun mk_spvec_entry i f =
haftmann@22964
    57
  HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, FloatArith.mk_float f);
obua@16784
    58
haftmann@22951
    59
fun mk_spmat_entry i e =
haftmann@22964
    60
  HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, e);
obua@16784
    61
obua@16784
    62
fun approx_vector prec pprt vector =
haftmann@22964
    63
  let
haftmann@22964
    64
    fun app (index, s) (lower, upper) =
haftmann@22964
    65
      let
haftmann@22964
    66
        val (flower, fupper) = approx_value prec pprt s
haftmann@23261
    67
        val index = HOLogic.mk_number HOLogic.natT (Integer.int index)
haftmann@22964
    68
        val elower = HOLogic.mk_prod (index, flower)
haftmann@22964
    69
        val eupper = HOLogic.mk_prod (index, fupper)
haftmann@22964
    70
      in (elower :: lower, eupper :: upper) end;
haftmann@22964
    71
  in
obua@23665
    72
    pairself (HOLogic.mk_list spvec_elemT) (Inttab.fold app vector ([], []))
haftmann@22964
    73
  end;
obua@16784
    74
haftmann@22964
    75
fun approx_matrix prec pprt vector =
haftmann@22964
    76
  let
haftmann@22964
    77
    fun app (index, v) (lower, upper) =
haftmann@22964
    78
      let
haftmann@22964
    79
        val (flower, fupper) = approx_vector prec pprt v
haftmann@23261
    80
        val index = HOLogic.mk_number HOLogic.natT (Integer.int index)
haftmann@22964
    81
        val elower = HOLogic.mk_prod (index, flower)
haftmann@22964
    82
        val eupper = HOLogic.mk_prod (index, fupper)
haftmann@22964
    83
      in (elower :: lower, eupper :: upper) end;
haftmann@22964
    84
  in
obua@23665
    85
    pairself (HOLogic.mk_list spmat_elemT) (Inttab.fold app vector ([], []))
haftmann@22964
    86
  end;
obua@16784
    87
obua@16784
    88
exception Nat_expected of int;
obua@16784
    89
haftmann@22964
    90
val zero_interval = approx_value 1 I "0"
obua@16784
    91
haftmann@22951
    92
fun set_elem vector index str =
haftmann@22951
    93
    if index < 0 then
haftmann@22951
    94
        raise (Nat_expected index)
haftmann@22964
    95
    else if (approx_value 1 I str) = zero_interval then
haftmann@22951
    96
        vector
haftmann@22951
    97
    else
haftmann@22951
    98
        Inttab.update (index, str) vector
obua@16784
    99
haftmann@22951
   100
fun set_vector matrix index vector =
obua@16784
   101
    if index < 0 then
haftmann@22951
   102
        raise (Nat_expected index)
obua@16784
   103
    else if Inttab.is_empty vector then
haftmann@22951
   104
        matrix
obua@16784
   105
    else
haftmann@22951
   106
        Inttab.update (index, vector) matrix
obua@16784
   107
obua@16784
   108
val empty_matrix = Inttab.empty
obua@16784
   109
val empty_vector = Inttab.empty
obua@16784
   110
obua@16784
   111
(* dual stuff *)
obua@16784
   112
obua@16784
   113
structure cplex = Cplex
obua@16784
   114
haftmann@22951
   115
fun transpose_matrix matrix =
haftmann@21056
   116
  let
haftmann@21056
   117
    fun upd j (i, s) =
haftmann@21056
   118
      Inttab.map_default (i, Inttab.empty) (Inttab.update (j, s));
haftmann@21056
   119
    fun updm (j, v) = Inttab.fold (upd j) v;
haftmann@21056
   120
  in Inttab.fold updm matrix empty_matrix end;
obua@16784
   121
obua@16784
   122
exception No_name of string;
obua@16784
   123
obua@16784
   124
exception Superfluous_constr_right_hand_sides
obua@16784
   125
haftmann@22951
   126
fun cplexProg c A b =
obua@16784
   127
    let
haftmann@22951
   128
        val ytable = ref Inttab.empty
haftmann@22951
   129
        fun indexof s =
haftmann@22951
   130
            if String.size s = 0 then raise (No_name s)
haftmann@22951
   131
            else case Int.fromString (String.extract(s, 1, NONE)) of
haftmann@22951
   132
                     SOME i => i | NONE => raise (No_name s)
obua@16784
   133
haftmann@22951
   134
        fun nameof i =
haftmann@22951
   135
            let
haftmann@22951
   136
                val s = "x"^(Int.toString i)
haftmann@22951
   137
                val _ = change ytable (Inttab.update (i, s))
haftmann@22951
   138
            in
haftmann@22951
   139
                s
haftmann@22951
   140
            end
haftmann@22951
   141
haftmann@22951
   142
        fun split_numstr s =
haftmann@22951
   143
            if String.isPrefix "-" s then (false,String.extract(s, 1, NONE))
haftmann@22951
   144
            else if String.isPrefix "+" s then (true, String.extract(s, 1, NONE))
haftmann@22951
   145
            else (true, s)
haftmann@22951
   146
haftmann@22951
   147
        fun mk_term index s =
haftmann@22951
   148
            let
haftmann@22951
   149
                val (p, s) = split_numstr s
haftmann@22951
   150
                val prod = cplex.cplexProd (cplex.cplexNum s, cplex.cplexVar (nameof index))
haftmann@22951
   151
            in
haftmann@22951
   152
                if p then prod else cplex.cplexNeg prod
haftmann@22951
   153
            end
obua@16784
   154
haftmann@22951
   155
        fun vec2sum vector =
haftmann@22951
   156
            cplex.cplexSum (Inttab.fold (fn (index, s) => fn list => (mk_term index s) :: list) vector [])
obua@16784
   157
haftmann@22951
   158
        fun mk_constr index vector c =
haftmann@22951
   159
            let
haftmann@22951
   160
                val s = case Inttab.lookup c index of SOME s => s | NONE => "0"
haftmann@22951
   161
                val (p, s) = split_numstr s
haftmann@22951
   162
                val num = if p then cplex.cplexNum s else cplex.cplexNeg (cplex.cplexNum s)
haftmann@22951
   163
            in
haftmann@22951
   164
                (NONE, cplex.cplexConstr (cplex.cplexLeq, (vec2sum vector, num)))
haftmann@22951
   165
            end
haftmann@22951
   166
haftmann@22951
   167
        fun delete index c = Inttab.delete index c handle Inttab.UNDEF _ => c
obua@16784
   168
haftmann@22951
   169
        val (list, b) = Inttab.fold
haftmann@22951
   170
                            (fn (index, v) => fn (list, c) => ((mk_constr index v c)::list, delete index c))
haftmann@22951
   171
                            A ([], b)
haftmann@22951
   172
        val _ = if Inttab.is_empty b then () else raise Superfluous_constr_right_hand_sides
obua@16784
   173
haftmann@22951
   174
        fun mk_free y = cplex.cplexBounds (cplex.cplexNeg cplex.cplexInf, cplex.cplexLeq,
haftmann@22951
   175
                                           cplex.cplexVar y, cplex.cplexLeq,
haftmann@22951
   176
                                           cplex.cplexInf)
obua@16784
   177
haftmann@22951
   178
        val yvars = Inttab.fold (fn (i, y) => fn l => (mk_free y)::l) (!ytable) []
obua@16784
   179
haftmann@22951
   180
        val prog = cplex.cplexProg ("original", cplex.cplexMaximize (vec2sum c), list, yvars)
obua@16784
   181
    in
haftmann@22951
   182
        (prog, indexof)
obua@16784
   183
    end
obua@16784
   184
obua@16784
   185
haftmann@22951
   186
fun dual_cplexProg c A b =
obua@16784
   187
    let
haftmann@22951
   188
        fun indexof s =
haftmann@22951
   189
            if String.size s = 0 then raise (No_name s)
haftmann@22951
   190
            else case Int.fromString (String.extract(s, 1, NONE)) of
haftmann@22951
   191
                     SOME i => i | NONE => raise (No_name s)
haftmann@22951
   192
haftmann@22951
   193
        fun nameof i = "y"^(Int.toString i)
obua@16784
   194
haftmann@22951
   195
        fun split_numstr s =
haftmann@22951
   196
            if String.isPrefix "-" s then (false,String.extract(s, 1, NONE))
haftmann@22951
   197
            else if String.isPrefix "+" s then (true, String.extract(s, 1, NONE))
haftmann@22951
   198
            else (true, s)
haftmann@22951
   199
haftmann@22951
   200
        fun mk_term index s =
haftmann@22951
   201
            let
haftmann@22951
   202
                val (p, s) = split_numstr s
haftmann@22951
   203
                val prod = cplex.cplexProd (cplex.cplexNum s, cplex.cplexVar (nameof index))
haftmann@22951
   204
            in
haftmann@22951
   205
                if p then prod else cplex.cplexNeg prod
haftmann@22951
   206
            end
obua@16784
   207
haftmann@22951
   208
        fun vec2sum vector =
haftmann@22951
   209
            cplex.cplexSum (Inttab.fold (fn (index, s) => fn list => (mk_term index s)::list) vector [])
obua@16784
   210
haftmann@22951
   211
        fun mk_constr index vector c =
haftmann@22951
   212
            let
haftmann@22951
   213
                val s = case Inttab.lookup c index of SOME s => s | NONE => "0"
haftmann@22951
   214
                val (p, s) = split_numstr s
haftmann@22951
   215
                val num = if p then cplex.cplexNum s else cplex.cplexNeg (cplex.cplexNum s)
haftmann@22951
   216
            in
haftmann@22951
   217
                (NONE, cplex.cplexConstr (cplex.cplexEq, (vec2sum vector, num)))
haftmann@22951
   218
            end
obua@16784
   219
haftmann@22951
   220
        fun delete index c = Inttab.delete index c handle Inttab.UNDEF _ => c
obua@16784
   221
haftmann@22951
   222
        val (list, c) = Inttab.fold
haftmann@22951
   223
                            (fn (index, v) => fn (list, c) => ((mk_constr index v c)::list, delete index c))
haftmann@22951
   224
                            (transpose_matrix A) ([], c)
haftmann@22951
   225
        val _ = if Inttab.is_empty c then () else raise Superfluous_constr_right_hand_sides
haftmann@22951
   226
haftmann@22951
   227
        val prog = cplex.cplexProg ("dual", cplex.cplexMinimize (vec2sum b), list, [])
obua@16784
   228
    in
haftmann@22951
   229
        (prog, indexof)
obua@16784
   230
    end
obua@16784
   231
haftmann@22951
   232
fun cut_vector size v =
haftmann@21056
   233
  let
haftmann@21056
   234
    val count = ref 0;
haftmann@21056
   235
    fun app (i, s) =  if (!count < size) then
haftmann@21056
   236
        (count := !count +1 ; Inttab.update (i, s))
haftmann@21056
   237
      else I
haftmann@21056
   238
  in
haftmann@21056
   239
    Inttab.fold app v empty_vector
haftmann@21056
   240
  end
obua@16784
   241
haftmann@22951
   242
fun cut_matrix vfilter vsize m =
haftmann@21056
   243
  let
haftmann@22951
   244
    fun app (i, v) =
haftmann@21056
   245
      if is_none (Inttab.lookup vfilter i) then I
haftmann@21056
   246
      else case vsize
haftmann@21056
   247
       of NONE => Inttab.update (i, v)
haftmann@21056
   248
        | SOME s => Inttab.update (i, cut_vector s v)
haftmann@21056
   249
  in Inttab.fold app m empty_matrix end
haftmann@21056
   250
wenzelm@17412
   251
fun v_elem_at v i = Inttab.lookup v i
wenzelm@17412
   252
fun m_elem_at m i = Inttab.lookup m i
obua@16784
   253
haftmann@22951
   254
fun v_only_elem v =
obua@16784
   255
    case Inttab.min_key v of
haftmann@22951
   256
        NONE => NONE
obua@16784
   257
      | SOME vmin => (case Inttab.max_key v of
haftmann@22951
   258
                          NONE => SOME vmin
haftmann@22951
   259
                        | SOME vmax => if vmin = vmax then SOME vmin else NONE)
obua@16784
   260
haftmann@21056
   261
fun v_fold f = Inttab.fold f;
haftmann@21056
   262
fun m_fold f = Inttab.fold f;
obua@16784
   263
obua@16784
   264
end;