src/HOL/Matrix_LP/fspmlp.ML

Tuned Euclidean Ring instance for polynomials

(*  Title:      HOL/Matrix_LP/fspmlp.ML
    Author:     Steven Obua
*)

signature FSPMLP =
sig
    type linprog
    type vector = FloatSparseMatrixBuilder.vector
    type matrix = FloatSparseMatrixBuilder.matrix

    val y : linprog -> term
    val A : linprog -> term * term
    val b : linprog -> term
    val c : linprog -> term * term
    val r12 : linprog -> term * term

    exception Load of string

    val load : string -> int -> bool -> linprog
end

structure Fspmlp : FSPMLP =
struct

type vector = FloatSparseMatrixBuilder.vector
type matrix = FloatSparseMatrixBuilder.matrix

type linprog = term * (term * term) * term * (term * term) * (term * term)

fun y (c1, _, _, _, _) = c1
fun A (_, c2, _, _, _) = c2
fun b (_, _, c3, _, _) = c3
fun c (_, _, _, c4, _) = c4
fun r12 (_, _, _, _, c6) = c6

structure CplexFloatSparseMatrixConverter =
MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = FloatSparseMatrixBuilder);

datatype bound_type = LOWER | UPPER

fun intbound_ord ((i1: int, b1),(i2,b2)) =
    if i1 < i2 then LESS
    else if i1 = i2 then
        (if b1 = b2 then EQUAL else if b1=LOWER then LESS else GREATER)
    else GREATER

structure Inttab = Table(type key = int val ord = (rev_order o int_ord));

structure VarGraph = Table(type key = int*bound_type val ord = intbound_ord);
(* key -> (float option) * (int -> (float * (((float * float) * key) list)))) *)
(* dest_key -> (sure_bound * (row_index -> (row_bound * (((coeff_lower * coeff_upper) * src_key) list)))) *)

exception Internal of string;

fun add_row_bound g dest_key row_index row_bound =
    let
        val x =
            case VarGraph.lookup g dest_key of
                NONE => (NONE, Inttab.update (row_index, (row_bound, [])) Inttab.empty)
              | SOME (sure_bound, f) =>
                (sure_bound,
                 case Inttab.lookup f row_index of
                     NONE => Inttab.update (row_index, (row_bound, [])) f
                   | SOME _ => raise (Internal "add_row_bound"))
    in
        VarGraph.update (dest_key, x) g
    end

fun update_sure_bound g (key as (_, btype)) bound =
    let
        val x =
            case VarGraph.lookup g key of
                NONE => (SOME bound, Inttab.empty)
              | SOME (NONE, f) => (SOME bound, f)
              | SOME (SOME old_bound, f) =>
                (SOME ((case btype of
                            UPPER => Float.min
                          | LOWER => Float.max)
                           old_bound bound), f)
    in
        VarGraph.update (key, x) g
    end

fun get_sure_bound g key =
    case VarGraph.lookup g key of
        NONE => NONE
      | SOME (sure_bound, _) => sure_bound

(*fun get_row_bound g key row_index =
    case VarGraph.lookup g key of
        NONE => NONE
      | SOME (sure_bound, f) =>
        (case Inttab.lookup f row_index of
             NONE => NONE
           | SOME (row_bound, _) => (sure_bound, row_bound))*)

fun add_edge g src_key dest_key row_index coeff =
    case VarGraph.lookup g dest_key of
        NONE => raise (Internal "add_edge: dest_key not found")
      | SOME (sure_bound, f) =>
        (case Inttab.lookup f row_index of
             NONE => raise (Internal "add_edge: row_index not found")
           | SOME (row_bound, sources) =>
             VarGraph.update (dest_key, (sure_bound, Inttab.update (row_index, (row_bound, (coeff, src_key) :: sources)) f)) g)

