src/HOL/Matrix/cplex/Cplex_tools.ML
author haftmann
Fri, 20 Oct 2006 10:44:33 +0200
changeset 21056 2cfe839e8d58
parent 17521 0f1c48de39f5
child 21858 05f57309170c
permissions -rw-r--r--
Symtab.foldl replaced by Symtab.fold

(*  Title:      HOL/Matrix/cplex/Cplex_tools.ML
    ID:         $Id$
    Author:     Steven Obua
*)

signature CPLEX =
sig
    
    datatype cplexTerm = cplexVar of string | cplexNum of string | cplexInf 
                       | cplexNeg of cplexTerm 
                       | cplexProd of cplexTerm * cplexTerm 
                       | cplexSum of (cplexTerm list) 
    
    datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq 
								  
    datatype cplexGoal = cplexMinimize of cplexTerm 
		       | cplexMaximize of cplexTerm
					  
    datatype cplexConstr = cplexConstr of cplexComp * 
					  (cplexTerm * cplexTerm)
					  
    datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm 
					  * cplexComp * cplexTerm
			 | cplexBound of cplexTerm * cplexComp * cplexTerm
					 
    datatype cplexProg = cplexProg of string 
				      * cplexGoal 
				      * ((string option * cplexConstr) 
					     list) 
				      * cplexBounds list
				      
    datatype cplexResult = Unbounded 
			 | Infeasible 
			 | Undefined
			 | Optimal of string * 
				      (((* name *) string * 
					(* value *) string) list)

    datatype cplexSolver = SOLVER_DEFAULT | SOLVER_CPLEX | SOLVER_GLPK
						     
    exception Load_cplexFile of string
    exception Load_cplexResult of string
    exception Save_cplexFile of string
    exception Execute of string

    val load_cplexFile : string -> cplexProg

    val save_cplexFile : string -> cplexProg -> unit 

    val elim_nonfree_bounds : cplexProg -> cplexProg

    val relax_strict_ineqs : cplexProg -> cplexProg

    val is_normed_cplexProg : cplexProg -> bool

    val get_solver : unit -> cplexSolver
    val set_solver : cplexSolver -> unit
    val solve : cplexProg -> cplexResult
end;

structure Cplex  : CPLEX =
struct

datatype cplexSolver = SOLVER_DEFAULT | SOLVER_CPLEX | SOLVER_GLPK

val cplexsolver = ref SOLVER_DEFAULT;
fun get_solver () = !cplexsolver;
fun set_solver s = (cplexsolver := s);

exception Load_cplexFile of string;
exception Load_cplexResult of string;
exception Save_cplexFile of string;
	  
datatype cplexTerm = cplexVar of string 
		   | cplexNum of string 
		   | cplexInf 
                   | cplexNeg of cplexTerm 
                   | cplexProd of cplexTerm * cplexTerm 
                   | cplexSum of (cplexTerm list) 
datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq 
datatype cplexGoal = cplexMinimize of cplexTerm | cplexMaximize of cplexTerm
datatype cplexConstr = cplexConstr of cplexComp * (cplexTerm * cplexTerm)
datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm 
				      * cplexComp * cplexTerm
                     | cplexBound of cplexTerm * cplexComp * cplexTerm 
datatype cplexProg = cplexProg of string 
				  * cplexGoal 
				  * ((string option * cplexConstr) list) 
				  * cplexBounds list
				  
fun rev_cmp cplexLe = cplexGe
  | rev_cmp cplexLeq = cplexGeq
  | rev_cmp cplexGe = cplexLe
  | rev_cmp cplexGeq = cplexLeq
  | rev_cmp cplexEq = cplexEq

fun the NONE = raise (Load_cplexFile "SOME expected")
  | the (SOME x) = x; 
    
fun modulo_signed is_something (cplexNeg u) = is_something u
  | modulo_signed is_something u = is_something u

fun is_Num (cplexNum _) = true
  | is_Num _ = false

fun is_Inf cplexInf = true
  | is_Inf _ = false

fun is_Var (cplexVar _) = true
  | is_Var _ = false

fun is_Neg (cplexNeg x ) = true
  | is_Neg _ = false

fun is_normed_Prod (cplexProd (t1, t2)) = 
    (is_Num t1) andalso (is_Var t2)
  | is_normed_Prod x = is_Var x

fun is_normed_Sum (cplexSum ts) = 
    (ts <> []) andalso forall (modulo_signed is_normed_Prod) ts
  | is_normed_Sum x = modulo_signed is_normed_Prod x

fun is_normed_Constr (cplexConstr (c, (t1, t2))) = 
    (is_normed_Sum t1) andalso (modulo_signed is_Num t2)

fun is_Num_or_Inf x = is_Inf x orelse is_Num x

fun is_normed_Bounds (cplexBounds (t1, c1, t2, c2, t3)) =
    (c1 = cplexLe orelse c1 = cplexLeq) andalso 
    (c2 = cplexLe orelse c2 = cplexLeq) andalso
    is_Var t2 andalso
    modulo_signed is_Num_or_Inf t1 andalso
    modulo_signed is_Num_or_Inf t3
  | is_normed_Bounds (cplexBound (t1, c, t2)) =
    (is_Var t1 andalso (modulo_signed is_Num_or_Inf t2)) 
    orelse
    (c <> cplexEq andalso 
     is_Var t2 andalso (modulo_signed is_Num_or_Inf t1))
	             		
