;; Load functionality from modules
(use-modules ((coasim markers)        :select (insert-sorted-idx))
	     ((coasim io)             :select (sequences-printer
					       marker-positions-printer))
	     ((coasim SNP haplotypes) :select (split-in-cases-controls))
	     ((srfi srfi-1)           :select (take)))


;; For pop.size ~10,000, rho=400 is a mutation rate of 0.01, or about
;; a centi-Morgan.  Since we are simulating a region of twice the
;; size, to test adding markers on both sides of the "real" interval,
;; we simulate with rho=800
(define rho 800) 

(define (random-positions range-start range-stop n)
  (let* ((msec (cdr (gettimeofday)))
	 (random-state (seed->random-state msec))
	 (range-size (- range-stop range-start))
	 (mk-pos (lambda ()
		   (+ range-start (random range-size random-state))))
	 (pos (let loop ((i 0) (acc '()))
		(if (= i n) acc
		    (loop (+ i 1) (cons (mk-pos) acc))))))
    (sort pos <)))


(define snp-markers
  (let* ((first-positions (random-positions 0.0 0.25 10))
	 (middle-positions (random-positions 0.25 0.75 40))
	 (last-positions (random-positions 0.75 1 10))
	 (positions (append first-positions middle-positions last-positions))
	 (mk-marker (lambda (pos) (snp-marker pos 0.1 0.9))))
    (map mk-marker positions)))

(define disease-marker
  (let ((pos (car (random-positions 0.25 0.75 1))))
    (trait-marker pos 0.18 0.22)))
	

(define markers-and-index (insert-sorted-idx snp-markers disease-marker))
(define markers   (car  markers-and-index))
(define trait-idx (cadr markers-and-index))

;; simulating 600 haplotypes -- with at least 18% mutants, this should
;; give us >100 cases (600*0.18 == 108)
(define haplotypes (simulate-sequences markers 600 :rho rho))

;; Split in cases and controls, based on trait marker
(define cases-and-controls (split-in-cases-controls haplotypes trait-idx))
(define cases    (take (car  cases-and-controls) 100))
(define controls (take (cadr cases-and-controls) 100))


(call-with-output-file "snp-positions.txt" 
  (marker-positions-printer snp-markers))
(call-with-output-file "trait-position.txt" 
  (marker-positions-printer (list disease-marker)))
(call-with-output-file "cases.txt"    (sequences-printer cases))
(call-with-output-file "controls.txt" (sequences-printer controls))