fun split_graph g =
  let
    fun split (key, (sure_bound, _)) (r1, r2) = case sure_bound
     of NONE => (r1, r2)
      | SOME bound =>  (case key
         of (u, UPPER) => (r1, Inttab.update (u, bound) r2)
          | (u, LOWER) => (Inttab.update (u, bound) r1, r2))
  in VarGraph.fold split g (Inttab.empty, Inttab.empty) end

(* If safe is true, termination is guaranteed, but the sure bounds may be not optimal (relative to the algorithm).
   If safe is false, termination is not guaranteed, but on termination the sure bounds are optimal (relative to the algorithm) *)
fun propagate_sure_bounds safe names g =
    let
        (* returns NONE if no new sure bound could be calculated, otherwise the new sure bound is returned *)
        fun calc_sure_bound_from_sources g (key as (_, btype)) =
            let
                fun mult_upper x (lower, upper) =
                    if Float.sign x = LESS then
                        Float.mult x lower
                    else
                        Float.mult x upper

                fun mult_lower x (lower, upper) =
                    if Float.sign x = LESS then
                        Float.mult x upper
                    else
                        Float.mult x lower

                val mult_btype = case btype of UPPER => mult_upper | LOWER => mult_lower

                fun calc_sure_bound (_, (row_bound, sources)) sure_bound =
                    let
                        fun add_src_bound (coeff, src_key) sum =
                            case sum of
                                NONE => NONE
                              | SOME x =>
                                (case get_sure_bound g src_key of
                                     NONE => NONE
                                   | SOME src_sure_bound => SOME (Float.add x (mult_btype src_sure_bound coeff)))
                    in
                        case fold add_src_bound sources (SOME row_bound) of
                            NONE => sure_bound
                          | new_sure_bound as (SOME new_bound) =>
                            (case sure_bound of
                                 NONE => new_sure_bound
                               | SOME old_bound =>
                                 SOME (case btype of
                                           UPPER => Float.min old_bound new_bound
                                         | LOWER => Float.max old_bound new_bound))
                    end
            in
                case VarGraph.lookup g key of
                    NONE => NONE
                  | SOME (sure_bound, f) =>
                    let
                        val x = Inttab.fold calc_sure_bound f sure_bound
                    in
                        if x = sure_bound then NONE else x
                    end
                end

        fun propagate (key, _) (g, b) =
            case calc_sure_bound_from_sources g key of
                NONE => (g,b)
              | SOME bound => (update_sure_bound g key bound,
                               if safe then
                                   case get_sure_bound g key of
                                       NONE => true
                                     | _ => b
                               else
                                   true)

        val (g, b) = VarGraph.fold propagate g (g, false)
    in
        if b then propagate_sure_bounds safe names g else g
    end

exception Load of string;

val empty_spvec = @{term "Nil :: real spvec"};
fun cons_spvec x xs = @{term "Cons :: nat * real => real spvec => real spvec"} $ x $ xs;
val empty_spmat = @{term "Nil :: real spmat"};
fun cons_spmat x xs = @{term "Cons :: nat * real spvec => real spmat => real spmat"} $ x $ xs;

