Skip to content
Browse files

Faster branch and bound, fix some convergence in branch and cut

  • Loading branch information...
1 parent 490c746 commit b3c9f293eeca34d995bad34ed8620cbc7704748a @pkhuong committed Oct 7, 2012
View
6 demo/branch-and-cut-fit/branch-and-cut-fit.asd
@@ -3,9 +3,11 @@
:licence "BSD"
:description "Branch and cut solver for minimax polynomial approximation"
:serial t
- :depends-on ("rational-simplex" "computable-reals")
+ :depends-on ("rational-simplex" "computable-reals" "sb-concurrency"
+ "sb-md5")
:components ((:file "utility")
(:file "linf-fit")
(:file "newton")
(:file "find-extrema")
- (:file "driver")))
+ (:file "driver")
+ (:file "print-pareto")))
View
302 demo/branch-and-cut-fit/driver.lisp
@@ -13,30 +13,31 @@
:adjustable t :fill-pointer 0))
(seen (make-hash-table)))
(flet ((point (x)
- (let ((x (round-to-double x)))
+ (let ((x (round-to-float x)))
(unless (or (gethash x seen)
(< x from)
(> x to))
(setf (gethash x seen) t)
(vector-push-extend (make-point x) points)))))
(point from)
- (point (double-float-from-bits (double-float-bits from) 1))
+ (point (float-from-bits (float-bits from) 1))
(point to)
- (point (double-float-from-bits (double-float-bits to) -1))
+ (point (float-from-bits (float-bits to) -1))
(loop with stride = (/ (- to from) 2)
for i from 1 upto count
for k = (1+ (pull-bits
(computable-reals:cos-r (computable-reals:*r
computable-reals:+pi-r+
(/ (1- (* 2 i))
(* 2 count))))))
- for x = (round-to-double (+ from (* k stride)))
+ for x = (round-to-float (+ from (* k stride)))
do (point x)
- (point (double-float-from-bits (double-float-bits x) -1))
- (point (double-float-from-bits (double-float-bits x) 1))))
+ (point (float-from-bits (float-bits x) -1))
+ (point (float-from-bits (float-bits x) 1))))
(values points seen)))
(defun eval-coefs (points coefs &optional round)
+ (declare (optimize debug))
(let ((poly (make-poly coefs round)))
(reduce #'max points
:key (lambda (x)
@@ -70,9 +71,10 @@
points)
(values (make-array (length points)
:adjustable t :fill-pointer (length points)
- :initial-contents points)))
+ :initial-contents points)
+ seen))
(chebyshev-nodes from to *initial-approximation*))
- (loop
+ (loop for i upfrom 0 do
(multiple-value-bind (diff coefs basis)
(solve-fit points :bounds bounds :advanced-basis basis)
(multiple-value-bind (actual-diff new-extrema distance root-count)
@@ -83,15 +85,16 @@
(if (zerop distance) 0 (log distance 2d0))
root-count)
(when (or (>= (+ diff (* *lp-gap* (abs diff))) actual-diff)
- (>= actual-diff cutoff))
+ (>= actual-diff cutoff)
+ (> i 4))
(return (values (make-state :points points
:point-set seen
:primal-solutions (make-hash-table :test #'equalp)
:nodes '()
:cutoff cutoff)
coefs diff basis)))
(map nil (lambda (extremum)
- (let ((extremum (round-to-double extremum)))
+ (let ((extremum (round-to-float extremum)))
(unless (gethash extremum seen)
(setf (gethash extremum seen) t)
(vector-push-extend (make-point extremum)
@@ -133,7 +136,7 @@
(map 'simple-vector
(lambda (x)
(multiple-value-bind (lo hi)
- (split-double x)
+ (split-float x)
(if (null lo)
x
(let* ((lo (rational lo))
@@ -161,19 +164,19 @@
(solve-fit (state-points *state*) :bounds bounds :advanced-basis basis)))
(defun branch-on-node (coefs bounds value basis)
- (let ((i (position-if #'split-double coefs)))
+ (let ((i (position-if #'split-float coefs)))
(and i
(let ((coef (elt coefs i))
(nodes '()))
(multiple-value-bind (lo hi)
- (split-double coef)
+ (split-float coef)
(format *trace-output*
" branch on ~A ~A~%"
i (mapcar (lambda (x)
(destructuring-bind (lo . hi) x
(cond ((and lo hi)
- (abs (- (double-float-bits hi)
- (double-float-bits lo))))
+ (abs (- (float-bits hi)
+ (float-bits lo))))
(lo :lo)
(hi :hi)
(t '-))))
@@ -202,7 +205,7 @@
(eval-node bounds basis)
(adjoin-solution (round-coefs coefs))
(cond ((every (lambda (x)
- (eql x (round-to-double x)))
+ (eql x (round-to-float x)))
coefs)
;; generate cuts
(multiple-value-bind (actual-diff extrema distance root-count)
@@ -215,7 +218,7 @@
(cond ((>= (+ value (* *lp-gap* (abs value))) actual-diff)
(let ((delta nil))
(map nil (lambda (extremum)
- (let ((extremum (round-to-double extremum)))
+ (let ((extremum (round-to-float extremum)))
(unless (gethash extremum (state-point-set *state*))
(setf (gethash extremum (state-point-set *state*)) t
delta t)
@@ -277,7 +280,7 @@
(setf *state* state
root-basis basis
root-points (copy-seq (state-points state)))
- (adjoin-solution (map 'simple-vector #'round-to-double coefs))
+ (adjoin-solution (map 'simple-vector #'round-to-float coefs))
(adjoin-solution (round-coefs coefs))
(update-points)
(let ((primal (best-upper-bound)))
@@ -316,6 +319,7 @@
(1+ (or (position 0 coefs :test-not #'eql :from-end t)
-1)))
+#+nil
(defun enumerate-pareto-front (degree from to function df d2f
&optional (cutoff 1d-1)
(cutoff-scale 50d0))
@@ -462,3 +466,265 @@
(map nil #'sb-thread:join-thread threads))
(setf done t)))
primals)))
+
+
+(defstruct int-node
+ ;; misnomer: it's actually degree+1 (# coefficients...)
+ degree bounds value)
+
+(defstruct search-state
+ ;; per degree
+ nogoods
+ pre-exploreds
+ bounds
+ ;; global
+ primals
+ queue
+ lock
+ cvar
+ ;; read-only
+ points
+ cutoff-scale
+ from to
+ function df d2f)
+
+(defun prune-trailing-nil (x)
+ (let ((rev (reverse x)))
+ (loop while rev
+ do (if (null (car rev))
+ (pop rev)
+ (return)))
+ (reverse rev)))
+
+(defun explore-one-node (node state)
+ (declare (optimize debug)
+ (type search-state state)
+ (type int-node node))
+ (let* ((degree (int-node-degree node))
+ (bounds (int-node-bounds node))
+ (value (int-node-value node))
+ (cutoff (cond ((>= degree 3)
+ (aref (search-state-bounds state)
+ (- degree 3)))
+ ((plusp degree)
+ (aref (search-state-bounds state) 0))
+ (t 1d300)))
+ (points (aref (search-state-points state) degree)))
+ (assert (= (length bounds) degree))
+ (when (and (plusp degree)
+ (eql 0 (elt bounds (1- degree))))
+ (return-from explore-one-node))
+ (unless (some #'null bounds)
+ (let* ((key (coerce bounds 'simple-vector))
+ (value (eval-coefs points key)))
+ (sb-thread:with-mutex ((search-state-lock state))
+ (setf (gethash key (search-state-primals state)) value)
+ (when (null bounds)
+ (let ((cutoffs (search-state-bounds state)))
+ (loop for i from degree below (length cutoffs)
+ do (setf (aref cutoffs i)
+ (min (aref cutoffs i) value)))))))
+ (return-from explore-one-node))
+ (when (and value (> value cutoff))
+ (sb-thread:with-mutex ((search-state-lock state))
+ (push bounds (aref (search-state-nogoods state) degree)))
+ (return-from explore-one-node))
+ (let ((nogoods (aref (search-state-nogoods state) degree)))
+ (when (some (lambda (nogood)
+ (every (lambda (nogood fixed)
+ (or (null nogood)
+ (eql nogood fixed)))
+ nogood bounds))
+ nogoods)
+ (return-from explore-one-node)))
+ (assert (plusp degree))
+ (multiple-value-bind (value coefs results)
+ (branch-and-cut (1- degree)
+ (search-state-from state)
+ (search-state-to state)
+ (search-state-function state)
+ (search-state-df state)
+ (search-state-d2f state)
+ :bounds
+ (mapcar (lambda (x)
+ (cons x x))
+ bounds)
+ :cutoff cutoff
+ :points points)
+ (sb-thread:with-mutex ((search-state-lock state))
+ (when (every #'null bounds)
+ (let ((cutoffs (search-state-bounds state)))
+ (when value
+ (loop for i from degree below (length cutoffs)
+ do (setf (aref cutoffs i)
+ (min (aref cutoffs i) value))))))
+ (let ((dst (search-state-primals state)))
+ (maphash (lambda (k v)
+ (setf (gethash k dst) v))
+ (state-primal-solutions results))))
+ (when (or (null value)
+ (zerop (hash-table-count (state-primal-solutions results)))
+ (> value cutoff))
+ (sb-thread:with-mutex ((search-state-lock state))
+ (format t "~A -> ~E / pruned~%" bounds value)
+ (let ((bounds (prune-trailing-nil bounds))
+ (nogoods (search-state-nogoods state)))
+ (loop for i from (max 1 (length bounds)) upto degree do
+ (push bounds (aref nogoods i)))))
+ (return-from explore-one-node))
+ (sb-thread:with-mutex ((search-state-lock state))
+ (let ((bounds (prune-trailing-nil bounds))
+ #+nil(scale (search-state-cutoff-scale state))
+ (cutoffs (search-state-bounds state))
+ (nogoods (search-state-nogoods state)))
+ (loop for i from (max 3 (length bounds)) upto degree do
+ (if (> value (* #+nil scale (aref cutoffs (- i 3))))
+ (push bounds (aref nogoods i))
+ (return))))
+ (let ((queue (search-state-queue state))
+ (seen (aref (search-state-pre-exploreds state)
+ degree))
+ (count 0))
+ (flet ((enqueue (i x)
+ (let ((bounds (copy-seq bounds)))
+ (setf (elt bounds i) x)
+ (unless (gethash bounds seen)
+ (setf (gethash bounds seen) t)
+ (incf count)
+ (sb-concurrency:enqueue
+ (make-int-node :degree degree
+ :bounds bounds
+ :value value)
+ queue)))))
+ (loop for i upto degree
+ for fixed in bounds
+ for x across coefs do
+ (when (null fixed)
+ (when (and (/= x 0) (< i degree))
+ (enqueue i 0))
+ (let ((abs-x (abs x))
+ (sign (signum x)))
+ (cond ((eql abs-x 2)
+ (enqueue i sign))
+ ((member abs-x '(0 1 2)))
+ ((< 1 abs-x 3)
+ (enqueue i (* 2 sign))
+ (enqueue i sign))
+ #+nil
+ ((< 1 abs-x 2)
+ (enqueue i (* 2 sign))
+ (enqueue i sign))
+ ((< 0 abs-x 1)
+ (enqueue i sign)))))))
+ (sb-thread:condition-broadcast (search-state-cvar state))
+ (format t "~A -> ~E (~A / ~A)~%"
+ bounds value (sb-concurrency:queue-count queue) count))))))
+
+(defvar *nthread* 11)
+
+(defun enumerate-pareto-front (degree from to
+ *loc-value* *loc-dvalue* *loc-d2value*
+ &optional (cutoff 1d300) (cutoff-scale 50d0))
+ (declare (optimize debug))
+ (let ((done nil)
+ (waitcount 0)
+ (stem rational-simplex:*instance-stem*)
+ (state (make-search-state
+ :nogoods (make-array (+ 2 degree)
+ :initial-element nil)
+ :pre-exploreds (map-into (make-array (+ 2 degree))
+ (lambda ()
+ (make-hash-table :test #'equalp)))
+ :bounds (make-array (+ 2 degree)
+ :initial-element (/ cutoff cutoff-scale))
+
+ :primals (make-hash-table :test #'equalp)
+ :queue (sb-concurrency:make-queue)
+ :lock (sb-thread:make-mutex)
+ :cvar (sb-thread:make-waitqueue)
+
+ :points
+ (coerce (loop for i upto (1+ degree)
+ collect
+ (let ((*dimension* i))
+ (chebyshev-nodes from to (ash 1 8))))
+ 'simple-vector)
+ :from from :to to
+ :cutoff-scale cutoff-scale
+ :function *loc-value*
+ :df *loc-dvalue*
+ :d2f *loc-d2value*))
+ (stdout *standard-output*))
+ (flet ((worker (state)
+ (let ((rational-simplex:*instance-stem*
+ (concatenate 'string stem
+ (format nil "-~A"
+ (sb-thread::thread-os-thread
+ sb-thread:*current-thread*))))
+ (*standard-output* stdout)
+ (*trace-output* (make-broadcast-stream))
+ computable-reals::(+log2-r+ (+r (ash-r (log-r2 1/7) 1) (log-r2 1/17)))
+ computable-reals::(+PI-R+ (-r (ash-r (atan-r1 1/10) 5)
+ (ash-r (atan-r1 1/515) 4)
+ (ash-r (atan-r1 1/239) 2)))
+ computable-reals::(+PI/2-R+ (ash-r +pi-r+ -1))
+ computable-reals::(+PI/4-R+ (ash-r +pi-r+ -2)))
+ (loop until done do
+ (sb-thread:with-mutex ((search-state-lock state))
+ (loop while (and (sb-concurrency:queue-empty-p
+ (search-state-queue state))
+ (not done))
+ do (incf waitcount)
+ (sb-thread:condition-wait
+ (search-state-cvar state)
+ (search-state-lock state))
+ (decf waitcount)))
+ (when done (return))
+ (loop
+ (multiple-value-bind (node fullp)
+ (sb-concurrency:dequeue (search-state-queue state))
+ (unless fullp
+ (return))
+ (multiple-value-bind (ok error)
+ (ignore-errors
+ (explore-one-node node state)
+ t)
+ (unless ok
+ (sb-thread:with-mutex ((search-state-lock state))
+ (format t "error: ~A~%" error))))))))))
+ (let ((threads (loop repeat *nthread*
+ collect
+ (sb-thread:make-thread (lambda ()
+ (worker state))
+ :name "worker"))))
+ (sb-thread:with-mutex ((search-state-lock state))
+ (loop for i from 0 upto (1+ degree) do
+ (sb-concurrency:enqueue (make-int-node
+ :degree i
+ :bounds (make-list i))
+ (search-state-queue state)))
+ (sb-thread:condition-broadcast (search-state-cvar state)))
+ (unwind-protect
+ (progn
+ (loop until (sb-thread:with-mutex ((search-state-lock state))
+ (when (and (= waitcount (length threads))
+ (sb-concurrency:queue-empty-p
+ (search-state-queue state)))
+ (setf done t)
+ (sb-thread:condition-broadcast
+ (search-state-cvar state))
+ t))
+ do (sleep 15))
+ (map nil #'sb-thread:join-thread threads))
+ (setf done t)))
+ (values (search-state-primals state)
+ state))))
+
+(defun print-hash-table (hash file)
+ (with-open-file (s file :direction :output
+ :if-exists :supersede)
+ (with-standard-io-syntax
+ (write (alexandria:hash-table-alist hash)
+ :readably t
+ :stream s)))
+ t)
View
82 demo/branch-and-cut-fit/find-extrema.lisp
@@ -1,20 +1,40 @@
(defun make-poly (coefs &optional round)
(declare (type simple-vector coefs))
+ (when (zerop (length coefs))
+ (return-from make-poly (lambda (x) x
+ 0)))
(if round
- (let ((coefs (map '(simple-array double-float 1)
- (lambda (x)
- (float x 1d0))
- coefs)))
- (lambda (x)
- (let ((x (float x 1d0))
- (acc 0d0)
- (x^i 1d0))
- (declare (type double-float x acc x^i))
- (map nil (lambda (coef)
- (incf acc (* x^i coef))
- (setf x^i (* x x^i)))
- coefs)
- (rational acc))))
+ (ecase *float-mode*
+ (double-float
+ (let ((coefs (map '(simple-array double-float 1)
+ (lambda (x)
+ (float x 1d0))
+ coefs)))
+ (lambda (x)
+ (let ((x (float x 1d0))
+ (acc 0d0)
+ (x^i 1d0))
+ (declare (type double-float x acc x^i))
+ (map nil (lambda (coef)
+ (incf acc (* x^i coef))
+ (setf x^i (* x x^i)))
+ coefs)
+ (rational acc)))))
+ (single-float
+ (let ((coefs (map '(simple-array single-float 1)
+ (lambda (x)
+ (float x 1s0))
+ coefs)))
+ (lambda (x)
+ (let ((x (float x 1s0))
+ (acc 0s0)
+ (x^i 1s0))
+ (declare (type single-float x acc x^i))
+ (map nil (lambda (coef)
+ (incf acc (* x^i coef))
+ (setf x^i (* x x^i)))
+ coefs)
+ (rational acc))))))
(let* ((lcm (reduce #'lcm coefs :key #'denominator))
(coefs (map 'simple-vector (lambda (x)
(* x lcm))
@@ -35,6 +55,8 @@
(/ acc (ash lcm (* lb-denom (1- (length coefs))))))))))
(defun make-dpoly (coefs &optional (diff 1) round)
+ (when (> diff (length coefs))
+ (return-from make-dpoly (make-poly #() round)))
(flet ((diff (coefs)
(let ((coefs (subseq coefs 1)))
(loop for i upfrom 0
@@ -44,7 +66,7 @@
(loop repeat diff do (setf coefs (diff coefs)))
(make-poly coefs round)))
-(defun approx-error (f df d2f coefs &optional round)
+(defun poly-approx-error (f df d2f coefs &optional round)
(values (let ((f~ (make-poly coefs round)))
(lambda (x)
(computable-reals:-r (funcall f x)
@@ -71,35 +93,35 @@
(pull-bits (funcall *loc-value* x)))))
(new-root (root lb ub)
(incf roots)
- (let ((root (round-to-double root)))
+ (let ((root (round-to-float root)))
(setf maximin-distance (max maximin-distance
- (min (abs (- (double-float-bits lb)
- (double-float-bits root)))
- (abs (- (double-float-bits ub)
- (double-float-bits root))))))
- (flet ((maybe-push (x &aux (x (round-to-double x)))
+ (min (abs (- (float-bits lb)
+ (float-bits root)))
+ (abs (- (float-bits ub)
+ (float-bits root))))))
+ (flet ((maybe-push (x &aux (x (round-to-float x)))
(unless (or (<= x lb)
(>= x ub))
(setf worst-diff (max worst-diff (delta x)))
(pushnew x new-extrema))))
(maybe-push root)
- (maybe-push (double-float-from-bits
- (1+ (double-float-bits root))))
- (maybe-push (double-float-from-bits
- (1- (double-float-bits root))))))))
+ (maybe-push (float-from-bits
+ (1+ (float-bits root))))
+ (maybe-push (float-from-bits
+ (1- (float-bits root))))))))
(setf worst-diff (delta (point-loc (aref points 0))))
(loop for i from 0
for j from 1 below (length points)
for lb = (point-loc (aref points i))
for ub = (point-loc (aref points j))
do (setf worst-diff (max worst-diff (delta ub)))
- when (> (- (double-float-bits ub) (double-float-bits lb)) 1)
- do (let ((root (and (> (- (double-float-bits ub)
- (double-float-bits lb))
+ when (> (- (float-bits ub) (float-bits lb)) 1)
+ do (let ((root (and (> (- (float-bits ub)
+ (float-bits lb))
1)
(newton lb ub ddelta-f d2delta-f))))
(when root
- (new-root (round-to-double root)
+ (new-root (round-to-float root)
lb ub)))))
(values worst-diff new-extrema maximin-distance roots)))
@@ -110,6 +132,6 @@
(defun find-error-extrema (coefs points &optional round)
(multiple-value-bind (delta-f ddelta-f d2delta-f)
- (approx-error *loc-value* *loc-dvalue* *loc-d2value* coefs round)
+ (poly-approx-error *loc-value* *loc-dvalue* *loc-d2value* coefs round)
(declare (ignore delta-f))
(%approximation-error-extrema (make-poly coefs) ddelta-f d2delta-f points)))
View
6 demo/branch-and-cut-fit/linf-fit.lisp
@@ -1,10 +1,10 @@
-(declaim (type (integer 1) *dimension*))
+(declaim (type unsigned-byte *dimension*))
(defvar *dimension* 6)
(defun powers (x &optional (degree (1- *dimension*)))
- (let ((x (float x 1d0)))
+ (let ((x (floatify x)))
(coerce (loop for i upto degree
- collect (round-to-double (expt x i)))
+ collect (round-to-float (expt x i)))
'simple-vector)))
(defvar *loc-parameters* 'powers)
View
14 demo/branch-and-cut-fit/newton.lisp
@@ -9,8 +9,8 @@
(defun initial-newton (lo hi f df)
(declare (optimize debug))
(let* ((*precision* 64)
- (lo (round-to-double lo))
- (hi (round-to-double hi))
+ (lo (round-to-float lo))
+ (hi (round-to-float hi))
(slo (signum (pull-bits (funcall f lo))))
(shi (signum (pull-bits (funcall f hi))))
(x (/ (+ lo hi) 2)))
@@ -49,8 +49,8 @@
(return lo))
((zerop shi)
(return hi))
- ((<= (abs (- (double-float-bits lo)
- (double-float-bits hi)))
+ ((<= (abs (- (float-bits lo)
+ (float-bits hi)))
1)
(return x)))
(loop for i upfrom 1 do
@@ -77,13 +77,13 @@
(setf x (if newtonp
(pull-bits x2)
(/ (+ lo hi) 2)))
- (when (<= (abs (- (double-float-bits lo)
- (double-float-bits hi)))
+ (when (<= (abs (- (float-bits lo)
+ (float-bits hi)))
1)
(return x))))))))
(defun newton (lo hi f df)
(multiple-value-bind (lo hi)
- (initial-newton (round-to-double lo) (round-to-double hi) f df)
+ (initial-newton (round-to-float lo) (round-to-float hi) f df)
(and lo hi
(bracketed-newton lo hi f df))))
View
226 demo/branch-and-cut-fit/print-pareto.lisp
@@ -0,0 +1,226 @@
+(defvar *points*)
+
+(defun eval-poly (coefs)
+ (let ((max 0)
+ (poly (make-poly coefs t)))
+ (map nil (lambda (point)
+ (let ((error (abs (- (funcall poly (point-loc point))
+ (point-value point)))))
+ (setf max (max max error))))
+ *points*)
+ max))
+
+(defstruct approx
+ coefs
+ error
+ lg-error
+ degree
+ non-zero
+ non-one
+ non-two
+ id)
+
+(defun find-non-dominated (approxs)
+ (let ((good (make-array (length approxs)
+ :fill-pointer 0)))
+ (map nil (lambda (x)
+ (unless (or (null x)
+ (some (lambda (y)
+ (and (not (eql x y)) y
+ (let ((cmp
+ (list (- (approx-error y)
+ (approx-error x))
+ (- (approx-degree y)
+ (approx-degree x))
+ (- (approx-non-zero y)
+ (approx-non-zero x))
+ (- (approx-non-one y)
+ (approx-non-one x))
+ (- (approx-non-two y)
+ (approx-non-two x)))))
+ (or (and (notany #'plusp cmp)
+ (some #'minusp cmp))
+ (and (every #'zerop cmp)
+ (< (approx-id y)
+ (approx-id x)))))))
+ approxs))
+ (vector-push-extend x good)))
+ approxs)
+ (coerce good 'simple-vector)))
+
+(defun coefs-pareto (coefs from to npoints function)
+ (let* ((*loc-value* function)
+ (*points* (chebyshev-nodes from to npoints))
+ (seen (make-hash-table :test #'equalp)))
+ (find-non-dominated
+ (map 'simple-vector
+ (lambda (coefs)
+ (let* ((coefs (car coefs))
+ (key (map '(simple-array double-float 1)
+ (lambda (x) (float x 1d0)) coefs)))
+ (unless (gethash key seen)
+ (setf (gethash key seen) t)
+ (let ((error (float (eval-poly coefs) 1d0)))
+ (make-approx
+ :coefs coefs
+ :error error
+ :lg-error (if (zerop error)
+ sb-ext:double-float-positive-infinity
+ (floor (- (log error 2d0))))
+ :degree (coefs-degree coefs)
+ :non-zero (count 0 coefs :test-not #'eql)
+ :non-one (count-if-not (lambda (x)
+ (member x '(-1 0 1)))
+ coefs)
+ :non-two (count-if-not (lambda (x)
+ (member x '(-2 -1 0 1 2)))
+ coefs)
+ :id (hash-table-count seen))))))
+ coefs))))
+
+(defun sort-by (values &rest accessors)
+ (let ((values (if (simple-vector-p values)
+ (copy-seq values)
+ (coerce values 'simple-vector))))
+ (dolist (accessor (reverse accessors) values)
+ (setf values (stable-sort values #'< :key accessor)))))
+
+(defvar *accessors* (list (cons #'approx-error "error")
+ (cons (lambda (x)
+ (- (approx-lg-error x)))
+ "lb_error")
+ (cons (lambda (x)
+ (1- (approx-degree x)))
+ "degree")
+ (cons #'approx-non-zero "non_zero")
+ (cons #'approx-non-one "non_one")
+ (cons #'approx-non-two "non_two")))
+
+(defvar *default-accessors* (list (cons (lambda (x)
+ (- (approx-lg-error x)))
+ "lb_error")
+ (cons (lambda (x)
+ (1- (approx-degree x)))
+ "degree")
+ (cons #'approx-error "error")
+ (cons #'approx-non-zero "non_zero")
+ (cons #'approx-non-one "non_one")
+ (cons #'approx-non-two "non_two")))
+
+(defvar *default-accessors2* (list (cons (lambda (x)
+ (1- (approx-degree x)))
+ "degree")
+ (cons (lambda (x)
+ (- (approx-lg-error x)))
+ "lb_error")
+ (cons #'approx-non-zero "non_zero")
+ (cons #'approx-non-one "non_one")
+ (cons #'approx-non-two "non_two")
+ (cons #'approx-error "error")))
+
+
+(defun md5 (x)
+ #+sbcl
+ (let ((sum (sb-md5:md5sum-string x)))
+ (format nil "~{~2,'0X~}" (coerce sum 'list)))
+ #-sbcl
+ (sxhash x))
+
+(defun md5-coefs (coefs)
+ #+sbcl
+ (ecase *float-mode*
+ (double-float
+ (let* ((coefs (map '(simple-array double-float 1) (lambda (x)
+ (float x 1d0))
+ coefs))
+ (nbytes (* 8 (length coefs)))
+ (bytes (make-array nbytes :element-type '(unsigned-byte 8))))
+ (sb-impl::ub8-bash-copy coefs 0 bytes 0 nbytes)
+ (let ((sum (sb-md5:md5sum-sequence bytes)))
+ (format nil "~{~2,'0X~}" (coerce sum 'list)))))
+ (single-float
+ (let* ((coefs (map '(simple-array single-float 1) (lambda (x)
+ (float x 1s0))
+ coefs))
+ (nbytes (* 4 (length coefs)))
+ (bytes (make-array nbytes :element-type '(unsigned-byte 8))))
+ (sb-impl::ub8-bash-copy coefs 0 bytes 0 nbytes)
+ (let ((sum (sb-md5:md5sum-sequence bytes)))
+ (format nil "~{~2,'0X~}" (coerce sum 'list))))))
+ #-sbcl
+ (sxhash coefs))
+
+(defun print-index (function long-description approx accessors)
+ (let ((approx (apply 'sort-by approx (mapcar #'car accessors)))
+ (accessors (mapcar #'car accessors))
+ (names (mapcar #'cdr accessors))
+ (*read-default-float-format* 'double-float))
+ (with-open-file (s (format nil "~A~{-~A~}"
+ function names)
+ :direction :output
+ :if-exists :supersede)
+ (format s "~A ~A~%" function long-description)
+ (dolist (name names)
+ (if (equal name "error")
+ (format s "~24,A " "error")
+ (format s "~8,A " name)))
+ (format s "| coefficients | rationals | hash~%")
+ (map nil (lambda (approx)
+ (loop for accessor in accessors
+ for name in names
+ do (let ((x (funcall accessor approx)))
+ (when (equal name "lb_error")
+ (setf x (- x)))
+ (cond ((eql x sb-ext:double-float-positive-infinity)
+ (format s "~24,A" "inf"))
+ ((floatp x)
+ (format s "~24,E " x))
+ (t
+ (format s "~8,A " x)))))
+ (let ((coefs (approx-coefs approx))
+ (*read-default-float-format* *float-mode*))
+ (setf coefs (subseq coefs 0 (coefs-degree coefs)))
+ (format s "| ~{~W ~}| ~{~A ~}| ~A-~A~%"
+ (map 'list (lambda (x)
+ (floatify x))
+ coefs)
+ (coerce coefs 'list)
+ function
+ (md5-coefs coefs))))
+ approx))))
+
+(defun print-all-indices (approx function long-description
+ &key (max-error sb-ext:double-float-positive-infinity))
+ (let ((approx (remove-if (lambda (x)
+ (> (approx-error x) max-error))
+ approx))
+ (seen (make-hash-table :test #'equalp)))
+ (let ((error (car *accessors*)))
+ (alexandria:map-permutations
+ (lambda (accessors)
+ (let* ((split (1+ (position error accessors)))
+ (first (subseq accessors 0 split))
+ (second (subseq accessors split)))
+ (setf second
+ (sort second #'<
+ :key (lambda (x)
+ (position x *accessors*))))
+ (let* ((accessors (append first second))
+ (key (mapcar #'cdr accessors)))
+ (unless (gethash key seen)
+ (setf (gethash key seen) t)
+ (print-index
+ function long-description
+ approx accessors)))))
+ *accessors*))))
+
+(defun print-default-index (approx function long-description
+ &key (max-error
+ sb-ext:double-float-positive-infinity))
+ (let ((approx (remove-if (lambda (x)
+ (> (approx-error x) max-error))
+ approx)))
+ (print-index function long-description
+ approx *default-accessors*)
+ (print-index function long-description
+ approx *default-accessors2*)))
View
278 demo/branch-and-cut-fit/query.lisp
@@ -0,0 +1,278 @@
+CL-USER> (defparameter *worker*
+ (sb-thread:make-thread
+ (lambda ()
+ (let ((*standard-output* sb-impl::*stdout*)
+ (*trace-output* sb-impl::*stdout*))
+ #+nil
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 -2 2
+ #'computable-reals:sin-r
+ #'computable-reals:cos-r
+ (lambda (x)
+ (computable-reals:-r (computable-reals:sin-r x)))))
+ "/tmp/sin-16")
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 0 1
+ (lambda (x)
+ (computable-reals:log-r (computable-reals:+r x 1)))
+ (lambda (x)
+ (let ((x (computable-reals:+r x 1)))
+ (computable-reals:/R x)))
+ (lambda (x)
+ (let ((x (computable-reals:+r x 1)))
+ (computable-reals:-r
+ (computable-reals:/R
+ (computable-reals:*r x x)))))))
+ "/tmp/log1px-16")
+ #+nil
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 1 2
+ (lambda (x)
+ (computable-reals:log-r x))
+ (lambda (x)
+ (computable-reals:/R x))
+ (lambda (x)
+ (computable-reals:-r
+ (computable-reals:/R
+ (computable-reals:*r x x))))))
+ "/tmp/log-16")
+ #+nil
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 0 1
+ #'computable-reals:exp-r
+ #'computable-reals:exp-r
+ #'computable-reals:exp-r))
+ "/tmp/exp-16")))))
+
+*WORKER*
+
+CL-USER> (defparameter *worker2*
+ (sb-thread:make-thread
+ (lambda ()
+ (sb-thread:join-thread *worker*)
+ (let ((*standard-output* sb-impl::*stdout*)
+ (*trace-output* sb-impl::*stdout*))
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16
+ (- (double-float-from-bits (double-float-bits (/ pi 2)) 2))
+ (double-float-from-bits (double-float-bits (/ pi 2)) 2)
+ #'computable-reals:sin-r
+ #'computable-reals:cos-r
+ (lambda (x)
+ (computable-reals:-r (computable-reals:sin-r x)))))
+ "/tmp/sin-narrow-16")
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 -1 1
+ #'computable-reals:exp-r
+ #'computable-reals:exp-r
+ #'computable-reals:exp-r))
+ "/tmp/exp-wide-16")))))
+*WORKER2*
+CL-USER> (defparameter *worker*
+ (sb-thread:make-thread
+ (lambda ()
+ (let ((*standard-output* sb-impl::*stdout*)
+ (*trace-output* sb-impl::*stdout*))
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 -1 1
+ #'computable-reals:atan-r
+ (lambda (x)
+ (computable-reals:/r (computable-reals:+r
+ 1
+ (computable-reals:*r x x))))
+ (lambda (x)
+ (let ((x^2+1 (computable-reals:+r
+ 1 (computable-reals:*r x x))))
+ (computable-reals:-r
+ (computable-reals:/r (computable-reals:*r 2 x)
+ x^2+1))))))
+ "/tmp/arctan-16")))))
+*WORKER*
+CL-USER> (defparameter *worker*
+ (sb-thread:make-thread
+ (lambda ()
+ (let ((*standard-output* sb-impl::*stdout*)
+ (*trace-output* sb-impl::*stdout*))
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 0 1
+ #'computable-reals:atan-r
+ (lambda (x)
+ (computable-reals:/r (computable-reals:+r
+ 1
+ (computable-reals:*r x x))))
+ (lambda (x)
+ (let ((x^2+1 (computable-reals:+r
+ 1 (computable-reals:*r x x))))
+ (computable-reals:-r
+ (computable-reals:/r (computable-reals:*r 2 x)
+ x^2+1))))))
+ "/tmp/arctan-2-16")))))
+*WORKER*
+CL-USER> (defparameter *worker*
+ (sb-thread:make-thread
+ (lambda ()
+ (let ((*standard-output* sb-impl::*stdout*)
+ (*trace-output* sb-impl::*stdout*))
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 0 1
+ (lambda (x)
+ (let ((x (computable-reals:+r 1 x)))
+ (computable-reals:log-r x 2)))
+ (lambda (x)
+ (let ((x (computable-reals:+r 1 x)))
+ (computable-reals:/r
+ (computable-reals:*r (computable-reals:log-r 2)
+ x))))
+ (lambda (x)
+ (let ((x (computable-reals:+r 1 x)))
+ (computable-reals:-r
+ (computable-reals:/r
+ (computable-reals:*r (computable-reals:log-r 2)
+ x x)))))))
+ "/tmp/lg1px-16")))))
+*WORKER*
+CL-USER> (Setf *float-mode* 'single-float)
+SINGLE-FLOAT
+CL-USER> (defparameter *worker*
+ (sb-thread:make-thread
+ (lambda ()
+ (let ((*standard-output* sb-impl::*stdout*)
+ (*trace-output* sb-impl::*stdout*))
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 0 1
+ (lambda (x)
+ (computable-reals:log-r (computable-reals:+r x 1)))
+ (lambda (x)
+ (let ((x (computable-reals:+r x 1)))
+ (computable-reals:/R x)))
+ (lambda (x)
+ (let ((x (computable-reals:+r x 1)))
+ (computable-reals:-r
+ (computable-reals:/R
+ (computable-reals:*r x x)))))))
+ "/tmp/single/log1px-16")
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 1 2
+ (lambda (x)
+ (computable-reals:log-r x))
+ (lambda (x)
+ (computable-reals:/R x))
+ (lambda (x)
+ (computable-reals:-r
+ (computable-reals:/R
+ (computable-reals:*r x x))))))
+ "/tmp/single/log1px-16")
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16
+ (- (float-from-bits (float-bits (/ pi 2)) 2))
+ (float-from-bits (float-bits (/ pi 2)) 2)
+ #'computable-reals:sin-r
+ #'computable-reals:cos-r
+ (lambda (x)
+ (computable-reals:-r (computable-reals:sin-r x)))))
+ "/tmp/single/sin-narrow-16")
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16
+ (- (float-from-bits (float-bits (/ pi 2)) 2))
+ (float-from-bits (float-bits (/ pi 2)) 2)
+ #'computable-reals:cos-r
+ (lambda (x)
+ (computable-reals:-r (computable-reals:sin-r x)))
+ (lambda (x)
+ (computable-reals:-r (computable-reals:cos-r x)))))
+ "/tmp/single/cos-narrow-16")
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 -1 1
+ #'computable-reals:exp-r
+ #'computable-reals:exp-r
+ #'computable-reals:exp-r))
+ "/tmp/single/exp-wide-16")
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 0 1
+ #'computable-reals:exp-r
+ #'computable-reals:exp-r
+ #'computable-reals:exp-r))
+ "/tmp/single/exp-16")
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 0 1
+ #'computable-reals:atan-r
+ (lambda (x)
+ (computable-reals:/r (computable-reals:+r
+ 1
+ (computable-reals:*r x x))))
+ (lambda (x)
+ (let ((x^2+1 (computable-reals:+r
+ 1 (computable-reals:*r x x))))
+ (computable-reals:-r
+ (computable-reals:/r (computable-reals:*r 2 x)
+ x^2+1))))))
+ "/tmp/single/arctan-16")
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 -1 1
+ #'computable-reals:atan-r
+ (lambda (x)
+ (computable-reals:/r (computable-reals:+r
+ 1
+ (computable-reals:*r x x))))
+ (lambda (x)
+ (let ((x^2+1 (computable-reals:+r
+ 1 (computable-reals:*r x x))))
+ (computable-reals:-r
+ (computable-reals:/r (computable-reals:*r 2 x)
+ x^2+1))))))
+ "/tmp/single/arctan-wide-16")
+
+ (print-hash-table
+ (time
+ (enumerate-pareto-front
+ 16 0 1
+ (lambda (x)
+ (let ((x (computable-reals:+r 1 x)))
+ (computable-reals:log-r x 2)))
+ (lambda (x)
+ (let ((x (computable-reals:+r 1 x)))
+ (computable-reals:/r
+ (computable-reals:*r (computable-reals:log-r 2)
+ x))))
+ (lambda (x)
+ (let ((x (computable-reals:+r 1 x)))
+ (computable-reals:-r
+ (computable-reals:/r
+ (computable-reals:*r (computable-reals:log-r 2)
+ x x)))))))
+ "/tmp/single/lg1px-16")))))
+*WORKER*
View
44 demo/branch-and-cut-fit/utility.lisp
@@ -1,4 +1,6 @@
-(defun double-float-bits (x)
+(defvar *float-mode* 'double-float)
+
+(defun %double-float-bits (x)
(let* ((x (float x 1d0))
(bits (logior (ash (sb-kernel:double-float-high-bits x) 32)
(sb-kernel:double-float-low-bits x))))
@@ -7,7 +9,7 @@
(- (ldb (byte 63 0) bits))
bits)))
-(defun double-float-from-bits (x &optional delta)
+(defun %double-float-from-bits (x &optional delta)
;; revert into sign/magnitude
(when delta (incf x delta))
(when (minusp x)
@@ -16,19 +18,45 @@
(sb-kernel:make-double-float (ash x -32)
(ldb (byte 32 0) x)))
-(defun round-to-double (x)
- (rational (float x 1d0)))
-(defun split-double (x)
- (let ((dx (round-to-double x)))
+(defun float-bits (x)
+ (ecase *float-mode*
+ (double-float
+ (%double-float-bits x))
+ (single-float
+ (let ((bits (sb-kernel:single-float-bits (float x 1s0))))
+ (if (logbitp 31 bits)
+ (- (ldb (byte 31 0) bits))
+ bits)))))
+
+(defun float-from-bits (x &optional delta)
+ (ecase *float-mode*
+ (double-float
+ (%double-float-from-bits x delta))
+ (single-float
+ (when delta
+ (incf x delta))
+ (when (minusp x)
+ (setf x (logior (ash -1 31)
+ (ldb (byte 63 0) (- x)))))
+ (sb-kernel:make-single-float x))))
+
+(defun round-to-float (x)
+ (rational (coerce x *float-mode*)))
+
+(defun floatify (x)
+ (coerce x *float-mode*))
+
+(defun split-float (x)
+ (let ((dx (round-to-float x)))
(cond ((= x dx)
nil)
((< x dx)
- (values (double-float-from-bits (1- (double-float-bits dx)))
+ (values (float-from-bits (1- (float-bits dx)))
dx))
(t
(values dx
- (double-float-from-bits (1+ (double-float-bits dx))))))))
+ (float-from-bits (1+ (float-bits dx))))))))
(defvar *precision* 64)
(defvar *max-accuracy* 4096)

0 comments on commit b3c9f29

Please sign in to comment.
Something went wrong with that request. Please try again.