(ns rlm.dna-melting

   (:refer-clojure :only [])

   (:require rlm.ns-rlm mobius.base))

 (rlm.ns-rlm/ns-clone mobius.base)

 (use :reload-all 'rlm.web-stuff)

 				

 

 (def pre "/home/r/MIT/6/6.122/DNA/samples/")

 

 

 

 (defrecord dna-data [temp glow name temp-units glow-units]

   clojure.lang.IFn

   (invoke [this k] (get this k))

   (invoke [this k not-found] (get this k not-found))

   (toString [this] "whatever"))

 

 (defmethod print-method dna-data [d w]

 	   (.write w (.toString d)))

 

 (defmethod print-dup dna-data [d w]

 	   (print-method d w))

 

 (defn load-data [source]

   (let [data (read-dataset source :delim \tab)

 	temp (vec ($ :col0 data))

 	glow (vec ($ :col1 data))]

     (dna-data. temp glow

 	        (last (re-split #"/" source))

 		"Raw Temperature (Volts)"

 		"Raw Glow (Volts)")))

 

 

 (def ?-1 (load-data (str pre "sample4-1.txt")))

 

 ;; I accidentally removed the cuvette while it was still recording!

 ;; drop last 5 seconds of invalid data!

 (def ?-2* (load-data (str pre "sample4-2.txt")))

 (def ?-2 (assoc ?-2*

 	   :temp (drop-last 500 (:temp ?-2*))

 	   :glow (drop-last 500 (:glow ?-2*))))

 

 (def ?-3 (load-data (str pre "sample4-3.txt")))

 (def A-1 (load-data (str pre "sampleA-1.txt")))

 (def A-2 (load-data (str pre "sampleA-2.txt")))

 (def A-3 (load-data (str pre "sampleA-3.txt")))

 (def B-2 (load-data (str pre "sampleB-2.txt")))

 (def B-3 (load-data (str pre "sampleB-3.txt")))

 (def C-1 (load-data (str pre "sampleC-1.txt")))

 (def C-2 (load-data (str pre "sampleC-2.txt")))

 (def C-3 (load-data (str pre "sampleC-3.txt")))

 

 

 (def pure-bleaching (load-data (str pre "sampleB-1(failed-to-heat).txt")))

 (def double-cycle (load-data (str pre "sampleC-3(double).txt")))

 (def triple-cycle (load-data (str pre "sampleC-3(triple).txt")))

 

 (def all-samples

      [?-1 ?-3 A-1 A-2 A-3 B-2 B-3 C-1 C-2 C-3])

 

 (def everything

      (concat [pure-bleaching double-cycle triple-cycle] all-samples))

 

 

 

 (defn dna-plot [d]

     (view

      (->

       (visual (:temp d))

       (set-y-label (:temp-units d))

       (set-x-label "Time (deci-seconds)")

       (set-title  (:name d))))

     (view

      (->

       (visual (:glow d))

       (set-y-label (:glow-units d))

       (set-x-label "Time (deci-seconds)")

       (set-title  (:name d))))

     d)

 

 (defn dna-melting [d]

   

    (->

     (scatter-plot (:temp d) (:glow d))

     (set-y-label (:glow-units d))

     (set-x-label (:temp-units d))

     (set-title (:name d))))

 

 

 (defmethod visual clojure.lang.LazySeq [ls]

 	   (visual (vec ls)))

 

 (defmethod visual clojure.lang.PersistentVector [v]

 	   (set-y-range (xy-plot (vec (range (count v)))  v)

 			      (apply min v) (apply max v)))

 

 (defmethod visual dna-data [d]

 	   (dna-plot d))

 

 

 

 

     

 

 ;;; Filter V_RTD and V_F to reduce noise

 

 

 

 ;(let [a (int 1.8e4) b (int 2.2e4)]

 ;  (visual (map (vec (temperature triple-cycle)) (range a b)))

 ;  (visual (map (vec (dna-fraction triple-cycle)) (range a b))))

 

 ;; Differentiate the resulting function with respect to temperature. 

 

 (def test-temp-range (vec (range 260 380 0.1)))

 

 (defn trans-print [& args] (apply println args) args)

 

 ;; (def search

 ;;      #(for [x (range -1000 1000 100) y (range -1000 1000 100)]

 ;; 	(do

 ;; 	  (print [x y])

 ;; 	  (try (trans-print (non-linear-model

 ;; 			     (ideal-dna-fraction 30e-6)

 ;; 			      test-dna-fraction test-temperature [x y]

 ;; 			     :method :newton-raphson))

 	       

 ;; 	       (catch Exception e (trans-print nil))))))

 

 

 ;;(non-linear-model (ideal-dna-fraction 30e-6) test-dna-fraction test-temperature [-184  -71e3])

 

 ;;(time (def answer (non-linear-model (ideal-dna-fraction 30e-6) test-dna-fraction test-temperature [-200  -90e3])))

 

 ;;we're assuming that the SYBR Green 1 doesn't glow at all when there is only

 ;;single stranded DNA, but it actually does. V_f only ever gets down to ~0.1 V

 ;;even with everything melted.

 

 

 ;;correct for thermal capacatance

 

 (defn find-max [v]

   (first (reduce (fn [a b] (if (> (last b) (last a)) b a)) (indexed v))))

 

 (defn find-min [v]

   (first (reduce (fn [a b] (if (< (last b) (last a)) b a)) (indexed v))))

 

 (defn plot-thermal-delay []

    (let [a (int 2.5e3) b (int 5e3)]

      (dna-plot (assoc A-1

 		 :temp (map (vec (:temp A-1)) (range a b))

 		 :glow (map (vec (:glow A-1)) (range a b))))))

      

 (defvar thermal-delay

      (let [a (int 2.5e3) b (int 5e3)]

        (- 

 	(find-min (map (vec (:glow A-1)) (range a b)))

 	(find-max (map (vec (:temp A-1)) (range a b)))))

      "There are two sources of thermal delay: thermal capacatance

       of the cuvette holder and thermal capacatance from the time

       it takes DNA to unwind.

 

       Use A-1 to find the thermal delay. There was almost no delay

       between heating and cooling for A-1, and thus the extremum of

       temperature and glow should temporally coincide if there were

       no delays. It turns out that there is about a 30 second lag!" )

 

 (defn thermal-delay-correction*

   "There is a certain delay between the RTD reporting changes in

    temperature and the temperature of the DNA actually changing.

    This shifts the temperature values so that the delay is eliminated."

   [thermal-delay #^dna-data d]

   (assoc d

     :name (str (:name d) " |thermal-correction| ")

     :temp (drop-last thermal-delay (:temp d))

     :glow (drop thermal-delay (:glow d))))

   

 (def thermal-delay-correction (partial thermal-delay-correction* thermal-delay))

 

 ;;correct for photo-bleaching

 

 (defn general-linear [[A B] x] (+ (* A x) B))

 

 (defvar bleaching-coef nil

    ;; (let [[A B]

    ;; 	 (:coefs				   

    ;; 	  (let [a (int 3e4) b (length (:glow pure-bleaching))]

    ;; 	    (non-linear-model

    ;; 	     general-linear

    ;; 	     (map (vec (:glow pure-bleaching)) (range a b))

    ;; 	     (range a b)

    ;; 	     [1 1])))]

    ;;   A)

   "This function corrects V_f for photo-bleaching over a period of time.

      Photo-bleaching was found to be linear(!) over quite a long period of time.

      Data is taken from sample B-1, which I left in the holder for an hour

      with no heating")

 

 (defn bleaching-correction*

   [bleaching-coef #^dna-data d]

   ;;the correction is clibrated to apply to the raw data only

   (assert (= (:glow-units d) "Raw Glow (Volts)"))

   ;;we need to modulate this correction factor

 

   (let [baseline (first (drop 3e4 (:temp pure-bleaching)))]

     (assoc d

       :name (str (:name d) " |bleaching-correction| ")

       :glow (map (fn [[i n]] (- n (*

 				   (/ n baseline)

 				   ((partial general-linear [bleaching-coef 0]) i))))

 		 (indexed (:glow d))))))

 

 (def bleaching-correction (partial bleaching-correction* bleaching-coef))

 

 ;; compensate for residual fluorescence when everything is melted

 

 (defn residual-fluorescence

   "the SYBR1 green dye doesn't completly stop glowing once everything

    is melted."

   [#^dna-data d]

   (let [low (apply min (:glow d))

 	high (apply max (:glow d))]

     (assoc d

       :name (str (:name d) " |residual-fluorescence| ")

       :glow (map (fn [f]

 		   (- f (* (- 1 (/ f high)) low))) (:glow d)))))

 

 

 

 

 ;; extract the cooling phase

 

 (defn cooling-phase

   "first stretch where the temperature is going down for a while

    equivalent to where the f curve starts to recover"

   [#^dna-data d]

   (let [t (:temp d)

 	f (:glow d)

 	dividing-line (reduce (fn [a b] (if (< (a 1) (b 1)) a b)) (indexed f))

 	t-new (drop (dividing-line 0) t)

 	f-new (drop (dividing-line 0) f)]

     (assoc d

       :name (str (:name d) " |extract cooling| ")

       :temp (vec t-new) :glow (vec f-new))))

 

 

 

 ;; filter noise

 ;;(http://clojure-log.n01se.net/date/2009-10-13.html by ambient

 ;;this is such a cool transform!!!

 

 

 (defn lo-pass

   [d coll]

   (reductions #(+ (* %2 d) (* %1 (- 1.0 d))) coll))

 

 

 ;;)

 

 (def reasonable-low-pass (partial lo-pass 1e-2))

 

 (defn low-pass-correction

   "filter the raw data through a low pass filter to remove

    electronic noise"

   [#^dna-data d]

   (assert (= (:glow-units d) "Raw Glow (Volts)"))

   (assert (= (:temp-units d) "Raw Temperature (Volts)"))

   (assoc d

     :name (str (:name d) " |low-pass| ")

     :temp (reasonable-low-pass (:temp d))

     :glow (reasonable-low-pass (:glow d))))

   

 

 

 

 

 ;; correct for units

 

 (defn RTD-voltage->temperature

   "transform RTD voltage to temperature.

   R_RTD = 1000 + 3.85 * (T - 273)

   T = ((R_RTD - 1000) / 3.85 ) + 273

   V_RTD = 15 * R_RTD / (R + R_RTD)

   R_RTD = R / (15/V_RTD - 1)"

   [coll]

   (let [R 14.97e3

 	V 15]

     (letfn [(T [R_RTD] (+ 273 (/ (- R_RTD 1000)  3.85))) ;; from lab notrd

 	    (R_RTD [V_RTD] (/ R (- (/ V V_RTD) 1)))]

       (map (comp T R_RTD) coll))))

 

 (defn normalize

   "Transform photodiode current to relative fluorescence."

   [coll]

   (let [low (apply min coll)

 	high (apply max coll)]

     (map #(/ (- % low) (- high low)) coll)))

 

 (defn unit-transform

   "turn raw voltage into temperature and dna-concentration"

   [#^dna-data d]

   (assert (= (:glow-units d) "Raw Glow (Volts)"))

   (assert (= (:temp-units d) "Raw Temperature (Volts)"))

   (assoc d

     :name (str (:name d) " |units| ")

     :temp-units "Temperature (Kelvin)"

     :glow-units "DNA-Fracton"

     :temp (RTD-voltage->temperature (:temp d))

     :glow (normalize (:glow d))))

 

 

 (defn temp-unit-transform

   "turn raw voltage into temperature and dna-concentration"

   [#^dna-data d]

   (assert (= (:temp-units d) "Raw Temperature (Volts)"))

   (assoc d

     :name (str (:name d) " |temp-units| ")

     :temp-units "Temperature (Kelvin)"

     :temp (RTD-voltage->temperature (:temp d))))

 

 (defn glow-unit-transform

   "turn raw voltage into temperature and dna-concentration"

   [#^dna-data d]

   (assert (= (:glow-units d) "Raw Glow (Volts)"))

   (assoc d

     :name (str (:name d) " |glow-units| ")

     :glow-units "DNA-Fracton"

     :glow (normalize (:glow d))))

 

 

 (defn thermal-efficiency

   [#^dna-data d]

   (assert (= (:temp-units d) "Temperature (Kelvin)"))

   (assert (= (:glow-units d) "Raw Glow (Volts)"))

   (assoc d

     :name (str (:name d) " |thermal-efficiency| ")

     :glow nil))

 

 

 

 (defn remove-first-values

   "remove the first minute or so of data so that the initial temperature of

    the cuvette itself doesn't matter"

   [#^dna-data d]

   (assoc d

     :name (str (:name d) " |remove-first| ")

     :temp (drop 600 (:temp d)) 

     :glow (drop 600 (:glow d))))

 

 (defn diff [coll]  (map - (rest coll) coll ))

 (defn temp-transform

   "my assumption will be that at the very end of the temperature data, the

    cuvette is exactly the same temperature as the RTD.  Then, I use newton's

    heat transfer and numerical integration to derive the heat of the cuvette

    for previous times"

   [k-value temperature-data]

   (let [T1 (vec (reverse temperature-data))

     ;; want fast indexed access, and instead of doing the problem in reverse

     ;; just reverse the vector at the start, treat it normally, and

     ;; then reverse at the end

 	initial-cuvette-temp (first T1)

 	initial-RTD-temp     (first T1)

 	step (fn [prev-T2 current-T1]

 	  (let [new-T2 (+ prev-T2 (* k-value (- current-T1 prev-T2)))]

 	    new-T2))]

     (reverse (reductions step initial-cuvette-temp T1))))

 

 (def k-value 0.1)

 

 (defn temp-corr-2 [k-value T1*]

   (let [T1 (reverse T1*)]

     (reverse  (cons (first T1)  (map + (map (partial * k-value) (diff T1)) (rest T1))))))

 

 (defn temp-corr-3 [initial-cuvette-temp k-value T1]

   (let [step

 	(fn [prev-T2 current-T1]

 	  (+ prev-T2 (* k-value (- current-T1 prev-T2))))]

     (rest (reductions step initial-cuvette-temp T1))))

 

 (defn newton-temp-correction* [k-value #^dna-data d]

   (assoc d

     :name (str (:name d) " |newton-temp| ")

     :temp (temp-transform k-value (:temp d))))

 (defn newton-temp-correction2* [k-value #^dna-data d]

   (assoc d

     :name (str (:name d) " |newton-temp-reverse| ")

     :temp (temp-corr-2 k-value (:temp d))))

 (defn newton-temp-correction3* [init-temp k-value #^dna-data d]

   (assoc d

     :name (str (:name d) " |newton-temp-reverse| ")

     :temp (temp-corr-3 init-temp k-value (:temp d))))

 

 (defn turn-e [n] (format "%5.2e" n))

 

 (defn newton-temp-correction4*

   "now with AIR!!!"

   [k-block-sample k-sample-air initial-sample-temp room-temp #^dna-data d]

   (let [old-temp (:temp d)

 	step (fn [prev-T2 current-T1]

 	       (+ prev-T2

 		  (* k-block-sample (- current-T1 prev-T2))

 		  (* k-sample-air   (- room-temp  prev-T2))))

 	new-temp (rest (reductions step initial-sample-temp old-temp))]

     (assoc d

       :name (str (:name d)

 		 " |temp-adjusted {" 

 		 "initial-sample:" (turn initial-sample-temp) " "

 		 "room-temp:" (turn room-temp) " "

 		 "k-block-sample:" (turn-e k-block-sample) " "

 		 "k-sample-air:" (turn-e k-sample-air) "}")

       :temp new-temp)))

 

 (defn plot-adjusted-temp

   "visualize the temperature correction"

   [adjuster #^dna-data d]

   (view

    (add-lines

     (visual (:temp d))

     (vec (range (count (:temp d))))

     (vec (:temp (adjuster d))))))

   

     

 

 

 

 ;; process data

 

 (defn process [#^dna-data d]

   (cooling-phase

    (unit-transform

     (residual-fluorescence

      (bleaching-correction

       (low-pass-correction

        ;(thermal-delay-correction

 	d))))))

 

 

 (defn first-pass

   "after this, the only thing stopping the melting curve from

    being a straight line is thermal capacatance"

   [#^dna-data data]

   (low-pass-correction

    ;(remove-first-values

     (bleaching-correction

      data)))

 

 

 

 (defn explore [k-value]

   (dna-melting (newton-temp-correction* k-value (first-pass A-2))))

 

 

 (defn explore2 [k-value]

   (dna-melting (newton-temp-correction2* k-value (first-pass A-2))))

   

 

 (defn explore3 [init-temp k-value]

   (dna-melting (newton-temp-correction3* init-temp  k-value (unit-transform (first-pass A-2)))))

   

 

 

 ;; get the three characteristic values we want.

 

 (defn ideal-dna-fraction* [concentration [delta-S delta-H] temperature]

     (let [R 8.31447215

 	  t temperature

 	  Ct concentration

 	  K_eq (exp (- (/ delta-S R) (/ delta-H (* R t))))

 	  f (/

 	     (+ 1 (* Ct K_eq)

 		(- (sqrt (+ 1 (* 2 Ct K_eq)))))

 	     (* Ct K_eq))]

   f))

 

 (def ideal-dna-fraction (partial ideal-dna-fraction* 30e-6))

 

 (defvar starting-values

   ;;  [-184 -71e3]

   [-904 -32e4]

   "chosen from theoritical estimation") 

 

 

 (defn analyze* [iterations #^dna-data d]

   (let [best-fit-results

 	(non-linear-model

 	 ideal-dna-fraction

 	 (:glow d)

 	 (:temp d)

 	 starting-values

 	:max-iter iterations)]

     best-fit-results))

 		    

 

 (def analyze (partial analyze* 200))

 

 (defn plot-results [d]

   (let [ans (analyze d)]

     (view (doto (xy-plot (:x ans) (:y ans))

       (set-title (:name d))

       (add-lines (:x ans) (:fitted ans))))

     ans))

 

 (defn zzz [d [A B]]

   (let [gar (vec (map (partial ideal-dna-fraction [A B]) (:temp d)))]

     (view (doto (xy-plot (:temp d) (:glow d))

       (set-title (:name d))

       (add-lines (:temp d) gar)))))

   

 

 

 

 (defn index-filter [pred coll]

   (when pred

     (for [[idx elt] (indexed coll) :when (pred elt)] idx)))

 

 

 

 

 ;; we're going to estimate two correction factors,

 ;; a constant capacative temperature delay,

 ;; and a linear decrease in brightness.  We'll do this

 ;; before unit-transforms.

 

 (defn zipmap-collect

   "Returns a map with the keys mapped to the corresponding vals.

    adds extra values to a list"

   [keys vals]

     (loop [map {}

            ks (seq keys)

            vs (seq vals)]

       (if (and ks vs)

         (recur (if (contains? map (first ks))

 		 (assoc map (first ks) (cons (first vs) (map (first ks))))

 		 (assoc map (first ks) (list (first vs))))

                (next ks)

                (next vs))

         map)))

 

 

 (defn little-process [thermal-cap bleaching #^dna-data d]

   ((partial thermal-delay-correction thermal-cap)

    ((partial bleaching-correction bleaching) d)))

 

 (def little-process (memoize little-process))

  

 ;(defn lookup-map [#^dna-data d]

   

 

 (defn converge-function

   "this is the vertical discrepancy after applying the two transforms"

   [#^dna-data d [thermal-cap bleaching] x]

   (let [new-d (little-process thermal-cap bleaching d)

 	spreads

 	(map (vec (:glow new-d))

 	     (index-filter (partial = x) (vec (:temp new-d))))

 	spread (- (apply max spreads) (apply min spreads))]

     spread))

 

 

 

 

 

 

 

 

 

 

 (defn  blah  [k] (comp 

 	   ( partial newton-temp-correction3* 280 k ) temp-unit-transform first-pass))

 

 ;(def ggg (comp 

 ;	  dna-melting ( partial thermal-delay-correction* 120 ) (blah 0.0018)))

 

 

 

 

 

 ;; Ensure that the melting function is uniquely valued by

 ;;   combining samples with identical temperature values.

 (defn cooling-function

   "just put it in a map for this!"

   [source]

   (let [data (cooling-phase source)

 	lookup (apply sorted-map

 		    (interleave

 		     (:temperature data) (:dna-fraction data)))]

     {:temperature (vec (keys lookup))

      :dna-fraction (vec (map lookup (keys lookup)))}))

 

 

 (defn turn [n]

   (format "%5.2f" n))

 

 (defn bleach-correction-2 [bleaching-constant #^dna-data d] 

   (assert (= (:glow-units d) "Raw Glow (Volts)"))

   (let [old-glow (vec (:glow d))

 	integral-old-glow (vec (reductions + old-glow))

 	correction (fn [t] (+

 			    (old-glow t)

 			    (* bleaching-constant (integral-old-glow t))))

 	new-glow (map + (map correction (range 0 (count old-glow))) old-glow)]

     ;;glow is related to the number of particles,

     ;;and the number destroyed is some fraction of

     ;;the number present and glowing.

     (assoc d

       :name (str (:name d) " |bleaching-corrected constant = " (turn bleaching-constant) "| ")

       :glow new-glow)))

     

     

 

 (defn ultimate-correction

   ([k-block-sample k-sample-air initial-sample-temp room-temp bleach-constant

     #^dna-data d

     ]

      ((comp

       (partial newton-temp-correction4*

 	       k-block-sample k-sample-air initial-sample-temp room-temp)

 

       temp-unit-transform

       (partial bleach-correction-2 bleach-constant)

       low-pass-correction) d))

 

   ([{:keys [test-k-sample-air

 	    test-k-block-sample

 	    test-sample-temp

 	    test-room-temp

 	    test-bleach]}

 		#^dna-data d]

   (ultimate-correction

     (mean test-k-block-sample) (mean test-k-sample-air)

     (mean test-sample-temp) (mean test-room-temp) (mean test-bleach) d)))

 

 

 ;;((uui 0.0020 0.0001 280 293) (bleach-correction-2 1.1e-5 ?-3))

 

 

 ;;; this is for ?-3 after 35 iterations

 (def opt {:test-k-sample-air '(9.375630358472117E-4 9.375630359053611E-4), :test-k-block-sample '(0.0010343754438043107 0.0010343754455447197), :test-room-temp '(294.9804782861611 294.98047828674316), :test-sample-temp '(293.36120605294127 293.3612060546875), :test-bleach '(1.771081481128931E-5 1.7710814863094127E-5), :align 1.5796052678842228E-4})

 

 

 ;;; C-1 after 20 iterations

 (def opt2

      {:test-k-sample-air '(9.936472587585448E-4 9.936491641998292E-4), :test-k-block-sample '(2.9390106201171873E-4 2.939580917358398E-4), :test-room-temp '(294.76531982421875 294.7653388977051), :test-sample-temp '(294.765625 294.765682220459), :test-bleach '(2.39045524597168E-5 2.3906250000000002E-5), :align 1.7562642331021098E-4})

 

 

 ;;(time (def result-?-3 (doall (take 10 (iterate (partial improve ?-3) start-values-overestimates)))))

 

 

 ;;(doall (map (comp dna-melting #(ultimate-correction % ?-3)) result-?-3))

 

 

 (defn subsample

   "subsample v every n points"

   [n v]

 

     (vec (map (vec v) (vec (range 0 (count (vec v)) n)))))

 

 

 (defn max-difference [v] (- (apply max v) (apply min v)))

    

 (defn alignment [#^dna-data d]

   (let [distribution (vec (vals

 			  (apply sorted-map

 				 (interleave (:temp d) (:glow d)))))]

     (mean (map variance (partition 10 distribution)))))

 

 (defn garigh [d]

   (let [distribution (vec (vals

 			   (apply sorted-map

 				  (interleave (:temp d) (:glow d)))))]

     (view (visual (vec  (map max-difference (partition 100 distribution)))))))

 

 

 (defn alignment-4 [#^dna-data d]

   (let [distribution (vec (vals

 			  (apply sorted-map

 				 (interleave (:temp d) (:glow d)))))]

     (reduce + (map sq (map max-difference (partition 100 distribution))))))

 

 

 (defn alignment-2 [#^dna-data d]

   (let [distribution

 	(vec

 	 (vals

 	  (apply sorted-map

 		 (interleave (subsample 2 (:temp d)) (subsample 2(:glow d))))))]

     (reduce +     (map max-difference (partition 30 distribution)))))

 (defn alignment-3 [#^dna-data d]

   (let [distribution

 	(vec

 	 (vals

 	  (apply sorted-map

 		 (interleave (subsample 2 (:temp d)) (subsample 2(:glow d))))))]

     (reduce max     (map max-difference (partition 30 distribution)))))

 

 

 (comment

   ;; this takes a while

   (def opts (zipmap all-samples (map #(last (take 20 (iterate (partial improve %) start-values))) all-samples))))

 

 ;(def opts (zipmap all-samples opts-vals))

 

 (defn mean-map [v s]

   (map mean (map (partial vector v) s)))

 

   

 (def start-values

      {:test-k-sample-air [0.000001 0.001]

       :test-k-block-sample   [0.001 0.03]

       :test-room-temp  [290 300]

       :test-sample-temp  [285 310]

       :test-bleach [0.00001 0.0009]})

 

 (def start-values-2 

      {:test-k-sample-air [1e-15 0.0005]

       :test-k-block-sample   [0.001 0.002]

       :test-room-temp  [290 300]

       :test-sample-temp  [290 300]

       :test-bleach [1e-7 0.0005]})

 

     

 

 (def start-values-3

      {:test-k-block-sample   [0.0005  0.00125 0.002]

       :test-k-sample-air [0.0000001 0.000001 0.0000019]

       :test-sample-temp  [280 300]

       :test-room-temp  [300 306]

       :test-bleach [0.000001 0.00001]})

          

 (defn condense [[a b :as v]]

   (if (< (/ (abs (- (apply max v) (apply min v))) (apply max v)) 1e-2) (vector (mean v)) v))

 

 (defn mass-shift

   [value v]

   (let [frac 2/3]

     (cond (>= (count v) 2)

 	  (let [[a b] v]

 	    (condense (vec (map #(float (+ (* frac %) (* (- 1 frac) value))) v))))

 	  (= (count v) 1)

 	  v)))

 

 

 (defn improve* [#^dna-data d

 	       {:keys [test-k-sample-air

 		       test-k-block-sample

 		       test-room-temp

 		       test-sample-temp

 		       test-bleach] :as domain}]

   (do (print \newline)

       (reduce (fn [m1 m2] (if (< (:align m1) (:align m2)) m1 m2))

 

 	      

 	      (for [k-sample-air (condense test-k-sample-air)

 		    k-block-sample (condense test-k-block-sample)

 		    room-temp (condense test-room-temp)

 		    initial-sample-temp (condense test-sample-temp)

 		    bleach-constant (condense test-bleach)]

 		

 		(do

 		  (print ".")

 		  (let [align

 			(alignment-2

 			 (ultimate-correction

 			  k-block-sample k-sample-air initial-sample-temp room-temp

 			  bleach-constant d))

 			align (if (Double/isNaN align) Double/POSITIVE_INFINITY  align)]

 		    {:test-k-sample-air (mass-shift k-sample-air test-k-sample-air)

 		     :test-k-block-sample(mass-shift k-block-sample test-k-block-sample)

 		     :test-room-temp (mass-shift room-temp test-room-temp)

 		     :test-sample-temp (mass-shift initial-sample-temp test-sample-temp)

 		     :test-bleach (mass-shift bleach-constant test-bleach)

 		     :align align}))))))

 

 

 

 (def println-repl (bound-fn [& args] (apply println args)))

 (def print-repl (bound-fn [& args] (apply print args)))

 

 

 

 (defn improve [#^dna-data d

 	       {:keys [test-k-block-sample

 		       test-k-sample-air

 		       test-sample-temp

 		       test-room-temp

 		       test-bleach] :as domain}]

        (reduce

        (fn [m1 m2] (if (< (:align m1) (:align m2)) m1 m2))

        (let [return-hash

 	     (fn

 	       [[k-block-sample k-sample-air initial-sample-temp

 		room-temp bleach-constant]]

 	       (do

 		 (print-repl ".")

 		 (let [align

 		       (alignment-4

 			(ultimate-correction

 			 k-block-sample k-sample-air initial-sample-temp room-temp

 			 bleach-constant d))

 		       align (if (Double/isNaN align) Double/POSITIVE_INFINITY  align)]

 		   {:test-k-sample-air (mass-shift k-sample-air test-k-sample-air)

 		    :test-k-block-sample(mass-shift k-block-sample test-k-block-sample)

 		    :test-room-temp (mass-shift room-temp test-room-temp)

 		    :test-sample-temp (mass-shift initial-sample-temp test-sample-temp)

 		    :test-bleach (mass-shift bleach-constant test-bleach)

 		    :align align})))]

 	 (pmap return-hash 

 	       (cartesian-product

 		 test-k-block-sample

 		 test-k-sample-air

 		 test-sample-temp

 		 test-room-temp

 		 test-bleach)))))

 

 (def improve (memoize improve))

 

 (defn iterate-test [n #^dna-data d]

   (last (take n (iterate (partial improve d) start-values-3))))

 

 

 (defn nested-vals [key coll]

   (for [x (tree-seq coll? seq coll) :when (contains? x key)]

        (get x key)))

 

 

 (defn quantum-efficiency* [q-constant #^dna-data d]

     ;; 0.8 at room temperature

     ;; quantum efficiency is somehow

     ;; inversely proportional to temperature

   (let [correction

 

 

 

 	(fn [t g] (/ g   (* 0.8 (exp (* q-constant (- 292 t))))))]

       (assoc d

 	:name (str (:name d) " |quantum-correction| ")

 	:glow (map correction (:temp d) (:glow d)))))

 

 (def quantum-efficiency (partial quantum-efficiency* 1.2e-2))

 

 (defn exp-decay [[A] x]

   (* 0.8 (exp (* A (- 298 x)))))

 

 ;(defn final-join [#^dna-data d]

   

   

 

 (def join-opts

      {

       A-1 {:test-k-sample-air [1.760184242508937E-6], :test-k-block-sample [0.0012930832259977858], :test-room-temp [304.41977945963544], :test-sample-temp [281.31689453125], :test-bleach [9.956068424799014E-6], :align 0.800815268639094}

       A-2 {:test-k-sample-air [3.017639376897326E-7], :test-k-block-sample [0.0011550507430608075], :test-room-temp [304.6666717529297], :test-sample-temp [282.6337890625], :test-bleach [9.5610075732111E-6], :align 0.8530570499185829}

       A-3 {:test-k-sample-air [2.6829059556281816E-6], :test-k-block-sample [0.002413294818252325], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6]}

 

       B-2 {:test-k-sample-air [1.6274943088016396E-6], :test-k-block-sample [0.0013714926317334175], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6], :align 0.7631166194402621}

       B-3

       {:test-k-sample-air [1.455345833771086E-6], :test-k-block-sample [0.0014841814991086721], :test-room-temp [304.6666717529297], :test-sample-temp [287.9835662841797], :test-bleach [9.965317076421343E-6], :align 0.5358800697386289}

 

       C-1

       {:test-k-sample-air [1.6907095717518434E-6], :test-k-block-sample [0.001278722076676786], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6], :align 0.7249279378017799}

       C-2

       {:test-k-sample-air [1.77524979487013E-6], :test-k-block-sample [0.001278722076676786], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6], :align 2.0474987460209793}

       C-3

       {:test-k-sample-air [1.010541116860016E-6], :test-k-block-sample [0.001601272417853276], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6], :align 0.7191738524207253}

 

       ?-1

       {:test-k-sample-air [1.5533593114014366E-6], :test-k-block-sample [0.0013761443551629782], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6], :align 1.3838506291264723}

       ?-3

       {:test-k-sample-air [2.940294715851148E-7], :test-k-block-sample [0.0018434864080821474], :test-room-temp [304.6666717529297], :test-sample-temp [281.31689453125], :test-bleach [9.965317076421343E-6], :align 0.48549623932134006}

 

       })

 

 

 (defn final-join [#^dna-data d]

   

   (let [

 	m (apply sorted-map (interleave (:temp d) (:glow d)))

 	distribution (vec (vals m))

 	ordinates (vec (keys m))

 			   ]

     (assoc d

       :name (str (:name d) " |joined| ")

       :temp (map (fn [region] (mean [(apply min region) (apply max region)]))

 		 (partition 30 ordinates))

 

       :glow (map (fn [region] (mean [(apply min region) (apply max region)]))

 		 (partition 30 distribution)))))

 

 

 (defn melting-point [#^dna-data d]

   (let [t (drop 150 (:temp d))

 	g (drop 150 (:glow d))]

     (nth t (find-min (map / (diff g) (diff t) )))))

 

 

 (def target-dir (file-str "/home/r/MIT/6/6.122/DNA/output"))

 

 (defn x-form-dna [#^dna-data d]

   (final-join (ultimate-correction (join-opts d) d)))

 

 (defn save-raw-dna-curve [#^dna-data d]

   (save (dna-melting d) (str target-dir File/separator (:name d) ".png")))

 

 (defn save-dna-curve-with-temp [#^dna-data d*]

   (let [

 

 	melt (melting-point d*)

 	d (assoc d* :name (str (:name d*) " |M.P. = " (turn melt) " |"))

 	]

     (save (add-lines (dna-melting d) (repeat (count (:glow d)) melt) (:glow d))

 	  (str target-dir File/separator (:name d) ".png"))))

 

   

   

 (defn process-all []

   (doall (pmap  save-dna-curve all-samples))

   (doall (pmap (comp save-dna-curve-with-temp x-form-dna) all-samples)))

 

 

 (def ?-temps (vec (map (comp melting-point x-form-dna) [?-1 ?-3])))

 (def A-temps (vec (map (comp melting-point x-form-dna) [A-1 A-2 A-3])))

 (def B-temps (vec (map (comp melting-point x-form-dna) [B-2 B-3])))

 (def C-temps (vec (map (comp melting-point x-form-dna) [C-1 C-2 C-3])))

 

 (defn report [temps]

   

   (sorted-map :data temps :std (sqrt (variance temps)) :mean (mean temps)))

 

 

 ;;; now it's time to start data reporting

 

 ;(defn process [#^dna-data d]

 ;  (ultimate-correction (opts d) d)

 

 

 

 (report A-temps)

 {:data [348.8771481779248 349.3989434763391 348.76079920968573],

  :mean 349.0122969546499,

  :std 0.33986161910548834}

 

 (report B-temps)

 {:data [349.0655539116902 349.92155530431506],

  :mean 349.4935546080026,

  :std 0.6052843894485812}

 

 (report C-temps)

 {:data [349.92341967511976 349.68832687337056 348.2872189461977],

  :mean 349.299655164896,

  :std 0.884639745344559}

 

 (report ?-temps)

 {:data [350.4525305192239 349.87377093318975],

  :mean 350.1631507262068,

  :std 0.4092448279511269}