(use-modules ((ice-9 format) :select (format))
	     ((coasim batch) :select (fold)))

(define no-leaves 20)
(define no-runs 100000)

(define (calc-scale-factor beta)
  (let ((no-growth-tree-length 
	 (let loop ((j 1) (sum 0))
	   (if (< j no-leaves) 
	       (loop (+ j 1) (+ sum (/ 1 j)))
	       (* 2 sum))))
	(growht-tree-length 
	 (let* ((no-iterations 10000)
		(branch-sum
		 (fold no-iterations (lambda (val sum) (+ val sum)) 0
		       (let* ((m (snp-marker 0.5 0 1))
			      (arg (simulate (list m) no-leaves :beta beta))
			      (tree (car (local-trees arg)))
			      (branch-length (total-branch-length tree)))
			 branch-length))))
	   (/ branch-sum no-iterations))))
    (/ no-growth-tree-length growht-tree-length)))

(define rho-scale-factor
  ;; map from beta to scale factor
  ;;          beta scale-factor
  (list (cons   0  1)
	(cons  10  (calc-scale-factor 10))
	(cons 100  (calc-scale-factor 100))))


(define (inc-counter counts val)
  (let ((cell (assq val counts)))
    (if cell
	(begin (assq-set! counts val (+ 1 (cdr cell)))
	       counts)
	(acons val 1 counts))))

(define (count-recombinations rho beta)
  (let loop ((n 0) (counts '()))
    (if (= n no-runs) counts
	(let* ((arg (simulate '() no-leaves 
			      :rho rho :beta beta
			      :keep-empty-intervals #t))
	       (recombs (no-recombinations arg)))
	  (loop (+ n 1) (inc-counter counts recombs))))))
       
(define (print-recomb-counts beta counts)
  (for-each 
   (lambda (c) (format #t "~d\t~d\t~d\n" beta (car c) (cdr c))) 
   counts))

(define (run beta)
  (let ((rho (* 0.3 (assoc-ref rho-scale-factor beta))))
    (print-recomb-counts beta (count-recombinations rho beta))))

(display "beta\tno\tcount\n")
(for-each run '(0 10 100))