Skip to content

Instantly share code, notes, and snippets.

# pkhuong/fractional-set-cover.lisp Created Apr 23, 2019

 ;;; Fractional set cover solver with adahedge for the dual surrogate multiplers. (defpackage "FRACTIONAL-SET-COVER" (:use "CL")) (in-package "FRACTIONAL-SET-COVER") (deftype dvec () `(simple-array double-float 1)) (defstruct tour (index (error "index must be provided") :type (and unsigned-byte fixnum) :read-only t) (cost (error "cost must be provided") :type double-float :read-only t) ;; offset by 2, so loc #2 = 0 (locations (error "locations must be provided") :type (simple-array (unsigned-byte 16) 1) :read-only t) (weight-in-average-constraint 0d0 :type double-float)) (defun sort-tours (tours multipliers) (declare (type (simple-array tour 1) tours) (type dvec multipliers) (optimize speed)) (map nil (lambda (tour) (let ((weight 0d0)) (declare (type double-float weight)) (map nil (lambda (loc) (incf weight (aref multipliers loc))) (tour-locations tour)) (setf (tour-weight-in-average-constraint tour) weight))) tours) (flet ((comparator (x y) (declare (type tour x y)) ;; We want to return T when ;; x-weight / x-cost > y-weight / y-cost. ;; Avoid the division with ;; x-weight * y-cost > y-weight * x-cost (> (* (tour-weight-in-average-constraint x) (tour-cost y)) (* (tour-weight-in-average-constraint y) (tour-cost x))))) (declare (inline sort)) (sort tours #'comparator))) (defun find-lower-bound (sorted-tours) (declare (type (simple-array tour 1) sorted-tours)) (let ((total-cost 0d0) (coverage 0d0)) (declare (type double-float total-cost coverage)) (loop for tour across sorted-tours do (let ((cost (tour-cost tour)) (weight (tour-weight-in-average-constraint tour))) (when (>= (+ coverage weight) 1d0) (incf total-cost (* cost (/ (- 1d0 coverage) (+ weight 1d-8)))) ;; The right-hand side is 1, and, when multiplied by the ;; lagrangian multiplier, the result must be equal to this ;; lower bound. (return-from find-lower-bound (values total-cost total-cost))) (incf total-cost cost) (incf coverage weight))) (format t "infeasible?~%") (values total-cost 10d0))) (defun find-solution (sorted-tours max-cost num-locs) (declare (type (simple-array tour 1) sorted-tours) (type double-float max-cost)) (let ((solution (make-array (length sorted-tours) :element-type 'double-float :initial-element 0d0)) (satisfaction (make-array num-locs :element-type 'double-float :initial-element -1d0)) (total-cost 0d0) (coverage -1d0)) (declare (type double-float total-cost coverage)) (flet ((observe-value (tour coefficient) (declare (type tour tour) (type (double-float 0d0 1d0) coefficient)) (incf (aref solution (tour-index tour)) coefficient) (incf total-cost (* (tour-cost tour) coefficient)) (incf coverage (* (tour-weight-in-average-constraint tour) coefficient)) (map nil (lambda (loc) (incf (aref satisfaction loc) coefficient)) (tour-locations tour)))) (loop for tour across sorted-tours do (let ((cost (tour-cost tour))) (when (>= (+ total-cost cost) max-cost) (observe-value tour (/ (- max-cost total-cost) cost)) (return)) (observe-value tour 1d0))) (values solution total-cost satisfaction coverage)))) (defstruct adahedge (log-num-variables 30 :type double-float :read-only t) (total-loss (error "total-loss must be provided") :type dvec :read-only t) (mix-gap 0d0 :type double-float)) (defun make-adahedge-state (num-locations num-tours) (make-adahedge :log-num-variables (log (float num-tours 0d0)) :total-loss (make-array num-locations :element-type 'double-float :initial-element 0d0))) (declaim (type (function (dvec) (values double-float &optional)) dmin dmax dsum)) (defun dmin (x) (declare (type dvec x) (optimize speed)) (let ((min sb-ext:double-float-positive-infinity)) (declare (type double-float min)) (map nil (lambda (x) (setf min (min x min))) x) (+ 0d0 min))) (defun dmax (x) (declare (type dvec x) (optimize speed)) (let ((max sb-ext:double-float-negative-infinity)) (declare (type double-float max)) (map nil (lambda (x) (setf max (max x max))) x) (+ 0d0 max))) (defun dsum (x) (declare (type dvec x) (optimize speed)) (let ((sum 0d0)) (declare (type double-float sum)) (map nil (lambda (x) (incf sum x)) x) (+ 0d0 sum))) (declaim (ftype (function (dvec double-float) (values dvec &optional)) dscalen)) (defun dscalen (x scale) (declare (type dvec x) (type double-float scale) (optimize speed)) (map-into x (lambda (x) (* scale x)) x)) (declaim (ftype (function (dvec dvec) (values dvec &optional)) dincn)) (defun dincn (dst increment) (declare (type dvec dst increment) (optimize speed)) (map-into dst #'+ dst increment)) ;; https://arxiv.org/pdf/1301.0534.pdf (defun mix (step-size total-loss) (declare (type double-float step-size) (type dvec total-loss)) (let* ((min-loss (dmin total-loss)) (weights (if (>= step-size sb-ext:double-float-positive-infinity) (map 'dvec (lambda (loss) (if (= loss min-loss) 1d0 0d0)) total-loss) (map 'dvec (lambda (loss) (exp (- (* step-size (- loss min-loss))))) total-loss))) (norm (dsum weights))) (values (dscalen weights (/ norm)) (- min-loss (/ (log (/ norm (length total-loss))) step-size))))) (defun adahedge-step (state loss-function) (declare (type adahedge state) (type (function (dvec) (values double-float dvec &optional)) loss-function)) (let ((step-size (if (zerop (adahedge-mix-gap state)) sb-ext:double-float-positive-infinity (/ (adahedge-log-num-variables state) (adahedge-mix-gap state))))) (multiple-value-bind (weights prev-mix-loss) (mix step-size (adahedge-total-loss state)) (multiple-value-bind (observed-loss observed-component-loss) (funcall loss-function weights) (dincn (adahedge-total-loss state) observed-component-loss) ;; XXX don't cons the weight vector (let ((mix-loss (nth-value 1 (mix step-size (adahedge-total-loss state))))) (incf (adahedge-mix-gap state) (max 0d0 (- observed-loss (- mix-loss prev-mix-loss))))) ;; XXX: also return count at dmin. (dmin (adahedge-total-loss state)))))) (defstruct set-cover (adahedge (error "adahedge must be provided") :type adahedge :read-only t) (tours (error "tours must be provided") :type (simple-array tour 1) :read-only t) (best-multipliers nil :type (or null dvec)) (num-solutions 0 :type (and fixnum unsigned-byte)) (sum-solutions (error "sum-solutions must be provided") :type dvec :read-only t) (sum-cost 0d0 :type double-float) (lower-bound 0d0 :type double-float)) (defun make-set-cover-state (tours num-locations) (make-set-cover :adahedge (make-adahedge-state num-locations (1+ (length tours))) :tours (copy-seq tours) :sum-solutions (make-array (length tours) :element-type 'double-float :initial-element 0d0))) (declaim (ftype (function (set-cover dvec) (values double-float dvec &optional)) set-cover-loss-function)) (defun set-cover-loss-function (state multipliers) (declare (type set-cover state) (type dvec multipliers) (values double-float dvec &optional)) (let ((tours (sort-tours (set-cover-tours state) multipliers))) (multiple-value-bind (new-lower-bound dual-scale) (find-lower-bound tours) (let* ((lower-bound (max (set-cover-lower-bound state) new-lower-bound)) (max-cost (max lower-bound (- (* (1+ (set-cover-num-solutions state)) lower-bound) (* (set-cover-sum-cost state) (1+ 1d-8)))))) (multiple-value-bind (solution cost satisfaction average-satisfaction) (find-solution tours max-cost (length multipliers)) (declare (type dvec solution satisfaction) (type double-float cost average-satisfaction)) (incf (set-cover-num-solutions state)) (dincn (set-cover-sum-solutions state) solution) (incf (set-cover-sum-cost state) cost) (assert (<= (/ (set-cover-sum-cost state) (set-cover-num-solutions state)) lower-bound)) (when (>= new-lower-bound (set-cover-lower-bound state)) (setf (set-cover-best-multipliers state) (dscalen (copy-seq multipliers) dual-scale)) (when (> new-lower-bound (set-cover-lower-bound state)) (format t "LB: ~A ~A~%" lower-bound dual-scale))) (setf (set-cover-lower-bound state) lower-bound) (values average-satisfaction satisfaction)))))) (defun set-cover-iteration (state) (declare (type set-cover state)) (values (/ (adahedge-step (set-cover-adahedge state) (lambda (multipliers) (set-cover-loss-function state multipliers))) (set-cover-num-solutions state)) (/ (set-cover-sum-cost state) (set-cover-num-solutions state)) (set-cover-lower-bound state))) (defun set-cover-solution (state) (declare (type set-cover state)) (values (dscalen (copy-seq (set-cover-sum-solutions state)) (/ 1d0 (set-cover-num-solutions state))) (/ (set-cover-sum-cost state) (set-cover-num-solutions state)) (set-cover-lower-bound state) (set-cover-best-multipliers state))) (defun print-set-cover-solution (file state) (multiple-value-bind (solution cost lower-bound multipliers) (set-cover-solution state) (with-open-file (s file :direction :output :if-does-not-exist :create :if-exists :supersede) (write (list solution cost lower-bound multipliers) :stream s :readably t)))) (defvar *default-tours* (vector (make-tour :index 0 :cost 10d0 :locations (coerce '(0 1 2) '(simple-array (unsigned-byte 16) 1))) (make-tour :index 1 :cost 1d0 :locations (coerce '(0 1) '(simple-array (unsigned-byte 16) 1))) (make-tour :index 2 :cost 1d0 :locations (coerce '(1 2) '(simple-array (unsigned-byte 16) 1))))) (defun read-instance (file) (let* ((*read-default-float-format* 'double-float) (sexp (with-open-file (s file) (read s))) (dst (make-array 16 :adjustable t :fill-pointer 0))) (destructuring-bind (capacity weights tours) sexp (format t "capacity: ~A~%" capacity) (loop for i upfrom 0 for (cost locs) across tours do (vector-push-extend (make-tour :index i :cost cost :locations (map '(simple-array (unsigned-byte 16) 1) (lambda (loc) (- loc 2)) locs)) dst)) (values (1- (length weights)) (coerce dst 'simple-vector))))) (defun read-weights (file) (let* ((*read-default-float-format* 'double-float) (sexp (with-open-file (s file) (read s)))) (destructuring-bind (capacity weights tours) sexp (declare (ignore capacity tours)) (map 'simple-vector #'ceiling (subseq weights 1))))) (defun generate-solution (fractional-solution tours weights) (declare (type dvec fractional-solution) (type (simple-array tour 1) tours) (type simple-vector weights)) (let ((capacity (* 10 1000 1000)) (total-weight (reduce #'+ weights)) (current-weight 0) (covered (make-array (length weights) :element-type 'bit :initial-element 0)) (to-use (remove-if #'zerop (map 'simple-vector #'cons fractional-solution tours) :key #'car)) (random-solution '())) (labels ((collect-unused () (loop for i upfrom 0 for coverage across covered when (zerop coverage) collect i)) (consider-item (probability tour) "Returns T if tour has been taken or should otherwise be removed from to-use." (cond ((zerop probability) t) ((every (lambda (loc) (plusp (aref covered loc))) (tour-locations tour)) t) ((> (random 1d0) probability) nil) (t (let ((new-tour '())) (loop for loc across (tour-locations tour) when (zerop (aref covered loc)) do (setf (aref covered loc) 1) (push loc new-tour) (incf current-weight (aref weights loc))) (push new-tour random-solution)) t))) (one-iter (to-use) (let ((to-reuse (make-array (length to-use) :fill-pointer 0))) (loop for entry across to-use do (when (<= (- total-weight current-weight) capacity) (push (collect-unused) random-solution) (return-from generate-solution random-solution)) (unless (consider-item (car entry) (cdr entry)) (vector-push entry to-reuse)) finally (return to-reuse))))) (loop (setf to-use (one-iter to-use)))))) (defun print-random-solution-as-tour (solution file) (with-open-file (s file :direction :output :if-exists :supersede :if-does-not-exist :create) (loop for tour in solution do (format s "1~%~{~A~%~}" (mapcar (lambda (loc) (+ loc 2)) tour))) (format s "-1~%"))) #|| FRACTIONAL-SET-COVER> (prog1 nil (setf (values *num-locs* *tours*) (read-instance "/Users/pkhuong/cvrp-santa/tours-2018-12-24T11:16:14-05:00.sexp"))) (defparameter *state* (make-set-cover-state *tours* *num-locs*)) FRACTIONAL-SET-COVER> (dotimes (i 200) (print (multiple-value-list (time (dotimes (i 100 (progn (sb-ext:gc :full t) (set-cover-iteration *state*))) (set-cover-iteration *state*)))))) FRACTIONAL-SET-COVER> (prog1 nil (setf *weights* (read-weights "/Users/pkhuong/cvrp-santa/tours-2018-12-24T11:16:14-05:00.sexp"))) ; in: PROG1 () ; (SETF FRACTIONAL-SET-COVER::*WEIGHTS* ; (FRACTIONAL-SET-COVER::READ-WEIGHTS ; "/Users/pkhuong/cvrp-santa/tours-2018-12-24T11:16:14-05:00.sexp")) ; ==> ; (SETQ FRACTIONAL-SET-COVER::*WEIGHTS* ; (FRACTIONAL-SET-COVER::READ-WEIGHTS ; "/Users/pkhuong/cvrp-santa/tours-2018-12-24T11:16:14-05:00.sexp")) ; ; caught WARNING: ; undefined variable: *WEIGHTS* ; ; compilation unit finished ; Undefined variable: ; *WEIGHTS* ; caught 1 WARNING condition NIL FRACTIONAL-SET-COVER> (defparameter *fractional-solution* (set-cover-solution *state*)) *FRACTIONAL-SOLUTION* ; compiling (DEFUN PRINT-RANDOM-SOLUTION-AS-TOUR ...) FRACTIONAL-SET-COVER> (print-random-solution-as-tour *random-solution* "/tmp/random-solution.tour") NIL FRACTIONAL-SET-COVER> (defparameter *random-solution* (generate-solution (set-cover-solution *state*) (set-cover-tours *state*) *weights*)) *RANDOM-SOLUTION* FRACTIONAL-SET-COVER> (print-random-solution-as-tour *random-solution* "/tmp/random-solution.tour") NIL ; compiling (DEFUN PRINT-RANDOM-SOLUTION-AS-TOUR ...) FRACTIONAL-SET-COVER> (print-random-solution-as-tour *random-solution* "/tmp/random-solution.tour") NIL FRACTIONAL-SET-COVER> (dotimes (i 100) (print-random-solution-as-tour *random-solution* (format nil "/tmp/random-solution-~A.tour" i))) NIL FRACTIONAL-SET-COVER> (dotimes (i 100) (print-random-solution-as-tour (generate-solution *fractional-solution* (set-cover-tours *state*) *weights*) (format nil "/tmp/random-solution-~A.tour" i))) ||#

### borgauf commented May 1, 2019

 ``````FRACTIONAL-SET-COVER> (prog1 nil (setf (values *num-locs* *tours*) (read-instance "/Users/pkhuong/cvrp-santa/tours-2018-12-24T11:16:14-05:00.sexp"))) `````` My being a beginner, what is happening here? What is the .sexp file? I take it this is inside a Lisp REPL?
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.