(use-modules (srfi srfi-1))

;; beta->rho map
(define rhos '((0 . 0.3) (10 . 1) (100 . 4)))

(define no-haplotypes 20)
(define markers
  (map (lambda (pos) (snp-marker pos 0.1 0.9))
       '(0.000 0.0000001 0.0000002 0.0000003 0.0000004 
	       0.0000005 0.0000006 0.0000007 0.0000008 0.0000009 
	       0.9999990 0.9999991 0.9999992 0.9999993 0.9999994
	       0.9999995 0.9999996 0.9999997 0.9999998 0.9999999)))

(define (bind2nd fun s) (lambda (f) (fun f s)))
(define (print-sequences port sequences)
  (let ((print-sequence (lambda (seq)
			  (map (lambda (x) 
				 (display x port)(display " " port)) 
			       seq)
			  (newline port))))
    (map print-sequence sequences)))
			
(define (print-recombination-times port times-and-ks)
  (let ((times (car  times-and-ks))
	(ks    (cadr times-and-ks))
	(print (lambda (t k)
		 (display t port)
		 (display " " port)
		 (display k port)
		 (newline port))))
    (map print times ks)))

(define (print-trees port trees)
  (map (bind2nd display port) trees))

(define (run-simulation beta)
  (let* ((rho (assoc-ref rhos beta))
	 (rc-times '())
	 (ks '())
	 (rc-cb (lambda (n1 n2 k)
		  (set! ks (cons k ks))
		  (set! rc-times (cons (event-time n1) rc-times))))
	 (arg (simulate markers no-haplotypes
			:rho rho :beta beta :keep-empty-intervals #t
			:recombination-callback rc-cb))
	 (seq (sequences arg))
	 (trees (local-trees arg))
	 
	 (seq-file  "seq.txt")
	 (time-file "recomb-times.txt")
	 (tree-file "trees.txt")
	 
	 (seq-writer  (bind2nd print-sequences seq))
	 (time-writer (bind2nd print-recombination-times (list rc-times ks)))
	 (tree-writer (bind2nd print-trees trees)))

    (map (lambda (l) (display l)(newline)) (map total-branch-length trees))

    (call-with-output-file "tree-lengths.txt"
      (lambda (port)
	(for-each (lambda (th) (display th port)(newline port))
		  (map total-branch-length trees))))

    (call-with-output-file seq-file  seq-writer)
    (call-with-output-file time-file time-writer)
    (call-with-output-file tree-file tree-writer)))

;; read the beta parameter from the command line options
(let* ((options (command-line))
       (beta (string->number (cadr options)))) ; option 0 is coasim_guile
					       ; option 1 is beta
  (run-simulation beta))