challenge: power set and set partitions

[adorf@eso.org (Hans-Martin Adorf)]

Dear MCLers,

after all the sad news about the not-so-bright future of MCL and all your
encouraging comments to change Apple's mind, you might be prepared for a
real (small) Lisp challenge:

I need to compute all the partitions of a given set (i.e. all subdivisions
of a given set into non-empty sets which cover the original set). Below I
append my quick hack. It does the job; however, it is not very Lispish in
style. 

The challenge consists in implementing the partitions function more
elegantly, preferrably using recursion. As you see my code uses a function
which generates the power set of a given set, a subproblem, for which an
elegant (and efficient) solution is also seeked.

I will summarize the solutions for partition and power-set offered to me.

In case you are interested in the original problem that stimulated the
partitions/power-set problem above: I am engaged in calculations for a
planned  astronomical hybrid interferometer within ESO's Very Large
Telescope project on Mt. Paranal, Chile.

If we will make use of any code forwarded to me, the author will be
properly acknowledged in the planned scientific paper.

Happy MCLing.

Hans-Martin Adorf

----------------- 
Hans-Martin Adorf
ST-ECF/ESO
Karl-Schwarzschild-Str. 2
D-85748 Garching b. Muenchen
Germany
Tel: +49-89-32006-261
Fax: +49-89-32006-480
Internet: adorf@eso.org


================================================================================
;;;;
;;;; partitions.lisp
;;;;
;;;; Compute the power-set and all partitions of an arbitrary set
;;;; 
;;;; Hans-Martin Adorf, ST-ECF, 26-02-94
;;;;

