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