Permalink
Browse files

Initial commit

  • Loading branch information...
0 parents commit 5e09578216e9f8423111fca1e99ecf467594b7c6 @pkhuong committed May 24, 2012
Showing with 761 additions and 0 deletions.
  1. +67 −0 EGlib-long-lines.patch
  2. +7 −0 rational-simplex.asd
  3. +531 −0 rational-simplex.lisp
  4. +156 −0 readme.md
67 EGlib-long-lines.patch
@@ -0,0 +1,67 @@
+--- src/eg_io.c 2011-09-27 13:05:55.000000000 +0200
++++ src/eg_io.c 2012-05-23 18:08:53.000000000 +0200
+@@ -430,26 +430,26 @@ struct EGioFile_st {int type; void*file;
+ /* ========================================================================= */
+ int EGioWrite(EGioFile_t*file,const char*const string)
+ {
+- char buf[EGio_BUFSIZE];
++ /* char buf[EGio_BUFSIZE]; */
+ int len;
+- buf[EGio_BUFSIZE-1] = 0;
+- snprintf(buf,EGio_BUFSIZE,"%s",string);
+- len = strlen(buf);
+- if(len<=0 || len >= EGio_BUFSIZE || buf[EGio_BUFSIZE-1]!=0) return 0;
++ /* buf[EGio_BUFSIZE-1] = 0; */
++ /* snprintf(buf,EGio_BUFSIZE,"%s",string); */
++ len = strlen(string);
++ /* if(len<=0 || len >= EGio_BUFSIZE || buf[EGio_BUFSIZE-1]!=0) return 0; */
+ switch(file->type)
+ {
+ case EGIO_PLAIN:
+- return fwrite(buf, (size_t)1, (size_t)len, (FILE*)(file->file));
++ return fwrite((void*)string, (size_t)1, (size_t)len, (FILE*)(file->file));
+ case EGIO_ZLIB:
+ #ifdef HAVE_LIBZ
+- return gzwrite((gzFile)(file->file),buf,(unsigned)len);
++ return gzwrite((gzFile)(file->file),(void*)string,(unsigned)len);
+ #else
+ fprintf(stderr,"no zlib support\n");
+ return 0;
+ #endif
+ case EGIO_BZLIB:
+ #ifdef HAVE_LIBBZ2
+- return BZ2_bzwrite((BZFILE*)(file->file),buf,len);
++ return BZ2_bzwrite((BZFILE*)(file->file),(void*)string,len);
+ #else
+ fprintf(stderr,"no bzip2 support\n");
+ return 0;
+@@ -462,13 +462,23 @@ int EGioWrite(EGioFile_t*file,const char
+ /* ========================================================================= */
+ int EGioPrintf(EGioFile_t*file,const char* format, ...)
+ {
+- char buf[EGio_BUFSIZE];
+- va_list va;
+- buf[EGio_BUFSIZE-1]=0;
++ char buf[EGio_BUFSIZE], * large_buf = NULL;
++ char * to_print = buf;
++ int count, ret;
++ va_list va, copy;
+ va_start(va,format);
+- vsnprintf(buf,EGio_BUFSIZE,format,va);
+- va_end(va);
+- return EGioWrite(file,buf);
++ va_copy(copy, va);
++ count = vsnprintf(buf,EGio_BUFSIZE,format,va);
++ if (count < EGio_BUFSIZE)
++ goto done;
++ to_print = large_buf = (char*)malloc(count+1);
++ vsnprintf(large_buf,count+1,format,copy);
++done:
++ ret = EGioWrite(file,to_print);
++ free(large_buf);
++ va_end(copy);
++ va_end(va);
++ return ret;
+ }
+ /* ========================================================================= */
+ EGioFile_t* EGioOpenFILE(FILE*ifile)
7 rational-simplex.asd
@@ -0,0 +1,7 @@
+(asdf:defsystem "rational-simplex"
+ :version "0.0.1"
+ :licence "BSD"
+ :description "Tiny modeling language and wrapper around QSopt-Exact"
+ :serial t
+ :depends-on ("cl-ppcre" #-sbcl "trivial-shell")
+ :components ((:file "rational-simplex")))
531 rational-simplex.lisp
@@ -0,0 +1,531 @@
+(defpackage "RATIONAL-SIMPLEX"
+ (:use "CL" "CL-PPCRE")
+ (:nicknames "LP")
+ (:export "MODEL" "*MODEL*" "NAME" "SENSE" "OBJ" "CONSTRAINTS"
+ "WITH-MODEL"
+ "VAR" "LOWER-BOUND" "UPPER-BOUND"
+ "LINEAR-EXPRESSION" "LINEXPR" "COEFS" "CONSTANT"
+ "CONSTRAINT" "LHS" "CMP" "RHS"
+
+ "ADD" "ADDF" "SUB" "SUBF" "SCALE" "REMOVE-VAR"
+ "ADD-CONSTRAINT" "DEL-CONSTRAINT" "CONSTRAIN" "POST"
+
+ "%PRINT-MODEL" "PRINT-FLOAT-MODEL" "PRINT-RATIONAL-MODEL"
+
+ "SOLVE" "SOLVE-WITH-PRINTER"
+ "*SOLVERS-PATH*" "*INSTANCE-STEM*"))
+
+(in-package "RATIONAL-SIMPLEX")
+(defun maybe-rational (x)
+ (if (numberp x)
+ (rational x)
+ x))
+(defun copy-hash-table (x)
+ (let ((dst (make-hash-table :test (hash-table-test x))))
+ (maphash (lambda (k v)
+ (setf (gethash k dst) v))
+ x)
+ dst))
+
+(defclass model ()
+ ((name :reader name :initarg :name :initform "Unnamed model")
+ (sense :reader sense :initarg :sense :initform :minimize)
+ (obj :accessor obj :initarg :obj)
+ (constraints :accessor %constraints
+ :initarg :constraints
+ :initform (make-hash-table))
+ (vars :reader %vars
+ :initarg :vars
+ :initform (make-hash-table))
+ (id-var :reader id-var
+ :initform (make-array 32 :adjustable t
+ :fill-pointer 0))))
+
+(defgeneric (setf sense) (sense model)
+ (:method (sense (model model))
+ (check-type sense (member :maximize :minimize))
+ (setf (slot-value model 'sense) sense)))
+
+(defgeneric constraints (model))
+(defgeneric (setf constraints) (value model))
+(defmethod constraints ((model model))
+ (copy-hash-table (%constraints model)))
+(defmethod (setf constraints) (value (model model))
+ (check-type value hash-table)
+ (setf (%constraints model) (copy-hash-table value)))
+
+(declaim (type model *model*))
+(defvar *model*)
+
+(defclass var ()
+ ((name :reader name :initarg :name)
+ (id :reader id :initform (length (id-var *model*)))
+ (model :reader %model :initform *model*)
+ (lower-bound :reader lower-bound
+ :initarg :lower-bound
+ :initform nil)
+ (upper-bound :reader upper-bound
+ :initarg :upper-bound
+ :initform nil)))
+
+(defclass linear-expression ()
+ ((coefs :reader %coefs
+ :initarg :coefs
+ :initform (make-hash-table))
+ (constant :reader constant
+ :initarg :constant
+ :initform 0)))
+
+(declaim (type unsigned-byte *constraint-id*))
+(defvar *constraint-id* 0)
+
+(defclass constraint ()
+ ((lhs :reader %lhs :initarg :lhs)
+ (cmp :reader cmp :initarg :cmp)
+ (rhs :reader rhs :initarg :rhs)
+ (model :reader %model :initform *model*)
+ (id :reader id :initform (incf *constraint-id*))))
+
+(defun adjust-linearexpr (dst delta-coefs)
+ (loop with coefs = (%coefs dst)
+ for (var . scale) in delta-coefs do
+ (check-type var var)
+ (check-type scale real)
+ (let ((coef (incf (gethash var coefs 0)
+ (rational scale))))
+ (when (zerop coef)
+ (remhash var coefs)))
+ finally (return dst)))
+
+(defun make-linexpr (&key constant coefs)
+ (let ((dst (make-instance 'linear-expression
+ :constant (rational (or constant 0)))))
+ (adjust-linearexpr dst coefs)))
+
+(defun copy-linexpr (src &key
+ constant
+ delta-constant
+ delta-coefs)
+ (let ((new (make-instance
+ 'linear-expression
+ :coefs (copy-hash-table
+ (%coefs src))
+ :constant (rational (+ (or constant
+ (constant src))
+ (or delta-constant
+ 0))))))
+ (adjust-linearexpr new delta-coefs)
+ new))
+
+(defun linexpr (&optional constant &rest coef-var)
+ (make-linexpr
+ :constant (or constant 0)
+ :coefs (loop for (coef var) on coef-var by #'cddr
+ do (when (typep coef 'var)
+ (rotatef coef var))
+ (check-type coef real)
+ (check-type var var)
+ collect (cons var coef))))
+
+(defgeneric %add (x y z))
+(defun add (x y &optional z)
+ (setf z (or z 1))
+ (when (realp y)
+ (rotatef y z))
+ (assert (realp z))
+ (%add (maybe-rational x)
+ (maybe-rational y)
+ (maybe-rational z)))
+
+(defun sub (x y &optional z)
+ (setf z (or z 1))
+ (add x y (- z)))
+
+(defun scale (x scale)
+ (add 0 x scale))
+
+(defun remove-var (linexpr var)
+ (check-type linexpr linear-expression)
+ (check-type var var)
+ (let ((coef (gethash var (%coefs linexpr))))
+ (if coef
+ (sub linexpr var)
+ linexpr)))
+
+(define-modify-macro addf (y &optional (z 1)) add)
+(define-modify-macro subf (y &optional (z 1)) sub)
+(define-modify-macro scalef (scale) scale)
+(define-modify-macro remove-varf (var) remove-var)
+
+(defmethod %add ((x real) (y real) z)
+ (+ x (* y z)))
+
+(defmethod %add ((x real) (y var) z)
+ (make-linexpr :constant x
+ :coefs `((,y . ,z))))
+
+(defmethod %add ((x var) (y real) z)
+ (make-linexpr :constant (* y z)
+ :coefs `((,x . 1))))
+
+(defmethod %add ((x var) (y var) z)
+ (make-linexpr :coefs `((,x . 1)
+ (,y . ,z))))
+
+(defmethod %add ((x linear-expression) (y real)
+ z)
+ (copy-linexpr x :delta-constant (* y z)))
+
+(defun linexpr-coefs (linexpr &optional scale)
+ (setf scale (or scale 1))
+ (let ((coefs '()))
+ (maphash (lambda (k v)
+ (push (cons k (* v scale))
+ coefs))
+ (%coefs linexpr))
+ (sort coefs #'< :key (lambda (x)
+ (id (car x))))))
+
+(defmethod %add ((x real) (y linear-expression)
+ z)
+ (make-linexpr :constant (+ x (* (constant y) z))
+ :coefs (linexpr-coefs y z)))
+
+(defmethod %add ((x linear-expression) (y var)
+ z)
+ (copy-linexpr x :delta-coefs `((,y . ,z))))
+
+(defmethod %add ((x var) (y linear-expression)
+ z)
+ (make-linexpr :constant (* z (constant y))
+ :coefs `(acons x 1
+ (linexpr-coefs y z))))
+
+(defmethod %add ((x linear-expression) (y linear-expression)
+ z)
+ (copy-linexpr
+ x
+ :delta-constant (* z (constant y))
+ :delta-coefs (linexpr-coefs y z)))
+
+(defun model (&key name sense)
+ (check-type name (or null string))
+ (check-type sense (member nil :minimize :maximize))
+ (make-instance 'model :name (or name "Unnamed model")
+ :obj (make-linexpr)
+ :sense (or sense :minimize)))
+
+(defmacro with-model ((&key name sense) &body body)
+ `(let ((*model* (model :name ,name :sense ,sense))
+ (*constraint-id* 0))
+ ,@body))
+
+(defun var (&key name (lower 0) upper obj)
+ (check-type name (or null string))
+ (check-type lower (or null real))
+ (check-type upper (or null real))
+ (check-type obj (or null real))
+ (let ((var (make-instance 'var
+ :name name
+ :lower-bound (maybe-rational
+ lower)
+ :upper-bound (maybe-rational
+ upper))))
+ (when obj
+ (addf (obj *model*) var (maybe-rational obj)))
+ (setf (gethash var (%vars *model*)) t)
+ (let ((id (id var))
+ (id-var (id-var *model*)))
+ (assert (= id (length id-var)))
+ (vector-push-extend var (id-var *model*)))
+ var))
+
+(defmethod print-object ((var var) stream)
+ (print-unreadable-object (var stream :type t)
+ (format stream "~A [~F:~F]"
+ (or (name var)
+ (format nil "v~A" (id var)))
+ (let ((lower (lower-bound var)))
+ (if lower (float lower 1d0) "-inf"))
+ (let ((upper (upper-bound var)))
+ (if upper (float upper 1d0) "+inf")))))
+
+(defun constraint (lhs cmp rhs)
+ (check-type lhs (or var linear-expression real))
+ (check-type rhs (or var linear-expression real))
+ (check-type cmp (member <= >= =))
+ (setf lhs (maybe-rational lhs)
+ rhs (maybe-rational rhs))
+ (when (and (realp lhs)
+ (realp rhs))
+ (return-from constraint
+ (make-instance 'constraint
+ :lhs '()
+ :cmp cmp
+ :rhs (- rhs lhs))))
+ (let* ((linexpr (add lhs rhs -1))
+ (coefs (linexpr-coefs linexpr)))
+ (every (lambda (coef)
+ (eql *model* (%model (car coef))))
+ coefs)
+ (make-instance 'constraint
+ :lhs coefs
+ :cmp cmp
+ :rhs (- (constant linexpr)))))
+
+(defgeneric lhs (object))
+(defmethod lhs ((constraint constraint))
+ (make-linexpr :coefs (%lhs constraint)))
+
+(defun add-constraint (model constraint)
+ (check-type constraint constraint)
+ (check-type model model)
+ (assert (eql model (%model constraint)))
+ (setf (gethash constraint (%constraints model)) t)
+ constraint)
+
+(defun del-constraint (model constraint)
+ (check-type constraint constraint)
+ (check-type model model)
+ (remhash constraint (%constraints model)))
+
+(defun constrain (lhs cmp rhs &optional model)
+ (add-constraint (or model *model*)
+ (constraint lhs cmp rhs)))
+
+(defun post (lhs cmp rhs &optional model)
+ (constrain lhs cmp rhs model))
+
+(defun print-coefs (coefs numberify stream)
+ (format stream "~{~A v~A~^ + ~}"
+ (loop for (var . scale) in coefs
+ collect (funcall numberify scale)
+ collect (id var))))
+
+(defun print-constraints (constraints numberify stream)
+ (let ((rows '()))
+ (maphash (lambda (k v) v
+ (push k rows))
+ constraints)
+
+ (setf rows (sort rows #'< :key #'id))
+ (dolist (row rows)
+ (format stream " C~A: " (id row))
+ (print-coefs (%lhs row) numberify stream)
+ (format stream " ~A ~A~%"
+ (ecase (cmp row)
+ (<= "<=")
+ (>= ">=")
+ (= "="))
+ (funcall numberify (rhs row))))))
+
+(defun print-bounds (vars numberify stream)
+ (let ((columns '()))
+ (maphash (lambda (k v) v
+ (push k columns))
+ vars)
+ (setf columns (sort columns #'< :key #'id))
+ (dolist (var columns)
+ (let ((lb (lower-bound var))
+ (ub (upper-bound var))
+ (id (id var)))
+ (cond ((or lb ub)
+ (format stream " ")
+ (when lb
+ (format stream "~A <= "
+ (funcall numberify lb)))
+ (format stream "v~A" id)
+ (when ub
+ (format stream " <= ~A"
+ (funcall numberify ub)))
+ (terpri stream))
+ (t
+ (format stream " v~A free~%"
+ id)))))))
+
+(defun %print-model (model numberify stream
+ &aux (*read-default-float-format* 'double-float))
+ (format stream "Problem~%")
+ (format stream " ~A~%" (or (name model)
+ "Unnamed model"))
+ (format stream "~A~%" (ecase (sense model)
+ (:minimize "Minimize")
+ (:maximize "Maximize")))
+ (format stream " obj: obj_var~%")
+ (format stream "Subject To~%")
+ (format stream " obj_con: ")
+ (let* ((obj (obj model))
+ (coefs (linexpr-coefs obj))
+ (constant (constant obj)))
+ (when coefs
+ (print-coefs coefs numberify stream))
+ (format stream " - obj_var = ~A~%"
+ (funcall numberify constant)))
+ (print-constraints (constraints model) numberify stream)
+ (format stream "Bounds~%")
+ (print-bounds (%vars *model*) numberify stream)
+ (format stream "End~%"))
+
+(defun print-float-model (model stream)
+ (let ((*read-default-float-format* 'double-float))
+ (%print-model model (lambda (x)
+ (float x 1d0))
+ stream)))
+
+(defun print-rational-model (model stream)
+ (let ((*read-default-float-format* 'double-float))
+ (%print-model model (lambda (x)
+ (rational x))
+ stream)))
+
+(defvar *solvers-path* (format nil "~A/bin/"
+ (directory-namestring *load-pathname*)))
+(defvar *instance-stem* "/tmp/rational-simplex-instance")
+(defvar *status-codes* '((1 . :optimal)
+ (2 . :infeasible)
+ (3 . :unbounded)
+ (4 . :unsolved)))
+
+(defun parse-output (s)
+ (let ((status nil)
+ (obj 0)
+ (var-values (make-hash-table :test #'equal))
+ (time nil)
+ (*read-default-float-format* 'double-float))
+ (labels ((parse-status (line)
+ (register-groups-bind ((#'parse-integer value))
+ ("LP Value: [-+0-9.]+, status ([0-9])" line)
+ (assert (not status))
+ (setf status value))
+ (register-groups-bind ((#'read-from-string sec))
+ ("Time for SOLVER: (.*) seconds." line)
+ (assert (not time))
+ (setf time sec)))
+ (parse-value (line)
+ (register-groups-bind (name (#'read-from-string value))
+ ("(.*) = ([-+0-9.e/]*)" line)
+ (assert (not (gethash name var-values)))
+ (setf (gethash name var-values) value)))
+ (parse-values (stream)
+ (loop for line = (read-line stream nil)
+ while line
+ do (parse-value line))))
+ (with-input-from-string (s s)
+ (loop for line = (read-line s nil)
+ while line
+ do (parse-status line)
+ when (equal line "Solution Values")
+ do (return (parse-values s)))))
+ (assert (and status time))
+ (values obj
+ (or (cdr (assoc status *status-codes*))
+ (error "Unknown status code ~A" status))
+ var-values
+ time)))
+
+(defun %solve-with-solver (solver printer basis-p
+ &aux
+ (lp-file (format nil "~A.lp"
+ *instance-stem*))
+ (bas-file (format nil "~A.bas"
+ *instance-stem*)))
+ (when printer
+ (with-open-file (s lp-file
+ :direction :output
+ :if-exists :supersede)
+ (funcall printer s)))
+ (let ((s #+sbcl (with-output-to-string (s)
+ (sb-ext:run-program (format nil "~A/~A"
+ *solvers-path* solver)
+ `("-O"
+ "-b" ,bas-file
+ ,@(and basis-p
+ `("-B" ,bas-file))
+ "-L" ,lp-file)
+ :output s
+ :error nil
+ :wait t))
+ #-sbcl
+ (trivial-shell:shell-command
+ (format nil "~A/~A -O -b ~A ~A -L ~A"
+ *solvers-path* solver
+ bas-file (if basis-p
+ (format nil "-B ~A" bas-file)
+ "")
+ lp-file))))
+ (parse-output s)))
+
+(defun double (x)
+ (float x 1d0))
+
+(defvar *solvers* '(("dbl_solver" double)
+ ("ldbl_solver" double)
+ ("float128_solver" double)
+ ("mpf_solver" double)
+ ("mpq_solver" rational)))
+
+(defun %solve (printer trace)
+ (let ((success nil)
+ (total-time 0)
+ prev-numberify)
+ (unwind-protect
+ (loop with count = (length *solvers*)
+ for i upfrom 1
+ for (solver numberify) in *solvers* do
+ (multiple-value-bind (obj status values time)
+ (%solve-with-solver solver
+ (and (not (eq numberify prev-numberify))
+ (lambda (stream)
+ (funcall printer numberify stream)))
+ success)
+ (setf prev-numberify numberify)
+ (format trace "~A: ~,3F sec ~F (~(~A~))~%"
+ solver time (double obj) status)
+ ;; FIXME: are basis file sane when !optimal/unbounded?
+ (setf success (member status '(:optimal :infeasible
+ :unbounded)))
+ (incf total-time time)
+ (when (= i count)
+ (return (values status obj values total-time)))))
+ (flet ((try-delete (x)
+ (when (probe-file x)
+ (delete-file x))))
+ (try-delete (format nil "~A.lp"
+ *instance-stem*))
+ (try-delete (format nil "~A.bas"
+ *instance-stem*))))))
+
+(defun solve-with-printer (printer &key (trace *trace-output*)
+ inexact double-only)
+ (let ((*solvers* (cond (double-only
+ (subseq *solvers* 0 1))
+ (inexact
+ (butlast *solvers*))
+ (t *solvers*))))
+ (%solve printer trace)))
+
+(defun solve (&key model (trace *trace-output*)
+ inexact double-only)
+ (setf model (or model *model*))
+ (multiple-value-bind (status obj values time)
+ (solve-with-printer (lambda (numberify stream)
+ (%print-model model numberify stream))
+ :trace trace
+ :inexact inexact
+ :double-only double-only)
+ (declare (ignore obj))
+ (let ((table (make-hash-table))
+ (obj nil)
+ (id-var (id-var model))
+ (*read-default-float-format* 'double-float))
+ (maphash (lambda (name value)
+ (cond ((equal name "obj_var")
+ (assert (not obj))
+ (setf obj value))
+ (t
+ (let* ((idx (parse-integer name :start 1))
+ (var (aref id-var idx)))
+ (assert (not (gethash var table)))
+ (setf (gethash var table) value)))))
+ values)
+ (values status obj table time))))
156 readme.md
@@ -0,0 +1,156 @@
+RATIONAL-SIMPLEX
+================
+
+This one-file system provides two things: a tiny modeling framework
+for linear programs, and a wrapper around Daniel Espinoza et al's
+QSopt-Exact, an exact (with rational arithmetic) solver for linear
+optimization problems.
+
+Installation
+------------
+
+The solver requires an installation of
+[QSopt-Exact](http://www.dii.uchile.cl/~daespino/ESolver_doc/main.html).
+The program depends on
+[EGlib](http://www.dii.uchile.cl/~daespino/EGlib_doc/main.html), an
+utility library, which has a bug that results in erroneous variable
+value output for very large fractions or floats. Both programs are
+GPL.
+
+If you intend to work with large denominators, you'll need to build
+both programs from sources (it only takes a few minutes). Applying
+`EGlib-long-lines.patch` with `patch -p0 < EGlib-long-lines.patch`
+from `EGlib-2.6.20/` before building EGlib (and QSopt_ex) should
+work... It's an ugly workaround, you probably don't want to look too
+closely (:
+
+`*solvers-path*` should point to the directory holding the QSopt_ex
+executables. It defaults to
+`[directory from which simplex.lisp is loaded]/bin`, which is probably
+not what you want if you're loading it with asdf or quicklisp.
+
+`*instance-stem*` defaults to `"/tmp/rational-simplex-instance"`;
+temporary files `/tmp/rational-simplex-instance.lp` and
+`/tmp/rational-simplex-instance.bas` will be created (and
+overwritten!). You probably want to change that to something more
+unique (and private) if you work on a shared machine.
+
+Other dependencies: `cl-ppcre` and `trivial-shell` on `#-sbcl`
+platforms (in which case, beware spaces in paths...).
+
+Modeling language
+-----------------
+
+The base object is the `model`. It stores the current objective,
+sense (should the objective function be minimised or maximised),
+variables and constraints. `with-model (&key name sense)` should be
+used in most cases (`name` is a string or nil, and `sense` is
+`:minimize`, the default, or `:maximise`).
+
+In the scope of `with-model`, `var &key name lower upper obj`
+instantiates a new variables associated with the current model. The
+lower bound defaults to 0, and the upper bound to none (`nil`). `obj`
+is the variable's coefficient in the objective function, and defaults
+to 0.
+
+Variables and real values can be composed together into immutable
+`linear-expression`s. `linexpr &optional constant &rest {coef var}*`
+creates a fresh linear expression corresponding to `constant +
+coef1*var1 + ...`. `add` or `sub` can be used to add or subtract
+reals, variables, or linear expressions; the optional third argument
+specifies a scaling value for the second argument. `scale` will scale
+it argument by a real, while `remove-var` can be used to completely
+remove a variable from a linear expression. `addf`, `subf`, `scalef`,
+`remove-varf` are convenience modify macros around the functions.
+
+The function `constraint lhs cmp rhs` returns a new constraints that
+asserts a relationship between two linear expressions (or variables or
+reals). `cmp` can be `<=`, `>=` or `=`. `add-constraint` and
+`del-constraint` add/delete a constraint to/from a model. `constrain
+lhs cmp rhs &optional model` or its synonym `post` create a constraint
+and add it to the model (defaults to `*model*`) in a single call. The
+return value is the new constraint itself, to easily `del-constraint`
+it.
+
+`constraints` is a place that ensures that a model always has unique
+ownership over the constraint store (copies are stored/returned).
+This can be used to simplify branch-and-bound type algorithms.
+
+`print-float-model` and `print-rational-model` can be used to print
+the current model in lp format. Values in the model, linear
+expressions, etc. are converted to rationals as early as possible, so
+`print-rational-model` does not lose precision. Most solvers work
+with floating point values; in that case `print-float-model` is more
+appropriate.
+
+Solver
+------
+
+`solve &key model trace inexact double-only` is used to solve a model
+specified in the framework described above. `model` defaults to
+`*model*`. `trace` specifies where short (a few lines line per solve)
+logging should be printed and defaults to `*trace-output*`. `inexact`
+defaults to `nil`; if true, the final solve in rational is skipped.
+`double-only` defaults to `nil`; if true, only the initial solve with
+double arithmetic is performed. The return values are a status
+(`:optimal`, `:infeasible`, `:unbounded` or `unsolved`), the objective
+value, a hash table from variables to values, and the total solution
+time.
+
+This is a wrapper around `solve-with-printer printer &key trace
+inexact double-only`. Instead of a model, this function receives a
+printer function which, givena function with which to canonicalize
+numbers (e.g. one that converts any real to a float) and a stream,
+prints an lp model to the stream. The return values are the last
+executed solver's status, the reported objective value, a hash table
+mapping variable names to values, and the total solution time. Note
+that the objective value is reported with little precision;
+associating the expression with a bogus variable and optimising on
+that variable will give full precision in the objective value.
+
+Example usage
+-------------
+
+This solves a silly 3-variable linear program, and reports the
+objective value and decision variables' values at an optimal
+solution.
+
+ CL-USER> (lp:with-model (:name "small" :sense :maximize)
+ (let ((x (lp:var :name "x" :lower 2 :obj 3))
+ (y (lp:var :name "y" :lower nil :obj 2))
+ (z (lp:var :name "z" :lower 1 :upper 10 :obj 4)))
+ (lp:post (lp:linexpr 0 3 x 2 y 1 z) '<= 12)
+ (lp:post (lp:add y x 5) '<= 10)
+ (multiple-value-bind (status obj values)
+ (lp:solve :double-only nil :inexact nil)
+ (assert (eql status :optimal))
+ (values obj (mapcar (lambda (var)
+ (cons var (gethash var values)))
+ (list x y z))))))
+ dbl_solver: 0.000 sec 0.0 (optimal)
+ ldbl_solver: 0.000 sec 0.0 (optimal)
+ float128_solver: 0.000 sec 0.0 (optimal)
+ mpf_solver: 0.010 sec 0.0 (optimal)
+ mpq_solver: 0.010 sec 0.0 (optimal)
+ 42
+ ((#<RATIONAL-SIMPLEX:VAR x [2.0:+inf]> . 2)
+ (#<RATIONAL-SIMPLEX:VAR y [-inf:+inf]> . -2)
+ (#<RATIONAL-SIMPLEX:VAR z [1.0:10.0]> . 10))
+
+Implementation
+--------------
+
+QSopt-Exact comes with a full adaptive-precision MIP solver... but it
+solution output routines seem broken on my linux, and I didn't
+particularly feel like going down that rabbit hole. Instead, I call
+LP solvers of increasing precision in sequence: double, long doubles,
+128 bit quads, multi-precision floats and finally rationals. The
+reason this works (and is interesting) is that the simplex is a
+combinatorial algorithm: the correctness of each step only depends on
+the *comparison* between two values (i.e. is `x >= y`). Better: the
+state determined at each step isn't numeric, but rather a basis (a set
+of variable or constraints). Thus, we can save the basis found by the
+double solver, and use it to warm start the more precise long double
+solver, etc. Our hope is that the float-based solvers will converge
+to the optimal solution, or very nearly, leaving only a few iterations
+until the final exact, slow, solver itself declares victory.

0 comments on commit 5e09578

Please sign in to comment.