Model Repository/Kingman's Coalescent

From Church Wiki
Jump to: navigation, search


(load/cd "../../../core/church-eval.ss")

;; coalescent
(church-top-level-eval
 '(

   (define (make-coalescent)
     (let ((U (mem (lambda (sym) (beta 1 1))))
           (V (mem (lambda (k) (beta 1 1)))))
       (lambda (p1 p2)
         (let ((U1 (U p1))
               (U2 (U p2)))
           (let loop ((i 1))
             (if (between? (V i) U1 U2)
                 i
                 (loop (+ i 1)) ))))))

   (define (between? v u1 u2)
     (or (and (> u1 v)
              (> v u2) )
         (and (> u2 v)
              (> v u1) )))
   
   (define (make-pairs a lst)
     (if (null? lst)
         '()
         (pair (list a (first lst))
               (make-pairs a (rest lst)) )))
   
   (define (all-pairs lst)
     (if (null? lst)
         '()
         (append (make-pairs (first lst) (rest lst))
                 (all-pairs (rest lst)) )))

   (define c (make-coalescent))
   
   ; tree => (0 tree1 tree2),tree => leaf
   (define leaf? symbol?)
   (define tree-root first)
   (define tree< second)
   (define tree> third)
   (define (any-leaf tree)
     (if (leaf? tree)
         tree 
         (any-leaf (tree< tree)) ))
   (define (leaves tree)
     (if (leaf? tree)
         (list tree)
         (append (leaves (tree< tree))
                 (leaves (tree> tree)) )))
   (define (extend-tree el tree)
     (if (leaf? tree)
         (list 0 tree el)
         (let* ((leaf< (any-leaf (tree< tree)))
                (leaf> (any-leaf (tree> tree)))
                (max-c< (c leaf< el))
                (max-c> (c leaf> el)) )
           (cond ((= max-c< max-c>)
                  (list 0 el tree))
                 ((< max-c< max-c>)
                  (list 0
                        (tree< tree)
                        (extend-tree el (tree> tree)) ))
                 (else
                  (list 0 
                        (extend-tree el (tree< tree))
                        (tree> tree) )))
           )))
                 
   (rest (list (repeat 14 (lambda () (beta 1 1)))
               ;; all pairs of coalescent merge points
               'all-pairs
               (map (lambda (a-pair) (pair a-pair
                                           (c (first a-pair)
                                              (second a-pair) )))
                    (all-pairs '(a b e d c)) )
               ;; autogenerated tree
               'auto-gen-tree
               (fold extend-tree 'a '(b c d e))

                 ))

   
   ))
Personal tools