fun term_of_goal (cplexMinimize x) = x
  | term_of_goal (cplexMaximize x) = x

fun is_normed_cplexProg (cplexProg (name, goal, constraints, bounds)) = 
    is_normed_Sum (term_of_goal goal) andalso
    forall (fn (_,x) => is_normed_Constr x) constraints andalso
    forall is_normed_Bounds bounds

fun is_NL s = s = "\n"

fun is_blank s = forall (fn c => c <> #"\n" andalso Char.isSpace c) (String.explode s)
		 
fun is_num a = 
    let 
	val b = String.explode a 
	fun num4 cs = forall Char.isDigit cs		      
	fun num3 [] = true
	  | num3 (ds as (c::cs)) = 
	    if c = #"+" orelse c = #"-" then
		num4 cs
	    else 
		num4 ds		    
	fun num2 [] = true
	  | num2 (c::cs) =
	    if c = #"e" orelse c = #"E" then num3 cs
	    else (Char.isDigit c) andalso num2 cs			 
	fun num1 [] = true
	  | num1 (c::cs) = 
	    if c = #"." then num2 cs 
	    else if c = #"e" orelse c = #"E" then num3 cs 
	    else (Char.isDigit c) andalso num1 cs
	fun num [] = true
	  | num (c::cs) = 
	    if c = #"." then num2 cs 
	    else (Char.isDigit c) andalso num1 cs
    in
	num b
    end  

fun is_delimiter s = s = "+" orelse s = "-" orelse s = ":"
		    
fun is_cmp s = s = "<" orelse s = ">" orelse s = "<=" 
		     orelse s = ">=" orelse s = "="

fun is_symbol a = 
    let
	val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~"
	fun is_symbol_char c = Char.isAlphaNum c orelse 
			       exists (fn d => d=c) symbol_char
	fun is_symbol_start c = is_symbol_char c andalso 
				not (Char.isDigit c) andalso 
				not (c= #".")
	val b = String.explode a
    in
	b <> [] andalso is_symbol_start (hd b) andalso 
	forall is_symbol_char b
    end

fun to_upper s = String.implode (map Char.toUpper (String.explode s))

fun keyword x = 
    let 	
	val a = to_upper x 
    in
	if a = "BOUNDS" orelse a = "BOUND" then
	    SOME "BOUNDS"
	else if a = "MINIMIZE" orelse a = "MINIMUM" orelse a = "MIN" then
	    SOME "MINIMIZE"
	else if a = "MAXIMIZE" orelse a = "MAXIMUM" orelse a = "MAX" then
	    SOME "MAXIMIZE"
	else if a = "ST" orelse a = "S.T." orelse a = "ST." then
	    SOME "ST"
	else if a = "FREE" orelse a = "END" then
	    SOME a
	else if a = "GENERAL" orelse a = "GENERALS" orelse a = "GEN" then
	    SOME "GENERAL"
	else if a = "INTEGER" orelse a = "INTEGERS" orelse a = "INT" then
	    SOME "INTEGER"
	else if a = "BINARY" orelse a = "BINARIES" orelse a = "BIN" then
	    SOME "BINARY"
	else if a = "INF" orelse a = "INFINITY" then
	    SOME "INF"
	else								   
	    NONE
    end
	        
val TOKEN_ERROR = ~1
val TOKEN_BLANK = 0
val TOKEN_NUM = 1
val TOKEN_DELIMITER = 2
val TOKEN_SYMBOL = 3
val TOKEN_LABEL = 4
val TOKEN_CMP = 5
val TOKEN_KEYWORD = 6
val TOKEN_NL = 7
		      
(* tokenize takes a list of chars as argument and returns a list of
   int * string pairs, each string representing a "cplex token", 
   and each int being one of TOKEN_NUM, TOKEN_DELIMITER, TOKEN_CMP
   or TOKEN_SYMBOL *)
fun tokenize s = 
    let
	val flist = [(is_NL, TOKEN_NL), 
		     (is_blank, TOKEN_BLANK), 
		     (is_num, TOKEN_NUM),  
                     (is_delimiter, TOKEN_DELIMITER), 
		     (is_cmp, TOKEN_CMP), 
		     (is_symbol, TOKEN_SYMBOL)]  
	fun match_helper [] s = (fn x => false, TOKEN_ERROR)
	  | match_helper (f::fs) s = 
	    if ((fst f) s) then f else match_helper fs s   
	fun match s = match_helper flist s
	fun tok s = 
	    if s = "" then [] else       
	    let 
		val h = String.substring (s,0,1) 
		val (f, j) = match h
		fun len i = 
		    if size s = i then i 
		    else if f (String.substring (s,0,i+1)) then 
			len (i+1) 
		    else i 
	    in
		if j < 0 then 
		    (if h = "\\" then [] 
		     else raise (Load_cplexFile ("token expected, found: "
						 ^s)))
		else 
		    let 
			val l = len 1 
			val u = String.substring (s,0,l)
			val v = String.extract (s,l,NONE)
		    in
			if j = 0 then tok v else (j, u) :: tok v
		    end
	    end                
    in
	tok s     
    end

exception Tokenize of string;

fun tokenize_general flist s = 
    let
	fun match_helper [] s = raise (Tokenize s)
	  | match_helper (f::fs) s = 
	    if ((fst f) s) then f else match_helper fs s   
	fun match s = match_helper flist s
	fun tok s = 
	    if s = "" then [] else       
	    let 
		val h = String.substring (s,0,1) 
		val (f, j) = match h
		fun len i = 
		    if size s = i then i 
		    else if f (String.substring (s,0,i+1)) then 
			len (i+1) 
		    else i 
		val l = len 1
	    in
		(j, String.substring (s,0,l)) :: tok (String.extract (s,l,NONE))
	    end
    in
	tok s     
    end

fun load_cplexFile name = 
    let 
	val f = TextIO.openIn name  
        val ignore_NL = ref true  
	val rest = ref []

	fun is_symbol s c = (fst c) = TOKEN_SYMBOL andalso (to_upper (snd c)) = s 
		   
	fun readToken_helper () =
	    if length (!rest) > 0 then  
		let val u = hd (!rest) in 
		    (
		     rest := tl (!rest); 
		     SOME u
		    ) 
		end
	    else 
		let val s = TextIO.inputLine f in
		    if s = "" then NONE else
		    let val t = tokenize s in
			if (length t >= 2 andalso 
			    snd(hd (tl t)) = ":") 
			then 
			    rest := (TOKEN_LABEL, snd (hd t)) :: (tl (tl t))
			else if (length t >= 2) andalso is_symbol "SUBJECT" (hd (t)) 
				andalso is_symbol "TO" (hd (tl t)) 
			then
			    rest := (TOKEN_SYMBOL, "ST") :: (tl (tl t))
			else
			    rest := t;
			readToken_helper ()  
		    end
		end           
	
	fun readToken_helper2 () = 
		let val c = readToken_helper () in
		    if c = NONE then NONE
                    else if !ignore_NL andalso fst (the c) = TOKEN_NL then
			readToken_helper2 ()
		    else if fst (the c) = TOKEN_SYMBOL 
			    andalso keyword (snd (the c)) <> NONE
		    then SOME (TOKEN_KEYWORD, the (keyword (snd (the c))))
		    else c
		end
		
	fun readToken () = readToken_helper2 ()
	    			   
	fun pushToken a = rest := (a::(!rest))
	    
	fun is_value token = 
	    fst token = TOKEN_NUM orelse (fst token = TOKEN_KEYWORD 
					  andalso snd token = "INF")	

        fun get_value token = 
	    if fst token = TOKEN_NUM then
		cplexNum (snd token)
	    else if fst token = TOKEN_KEYWORD andalso snd token = "INF" 
	    then
		cplexInf
	    else
		raise (Load_cplexFile "num expected")
					    
	fun readTerm_Product only_num = 
	    let val c = readToken () in
		if c = NONE then NONE
		else if fst (the c) = TOKEN_SYMBOL 
		then (
		    if only_num then (pushToken (the c); NONE) 
		    else SOME (cplexVar (snd (the c)))
		    ) 
		else if only_num andalso is_value (the c) then
		    SOME (get_value (the c))
		else if is_value (the c) then  
		    let val t1 = get_value (the c) 
			val d = readToken () 
		    in
			if d = NONE then SOME t1 
			else if fst (the d) = TOKEN_SYMBOL then 
			    SOME (cplexProd (t1, cplexVar (snd (the d))))
			else
			    (pushToken (the d); SOME t1)
		    end
		else (pushToken (the c); NONE)
	    end
	    
	fun readTerm_Signed only_signed only_num = 
	    let
		val c = readToken ()
	    in
		if c = NONE then NONE 
		else 
		    let val d = the c in 
			if d = (TOKEN_DELIMITER, "+") then 
			    readTerm_Product only_num
			 else if d = (TOKEN_DELIMITER, "-") then 
			     SOME (cplexNeg (the (readTerm_Product 
						      only_num)))
			 else (pushToken d; 
			       if only_signed then NONE 
			       else readTerm_Product only_num)
		    end
	    end
    
	fun readTerm_Sum first_signed = 
	    let val c = readTerm_Signed first_signed false in
		if c = NONE then [] else (the c)::(readTerm_Sum true)
	    end
	    
	fun readTerm () = 
	    let val c = readTerm_Sum false in
		if c = [] then NONE 
		else if tl c = [] then SOME (hd c)
		else SOME (cplexSum c) 
	    end
	    
	fun readLabeledTerm () = 
	    let val c = readToken () in
		if c = NONE then (NONE, NONE) 
		else if fst (the c) = TOKEN_LABEL then 
		    let val t = readTerm () in
			if t = NONE then 
			    raise (Load_cplexFile ("term after label "^
						   (snd (the c))^
						   " expected"))
			else (SOME (snd (the c)), t)
		    end
		else (pushToken (the c); (NONE, readTerm ()))
	    end
            
	fun readGoal () = 
	    let 
		val g = readToken ()
	    in
    		if g = SOME (TOKEN_KEYWORD, "MAXIMIZE") then 
		    cplexMaximize (the (snd (readLabeledTerm ()))) 
		else if g = SOME (TOKEN_KEYWORD, "MINIMIZE") then 
		    cplexMinimize (the (snd (readLabeledTerm ())))
		else raise (Load_cplexFile "MAXIMIZE or MINIMIZE expected")
	    end
	    
	fun str2cmp b = 
	    (case b of 
		 "<" => cplexLe 
	       | "<=" => cplexLeq 
	       | ">" => cplexGe 
	       | ">=" => cplexGeq 
               | "=" => cplexEq
	       | _ => raise (Load_cplexFile (b^" is no TOKEN_CMP")))
			    
	fun readConstraint () = 
            let 
		val t = readLabeledTerm () 
		fun make_constraint b t1 t2 =
                    cplexConstr 
			(str2cmp b,
			 (t1, t2))		
	    in
		if snd t = NONE then NONE
		else 
		    let val c = readToken () in
			if c = NONE orelse fst (the c) <> TOKEN_CMP 
			then raise (Load_cplexFile "TOKEN_CMP expected")
			else 
			    let val n = readTerm_Signed false true in
				if n = NONE then 
				    raise (Load_cplexFile "num expected")
				else
				    SOME (fst t,
					  make_constraint (snd (the c)) 
							  (the (snd t)) 
							  (the n))
			    end
		    end            
	    end
        
        fun readST () = 
	    let 
		fun readbody () = 
		    let val t = readConstraint () in
			if t = NONE then []
			else if (is_normed_Constr (snd (the t))) then
			    (the t)::(readbody ())
			else if (fst (the t) <> NONE) then 
			    raise (Load_cplexFile 
				       ("constraint '"^(the (fst (the t)))^ 
					"'is not normed"))
			else
			    raise (Load_cplexFile 
				       "constraint is not normed")
		    end
	    in    		
		if readToken () = SOME (TOKEN_KEYWORD, "ST") 
		then 
		    readbody ()
		else 
		    raise (Load_cplexFile "ST expected")
	    end
	    
	fun readCmp () = 
	    let val c = readToken () in
		if c = NONE then NONE
		else if fst (the c) = TOKEN_CMP then 
		    SOME (str2cmp (snd (the c)))
		else (pushToken (the c); NONE)
	    end
	    
	fun skip_NL () =
	    let val c = readToken () in
		if c <> NONE andalso fst (the c) = TOKEN_NL then
		    skip_NL ()
		else
		    (pushToken (the c); ())
	    end

	fun is_var (cplexVar _) = true
	  | is_var _ = false
	    
	fun make_bounds c t1 t2 = 
	    cplexBound (t1, c, t2)

	fun readBound () = 	    
	    let 		
		val _ = skip_NL ()
		val t1 = readTerm ()
	    in
		if t1 = NONE then NONE 
		else 
		    let 
			val c1 = readCmp () 
		    in
			if c1 = NONE then
			    let 
				val c = readToken () 
			    in
				if c = SOME (TOKEN_KEYWORD, "FREE") then
				    SOME (
				    cplexBounds (cplexNeg cplexInf,
						 cplexLeq,
						 the t1,
						 cplexLeq,
						 cplexInf))
				else
				    raise (Load_cplexFile "FREE expected")
			    end
			else
			    let 
				val t2 = readTerm () 
			    in
				if t2 = NONE then
				    raise (Load_cplexFile "term expected")
				else
				    let val c2 = readCmp () in
					if c2 = NONE then
					    SOME (make_bounds (the c1) 
							      (the t1)
							      (the t2))
					else
					    SOME (
					    cplexBounds (the t1,
							 the c1,
							 the t2,
							 the c2,
							 the (readTerm())))
				    end
			    end
		    end
	    end
	    
	fun readBounds () =
	    let 
		fun makestring b = "?"
		fun readbody () = 
		    let 
			val b = readBound ()
		    in
			if b = NONE then []
			else if (is_normed_Bounds (the b)) then
			    (the b)::(readbody())
			else (
			    raise (Load_cplexFile 
				       ("bounds are not normed in: "^
					(makestring (the b)))))
		    end
	    in
		if readToken () = SOME (TOKEN_KEYWORD, "BOUNDS") then 
		    readbody ()
		else raise (Load_cplexFile "BOUNDS expected")
	    end

        fun readEnd () = 
	    if readToken () = SOME (TOKEN_KEYWORD, "END") then ()
	    else raise (Load_cplexFile "END expected")
		    
	val result_Goal = readGoal ()  
	val result_ST = readST ()  
	val _ =	ignore_NL := false 
        val result_Bounds = readBounds ()
        val _ = ignore_NL := true
        val _ = readEnd ()
	val _ = TextIO.closeIn f
    in
	cplexProg (name, result_Goal, result_ST, result_Bounds)
    end

fun save_cplexFile filename (cplexProg (name, goal, constraints, bounds)) =
    let
	val f = TextIO.openOut filename
	
	fun basic_write s = TextIO.output(f, s)

	val linebuf = ref ""
	fun buf_flushline s = 
	    (basic_write (!linebuf); 
	     basic_write "\n";
	     linebuf := s)
	fun buf_add s = linebuf := (!linebuf) ^ s

	fun write s = 
	    if (String.size s) + (String.size (!linebuf)) >= 250 then 
		buf_flushline ("    "^s)
	    else 
		buf_add s

        fun writeln s = (buf_add s; buf_flushline "")
        
	fun write_term (cplexVar x) = write x
	  | write_term (cplexNum x) = write x
	  | write_term cplexInf = write "inf"
	  | write_term (cplexProd (cplexNum "1", b)) = write_term b
	  | write_term (cplexProd (a, b)) = 
	    (write_term a; write " "; write_term b)
          | write_term (cplexNeg x) = (write " - "; write_term x)
          | write_term (cplexSum ts) = write_terms ts
	and write_terms [] = ()
	  | write_terms (t::ts) = 
	    (if (not (is_Neg t)) then write " + " else (); 
	     write_term t; write_terms ts)

	fun write_goal (cplexMaximize term) = 
	    (writeln "MAXIMIZE"; write_term term; writeln "")
	  | write_goal (cplexMinimize term) = 
	    (writeln "MINIMIZE"; write_term term; writeln "")

	fun write_cmp cplexLe = write "<"
	  | write_cmp cplexLeq = write "<="
	  | write_cmp cplexEq = write "="
	  | write_cmp cplexGe = write ">"
	  | write_cmp cplexGeq = write ">="
        
	fun write_constr (cplexConstr (cmp, (a,b))) = 
	    (write_term a;
	     write " ";
	     write_cmp cmp;
	     write " ";
	     write_term b)
        
	fun write_constraints [] = ()
	  | write_constraints (c::cs) = 
	    (if (fst c <> NONE) 
	     then 
		 (write (the (fst c)); write ": ") 
	     else 
		 ();
	     write_constr (snd c);
	     writeln "";
	     write_constraints cs)

	fun write_bounds [] = ()
	  | write_bounds ((cplexBounds (t1,c1,t2,c2,t3))::bs) =
	    ((if t1 = cplexNeg cplexInf andalso t3 = cplexInf
		 andalso (c1 = cplexLeq orelse c1 = cplexLe) 
		 andalso (c2 = cplexLeq orelse c2 = cplexLe) 
	      then
		  (write_term t2; write " free")
	      else		 		 
		  (write_term t1; write " "; write_cmp c1; write " ";
		   write_term t2; write " "; write_cmp c2; write " ";
		   write_term t3)
	     ); writeln ""; write_bounds bs)
	  | write_bounds ((cplexBound (t1, c, t2)) :: bs) = 
	    (write_term t1; write " "; 
	     write_cmp c; write " ";
	     write_term t2; writeln ""; write_bounds bs)
	     
	val _ = write_goal goal
        val _ = (writeln ""; writeln "ST")
	val _ = write_constraints constraints
        val _ = (writeln ""; writeln "BOUNDS")
	val _ = write_bounds bounds
        val _ = (writeln ""; writeln "END") 
        val _ = TextIO.closeOut f
    in
	()
    end

fun norm_Constr (constr as cplexConstr (c, (t1, t2))) = 
    if not (modulo_signed is_Num t2) andalso  
       modulo_signed is_Num t1 
    then
	[cplexConstr (rev_cmp c, (t2, t1))]
    else if (c = cplexLe orelse c = cplexLeq) andalso
	    (t1 = (cplexNeg cplexInf) orelse t2 = cplexInf)	    
    then 
	[]
    else if (c = cplexGe orelse c = cplexGeq) andalso
	    (t1 = cplexInf orelse t2 = cplexNeg cplexInf)
    then
	[]
    else
	[constr]

fun bound2constr (cplexBounds (t1,c1,t2,c2,t3)) =
    (norm_Constr(cplexConstr (c1, (t1, t2))))
    @ (norm_Constr(cplexConstr (c2, (t2, t3))))
  | bound2constr (cplexBound (t1, cplexEq, t2)) =
    (norm_Constr(cplexConstr (cplexLeq, (t1, t2))))
    @ (norm_Constr(cplexConstr (cplexLeq, (t2, t1))))
  | bound2constr (cplexBound (t1, c1, t2)) = 
    norm_Constr(cplexConstr (c1, (t1,t2)))

val emptyset = Symtab.empty

fun singleton v = Symtab.update (v, ()) emptyset

fun merge a b = Symtab.merge (op =) (a, b)

fun mergemap f ts = fold (fn x => fn table => merge table (f x)) ts Symtab.empty

fun diff a b = Symtab.fold (fn (k, v) => fn a => 
			 (Symtab.delete k a) handle UNDEF => a) 
		     b a
		    
fun collect_vars (cplexVar v) = singleton v
  | collect_vars (cplexNeg t) = collect_vars t
  | collect_vars (cplexProd (t1, t2)) = 
    merge (collect_vars t1) (collect_vars t2)
  | collect_vars (cplexSum ts) = mergemap collect_vars ts
  | collect_vars _ = emptyset

(* Eliminates all nonfree bounds from the linear program and produces an
   equivalent program with only free bounds 
   IF for the input program P holds: is_normed_cplexProg P *)
fun elim_nonfree_bounds (cplexProg (name, goal, constraints, bounds)) =
    let
	fun collect_constr_vars (_, cplexConstr (c, (t1,_))) =
	    (collect_vars t1)   
	
	val cvars = merge (collect_vars (term_of_goal goal)) 
			  (mergemap collect_constr_vars constraints)         
		    
	fun collect_lower_bounded_vars 
	    (cplexBounds (t1, c1, cplexVar v, c2, t3)) =
	    singleton v
	  |  collect_lower_bounded_vars 
		 (cplexBound (_, cplexLe, cplexVar v)) =
	     singleton v
	  |  collect_lower_bounded_vars 
		 (cplexBound (_, cplexLeq, cplexVar v)) = 
	     singleton v
	  |  collect_lower_bounded_vars 
		 (cplexBound (cplexVar v, cplexGe,_)) = 
	     singleton v
	  |  collect_lower_bounded_vars 
		 (cplexBound (cplexVar v, cplexGeq, _)) = 
	     singleton v
	  | collect_lower_bounded_vars
	    (cplexBound (cplexVar v, cplexEq, _)) = 
	    singleton v
	  |  collect_lower_bounded_vars _ = emptyset
	
	val lvars = mergemap collect_lower_bounded_vars bounds        
	val positive_vars = diff cvars lvars			   
	val zero = cplexNum "0"
	
	fun make_pos_constr v = 
	    (NONE, cplexConstr (cplexGeq, ((cplexVar v), zero)))
	
	fun make_free_bound v = 
	    cplexBounds (cplexNeg cplexInf, cplexLeq, 
			 cplexVar v, cplexLeq,
			 cplexInf)
	
	val pos_constrs = rev (Symtab.fold
			      (fn (k, v) => cons (make_pos_constr k))
			      positive_vars [])
        val bound_constrs = map (pair NONE)
				(maps bound2constr bounds)
	val constraints' = constraints @ pos_constrs @ bound_constrs
	val bounds' = rev (Symtab.fold (fn (v, _) => cons (make_free_bound v)) cvars []);
    in
	cplexProg (name, goal, constraints', bounds')
    end

fun relax_strict_ineqs (cplexProg (name, goals, constrs, bounds)) = 
    let
	fun relax cplexLe = cplexLeq
	  | relax cplexGe = cplexGeq
	  | relax x = x
	
	fun relax_constr (n, cplexConstr(c, (t1, t2))) = 
	    (n, cplexConstr(relax c, (t1, t2)))
	    
	fun relax_bounds (cplexBounds (t1, c1, t2, c2, t3)) = 
	    cplexBounds (t1, relax c1, t2, relax c2, t3)
	  | relax_bounds (cplexBound (t1, c, t2)) = 
	    cplexBound (t1, relax c, t2)
    in
	cplexProg (name, 
		   goals, 
		   map relax_constr constrs, 
		   map relax_bounds bounds)
    end

datatype cplexResult = Unbounded 
		     | Infeasible 
		     | Undefined
		     | Optimal of string * ((string * string) list)

fun is_separator x = forall (fn c => c = #"-") (String.explode x)

fun is_sign x = (x = "+" orelse x = "-")

fun is_colon x = (x = ":")

fun is_resultsymbol a = 
    let
	val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~-"
	fun is_symbol_char c = Char.isAlphaNum c orelse 
			       exists (fn d => d=c) symbol_char
	fun is_symbol_start c = is_symbol_char c andalso 
				not (Char.isDigit c) andalso 
				not (c= #".") andalso
				not (c= #"-")
	val b = String.explode a
    in
	b <> [] andalso is_symbol_start (hd b) andalso 
	forall is_symbol_char b
    end

val TOKEN_SIGN = 100
val TOKEN_COLON = 101
val TOKEN_SEPARATOR = 102

fun load_glpkResult name = 
    let
	val flist = [(is_NL, TOKEN_NL), 
		     (is_blank, TOKEN_BLANK),
		     (is_num, TOKEN_NUM),
		     (is_sign, TOKEN_SIGN),
                     (is_colon, TOKEN_COLON),
		     (is_cmp, TOKEN_CMP),
		     (is_resultsymbol, TOKEN_SYMBOL),
		     (is_separator, TOKEN_SEPARATOR)]

	val tokenize = tokenize_general flist

	val f = TextIO.openIn name  

	val rest = ref []
		   
	fun readToken_helper () =
	    if length (!rest) > 0 then  
		let val u = hd (!rest) in 
		    (
		     rest := tl (!rest); 
		     SOME u
		    ) 
		end
	    else 
		let val s = TextIO.inputLine f in
		    if s = "" then NONE else (rest := tokenize s; readToken_helper())					     
		end           
		
	fun is_tt tok ty = (tok <> NONE andalso (fst (the tok)) = ty)

	fun pushToken a = if a = NONE then () else (rest := ((the a)::(!rest)))
	
	fun readToken () =
	    let val t = readToken_helper () in
		if is_tt t TOKEN_BLANK then 
		    readToken ()
		else if is_tt t TOKEN_NL then 
		    let val t2 = readToken_helper () in
			if is_tt t2 TOKEN_SIGN then 
			    (pushToken (SOME (TOKEN_SEPARATOR, "-")); t)
			else 
			    (pushToken t2; t)
		    end
		else if is_tt t TOKEN_SIGN then
		    let val t2 = readToken_helper () in
			if is_tt t2 TOKEN_NUM then
			    (SOME (TOKEN_NUM, (snd (the t))^(snd (the t2))))
			else
			    (pushToken t2; t)
		    end
		else
		    t
	    end	    		

        fun readRestOfLine P = 
	    let 
		val t = readToken ()
	    in
		if is_tt t TOKEN_NL orelse t = NONE 
		then P
		else readRestOfLine P
	    end			
							 
	fun readHeader () = 
	    let
		fun readStatus () = readRestOfLine ("STATUS", snd (the (readToken ())))
		fun readObjective () = readRestOfLine ("OBJECTIVE", snd (the (readToken (); readToken (); readToken ())))
		val t1 = readToken ()
		val t2 = readToken ()
	    in
		if is_tt t1 TOKEN_SYMBOL andalso is_tt t2 TOKEN_COLON 
		then
		    case to_upper (snd (the t1)) of
			"STATUS" => (readStatus ())::(readHeader ())
		      | "OBJECTIVE" => (readObjective())::(readHeader ())
		      | _ => (readRestOfLine (); readHeader ())
		else
		    (pushToken t2; pushToken t1; [])
	    end

	fun skip_until_sep () = 
	    let val x = readToken () in
		if is_tt x TOKEN_SEPARATOR then
		    readRestOfLine ()
		else 
		    skip_until_sep ()
	    end

	fun load_value () =
	    let 
		val t1 = readToken ()
		val t2 = readToken ()
	    in
		if is_tt t1 TOKEN_NUM andalso is_tt t2 TOKEN_SYMBOL then
		    let 
			val t = readToken ()
			val state = if is_tt t TOKEN_NL then readToken () else t
			val _ = if is_tt state TOKEN_SYMBOL then () else raise (Load_cplexResult "state expected")
			val k = readToken ()
		    in
			if is_tt k TOKEN_NUM then
			    readRestOfLine (SOME (snd (the t2), snd (the k)))
			else
			    raise (Load_cplexResult "number expected")
		    end				
		else
		    (pushToken t2; pushToken t1; NONE)
	    end

	fun load_values () = 
	    let val v = load_value () in
		if v = NONE then [] else (the v)::(load_values ())
	    end
	
	val header = readHeader () 

	val result = 
	    case AList.lookup (op =) header "STATUS" of
		SOME "INFEASIBLE" => Infeasible
	      | SOME "UNBOUNDED" => Unbounded
	      | SOME "OPTIMAL" => Optimal (the (AList.lookup (op =) header "OBJECTIVE"), 
					   (skip_until_sep (); 
					    skip_until_sep ();
					    load_values ()))
	      | _ => Undefined

	val _ = TextIO.closeIn f
    in
	result
    end
    handle (Tokenize s) => raise (Load_cplexResult ("Tokenize: "^s))
	 | Option => raise (Load_cplexResult "Option")
	 | x => raise x

fun load_cplexResult name = 
    let
	val flist = [(is_NL, TOKEN_NL), 
		     (is_blank, TOKEN_BLANK),
		     (is_num, TOKEN_NUM),
		     (is_sign, TOKEN_SIGN),
                     (is_colon, TOKEN_COLON),
		     (is_cmp, TOKEN_CMP),
		     (is_resultsymbol, TOKEN_SYMBOL)]

	val tokenize = tokenize_general flist

	val f = TextIO.openIn name  

	val rest = ref []
		   
	fun readToken_helper () =
	    if length (!rest) > 0 then  
		let val u = hd (!rest) in 
		    (
		     rest := tl (!rest); 
		     SOME u
		    ) 
		end
	    else 
		let 
		    val s = TextIO.inputLine f 
		in		    
		    if s = "" then NONE else (rest := tokenize s; readToken_helper())					     
		end           
		
	fun is_tt tok ty = (tok <> NONE andalso (fst (the tok)) = ty)

	fun pushToken a = if a = NONE then () else (rest := ((the a)::(!rest)))
	
	fun readToken () =
	    let val t = readToken_helper () in
		if is_tt t TOKEN_BLANK then 
		    readToken ()
		else if is_tt t TOKEN_SIGN then
		    let val t2 = readToken_helper () in
			if is_tt t2 TOKEN_NUM then
			    (SOME (TOKEN_NUM, (snd (the t))^(snd (the t2))))
			else
			    (pushToken t2; t)
		    end
		else
		    t
	    end	    		

        fun readRestOfLine P = 
	    let 
		val t = readToken ()
	    in
		if is_tt t TOKEN_NL orelse t = NONE 
		then P
		else readRestOfLine P
	    end			

	fun readHeader () = 
	    let
		fun readStatus () = readRestOfLine ("STATUS", snd (the (readToken ())))
		fun readObjective () = 
		    let
			val t = readToken ()
		    in
			if is_tt t TOKEN_SYMBOL andalso to_upper (snd (the t)) = "VALUE" then 
			    readRestOfLine ("OBJECTIVE", snd (the (readToken())))
			else
			    readRestOfLine ("OBJECTIVE_NAME", snd (the t))
		    end
		  
		val t = readToken ()
	    in
		if is_tt t TOKEN_SYMBOL then
		    case to_upper (snd (the t)) of 
			"STATUS" => (readStatus ())::(readHeader ())
		      | "OBJECTIVE" => (readObjective ())::(readHeader ())
		      | "SECTION" => (pushToken t; [])
		      | _ => (readRestOfLine (); readHeader ())
		else
		    (readRestOfLine (); readHeader ())
	    end

	fun skip_nls () = 
	    let val x = readToken () in
		if is_tt x TOKEN_NL then 
		    skip_nls () 
		else 
		    (pushToken x; ())
	    end

	fun skip_paragraph () = 
	    if is_tt (readToken ()) TOKEN_NL then
		(if is_tt (readToken ()) TOKEN_NL then
		     skip_nls ()
		 else
		     skip_paragraph ())
	    else			
		skip_paragraph ()	

	fun load_value () =
	    let 
		val t1 = readToken ()
		val t1 = if is_tt t1 TOKEN_SYMBOL andalso snd (the t1) = "A" then readToken () else t1
	    in
		if is_tt t1 TOKEN_NUM then
		    let 
			val name = readToken ()
			val status = readToken ()
			val value = readToken ()
		    in
			if is_tt name TOKEN_SYMBOL andalso
			   is_tt status TOKEN_SYMBOL andalso
			   is_tt value TOKEN_NUM 
			then
			    readRestOfLine (SOME (snd (the name), snd (the value)))
			else
			    raise (Load_cplexResult "column line expected")
		    end
		else 
		    (pushToken t1; NONE)
	    end

	fun load_values () = 
	    let val v = load_value () in
		if v = NONE then [] else (the v)::(load_values ())
	    end
	
	val header = readHeader ()

	val result = 
	    case AList.lookup (op =) header "STATUS" of
		SOME "INFEASIBLE" => Infeasible
	      | SOME "NONOPTIMAL" => Unbounded
	      | SOME "OPTIMAL" => Optimal (the (AList.lookup (op =) header "OBJECTIVE"), 
					   (skip_paragraph (); 
					    skip_paragraph (); 
					    skip_paragraph (); 
					    skip_paragraph (); 
					    skip_paragraph ();
					    load_values ()))
	      | _ => Undefined

	val _ = TextIO.closeIn f
    in
	result
    end
    handle (Tokenize s) => raise (Load_cplexResult ("Tokenize: "^s))
	 | Option => raise (Load_cplexResult "Option")
	 | x => raise x

exception Execute of string;

fun tmp_file s = Path.pack (Path.expand (File.tmp_path (Path.make [s])));
fun wrap s = "\""^s^"\"";

fun solve_glpk prog =
    let 
	val name = LargeInt.toString (Time.toMicroseconds (Time.now ()))	
	val lpname = tmp_file (name^".lp")
	val resultname = tmp_file (name^".txt")
	val _ = save_cplexFile lpname prog		
	val cplex_path = getenv "GLPK_PATH"
	val cplex = if cplex_path = "" then "glpsol" else cplex_path
	val command = (wrap cplex)^" --lpt "^(wrap lpname)^" --output "^(wrap resultname)
	val answer = execute command
    in
	let
	    val result = load_glpkResult resultname
	    val _ = OS.FileSys.remove lpname
	    val _ = OS.FileSys.remove resultname
	in
	    result
	end 
	handle (Load_cplexResult s) => raise (Execute ("Load_cplexResult: "^s^"\nExecute: "^answer))
	     | _ => raise (Execute answer)
    end		     
    
fun solve_cplex prog =
    let 
	fun write_script s lp r =
	    let
		val f = TextIO.openOut s
		val _ = TextIO.output (f, "read\n"^lp^"\noptimize\nwrite\n"^r^"\nquit")
		val _ = TextIO.closeOut f
	    in 
		() 
	    end 
	    
	val name = LargeInt.toString (Time.toMicroseconds (Time.now ()))	
	val lpname = tmp_file (name^".lp")
	val resultname = tmp_file (name^".txt")
	val scriptname = tmp_file (name^".script")
	val _ = save_cplexFile lpname prog		
	val cplex_path = getenv "CPLEX_PATH"
	val cplex = if cplex_path = "" then "cplex" else cplex_path
	val _ = write_script scriptname lpname resultname
	val command = (wrap cplex)^" < "^(wrap scriptname)^" > /dev/null"
	val answer = "return code "^(Int.toString (system command))
    in
	let
	    val result = load_cplexResult resultname
	    val _ = OS.FileSys.remove lpname
	    val _ = OS.FileSys.remove resultname
	    val _ = OS.FileSys.remove scriptname
	in
	    result
	end 
    end

fun solve prog = 
    case get_solver () of
      SOLVER_DEFAULT => 
        (case getenv "LP_SOLVER" of
	   "CPLEX" => solve_cplex prog
         | "GLPK" => solve_glpk prog
         | _ => raise (Execute ("LP_SOLVER must be set to CPLEX or to GLPK")))
    | SOLVER_CPLEX => solve_cplex prog
    | SOLVER_GLPK => solve_glpk prog
		   
end;

(*
val demofile = "/home/obua/flyspeck/kepler/LP/cplexPent2.lp45"
val demoout = "/home/obua/flyspeck/kepler/LP/test.out"
val demoresult = "/home/obua/flyspeck/kepler/LP/try/test2.sol"
	   
fun loadcplex () = Cplex.relax_strict_ineqs 
		   (Cplex.load_cplexFile demofile)

fun writecplex lp = Cplex.save_cplexFile demoout lp

fun test () = 
    let
	val lp = loadcplex ()
	val lp2 = Cplex.elim_nonfree_bounds lp
    in
	writecplex lp2
    end

fun loadresult () = Cplex.load_cplexResult demoresult;
*)

(*val prog = Cplex.load_cplexFile "/home/obua/tmp/pent/graph_0.lpt";
val _ = Cplex.solve prog;*)