rlm@0: (ns rlm.dna-melting
rlm@0:   (:refer-clojure :only [])
rlm@0:   (:require rlm.ns-rlm mobius.base))
rlm@0: (rlm.ns-rlm/ns-clone mobius.base)
rlm@0: (use :reload-all 'rlm.web-stuff)
rlm@0: 				
rlm@0: 
rlm@0: (def pre "/home/r/MIT/6/6.122/DNA/samples/")
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: (defrecord dna-data [temp glow name temp-units glow-units]
rlm@0:   clojure.lang.IFn
rlm@0:   (invoke [this k] (get this k))
rlm@0:   (invoke [this k not-found] (get this k not-found))
rlm@0:   (toString [this] "whatever"))
rlm@0: 
rlm@0: (defmethod print-method dna-data [d w]
rlm@0: 	   (.write w (.toString d)))
rlm@0: 
rlm@0: (defmethod print-dup dna-data [d w]
rlm@0: 	   (print-method d w))
rlm@0: 
rlm@0: (defn load-data [source]
rlm@0:   (let [data (read-dataset source :delim \tab)
rlm@0: 	temp (vec ($ :col0 data))
rlm@0: 	glow (vec ($ :col1 data))]
rlm@0:     (dna-data. temp glow
rlm@0: 	        (last (re-split #"/" source))
rlm@0: 		"Raw Temperature (Volts)"
rlm@0: 		"Raw Glow (Volts)")))
rlm@0: 
rlm@0: 
rlm@0: (def ?-1 (load-data (str pre "sample4-1.txt")))
rlm@0: 
rlm@0: ;; I accidentally removed the cuvette while it was still recording!
rlm@0: ;; drop last 5 seconds of invalid data!
rlm@0: (def ?-2* (load-data (str pre "sample4-2.txt")))
rlm@0: (def ?-2 (assoc ?-2*
rlm@0: 	   :temp (drop-last 500 (:temp ?-2*))
rlm@0: 	   :glow (drop-last 500 (:glow ?-2*))))
rlm@0: 
rlm@0: (def ?-3 (load-data (str pre "sample4-3.txt")))
rlm@0: (def A-1 (load-data (str pre "sampleA-1.txt")))
rlm@0: (def A-2 (load-data (str pre "sampleA-2.txt")))
rlm@0: (def A-3 (load-data (str pre "sampleA-3.txt")))
rlm@0: (def B-2 (load-data (str pre "sampleB-2.txt")))
rlm@0: (def B-3 (load-data (str pre "sampleB-3.txt")))
rlm@0: (def C-1 (load-data (str pre "sampleC-1.txt")))
rlm@0: (def C-2 (load-data (str pre "sampleC-2.txt")))
rlm@0: (def C-3 (load-data (str pre "sampleC-3.txt")))
rlm@0: 
rlm@0: 
rlm@0: (def pure-bleaching (load-data (str pre "sampleB-1(failed-to-heat).txt")))
rlm@0: (def double-cycle (load-data (str pre "sampleC-3(double).txt")))
rlm@0: (def triple-cycle (load-data (str pre "sampleC-3(triple).txt")))
rlm@0: 
rlm@0: (def all-samples
rlm@0:      [?-1 ?-3 A-1 A-2 A-3 B-2 B-3 C-1 C-2 C-3])
rlm@0: 
rlm@0: (def everything
rlm@0:      (concat [pure-bleaching double-cycle triple-cycle] all-samples))
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: (defn dna-plot [d]
rlm@0:     (view
rlm@0:      (->
rlm@0:       (visual (:temp d))
rlm@0:       (set-y-label (:temp-units d))
rlm@0:       (set-x-label "Time (deci-seconds)")
rlm@0:       (set-title  (:name d))))
rlm@0:     (view
rlm@0:      (->
rlm@0:       (visual (:glow d))
rlm@0:       (set-y-label (:glow-units d))
rlm@0:       (set-x-label "Time (deci-seconds)")
rlm@0:       (set-title  (:name d))))
rlm@0:     d)
rlm@0: 
rlm@0: (defn dna-melting [d]
rlm@0:   
rlm@0:    (->
rlm@0:     (scatter-plot (:temp d) (:glow d))
rlm@0:     (set-y-label (:glow-units d))
rlm@0:     (set-x-label (:temp-units d))
rlm@0:     (set-title (:name d))))
rlm@0: 
rlm@0: 
rlm@0: (defmethod visual clojure.lang.LazySeq [ls]
rlm@0: 	   (visual (vec ls)))
rlm@0: 
rlm@0: (defmethod visual clojure.lang.PersistentVector [v]
rlm@0: 	   (set-y-range (xy-plot (vec (range (count v)))  v)
rlm@0: 			      (apply min v) (apply max v)))
rlm@0: 
rlm@0: (defmethod visual dna-data [d]
rlm@0: 	   (dna-plot d))
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0:     
rlm@0: 
rlm@0: ;;; Filter V_RTD and V_F to reduce noise
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: ;(let [a (int 1.8e4) b (int 2.2e4)]
rlm@0: ;  (visual (map (vec (temperature triple-cycle)) (range a b)))
rlm@0: ;  (visual (map (vec (dna-fraction triple-cycle)) (range a b))))
rlm@0: 
rlm@0: ;; Differentiate the resulting function with respect to temperature. 
rlm@0: 
rlm@0: (def test-temp-range (vec (range 260 380 0.1)))
rlm@0: 
rlm@0: (defn trans-print [& args] (apply println args) args)
rlm@0: 
rlm@0: ;; (def search
rlm@0: ;;      #(for [x (range -1000 1000 100) y (range -1000 1000 100)]
rlm@0: ;; 	(do
rlm@0: ;; 	  (print [x y])
rlm@0: ;; 	  (try (trans-print (non-linear-model
rlm@0: ;; 			     (ideal-dna-fraction 30e-6)
rlm@0: ;; 			      test-dna-fraction test-temperature [x y]
rlm@0: ;; 			     :method :newton-raphson))
rlm@0: 	       
rlm@0: ;; 	       (catch Exception e (trans-print nil))))))
rlm@0: 
rlm@0: 
rlm@0: ;;(non-linear-model (ideal-dna-fraction 30e-6) test-dna-fraction test-temperature [-184  -71e3])
rlm@0: 
rlm@0: ;;(time (def answer (non-linear-model (ideal-dna-fraction 30e-6) test-dna-fraction test-temperature [-200  -90e3])))
rlm@0: 
rlm@0: ;;we're assuming that the SYBR Green 1 doesn't glow at all when there is only
rlm@0: ;;single stranded DNA, but it actually does. V_f only ever gets down to ~0.1 V
rlm@0: ;;even with everything melted.
rlm@0: 
rlm@0: 
rlm@0: ;;correct for thermal capacatance
rlm@0: 
rlm@0: (defn find-max [v]
rlm@0:   (first (reduce (fn [a b] (if (> (last b) (last a)) b a)) (indexed v))))
rlm@0: 
rlm@0: (defn find-min [v]
rlm@0:   (first (reduce (fn [a b] (if (< (last b) (last a)) b a)) (indexed v))))
rlm@0: 
rlm@0: (defn plot-thermal-delay []
rlm@0:    (let [a (int 2.5e3) b (int 5e3)]
rlm@0:      (dna-plot (assoc A-1
rlm@0: 		 :temp (map (vec (:temp A-1)) (range a b))
rlm@0: 		 :glow (map (vec (:glow A-1)) (range a b))))))
rlm@0:      
rlm@0: (defvar thermal-delay
rlm@0:      (let [a (int 2.5e3) b (int 5e3)]
rlm@0:        (- 
rlm@0: 	(find-min (map (vec (:glow A-1)) (range a b)))
rlm@0: 	(find-max (map (vec (:temp A-1)) (range a b)))))
rlm@0:      "There are two sources of thermal delay: thermal capacatance
rlm@0:       of the cuvette holder and thermal capacatance from the time
rlm@0:       it takes DNA to unwind.
rlm@0: 
rlm@0:       Use A-1 to find the thermal delay. There was almost no delay
rlm@0:       between heating and cooling for A-1, and thus the extremum of
rlm@0:       temperature and glow should temporally coincide if there were
rlm@0:       no delays. It turns out that there is about a 30 second lag!" )
rlm@0: 
rlm@0: (defn thermal-delay-correction*
rlm@0:   "There is a certain delay between the RTD reporting changes in
rlm@0:    temperature and the temperature of the DNA actually changing.
rlm@0:    This shifts the temperature values so that the delay is eliminated."
rlm@0:   [thermal-delay #^dna-data d]
rlm@0:   (assoc d
rlm@0:     :name (str (:name d) " |thermal-correction| ")
rlm@0:     :temp (drop-last thermal-delay (:temp d))
rlm@0:     :glow (drop thermal-delay (:glow d))))
rlm@0:   
rlm@0: (def thermal-delay-correction (partial thermal-delay-correction* thermal-delay))
rlm@0: 
rlm@0: ;;correct for photo-bleaching
rlm@0: 
rlm@0: (defn general-linear [[A B] x] (+ (* A x) B))
rlm@0: 
rlm@0: (defvar bleaching-coef nil
rlm@0:    ;; (let [[A B]
rlm@0:    ;; 	 (:coefs				   
rlm@0:    ;; 	  (let [a (int 3e4) b (length (:glow pure-bleaching))]
rlm@0:    ;; 	    (non-linear-model
rlm@0:    ;; 	     general-linear
rlm@0:    ;; 	     (map (vec (:glow pure-bleaching)) (range a b))
rlm@0:    ;; 	     (range a b)
rlm@0:    ;; 	     [1 1])))]
rlm@0:    ;;   A)
rlm@0:   "This function corrects V_f for photo-bleaching over a period of time.
rlm@0:      Photo-bleaching was found to be linear(!) over quite a long period of time.
rlm@0:      Data is taken from sample B-1, which I left in the holder for an hour
rlm@0:      with no heating")
rlm@0: 
rlm@0: (defn bleaching-correction*
rlm@0:   [bleaching-coef #^dna-data d]
rlm@0:   ;;the correction is clibrated to apply to the raw data only
rlm@0:   (assert (= (:glow-units d) "Raw Glow (Volts)"))
rlm@0:   ;;we need to modulate this correction factor
rlm@0: 
rlm@0:   (let [baseline (first (drop 3e4 (:temp pure-bleaching)))]
rlm@0:     (assoc d
rlm@0:       :name (str (:name d) " |bleaching-correction| ")
rlm@0:       :glow (map (fn [[i n]] (- n (*
rlm@0: 				   (/ n baseline)
rlm@0: 				   ((partial general-linear [bleaching-coef 0]) i))))
rlm@0: 		 (indexed (:glow d))))))
rlm@0: 
rlm@0: (def bleaching-correction (partial bleaching-correction* bleaching-coef))
rlm@0: 
rlm@0: ;; compensate for residual fluorescence when everything is melted
rlm@0: 
rlm@0: (defn residual-fluorescence
rlm@0:   "the SYBR1 green dye doesn't completly stop glowing once everything
rlm@0:    is melted."
rlm@0:   [#^dna-data d]
rlm@0:   (let [low (apply min (:glow d))
rlm@0: 	high (apply max (:glow d))]
rlm@0:     (assoc d
rlm@0:       :name (str (:name d) " |residual-fluorescence| ")
rlm@0:       :glow (map (fn [f]
rlm@0: 		   (- f (* (- 1 (/ f high)) low))) (:glow d)))))
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: ;; extract the cooling phase
rlm@0: 
rlm@0: (defn cooling-phase
rlm@0:   "first stretch where the temperature is going down for a while
rlm@0:    equivalent to where the f curve starts to recover"
rlm@0:   [#^dna-data d]
rlm@0:   (let [t (:temp d)
rlm@0: 	f (:glow d)
rlm@0: 	dividing-line (reduce (fn [a b] (if (< (a 1) (b 1)) a b)) (indexed f))
rlm@0: 	t-new (drop (dividing-line 0) t)
rlm@0: 	f-new (drop (dividing-line 0) f)]
rlm@0:     (assoc d
rlm@0:       :name (str (:name d) " |extract cooling| ")
rlm@0:       :temp (vec t-new) :glow (vec f-new))))
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: ;; filter noise
rlm@0: ;;(http://clojure-log.n01se.net/date/2009-10-13.html by ambient
rlm@0: ;;this is such a cool transform!!!
rlm@0: 
rlm@0: 
rlm@0: (defn lo-pass
rlm@0:   [d coll]
rlm@0:   (reductions #(+ (* %2 d) (* %1 (- 1.0 d))) coll))
rlm@0: 
rlm@0: 
rlm@0: ;;)
rlm@0: 
rlm@0: (def reasonable-low-pass (partial lo-pass 1e-2))
rlm@0: 
rlm@0: (defn low-pass-correction
rlm@0:   "filter the raw data through a low pass filter to remove
rlm@0:    electronic noise"
rlm@0:   [#^dna-data d]
rlm@0:   (assert (= (:glow-units d) "Raw Glow (Volts)"))
rlm@0:   (assert (= (:temp-units d) "Raw Temperature (Volts)"))
rlm@0:   (assoc d
rlm@0:     :name (str (:name d) " |low-pass| ")
rlm@0:     :temp (reasonable-low-pass (:temp d))
rlm@0:     :glow (reasonable-low-pass (:glow d))))
rlm@0:   
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: ;; correct for units
rlm@0: 
rlm@0: (defn RTD-voltage->temperature
rlm@0:   "transform RTD voltage to temperature.
rlm@0:   R_RTD = 1000 + 3.85 * (T - 273)
rlm@0:   T = ((R_RTD - 1000) / 3.85 ) + 273
rlm@0:   V_RTD = 15 * R_RTD / (R + R_RTD)
rlm@0:   R_RTD = R / (15/V_RTD - 1)"
rlm@0:   [coll]
rlm@0:   (let [R 14.97e3
rlm@0: 	V 15]
rlm@0:     (letfn [(T [R_RTD] (+ 273 (/ (- R_RTD 1000)  3.85))) ;; from lab notrd
rlm@0: 	    (R_RTD [V_RTD] (/ R (- (/ V V_RTD) 1)))]
rlm@0:       (map (comp T R_RTD) coll))))
rlm@0: 
rlm@0: (defn normalize
rlm@0:   "Transform photodiode current to relative fluorescence."
rlm@0:   [coll]
rlm@0:   (let [low (apply min coll)
rlm@0: 	high (apply max coll)]
rlm@0:     (map #(/ (- % low) (- high low)) coll)))
rlm@0: 
rlm@0: (defn unit-transform
rlm@0:   "turn raw voltage into temperature and dna-concentration"
rlm@0:   [#^dna-data d]
rlm@0:   (assert (= (:glow-units d) "Raw Glow (Volts)"))
rlm@0:   (assert (= (:temp-units d) "Raw Temperature (Volts)"))
rlm@0:   (assoc d
rlm@0:     :name (str (:name d) " |units| ")
rlm@0:     :temp-units "Temperature (Kelvin)"
rlm@0:     :glow-units "DNA-Fracton"
rlm@0:     :temp (RTD-voltage->temperature (:temp d))
rlm@0:     :glow (normalize (:glow d))))
rlm@0: 
rlm@0: 
rlm@0: (defn temp-unit-transform
rlm@0:   "turn raw voltage into temperature and dna-concentration"
rlm@0:   [#^dna-data d]
rlm@0:   (assert (= (:temp-units d) "Raw Temperature (Volts)"))
rlm@0:   (assoc d
rlm@0:     :name (str (:name d) " |temp-units| ")
rlm@0:     :temp-units "Temperature (Kelvin)"
rlm@0:     :temp (RTD-voltage->temperature (:temp d))))
rlm@0: 
rlm@0: (defn glow-unit-transform
rlm@0:   "turn raw voltage into temperature and dna-concentration"
rlm@0:   [#^dna-data d]
rlm@0:   (assert (= (:glow-units d) "Raw Glow (Volts)"))
rlm@0:   (assoc d
rlm@0:     :name (str (:name d) " |glow-units| ")
rlm@0:     :glow-units "DNA-Fracton"
rlm@0:     :glow (normalize (:glow d))))
rlm@0: 
rlm@0: 
rlm@0: (defn thermal-efficiency
rlm@0:   [#^dna-data d]
rlm@0:   (assert (= (:temp-units d) "Temperature (Kelvin)"))
rlm@0:   (assert (= (:glow-units d) "Raw Glow (Volts)"))
rlm@0:   (assoc d
rlm@0:     :name (str (:name d) " |thermal-efficiency| ")
rlm@0:     :glow nil))
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: (defn remove-first-values
rlm@0:   "remove the first minute or so of data so that the initial temperature of
rlm@0:    the cuvette itself doesn't matter"
rlm@0:   [#^dna-data d]
rlm@0:   (assoc d
rlm@0:     :name (str (:name d) " |remove-first| ")
rlm@0:     :temp (drop 600 (:temp d)) 
rlm@0:     :glow (drop 600 (:glow d))))
rlm@0: 
rlm@0: (defn diff [coll]  (map - (rest coll) coll ))
rlm@0: (defn temp-transform
rlm@0:   "my assumption will be that at the very end of the temperature data, the
rlm@0:    cuvette is exactly the same temperature as the RTD.  Then, I use newton's
rlm@0:    heat transfer and numerical integration to derive the heat of the cuvette
rlm@0:    for previous times"
rlm@0:   [k-value temperature-data]
rlm@0:   (let [T1 (vec (reverse temperature-data))
rlm@0:     ;; want fast indexed access, and instead of doing the problem in reverse
rlm@0:     ;; just reverse the vector at the start, treat it normally, and
rlm@0:     ;; then reverse at the end
rlm@0: 	initial-cuvette-temp (first T1)
rlm@0: 	initial-RTD-temp     (first T1)
rlm@0: 	step (fn [prev-T2 current-T1]
rlm@0: 	  (let [new-T2 (+ prev-T2 (* k-value (- current-T1 prev-T2)))]
rlm@0: 	    new-T2))]
rlm@0:     (reverse (reductions step initial-cuvette-temp T1))))
rlm@0: 
rlm@0: (def k-value 0.1)
rlm@0: 
rlm@0: (defn temp-corr-2 [k-value T1*]
rlm@0:   (let [T1 (reverse T1*)]
rlm@0:     (reverse  (cons (first T1)  (map + (map (partial * k-value) (diff T1)) (rest T1))))))
rlm@0: 
rlm@0: (defn temp-corr-3 [initial-cuvette-temp k-value T1]
rlm@0:   (let [step
rlm@0: 	(fn [prev-T2 current-T1]
rlm@0: 	  (+ prev-T2 (* k-value (- current-T1 prev-T2))))]
rlm@0:     (rest (reductions step initial-cuvette-temp T1))))
rlm@0: 
rlm@0: (defn newton-temp-correction* [k-value #^dna-data d]
rlm@0:   (assoc d
rlm@0:     :name (str (:name d) " |newton-temp| ")
rlm@0:     :temp (temp-transform k-value (:temp d))))
rlm@0: (defn newton-temp-correction2* [k-value #^dna-data d]
rlm@0:   (assoc d
rlm@0:     :name (str (:name d) " |newton-temp-reverse| ")
rlm@0:     :temp (temp-corr-2 k-value (:temp d))))
rlm@0: (defn newton-temp-correction3* [init-temp k-value #^dna-data d]
rlm@0:   (assoc d
rlm@0:     :name (str (:name d) " |newton-temp-reverse| ")
rlm@0:     :temp (temp-corr-3 init-temp k-value (:temp d))))
rlm@0: 
rlm@0: (defn turn-e [n] (format "%5.2e" n))
rlm@0: 
rlm@0: (defn newton-temp-correction4*
rlm@0:   "now with AIR!!!"
rlm@0:   [k-block-sample k-sample-air initial-sample-temp room-temp #^dna-data d]
rlm@0:   (let [old-temp (:temp d)
rlm@0: 	step (fn [prev-T2 current-T1]
rlm@0: 	       (+ prev-T2
rlm@0: 		  (* k-block-sample (- current-T1 prev-T2))
rlm@0: 		  (* k-sample-air   (- room-temp  prev-T2))))
rlm@0: 	new-temp (rest (reductions step initial-sample-temp old-temp))]
rlm@0:     (assoc d
rlm@0:       :name (str (:name d)
rlm@0: 		 " |temp-adjusted {" 
rlm@0: 		 "initial-sample:" (turn initial-sample-temp) " "
rlm@0: 		 "room-temp:" (turn room-temp) " "
rlm@0: 		 "k-block-sample:" (turn-e k-block-sample) " "
rlm@0: 		 "k-sample-air:" (turn-e k-sample-air) "}")
rlm@0:       :temp new-temp)))
rlm@0: 
rlm@0: (defn plot-adjusted-temp
rlm@0:   "visualize the temperature correction"
rlm@0:   [adjuster #^dna-data d]
rlm@0:   (view
rlm@0:    (add-lines
rlm@0:     (visual (:temp d))
rlm@0:     (vec (range (count (:temp d))))
rlm@0:     (vec (:temp (adjuster d))))))
rlm@0:   
rlm@0:     
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: ;; process data
rlm@0: 
rlm@0: (defn process [#^dna-data d]
rlm@0:   (cooling-phase
rlm@0:    (unit-transform
rlm@0:     (residual-fluorescence
rlm@0:      (bleaching-correction
rlm@0:       (low-pass-correction
rlm@0:        ;(thermal-delay-correction
rlm@0: 	d))))))
rlm@0: 
rlm@0: 
rlm@0: (defn first-pass
rlm@0:   "after this, the only thing stopping the melting curve from
rlm@0:    being a straight line is thermal capacatance"
rlm@0:   [#^dna-data data]
rlm@0:   (low-pass-correction
rlm@0:    ;(remove-first-values
rlm@0:     (bleaching-correction
rlm@0:      data)))
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: (defn explore [k-value]
rlm@0:   (dna-melting (newton-temp-correction* k-value (first-pass A-2))))
rlm@0: 
rlm@0: 
rlm@0: (defn explore2 [k-value]
rlm@0:   (dna-melting (newton-temp-correction2* k-value (first-pass A-2))))
rlm@0:   
rlm@0: 
rlm@0: (defn explore3 [init-temp k-value]
rlm@0:   (dna-melting (newton-temp-correction3* init-temp  k-value (unit-transform (first-pass A-2)))))
rlm@0:   
rlm@0: 
rlm@0: 
rlm@0: ;; get the three characteristic values we want.
rlm@0: 
rlm@0: (defn ideal-dna-fraction* [concentration [delta-S delta-H] temperature]
rlm@0:     (let [R 8.31447215
rlm@0: 	  t temperature
rlm@0: 	  Ct concentration
rlm@0: 	  K_eq (exp (- (/ delta-S R) (/ delta-H (* R t))))
rlm@0: 	  f (/
rlm@0: 	     (+ 1 (* Ct K_eq)
rlm@0: 		(- (sqrt (+ 1 (* 2 Ct K_eq)))))
rlm@0: 	     (* Ct K_eq))]
rlm@0:   f))
rlm@0: 
rlm@0: (def ideal-dna-fraction (partial ideal-dna-fraction* 30e-6))
rlm@0: 
rlm@0: (defvar starting-values
rlm@0:   ;;  [-184 -71e3]
rlm@0:   [-904 -32e4]
rlm@0:   "chosen from theoritical estimation") 
rlm@0: 
rlm@0: 
rlm@0: (defn analyze* [iterations #^dna-data d]
rlm@0:   (let [best-fit-results
rlm@0: 	(non-linear-model
rlm@0: 	 ideal-dna-fraction
rlm@0: 	 (:glow d)
rlm@0: 	 (:temp d)
rlm@0: 	 starting-values
rlm@0: 	:max-iter iterations)]
rlm@0:     best-fit-results))
rlm@0: 		    
rlm@0: 
rlm@0: (def analyze (partial analyze* 200))
rlm@0: 
rlm@0: (defn plot-results [d]
rlm@0:   (let [ans (analyze d)]
rlm@0:     (view (doto (xy-plot (:x ans) (:y ans))
rlm@0:       (set-title (:name d))
rlm@0:       (add-lines (:x ans) (:fitted ans))))
rlm@0:     ans))
rlm@0: 
rlm@0: (defn zzz [d [A B]]
rlm@0:   (let [gar (vec (map (partial ideal-dna-fraction [A B]) (:temp d)))]
rlm@0:     (view (doto (xy-plot (:temp d) (:glow d))
rlm@0:       (set-title (:name d))
rlm@0:       (add-lines (:temp d) gar)))))
rlm@0:   
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: (defn index-filter [pred coll]
rlm@0:   (when pred
rlm@0:     (for [[idx elt] (indexed coll) :when (pred elt)] idx)))
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: ;; we're going to estimate two correction factors,
rlm@0: ;; a constant capacative temperature delay,
rlm@0: ;; and a linear decrease in brightness.  We'll do this
rlm@0: ;; before unit-transforms.
rlm@0: 
rlm@0: (defn zipmap-collect
rlm@0:   "Returns a map with the keys mapped to the corresponding vals.
rlm@0:    adds extra values to a list"
rlm@0:   [keys vals]
rlm@0:     (loop [map {}
rlm@0:            ks (seq keys)
rlm@0:            vs (seq vals)]
rlm@0:       (if (and ks vs)
rlm@0:         (recur (if (contains? map (first ks))
rlm@0: 		 (assoc map (first ks) (cons (first vs) (map (first ks))))
rlm@0: 		 (assoc map (first ks) (list (first vs))))
rlm@0:                (next ks)
rlm@0:                (next vs))
rlm@0:         map)))
rlm@0: 
rlm@0: 
rlm@0: (defn little-process [thermal-cap bleaching #^dna-data d]
rlm@0:   ((partial thermal-delay-correction thermal-cap)
rlm@0:    ((partial bleaching-correction bleaching) d)))
rlm@0: 
rlm@0: (def little-process (memoize little-process))
rlm@0:  
rlm@0: ;(defn lookup-map [#^dna-data d]
rlm@0:   
rlm@0: 
rlm@0: (defn converge-function
rlm@0:   "this is the vertical discrepancy after applying the two transforms"
rlm@0:   [#^dna-data d [thermal-cap bleaching] x]
rlm@0:   (let [new-d (little-process thermal-cap bleaching d)
rlm@0: 	spreads
rlm@0: 	(map (vec (:glow new-d))
rlm@0: 	     (index-filter (partial = x) (vec (:temp new-d))))
rlm@0: 	spread (- (apply max spreads) (apply min spreads))]
rlm@0:     spread))
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: (defn  blah  [k] (comp 
rlm@0: 	   ( partial newton-temp-correction3* 280 k ) temp-unit-transform first-pass))
rlm@0: 
rlm@0: ;(def ggg (comp 
rlm@0: ;	  dna-melting ( partial thermal-delay-correction* 120 ) (blah 0.0018)))
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: ;; Ensure that the melting function is uniquely valued by
rlm@0: ;;   combining samples with identical temperature values.
rlm@0: (defn cooling-function
rlm@0:   "just put it in a map for this!"
rlm@0:   [source]
rlm@0:   (let [data (cooling-phase source)
rlm@0: 	lookup (apply sorted-map
rlm@0: 		    (interleave
rlm@0: 		     (:temperature data) (:dna-fraction data)))]
rlm@0:     {:temperature (vec (keys lookup))
rlm@0:      :dna-fraction (vec (map lookup (keys lookup)))}))
rlm@0: 
rlm@0: 
rlm@0: (defn turn [n]
rlm@0:   (format "%5.2f" n))
rlm@0: 
rlm@0: (defn bleach-correction-2 [bleaching-constant #^dna-data d] 
rlm@0:   (assert (= (:glow-units d) "Raw Glow (Volts)"))
rlm@0:   (let [old-glow (vec (:glow d))
rlm@0: 	integral-old-glow (vec (reductions + old-glow))
rlm@0: 	correction (fn [t] (+
rlm@0: 			    (old-glow t)
rlm@0: 			    (* bleaching-constant (integral-old-glow t))))
rlm@0: 	new-glow (map + (map correction (range 0 (count old-glow))) old-glow)]
rlm@0:     ;;glow is related to the number of particles,
rlm@0:     ;;and the number destroyed is some fraction of
rlm@0:     ;;the number present and glowing.
rlm@0:     (assoc d
rlm@0:       :name (str (:name d) " |bleaching-corrected constant = " (turn bleaching-constant) "| ")
rlm@0:       :glow new-glow)))
rlm@0:     
rlm@0:     
rlm@0: 
rlm@0: (defn ultimate-correction
rlm@0:   ([k-block-sample k-sample-air initial-sample-temp room-temp bleach-constant
rlm@0:     #^dna-data d
rlm@0:     ]
rlm@0:      ((comp
rlm@0:       (partial newton-temp-correction4*
rlm@0: 	       k-block-sample k-sample-air initial-sample-temp room-temp)
rlm@0: 
rlm@0:       temp-unit-transform
rlm@0:       (partial bleach-correction-2 bleach-constant)
rlm@0:       low-pass-correction) d))
rlm@0: 
rlm@0:   ([{:keys [test-k-sample-air
rlm@0: 	    test-k-block-sample
rlm@0: 	    test-sample-temp
rlm@0: 	    test-room-temp
rlm@0: 	    test-bleach]}
rlm@0: 		#^dna-data d]
rlm@0:   (ultimate-correction
rlm@0:     (mean test-k-block-sample) (mean test-k-sample-air)
rlm@0:     (mean test-sample-temp) (mean test-room-temp) (mean test-bleach) d)))
rlm@0: 
rlm@0: 
rlm@0: ;;((uui 0.0020 0.0001 280 293) (bleach-correction-2 1.1e-5 ?-3))
rlm@0: 
rlm@0: 
rlm@0: ;;; this is for ?-3 after 35 iterations
rlm@0: (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})
rlm@0: 
rlm@0: 
rlm@0: ;;; C-1 after 20 iterations
rlm@0: (def opt2
rlm@0:      {: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})
rlm@0: 
rlm@0: 
rlm@0: ;;(time (def result-?-3 (doall (take 10 (iterate (partial improve ?-3) start-values-overestimates)))))
rlm@0: 
rlm@0: 
rlm@0: ;;(doall (map (comp dna-melting #(ultimate-correction % ?-3)) result-?-3))
rlm@0: 
rlm@0: 
rlm@0: (defn subsample
rlm@0:   "subsample v every n points"
rlm@0:   [n v]
rlm@0: 
rlm@0:     (vec (map (vec v) (vec (range 0 (count (vec v)) n)))))
rlm@0: 
rlm@0: 
rlm@0: (defn max-difference [v] (- (apply max v) (apply min v)))
rlm@0:    
rlm@0: (defn alignment [#^dna-data d]
rlm@0:   (let [distribution (vec (vals
rlm@0: 			  (apply sorted-map
rlm@0: 				 (interleave (:temp d) (:glow d)))))]
rlm@0:     (mean (map variance (partition 10 distribution)))))
rlm@0: 
rlm@0: (defn garigh [d]
rlm@0:   (let [distribution (vec (vals
rlm@0: 			   (apply sorted-map
rlm@0: 				  (interleave (:temp d) (:glow d)))))]
rlm@0:     (view (visual (vec  (map max-difference (partition 100 distribution)))))))
rlm@0: 
rlm@0: 
rlm@0: (defn alignment-4 [#^dna-data d]
rlm@0:   (let [distribution (vec (vals
rlm@0: 			  (apply sorted-map
rlm@0: 				 (interleave (:temp d) (:glow d)))))]
rlm@0:     (reduce + (map sq (map max-difference (partition 100 distribution))))))
rlm@0: 
rlm@0: 
rlm@0: (defn alignment-2 [#^dna-data d]
rlm@0:   (let [distribution
rlm@0: 	(vec
rlm@0: 	 (vals
rlm@0: 	  (apply sorted-map
rlm@0: 		 (interleave (subsample 2 (:temp d)) (subsample 2(:glow d))))))]
rlm@0:     (reduce +     (map max-difference (partition 30 distribution)))))
rlm@0: (defn alignment-3 [#^dna-data d]
rlm@0:   (let [distribution
rlm@0: 	(vec
rlm@0: 	 (vals
rlm@0: 	  (apply sorted-map
rlm@0: 		 (interleave (subsample 2 (:temp d)) (subsample 2(:glow d))))))]
rlm@0:     (reduce max     (map max-difference (partition 30 distribution)))))
rlm@0: 
rlm@0: 
rlm@0: (comment
rlm@0:   ;; this takes a while
rlm@0:   (def opts (zipmap all-samples (map #(last (take 20 (iterate (partial improve %) start-values))) all-samples))))
rlm@0: 
rlm@0: ;(def opts (zipmap all-samples opts-vals))
rlm@0: 
rlm@0: (defn mean-map [v s]
rlm@0:   (map mean (map (partial vector v) s)))
rlm@0: 
rlm@0:   
rlm@0: (def start-values
rlm@0:      {:test-k-sample-air [0.000001 0.001]
rlm@0:       :test-k-block-sample   [0.001 0.03]
rlm@0:       :test-room-temp  [290 300]
rlm@0:       :test-sample-temp  [285 310]
rlm@0:       :test-bleach [0.00001 0.0009]})
rlm@0: 
rlm@0: (def start-values-2 
rlm@0:      {:test-k-sample-air [1e-15 0.0005]
rlm@0:       :test-k-block-sample   [0.001 0.002]
rlm@0:       :test-room-temp  [290 300]
rlm@0:       :test-sample-temp  [290 300]
rlm@0:       :test-bleach [1e-7 0.0005]})
rlm@0: 
rlm@0:     
rlm@0: 
rlm@0: (def start-values-3
rlm@0:      {:test-k-block-sample   [0.0005  0.00125 0.002]
rlm@0:       :test-k-sample-air [0.0000001 0.000001 0.0000019]
rlm@0:       :test-sample-temp  [280 300]
rlm@0:       :test-room-temp  [300 306]
rlm@0:       :test-bleach [0.000001 0.00001]})
rlm@0:          
rlm@0: (defn condense [[a b :as v]]
rlm@0:   (if (< (/ (abs (- (apply max v) (apply min v))) (apply max v)) 1e-2) (vector (mean v)) v))
rlm@0: 
rlm@0: (defn mass-shift
rlm@0:   [value v]
rlm@0:   (let [frac 2/3]
rlm@0:     (cond (>= (count v) 2)
rlm@0: 	  (let [[a b] v]
rlm@0: 	    (condense (vec (map #(float (+ (* frac %) (* (- 1 frac) value))) v))))
rlm@0: 	  (= (count v) 1)
rlm@0: 	  v)))
rlm@0: 
rlm@0: 
rlm@0: (defn improve* [#^dna-data d
rlm@0: 	       {:keys [test-k-sample-air
rlm@0: 		       test-k-block-sample
rlm@0: 		       test-room-temp
rlm@0: 		       test-sample-temp
rlm@0: 		       test-bleach] :as domain}]
rlm@0:   (do (print \newline)
rlm@0:       (reduce (fn [m1 m2] (if (< (:align m1) (:align m2)) m1 m2))
rlm@0: 
rlm@0: 	      
rlm@0: 	      (for [k-sample-air (condense test-k-sample-air)
rlm@0: 		    k-block-sample (condense test-k-block-sample)
rlm@0: 		    room-temp (condense test-room-temp)
rlm@0: 		    initial-sample-temp (condense test-sample-temp)
rlm@0: 		    bleach-constant (condense test-bleach)]
rlm@0: 		
rlm@0: 		(do
rlm@0: 		  (print ".")
rlm@0: 		  (let [align
rlm@0: 			(alignment-2
rlm@0: 			 (ultimate-correction
rlm@0: 			  k-block-sample k-sample-air initial-sample-temp room-temp
rlm@0: 			  bleach-constant d))
rlm@0: 			align (if (Double/isNaN align) Double/POSITIVE_INFINITY  align)]
rlm@0: 		    {:test-k-sample-air (mass-shift k-sample-air test-k-sample-air)
rlm@0: 		     :test-k-block-sample(mass-shift k-block-sample test-k-block-sample)
rlm@0: 		     :test-room-temp (mass-shift room-temp test-room-temp)
rlm@0: 		     :test-sample-temp (mass-shift initial-sample-temp test-sample-temp)
rlm@0: 		     :test-bleach (mass-shift bleach-constant test-bleach)
rlm@0: 		     :align align}))))))
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: (def println-repl (bound-fn [& args] (apply println args)))
rlm@0: (def print-repl (bound-fn [& args] (apply print args)))
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: (defn improve [#^dna-data d
rlm@0: 	       {:keys [test-k-block-sample
rlm@0: 		       test-k-sample-air
rlm@0: 		       test-sample-temp
rlm@0: 		       test-room-temp
rlm@0: 		       test-bleach] :as domain}]
rlm@0:        (reduce
rlm@0:        (fn [m1 m2] (if (< (:align m1) (:align m2)) m1 m2))
rlm@0:        (let [return-hash
rlm@0: 	     (fn
rlm@0: 	       [[k-block-sample k-sample-air initial-sample-temp
rlm@0: 		room-temp bleach-constant]]
rlm@0: 	       (do
rlm@0: 		 (print-repl ".")
rlm@0: 		 (let [align
rlm@0: 		       (alignment-4
rlm@0: 			(ultimate-correction
rlm@0: 			 k-block-sample k-sample-air initial-sample-temp room-temp
rlm@0: 			 bleach-constant d))
rlm@0: 		       align (if (Double/isNaN align) Double/POSITIVE_INFINITY  align)]
rlm@0: 		   {:test-k-sample-air (mass-shift k-sample-air test-k-sample-air)
rlm@0: 		    :test-k-block-sample(mass-shift k-block-sample test-k-block-sample)
rlm@0: 		    :test-room-temp (mass-shift room-temp test-room-temp)
rlm@0: 		    :test-sample-temp (mass-shift initial-sample-temp test-sample-temp)
rlm@0: 		    :test-bleach (mass-shift bleach-constant test-bleach)
rlm@0: 		    :align align})))]
rlm@0: 	 (pmap return-hash 
rlm@0: 	       (cartesian-product
rlm@0: 		 test-k-block-sample
rlm@0: 		 test-k-sample-air
rlm@0: 		 test-sample-temp
rlm@0: 		 test-room-temp
rlm@0: 		 test-bleach)))))
rlm@0: 
rlm@0: (def improve (memoize improve))
rlm@0: 
rlm@0: (defn iterate-test [n #^dna-data d]
rlm@0:   (last (take n (iterate (partial improve d) start-values-3))))
rlm@0: 
rlm@0: 
rlm@0: (defn nested-vals [key coll]
rlm@0:   (for [x (tree-seq coll? seq coll) :when (contains? x key)]
rlm@0:        (get x key)))
rlm@0: 
rlm@0: 
rlm@0: (defn quantum-efficiency* [q-constant #^dna-data d]
rlm@0:     ;; 0.8 at room temperature
rlm@0:     ;; quantum efficiency is somehow
rlm@0:     ;; inversely proportional to temperature
rlm@0:   (let [correction
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: 	(fn [t g] (/ g   (* 0.8 (exp (* q-constant (- 292 t))))))]
rlm@0:       (assoc d
rlm@0: 	:name (str (:name d) " |quantum-correction| ")
rlm@0: 	:glow (map correction (:temp d) (:glow d)))))
rlm@0: 
rlm@0: (def quantum-efficiency (partial quantum-efficiency* 1.2e-2))
rlm@0: 
rlm@0: (defn exp-decay [[A] x]
rlm@0:   (* 0.8 (exp (* A (- 298 x)))))
rlm@0: 
rlm@0: ;(defn final-join [#^dna-data d]
rlm@0:   
rlm@0:   
rlm@0: 
rlm@0: (def join-opts
rlm@0:      {
rlm@0:       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}
rlm@0:       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}
rlm@0:       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]}
rlm@0: 
rlm@0:       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}
rlm@0:       B-3
rlm@0:       {: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}
rlm@0: 
rlm@0:       C-1
rlm@0:       {: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}
rlm@0:       C-2
rlm@0:       {: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}
rlm@0:       C-3
rlm@0:       {: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}
rlm@0: 
rlm@0:       ?-1
rlm@0:       {: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}
rlm@0:       ?-3
rlm@0:       {: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}
rlm@0: 
rlm@0:       })
rlm@0: 
rlm@0: 
rlm@0: (defn final-join [#^dna-data d]
rlm@0:   
rlm@0:   (let [
rlm@0: 	m (apply sorted-map (interleave (:temp d) (:glow d)))
rlm@0: 	distribution (vec (vals m))
rlm@0: 	ordinates (vec (keys m))
rlm@0: 			   ]
rlm@0:     (assoc d
rlm@0:       :name (str (:name d) " |joined| ")
rlm@0:       :temp (map (fn [region] (mean [(apply min region) (apply max region)]))
rlm@0: 		 (partition 30 ordinates))
rlm@0: 
rlm@0:       :glow (map (fn [region] (mean [(apply min region) (apply max region)]))
rlm@0: 		 (partition 30 distribution)))))
rlm@0: 
rlm@0: 
rlm@0: (defn melting-point [#^dna-data d]
rlm@0:   (let [t (drop 150 (:temp d))
rlm@0: 	g (drop 150 (:glow d))]
rlm@0:     (nth t (find-min (map / (diff g) (diff t) )))))
rlm@0: 
rlm@0: 
rlm@0: (def target-dir (file-str "/home/r/MIT/6/6.122/DNA/output"))
rlm@0: 
rlm@0: (defn x-form-dna [#^dna-data d]
rlm@0:   (final-join (ultimate-correction (join-opts d) d)))
rlm@0: 
rlm@0: (defn save-raw-dna-curve [#^dna-data d]
rlm@0:   (save (dna-melting d) (str target-dir File/separator (:name d) ".png")))
rlm@0: 
rlm@0: (defn save-dna-curve-with-temp [#^dna-data d*]
rlm@0:   (let [
rlm@0: 
rlm@0: 	melt (melting-point d*)
rlm@0: 	d (assoc d* :name (str (:name d*) " |M.P. = " (turn melt) " |"))
rlm@0: 	]
rlm@0:     (save (add-lines (dna-melting d) (repeat (count (:glow d)) melt) (:glow d))
rlm@0: 	  (str target-dir File/separator (:name d) ".png"))))
rlm@0: 
rlm@0:   
rlm@0:   
rlm@0: (defn process-all []
rlm@0:   (doall (pmap  save-dna-curve all-samples))
rlm@0:   (doall (pmap (comp save-dna-curve-with-temp x-form-dna) all-samples)))
rlm@0: 
rlm@0: 
rlm@0: (def ?-temps (vec (map (comp melting-point x-form-dna) [?-1 ?-3])))
rlm@0: (def A-temps (vec (map (comp melting-point x-form-dna) [A-1 A-2 A-3])))
rlm@0: (def B-temps (vec (map (comp melting-point x-form-dna) [B-2 B-3])))
rlm@0: (def C-temps (vec (map (comp melting-point x-form-dna) [C-1 C-2 C-3])))
rlm@0: 
rlm@0: (defn report [temps]
rlm@0:   
rlm@0:   (sorted-map :data temps :std (sqrt (variance temps)) :mean (mean temps)))
rlm@0: 
rlm@0: 
rlm@0: ;;; now it's time to start data reporting
rlm@0: 
rlm@0: ;(defn process [#^dna-data d]
rlm@0: ;  (ultimate-correction (opts d) d)
rlm@0: 
rlm@0: 
rlm@0: 
rlm@0: (report A-temps)
rlm@0: {:data [348.8771481779248 349.3989434763391 348.76079920968573],
rlm@0:  :mean 349.0122969546499,
rlm@0:  :std 0.33986161910548834}
rlm@0: 
rlm@0: (report B-temps)
rlm@0: {:data [349.0655539116902 349.92155530431506],
rlm@0:  :mean 349.4935546080026,
rlm@0:  :std 0.6052843894485812}
rlm@0: 
rlm@0: (report C-temps)
rlm@0: {:data [349.92341967511976 349.68832687337056 348.2872189461977],
rlm@0:  :mean 349.299655164896,
rlm@0:  :std 0.884639745344559}
rlm@0: 
rlm@0: (report ?-temps)
rlm@0: {:data [350.4525305192239 349.87377093318975],
rlm@0:  :mean 350.1631507262068,
rlm@0:  :std 0.4092448279511269}
rlm@0: 
rlm@0: