fractional-set-cover.lisp

;;; 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)))
||#

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?

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?