fun calcr safe_propagation xlen names prec A b =
    let
        fun test_1 (lower, upper) =
            if lower = upper then
                (if Float.eq (lower, (~1, 0)) then ~1
                 else if Float.eq (lower, (1, 0)) then 1
                 else 0)
            else 0

        fun calcr (row_index, a) g =
            let
                val b =  FloatSparseMatrixBuilder.v_elem_at b row_index
                val (_, b2) = FloatArith.approx_decstr_by_bin prec (case b of NONE => "0" | SOME b => b)
                val approx_a = FloatSparseMatrixBuilder.v_fold (fn (i, s) => fn l =>
                                                                   (i, FloatArith.approx_decstr_by_bin prec s)::l) a []

                fun fold_dest_nodes (dest_index, dest_value) g =
                    let
                        val dest_test = test_1 dest_value
                    in
                        if dest_test = 0 then
                            g
                        else let
                                val (dest_key as (_, dest_btype), row_bound) =
                                    if dest_test = ~1 then
                                        ((dest_index, LOWER), Float.neg b2)
                                    else
                                        ((dest_index, UPPER), b2)

                                fun fold_src_nodes (src_index, src_value as (src_lower, src_upper)) g =
                                    if src_index = dest_index then g
                                    else
                                        let
                                            val coeff = case dest_btype of
                                                            UPPER => (Float.neg src_upper, Float.neg src_lower)
                                                          | LOWER => src_value
                                        in
                                            if Float.sign src_lower = LESS then
                                                add_edge g (src_index, UPPER) dest_key row_index coeff
                                            else
                                                add_edge g (src_index, LOWER) dest_key row_index coeff
                                        end
                            in
                                fold fold_src_nodes approx_a (add_row_bound g dest_key row_index row_bound)
                            end
                    end
            in
                case approx_a of
                    [] => g
                  | [(u, a)] =>
                    let
                        val atest = test_1 a
                    in
                        if atest = ~1 then
                            update_sure_bound g (u, LOWER) (Float.neg b2)
                        else if atest = 1 then
                            update_sure_bound g (u, UPPER) b2
                        else
                            g
                    end
                  | _ => fold fold_dest_nodes approx_a g
            end

        val g = FloatSparseMatrixBuilder.m_fold calcr A VarGraph.empty

        val g = propagate_sure_bounds safe_propagation names g

        val (r1, r2) = split_graph g

        fun add_row_entry m index f vname value =
            let
                val v = (case value of 
                             SOME value => FloatSparseMatrixBuilder.mk_spvec_entry 0 value
                           | NONE => FloatSparseMatrixBuilder.mk_spvec_entry' 0 (f $ (Var ((vname,0), HOLogic.realT))))
                val vec = cons_spvec v empty_spvec
            in
                cons_spmat (FloatSparseMatrixBuilder.mk_spmat_entry index vec) m
            end

        fun abs_estimate i r1 r2 =
            if i = 0 then
                let val e = empty_spmat in (e, e) end
            else
                let
                    val index = xlen-i
                    val (r12_1, r12_2) = abs_estimate (i-1) r1 r2
                    val b1 = Inttab.lookup r1 index
                    val b2 = Inttab.lookup r2 index
                in
                    (add_row_entry r12_1 index @{term "lbound :: real => real"} ((names index)^"l") b1, 
                     add_row_entry r12_2 index @{term "ubound :: real => real"} ((names index)^"u") b2)
                end

        val (r1, r2) = abs_estimate xlen r1 r2

    in
        (r1, r2)
    end

fun load filename prec safe_propagation =
    let
        val prog = Cplex.load_cplexFile filename
        val prog = Cplex.elim_nonfree_bounds prog
        val prog = Cplex.relax_strict_ineqs prog
        val (maximize, c, A, b, (xlen, names, _)) = CplexFloatSparseMatrixConverter.convert_prog prog                       
        val (r1, r2) = calcr safe_propagation xlen names prec A b
        val _ = if maximize then () else raise Load "sorry, cannot handle minimization problems"
        val (dualprog, indexof) = FloatSparseMatrixBuilder.dual_cplexProg c A b
        val results = Cplex.solve dualprog
        val (_, v) = CplexFloatSparseMatrixConverter.convert_results results indexof
        (*val A = FloatSparseMatrixBuilder.cut_matrix v NONE A*)
        fun id x = x
        val v = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 v
        val b = FloatSparseMatrixBuilder.transpose_matrix (FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 b)
        val c = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 c
        val (y1, _) = FloatSparseMatrixBuilder.approx_matrix prec Float.positive_part v
        val A = FloatSparseMatrixBuilder.approx_matrix prec id A
        val (_,b2) = FloatSparseMatrixBuilder.approx_matrix prec id b
        val c = FloatSparseMatrixBuilder.approx_matrix prec id c
    in
        (y1, A, b2, c, (r1, r2))
    end handle CplexFloatSparseMatrixConverter.Converter s => (raise (Load ("Converter: "^s)))

end