src/HOL/Matrix/fspmlp.ML
changeset 15178 5f621aa35c25
child 15531 08c8dad8e399
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Matrix/fspmlp.ML	Fri Sep 03 17:10:36 2004 +0200
@@ -0,0 +1,303 @@
+signature FSPMLP = 
+sig
+    type linprog
+
+    val y : linprog -> cterm
+    val A : linprog -> cterm * cterm
+    val b : linprog -> cterm
+    val c : linprog -> cterm * cterm
+    val r : linprog -> cterm
+
+    exception Load of string
+		       
+    val load : string -> int -> bool -> linprog
+end
+
+structure fspmlp : FSPMLP = 
+struct
+
+type linprog = cterm * (cterm * cterm) * cterm * (cterm * cterm) * cterm 
+
+fun y (c1, c2, c3, c4, c5) = c1
+fun A (c1, c2, c3, c4, c5) = c2
+fun b (c1, c2, c3, c4, c5) = c3
+fun c (c1, c2, c3, c4, c5) = c4
+fun r (c1, c2, c3, c4, c5) = c5
+
+structure CplexFloatSparseMatrixConverter = 
+MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = FloatSparseMatrixBuilder);
+
+datatype bound_type = LOWER | UPPER
+
+fun intbound_ord ((i1, 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 = TableFun(type key = int val ord = (rev_order o int_ord));
+
+structure VarGraph = TableFun(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 => FloatArith.min 
+			  | LOWER => FloatArith.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 ((r1, r2), (key, (sure_bound, _))) = 
+	    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.foldl split ((Inttab.empty, Inttab.empty), g)
+    end
+
+fun it2list t = 
+    let
+	fun tolist (l, a) = a::l
+    in
+	Inttab.foldl tolist ([], t)
+    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 FloatArith.is_negative x then
+			FloatArith.mul x lower
+		    else
+			FloatArith.mul x upper
+			
+		fun mult_lower x (lower, upper) = 
+		    if FloatArith.is_negative x then
+			FloatArith.mul x upper
+		    else
+			FloatArith.mul x lower
+
+		val mult_btype = case btype of UPPER => mult_upper | LOWER => mult_lower
+
+		fun calc_sure_bound (sure_bound, (row_index, (row_bound, sources))) = 
+		    let
+			fun add_src_bound (sum, (coeff, src_key)) = 
+			    case sum of 
+				None => None
+			      | Some x => 
+				(case get_sure_bound g src_key of
+				     None => None
+				   | Some src_sure_bound => Some (FloatArith.add x (mult_btype src_sure_bound coeff)))
+		    in
+			case foldl add_src_bound (Some row_bound, sources) 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 => FloatArith.min old_bound new_bound
+					 | LOWER => FloatArith.max old_bound new_bound))				 
+		    end		
+	    in
+		case VarGraph.lookup (g, key) of
+		    None => None
+		  | Some (sure_bound, f) =>
+		    let
+			val x = Inttab.foldl calc_sure_bound (sure_bound, f) 
+		    in
+			if x = sure_bound then None else x
+		    end		
+    	    end
+
+	fun propagate ((g, b), (key, _)) = 
+	    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.foldl propagate ((g, false), g) 
+    in
+	if b then propagate_sure_bounds safe names g else g	
+    end	    
+    		
+exception Load of string;
+
+fun calcr safe_propagation xlen names prec A b = 
+    let
+	val empty = Inttab.empty
+
+	fun instab t i x = Inttab.update ((i,x), t)
+
+	fun isnegstr x = String.isPrefix "-" x
+	fun negstr x = if isnegstr x then String.extract (x, 1, NONE) else "-"^x
+
+	fun test_1 (lower, upper) = 
+	    if lower = upper then
+		(if FloatArith.is_equal lower (IntInf.fromInt ~1, FloatArith.izero) then ~1 
+		 else if FloatArith.is_equal lower (IntInf.fromInt 1, FloatArith.izero) then 1
+		 else 0)
+	    else 0	
+
+	fun calcr (g, (row_index, a)) = 
+	    let				
+		val b =  FloatSparseMatrixBuilder.v_elem_at b row_index
+		val (_, b2) = ExactFloatingPoint.approx_decstr_by_bin prec (case b of None => "0" | Some b => b)
+		val approx_a = FloatSparseMatrixBuilder.v_fold (fn (l, (i,s)) => 
+								   (i, ExactFloatingPoint.approx_decstr_by_bin prec s)::l) [] a
+			       
+		fun fold_dest_nodes (g, (dest_index, dest_value)) = 
+		    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), FloatArith.neg b2)
+				    else
+					((dest_index, UPPER), b2)
+					
+				fun fold_src_nodes (g, (src_index, src_value as (src_lower, src_upper))) = 
+				    if src_index = dest_index then g
+				    else
+					let
+					    val coeff = case dest_btype of 
+							    UPPER => (FloatArith.neg src_upper, FloatArith.neg src_lower)
+							  | LOWER => src_value
+					in
+					    if FloatArith.is_negative src_lower 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	    
+				foldl fold_src_nodes ((add_row_bound g dest_key row_index row_bound), approx_a)
+			    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) (FloatArith.neg b2)
+			else if atest = 1 then
+			    update_sure_bound g (u, UPPER) b2
+			else
+			    g
+		    end
+		  | _ => foldl fold_dest_nodes (g, approx_a)
+	    end
+	
+	val g = FloatSparseMatrixBuilder.m_fold calcr VarGraph.empty A
+	val g = propagate_sure_bounds safe_propagation names g
+
+	val (r1, r2) = split_graph g
+	
+	fun abs_estimate i r1 r2 = 
+	    if i = 0 then FloatSparseMatrixBuilder.empty_spmat
+	    else
+		let
+		    val index = xlen-i
+		    val r = abs_estimate (i-1) r1 r2 
+		    val b1 = case Inttab.lookup (r1, index) of None => raise (Load ("x-value not bounded from below: "^(names index))) | Some x => x
+		    val b2 = case Inttab.lookup (r2, index) of None => raise (Load ("x-value not bounded from above: "^(names index))) | Some x => x
+		    val abs_max = FloatArith.max (FloatArith.neg (FloatArith.negative_part b1)) (FloatArith.positive_part b2)    
+		    val vec = FloatSparseMatrixBuilder.cons_spvec (FloatSparseMatrixBuilder.mk_spvec_entry 0 abs_max) FloatSparseMatrixBuilder.empty_spvec
+		in
+		    FloatSparseMatrixBuilder.cons_spmat (FloatSparseMatrixBuilder.mk_spmat_entry index vec) r
+		end		    		    
+    in
+	FloatSparseMatrixBuilder.sign_term (abs_estimate xlen 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 r = 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 (optimal,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 FloatArith.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, r)
+    end handle CplexFloatSparseMatrixConverter.Converter s => (raise (Load ("Converter: "^s)))
+
+end
\ No newline at end of file