(defvar *offset* 0 "offset value for generation of sets")
(defvar *unique* nil "flag indicating whether only unique partitions should
be returned")

;;;-----------------------------------------------------------------------------
;;; The following code is an application of partitions to the computation
of 
;;; the triple correlation function of an astronomical hybrid
interferometer.
;;; The idea is to generate Mathematica code which can further be evaluated

;;; in Mathematica.
;;;-----------------------------------------------------------------------------
(defun to-file (file-name expr)
  (with-open-file (stream file-name :direction :output :if-exists
:supersede)
    (princ expr stream)))
#|
(to-file #P"HMA_PB_i80:test.math" 
         (format-partitions (partitions '(\u \v \w \x \y \z))))
|#

(defun format-partitions (partitions)
  "Format partitions for Lisp or Mathematica evaluation"
  (plus (mapcar #'format-partition partitions)))
;; (format-partitions (partitions '(\u \v \w)))
;; (format-partitions (partitions '(\u \v \w \x \y \z)))

(defun format-partition (partition)
  (times (mapcar #'(lambda (set) (f (plus set))) partition)))
;; (format-partition '((a) (b c)))

#|
;;; for Lisp evaluation
(defun f (set)
  (list 'f set))

(defun plus (set)
  (if (< (length set) 2)
    set
    (cons '+ set)))

(defun times (set)
  (if (< (length set) 2)
    set
    (cons '* set)))
|#

;;; for Mathematica evaluation
(defun f (set)
  (format nil "f[~a]" set))
;; (f '(a b c))
;; (f (plus '(a b c)))

(defun plus (set)
  (if (< (length set) 2)
    (format nil "~{~a~}" set)
    (format nil "Plus[~{~a~^, ~}]" set)))
;; (plus '(u v))

(defun times (set)
  (if (< (length set) 2)
    (format nil "~{~a~}" set)
    (format nil "Times[~{~a~^, ~}]" set)))
;; (times '(u v w))
          
;;;-----------------------------------------------------------------------------
;;; Compute partitions
;;;-----------------------------------------------------------------------------
(defun partitions (set)
  "Compute all partitions of a general n-element set"
  (let ((result
         (cond ((null set) nil)
               ((= (length set) 1) (list (list set)))
               (t (partitions-aux (base-set set)))
               )))
    (values result (length result))))
;; (partitions '())
;; (partitions '(u))
;; (partitions '(u v))
;; (partitions '(u v w))
;; (partitions '(u v w x))
;; (partitions '(u v w x y))
;; (partitions '(u v w x y z))

(defun base-set (set)
  "Auxiliary function for partitions"
  (let* ((head (first set))
         (tail (rest set))
         (power-set (power-set tail))
         (base-set (mapcar #'(lambda (x) (cons head x)) power-set)))
    (mapcar #'(lambda (x) (list x (set-complement x set))) base-set)))
;; (base-set '(a b c))

(defun partitions-aux (base-set)
  (mapcan #'(lambda (pair)
              (new-partitions (first pair) (partitions (second pair))))
          base-set))
;; (partitions '(a b c))

(defun new-partitions (head partitions)
  (if (null partitions)
    (list (list head))
    (mapcar #'(lambda (x) (append (list head) x)) partitions)))

#|
;;;-----------------------------------------------------------------------------
;; dead code; works, but unused
;;;-----------------------------------------------------------------------------
(defun 2-set-partitions (n)
  "Generate all partitions of an n-element set (of intergers) into 2
subsets"
  (let ((result 
         (mapcan #'(lambda (k) (k-partitions k n))
                 (if *unique*
                   (integers 0 (ceiling (/ (1+ n) 2)))    ; don't compute
duplicates
                   (integers 0 (1+ n))
                   ))))
    (values result (length result))))
;; (2-set-partitions 2)
;; (setf *unique* nil)
;; (2-set-partitions 3)
;; (2-set-partitions 4)

(defun k-partitions (k n)
  "Generate all partitions of an n-element set (of integers) into 2 subsets

   with k elements and n-k elements, respectively"
  (let* ((set (make-set n))
         (result 
          (if (zerop k)
            (list (list nil set))
            (mapcar #'(lambda (subset) 
                        (list subset 
                              (set-complement subset set)))
                    (k-subsets k n)))))
    ;; for a proper partion one must remove duplicates from the result
    (when (and *unique* (= n (* 2 k))) 
      (setf result (half-seq result)))      ; cut half
    (values result (length result))))
;; (k-partitions 1 2)
;; (k-partitions 2 4)
;; (k-partitions 3 5)
;; (k-partitions 0 3)
;; (k-partitions 3 3)
|#

;;;-----------------------------------------------------------------------------
;;; Compute power set
;;;-----------------------------------------------------------------------------
(defun power-set (set)
  "Generate the power-set (i.e. the set of all subsets) for an arbitrary
n-element set"
  (let* ((n (length set))
         (*offset* 0)
         (power-set (power-set-n n))   ; generate power-set for integers
         (result (dotimes (i n power-set)
                   (nsubst (elt set i) i power-set)))) ; replace integers
by set-elements
    (values result (length result))))
;; (power-set '(a b c d))

(defun power-set-n (n)
  "Generate the power-set (i.e. the set of all subsets) for an n-element
set of integers"
  (let ((result (mapcan #'(lambda (k) (k-subsets k n))
                        (integers 0 (1+ n)))))
    (push nil result)                   ; the empty set is missing
    (values result (length result))))
;; (power-set-n 2)
;; (setf *offset* 0)
;; (power-set-n 3)
;; (power-set-n 4)

(defun k-subsets (k n &optional (subsets-so-far (append-elements nil n)))
  "Generate all k-element subsets for a set of n elements"
  (let ((result 
         (if (= (length (first subsets-so-far)) k) 
           subsets-so-far
           (k-subsets k n 
                      (mapcan #'(lambda (subset) (append-elements subset
n))
                              subsets-so-far)))))
    (values result (length result))))
;; (setf *offset* 1)
;; (k-subsets 0 3)
;; (k-subsets 1 3)
;; (k-subsets 2 3)
;; (k-subsets 3 3)
;; (k-subsets 4 7)
;; (k-subsets 2 3 '((0) (1) (2)))
;; (k-subsets 5 7 '((0 3 4) (0 3 5) (0 3 6)))

(defun append-elements (subset n)
  "Append indices up to n + *offset* (exclusively) to a given subset of
indices, 
  e.g. for n + *offset* = 6: (1 3) -> ((1 3 4) (1 3 5))"
  (if (null subset)
    (mapcar #'list (make-set n))      ; return list of parenthesized
integers
    (do ((i (1+ (last-elt subset)) (1+ i))
         (nn (+ n *offset*))
         (result))
        ((= i nn) (reverse result))
      (push (append subset (list i)) result))))
;; (append-elements nil 5)
;; (append-elements '(1 3) 6)
;; (append-elements '(1 4) 4)

;;;-----------------------------------------------------------------------------
;;; Auxiliary functions
;;;-----------------------------------------------------------------------------
(defun half-seq (seq)
  "Return the first half of a sequence, middle element inclusive"
  (subseq seq 0 (ceiling (/ (length seq) 2))))
;; (half-seq '(1 2 3 4))
;; (half-seq '(1 2 3 4 5))

(defun make-set (n)
  "Generate a set of integers"
  (integers *offset* (+ *offset* n)))

(defun integers (m n)
  "Generate a list of integers in the range of m (inclusively) to n
(exclusively)"
  (do ((i m (1+ i))
       result)
      ((= i n) (reverse result))
    (push i result)))
;; (integers 3 7)

(defun set-complement (subset set)
  "MCL returns the reverse of the set-difference"
  (reverse (set-difference set subset)))        ; revert result for
aesthetical reasons
;; (set-complement '(a) '(a b c))

(defun last-elt (seq)
  (first (last seq)))