Seidel's Linear Programming Algorithm

seidel-lp.rkt

#lang racket

; A non-random Racket implementation of the algorithm given in:
;   Small-Dimension Linear Programming and Convex Hulls Made Easy by R. Seidel

; There are a number of other documents floating around describing it,
; including another Seidel paper. I think all of them have at least typos,
; if not conceptual errors. The pseudo code at the end of the above paper
; is, I think, the least buggy. All the λ stuff is a bit odd; they seem
; to be mixing some numerical tolerance stuff with the variable bounds?
; All the logic seems to be correct, which is the important part.

; This code isn't quite right, either. Linear programming is hard!

; - Jay Kominek, September 2011

(define vref vector-ref)

(define (vector-drop-idx vec idx)
  (for/vector ([v vec]
               [i (in-naturals)]
               #:when (not (= i idx)))
              v))

(define (max-nonzero-idx vec limit)
  (let-values ([(idx _)
                (for/fold ([idx -1]
                           [max-v #f])
                  ([v vec]
                   [i (in-naturals)]
                   #:when (and (< i limit) (not (= v 0))))
                  (if (or (< idx 0) (> v max-v))
                      (values i v)
                      (values idx max-v)))])
    (if (>= idx 0)
        idx
        #f)))

(define (LP-base-case c l u A)
  (let-values ([(high low z)
                (for/fold ([high (vector-ref u 0)]
                           [low  (vector-ref l 0)]
                           [z    0])
                  ([a A])
                  (let ([a0 (vref a 0)])
                    (cond
                      [(> a0 0) (values (min high (/ (vref a 1) a0))
                                        low
                                        z)]
                      [(< a0 0) (values high
                                        (max low (/ (vref a 1) a0))
                                        z)]
                      [else     (values high
                                        low
                                        (min z (vref a 1)))])))])
    (if (and (>= z 0)
             (> high low))
        (if (>= (vref c 0) 0)
            (vector high)
            (vector low))
        #f)))

(define (LP-unconstrained-case c l u)
  (for/vector ([c_i c]
               [l_i l]
               [u_i u])
              (if (>= c_i 0) u_i l_i)))
                  
(define (LP-inductive d c l u A)
  ; try throwing away a constraint and see if the solution to
  ; the reduced case also works for us.
  (let* ([removed-constraint (car A)]
         [remaining-constraints (cdr A)]
         [x* (LP d c l u remaining-constraints)])
    (cond
      [(not x*) #f]
      [(<= (for/fold ([sum 0])
             ([x_i x*]
              [a_i removed-constraint])
             (+ sum (* x_i a_i)))
           (vref removed-constraint d))
       x*]
      [else
       ; reduced solution doesn't satisfy o ur one last constraint, fixup
       (LP-inductive-reduce-dimension x* d c l u removed-constraint remaining-constraints)])))

(define (LP-inductive-reduce-dimension x* d c l u A0 A)
  (let ([k (max-nonzero-idx A0 d)])
    (if (not k)
        #f
        (let* ([f (for/vector ([i (in-naturals)]
                               [a_i A0]
                               #:when (not (= i k)))
                    (cond [(= i d) (vref u k)]
                          [else (- (/ a_i (vref A0 k)))]))]
               [g (for/vector ([i (in-naturals)]
                               [a_i A0]
                               #:when (not (= i k)))
                    (cond [(= i d) (- (vref l k))]
                          [else (- (/ a_i (vref A0 k)))]))]
               [Ap (for/list ([b A])
                     (for/vector ([b_i b]
                                  [a_i A0]
                                  [i (in-naturals)]
                                  #:when (not (= i k)))
                       (- b_i (* (/ (vref b k) (vref A0 k)) a_i))))]
               [cp (for/vector ([c_i c]
                                [a_i A0]
                                [i (in-naturals)]
                                #:when (not (= i k)))
                     (- c_i (* (/ (vref c k) (vref A0 k)) a_i)))]
               [xp (LP (sub1 d)
                       cp
                       (for/vector ([l_i l] [i (in-naturals)] #:when (not (= i k))) l_i)
                       (for/vector ([u_i u] [i (in-naturals)] #:when (not (= i k))) u_i)
                       (cons f (cons g Ap)))])
          (if (not xp)
              #f
              (let ([x (build-vector d (lambda (i)
                                         (cond [(< i k) (vref xp i)]
                                               [(= i k) 0]
                                               [(> i k) (vref xp (sub1 i))])))])
                (vector-set! x k
                             (/ (- (vref A0 d)
                                   (for/fold ([sum 0])
                                     ([x_i x]
                                      [a_i A0])
                                     (+ sum (* a_i x_i))))
                                (vref A0 k)))
                x))))))

; entry point
;  d is dimensions
;  c is a vector representing the function to maximize
;  l, u are vectors representing the upper and lower bounds for each dimension
;  A is a list of d+1 element vectors representing the constraints
(define (LP d c l u A)
  (cond
    [(= d 1)          (LP-base-case c l u A)]
    [(= (length A) 0) (LP-unconstrained-case c l u)]
    [else             (LP-inductive d c l u A)]))

; random little thing i made up
(LP 2
    #(-2 1)
    #(-10 -10)
    #( 10  10)
    '(#(1 0 5.5)
      #(-1 0 5)
      #(0 -1 6)
      #(0 1 6)
      #(1 1 10)))

; sample problem from the lp_solve docs
(LP 2
    #(143 60)
    #(0 0)
    #(10000 10000)
    '(#(120 210 15000)
      #(110 30 4000)
      #(1 1 75)))

; is there any space inside of a cube?
(LP 3
    #(1 -1 1)
    #(-100 -100 -100)
    #(100 100 100)
    '(#( 1.0  0.0  0.0 0.5)
      #(-1.0  0.0  0.0 0.5)
      #( 0.0  1.0  0.0 0.5)
      #( 0.0 -1.0  0.0 0.5)
      #( 0.0  0.0  1.0 0.5)
      #( 0.0  0.0 -1.0 0.5)))

; shouldn't be feasible
(LP 3
    #(1 1 1)
    #(-100 -100 -100)
    #(100 100 100)
    '(#(-1.0 0.0 0.0 0.49)
      #(1.0 0.0 0.0 -0.51)
      ))