Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
;;; 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

This comment has been minimized.

Copy link

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?

@marcoxa

This comment has been minimized.

Copy link

marcoxa commented Sep 18, 2019

Hi borgauf

it is at the REPL (or so it seems). The FRACTIONAL-SET-COVER> is the prompt; some implementations change it to reflect the current package. The PROG1 is to ensure you get NIL as a result the rest is just invoking the solver.

MA

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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.