rlm@0
|
1 (ns rlm.dna-melting
|
rlm@0
|
2 (:refer-clojure :only [])
|
rlm@0
|
3 (:require rlm.ns-rlm mobius.base))
|
rlm@0
|
4 (rlm.ns-rlm/ns-clone mobius.base)
|
rlm@0
|
5 (use :reload-all 'rlm.web-stuff)
|
rlm@0
|
6
|
rlm@0
|
7
|
rlm@0
|
8 (def pre "/home/r/MIT/6/6.122/DNA/samples/")
|
rlm@0
|
9
|
rlm@0
|
10
|
rlm@0
|
11
|
rlm@0
|
12 (defrecord dna-data [temp glow name temp-units glow-units]
|
rlm@0
|
13 clojure.lang.IFn
|
rlm@0
|
14 (invoke [this k] (get this k))
|
rlm@0
|
15 (invoke [this k not-found] (get this k not-found))
|
rlm@0
|
16 (toString [this] "whatever"))
|
rlm@0
|
17
|
rlm@0
|
18 (defmethod print-method dna-data [d w]
|
rlm@0
|
19 (.write w (.toString d)))
|
rlm@0
|
20
|
rlm@0
|
21 (defmethod print-dup dna-data [d w]
|
rlm@0
|
22 (print-method d w))
|
rlm@0
|
23
|
rlm@0
|
24 (defn load-data [source]
|
rlm@0
|
25 (let [data (read-dataset source :delim \tab)
|
rlm@0
|
26 temp (vec ($ :col0 data))
|
rlm@0
|
27 glow (vec ($ :col1 data))]
|
rlm@0
|
28 (dna-data. temp glow
|
rlm@0
|
29 (last (re-split #"/" source))
|
rlm@0
|
30 "Raw Temperature (Volts)"
|
rlm@0
|
31 "Raw Glow (Volts)")))
|
rlm@0
|
32
|
rlm@0
|
33
|
rlm@0
|
34 (def ?-1 (load-data (str pre "sample4-1.txt")))
|
rlm@0
|
35
|
rlm@0
|
36 ;; I accidentally removed the cuvette while it was still recording!
|
rlm@0
|
37 ;; drop last 5 seconds of invalid data!
|
rlm@0
|
38 (def ?-2* (load-data (str pre "sample4-2.txt")))
|
rlm@0
|
39 (def ?-2 (assoc ?-2*
|
rlm@0
|
40 :temp (drop-last 500 (:temp ?-2*))
|
rlm@0
|
41 :glow (drop-last 500 (:glow ?-2*))))
|
rlm@0
|
42
|
rlm@0
|
43 (def ?-3 (load-data (str pre "sample4-3.txt")))
|
rlm@0
|
44 (def A-1 (load-data (str pre "sampleA-1.txt")))
|
rlm@0
|
45 (def A-2 (load-data (str pre "sampleA-2.txt")))
|
rlm@0
|
46 (def A-3 (load-data (str pre "sampleA-3.txt")))
|
rlm@0
|
47 (def B-2 (load-data (str pre "sampleB-2.txt")))
|
rlm@0
|
48 (def B-3 (load-data (str pre "sampleB-3.txt")))
|
rlm@0
|
49 (def C-1 (load-data (str pre "sampleC-1.txt")))
|
rlm@0
|
50 (def C-2 (load-data (str pre "sampleC-2.txt")))
|
rlm@0
|
51 (def C-3 (load-data (str pre "sampleC-3.txt")))
|
rlm@0
|
52
|
rlm@0
|
53
|
rlm@0
|
54 (def pure-bleaching (load-data (str pre "sampleB-1(failed-to-heat).txt")))
|
rlm@0
|
55 (def double-cycle (load-data (str pre "sampleC-3(double).txt")))
|
rlm@0
|
56 (def triple-cycle (load-data (str pre "sampleC-3(triple).txt")))
|
rlm@0
|
57
|
rlm@0
|
58 (def all-samples
|
rlm@0
|
59 [?-1 ?-3 A-1 A-2 A-3 B-2 B-3 C-1 C-2 C-3])
|
rlm@0
|
60
|
rlm@0
|
61 (def everything
|
rlm@0
|
62 (concat [pure-bleaching double-cycle triple-cycle] all-samples))
|
rlm@0
|
63
|
rlm@0
|
64
|
rlm@0
|
65
|
rlm@0
|
66 (defn dna-plot [d]
|
rlm@0
|
67 (view
|
rlm@0
|
68 (->
|
rlm@0
|
69 (visual (:temp d))
|
rlm@0
|
70 (set-y-label (:temp-units d))
|
rlm@0
|
71 (set-x-label "Time (deci-seconds)")
|
rlm@0
|
72 (set-title (:name d))))
|
rlm@0
|
73 (view
|
rlm@0
|
74 (->
|
rlm@0
|
75 (visual (:glow d))
|
rlm@0
|
76 (set-y-label (:glow-units d))
|
rlm@0
|
77 (set-x-label "Time (deci-seconds)")
|
rlm@0
|
78 (set-title (:name d))))
|
rlm@0
|
79 d)
|
rlm@0
|
80
|
rlm@0
|
81 (defn dna-melting [d]
|
rlm@0
|
82
|
rlm@0
|
83 (->
|
rlm@0
|
84 (scatter-plot (:temp d) (:glow d))
|
rlm@0
|
85 (set-y-label (:glow-units d))
|
rlm@0
|
86 (set-x-label (:temp-units d))
|
rlm@0
|
87 (set-title (:name d))))
|
rlm@0
|
88
|
rlm@0
|
89
|
rlm@0
|
90 (defmethod visual clojure.lang.LazySeq [ls]
|
rlm@0
|
91 (visual (vec ls)))
|
rlm@0
|
92
|
rlm@0
|
93 (defmethod visual clojure.lang.PersistentVector [v]
|
rlm@0
|
94 (set-y-range (xy-plot (vec (range (count v))) v)
|
rlm@0
|
95 (apply min v) (apply max v)))
|
rlm@0
|
96
|
rlm@0
|
97 (defmethod visual dna-data [d]
|
rlm@0
|
98 (dna-plot d))
|
rlm@0
|
99
|
rlm@0
|
100
|
rlm@0
|
101
|
rlm@0
|
102
|
rlm@0
|
103
|
rlm@0
|
104
|
rlm@0
|
105 ;;; Filter V_RTD and V_F to reduce noise
|
rlm@0
|
106
|
rlm@0
|
107
|
rlm@0
|
108
|
rlm@0
|
109 ;(let [a (int 1.8e4) b (int 2.2e4)]
|
rlm@0
|
110 ; (visual (map (vec (temperature triple-cycle)) (range a b)))
|
rlm@0
|
111 ; (visual (map (vec (dna-fraction triple-cycle)) (range a b))))
|
rlm@0
|
112
|
rlm@0
|
113 ;; Differentiate the resulting function with respect to temperature.
|
rlm@0
|
114
|
rlm@0
|
115 (def test-temp-range (vec (range 260 380 0.1)))
|
rlm@0
|
116
|
rlm@0
|
117 (defn trans-print [& args] (apply println args) args)
|
rlm@0
|
118
|
rlm@0
|
119 ;; (def search
|
rlm@0
|
120 ;; #(for [x (range -1000 1000 100) y (range -1000 1000 100)]
|
rlm@0
|
121 ;; (do
|
rlm@0
|
122 ;; (print [x y])
|
rlm@0
|
123 ;; (try (trans-print (non-linear-model
|
rlm@0
|
124 ;; (ideal-dna-fraction 30e-6)
|
rlm@0
|
125 ;; test-dna-fraction test-temperature [x y]
|
rlm@0
|
126 ;; :method :newton-raphson))
|
rlm@0
|
127
|
rlm@0
|
128 ;; (catch Exception e (trans-print nil))))))
|
rlm@0
|
129
|
rlm@0
|
130
|
rlm@0
|
131 ;;(non-linear-model (ideal-dna-fraction 30e-6) test-dna-fraction test-temperature [-184 -71e3])
|
rlm@0
|
132
|
rlm@0
|
133 ;;(time (def answer (non-linear-model (ideal-dna-fraction 30e-6) test-dna-fraction test-temperature [-200 -90e3])))
|
rlm@0
|
134
|
rlm@0
|
135 ;;we're assuming that the SYBR Green 1 doesn't glow at all when there is only
|
rlm@0
|
136 ;;single stranded DNA, but it actually does. V_f only ever gets down to ~0.1 V
|
rlm@0
|
137 ;;even with everything melted.
|
rlm@0
|
138
|
rlm@0
|
139
|
rlm@0
|
140 ;;correct for thermal capacatance
|
rlm@0
|
141
|
rlm@0
|
142 (defn find-max [v]
|
rlm@0
|
143 (first (reduce (fn [a b] (if (> (last b) (last a)) b a)) (indexed v))))
|
rlm@0
|
144
|
rlm@0
|
145 (defn find-min [v]
|
rlm@0
|
146 (first (reduce (fn [a b] (if (< (last b) (last a)) b a)) (indexed v))))
|
rlm@0
|
147
|
rlm@0
|
148 (defn plot-thermal-delay []
|
rlm@0
|
149 (let [a (int 2.5e3) b (int 5e3)]
|
rlm@0
|
150 (dna-plot (assoc A-1
|
rlm@0
|
151 :temp (map (vec (:temp A-1)) (range a b))
|
rlm@0
|
152 :glow (map (vec (:glow A-1)) (range a b))))))
|
rlm@0
|
153
|
rlm@0
|
154 (defvar thermal-delay
|
rlm@0
|
155 (let [a (int 2.5e3) b (int 5e3)]
|
rlm@0
|
156 (-
|
rlm@0
|
157 (find-min (map (vec (:glow A-1)) (range a b)))
|
rlm@0
|
158 (find-max (map (vec (:temp A-1)) (range a b)))))
|
rlm@0
|
159 "There are two sources of thermal delay: thermal capacatance
|
rlm@0
|
160 of the cuvette holder and thermal capacatance from the time
|
rlm@0
|
161 it takes DNA to unwind.
|
rlm@0
|
162
|
rlm@0
|
163 Use A-1 to find the thermal delay. There was almost no delay
|
rlm@0
|
164 between heating and cooling for A-1, and thus the extremum of
|
rlm@0
|
165 temperature and glow should temporally coincide if there were
|
rlm@0
|
166 no delays. It turns out that there is about a 30 second lag!" )
|
rlm@0
|
167
|
rlm@0
|
168 (defn thermal-delay-correction*
|
rlm@0
|
169 "There is a certain delay between the RTD reporting changes in
|
rlm@0
|
170 temperature and the temperature of the DNA actually changing.
|
rlm@0
|
171 This shifts the temperature values so that the delay is eliminated."
|
rlm@0
|
172 [thermal-delay #^dna-data d]
|
rlm@0
|
173 (assoc d
|
rlm@0
|
174 :name (str (:name d) " |thermal-correction| ")
|
rlm@0
|
175 :temp (drop-last thermal-delay (:temp d))
|
rlm@0
|
176 :glow (drop thermal-delay (:glow d))))
|
rlm@0
|
177
|
rlm@0
|
178 (def thermal-delay-correction (partial thermal-delay-correction* thermal-delay))
|
rlm@0
|
179
|
rlm@0
|
180 ;;correct for photo-bleaching
|
rlm@0
|
181
|
rlm@0
|
182 (defn general-linear [[A B] x] (+ (* A x) B))
|
rlm@0
|
183
|
rlm@0
|
184 (defvar bleaching-coef nil
|
rlm@0
|
185 ;; (let [[A B]
|
rlm@0
|
186 ;; (:coefs
|
rlm@0
|
187 ;; (let [a (int 3e4) b (length (:glow pure-bleaching))]
|
rlm@0
|
188 ;; (non-linear-model
|
rlm@0
|
189 ;; general-linear
|
rlm@0
|
190 ;; (map (vec (:glow pure-bleaching)) (range a b))
|
rlm@0
|
191 ;; (range a b)
|
rlm@0
|
192 ;; [1 1])))]
|
rlm@0
|
193 ;; A)
|
rlm@0
|
194 "This function corrects V_f for photo-bleaching over a period of time.
|
rlm@0
|
195 Photo-bleaching was found to be linear(!) over quite a long period of time.
|
rlm@0
|
196 Data is taken from sample B-1, which I left in the holder for an hour
|
rlm@0
|
197 with no heating")
|
rlm@0
|
198
|
rlm@0
|
199 (defn bleaching-correction*
|
rlm@0
|
200 [bleaching-coef #^dna-data d]
|
rlm@0
|
201 ;;the correction is clibrated to apply to the raw data only
|
rlm@0
|
202 (assert (= (:glow-units d) "Raw Glow (Volts)"))
|
rlm@0
|
203 ;;we need to modulate this correction factor
|
rlm@0
|
204
|
rlm@0
|
205 (let [baseline (first (drop 3e4 (:temp pure-bleaching)))]
|
rlm@0
|
206 (assoc d
|
rlm@0
|
207 :name (str (:name d) " |bleaching-correction| ")
|
rlm@0
|
208 :glow (map (fn [[i n]] (- n (*
|
rlm@0
|
209 (/ n baseline)
|
rlm@0
|
210 ((partial general-linear [bleaching-coef 0]) i))))
|
rlm@0
|
211 (indexed (:glow d))))))
|
rlm@0
|
212
|
rlm@0
|
213 (def bleaching-correction (partial bleaching-correction* bleaching-coef))
|
rlm@0
|
214
|
rlm@0
|
215 ;; compensate for residual fluorescence when everything is melted
|
rlm@0
|
216
|
rlm@0
|
217 (defn residual-fluorescence
|
rlm@0
|
218 "the SYBR1 green dye doesn't completly stop glowing once everything
|
rlm@0
|
219 is melted."
|
rlm@0
|
220 [#^dna-data d]
|
rlm@0
|
221 (let [low (apply min (:glow d))
|
rlm@0
|
222 high (apply max (:glow d))]
|
rlm@0
|
223 (assoc d
|
rlm@0
|
224 :name (str (:name d) " |residual-fluorescence| ")
|
rlm@0
|
225 :glow (map (fn [f]
|
rlm@0
|
226 (- f (* (- 1 (/ f high)) low))) (:glow d)))))
|
rlm@0
|
227
|
rlm@0
|
228
|
rlm@0
|
229
|
rlm@0
|
230
|
rlm@0
|
231 ;; extract the cooling phase
|
rlm@0
|
232
|
rlm@0
|
233 (defn cooling-phase
|
rlm@0
|
234 "first stretch where the temperature is going down for a while
|
rlm@0
|
235 equivalent to where the f curve starts to recover"
|
rlm@0
|
236 [#^dna-data d]
|
rlm@0
|
237 (let [t (:temp d)
|
rlm@0
|
238 f (:glow d)
|
rlm@0
|
239 dividing-line (reduce (fn [a b] (if (< (a 1) (b 1)) a b)) (indexed f))
|
rlm@0
|
240 t-new (drop (dividing-line 0) t)
|
rlm@0
|
241 f-new (drop (dividing-line 0) f)]
|
rlm@0
|
242 (assoc d
|
rlm@0
|
243 :name (str (:name d) " |extract cooling| ")
|
rlm@0
|
244 :temp (vec t-new) :glow (vec f-new))))
|
rlm@0
|
245
|
rlm@0
|
246
|
rlm@0
|
247
|
rlm@0
|
248 ;; filter noise
|
rlm@0
|
249 ;;(http://clojure-log.n01se.net/date/2009-10-13.html by ambient
|
rlm@0
|
250 ;;this is such a cool transform!!!
|
rlm@0
|
251
|
rlm@0
|
252
|
rlm@0
|
253 (defn lo-pass
|
rlm@0
|
254 [d coll]
|
rlm@0
|
255 (reductions #(+ (* %2 d) (* %1 (- 1.0 d))) coll))
|
rlm@0
|
256
|
rlm@0
|
257
|
rlm@0
|
258 ;;)
|
rlm@0
|
259
|
rlm@0
|
260 (def reasonable-low-pass (partial lo-pass 1e-2))
|
rlm@0
|
261
|
rlm@0
|
262 (defn low-pass-correction
|
rlm@0
|
263 "filter the raw data through a low pass filter to remove
|
rlm@0
|
264 electronic noise"
|
rlm@0
|
265 [#^dna-data d]
|
rlm@0
|
266 (assert (= (:glow-units d) "Raw Glow (Volts)"))
|
rlm@0
|
267 (assert (= (:temp-units d) "Raw Temperature (Volts)"))
|
rlm@0
|
268 (assoc d
|
rlm@0
|
269 :name (str (:name d) " |low-pass| ")
|
rlm@0
|
270 :temp (reasonable-low-pass (:temp d))
|
rlm@0
|
271 :glow (reasonable-low-pass (:glow d))))
|
rlm@0
|
272
|
rlm@0
|
273
|
rlm@0
|
274
|
rlm@0
|
275
|
rlm@0
|
276
|
rlm@0
|
277 ;; correct for units
|
rlm@0
|
278
|
rlm@0
|
279 (defn RTD-voltage->temperature
|
rlm@0
|
280 "transform RTD voltage to temperature.
|
rlm@0
|
281 R_RTD = 1000 + 3.85 * (T - 273)
|
rlm@0
|
282 T = ((R_RTD - 1000) / 3.85 ) + 273
|
rlm@0
|
283 V_RTD = 15 * R_RTD / (R + R_RTD)
|
rlm@0
|
284 R_RTD = R / (15/V_RTD - 1)"
|
rlm@0
|
285 [coll]
|
rlm@0
|
286 (let [R 14.97e3
|
rlm@0
|
287 V 15]
|
rlm@0
|
288 (letfn [(T [R_RTD] (+ 273 (/ (- R_RTD 1000) 3.85))) ;; from lab notrd
|
rlm@0
|
289 (R_RTD [V_RTD] (/ R (- (/ V V_RTD) 1)))]
|
rlm@0
|
290 (map (comp T R_RTD) coll))))
|
rlm@0
|
291
|
rlm@0
|
292 (defn normalize
|
rlm@0
|
293 "Transform photodiode current to relative fluorescence."
|
rlm@0
|
294 [coll]
|
rlm@0
|
295 (let [low (apply min coll)
|
rlm@0
|
296 high (apply max coll)]
|
rlm@0
|
297 (map #(/ (- % low) (- high low)) coll)))
|
rlm@0
|
298
|
rlm@0
|
299 (defn unit-transform
|
rlm@0
|
300 "turn raw voltage into temperature and dna-concentration"
|
rlm@0
|
301 [#^dna-data d]
|
rlm@0
|
302 (assert (= (:glow-units d) "Raw Glow (Volts)"))
|
rlm@0
|
303 (assert (= (:temp-units d) "Raw Temperature (Volts)"))
|
rlm@0
|
304 (assoc d
|
rlm@0
|
305 :name (str (:name d) " |units| ")
|
rlm@0
|
306 :temp-units "Temperature (Kelvin)"
|
rlm@0
|
307 :glow-units "DNA-Fracton"
|
rlm@0
|
308 :temp (RTD-voltage->temperature (:temp d))
|
rlm@0
|
309 :glow (normalize (:glow d))))
|
rlm@0
|
310
|
rlm@0
|
311
|
rlm@0
|
312 (defn temp-unit-transform
|
rlm@0
|
313 "turn raw voltage into temperature and dna-concentration"
|
rlm@0
|
314 [#^dna-data d]
|
rlm@0
|
315 (assert (= (:temp-units d) "Raw Temperature (Volts)"))
|
rlm@0
|
316 (assoc d
|
rlm@0
|
317 :name (str (:name d) " |temp-units| ")
|
rlm@0
|
318 :temp-units "Temperature (Kelvin)"
|
rlm@0
|
319 :temp (RTD-voltage->temperature (:temp d))))
|
rlm@0
|
320
|
rlm@0
|
321 (defn glow-unit-transform
|
rlm@0
|
322 "turn raw voltage into temperature and dna-concentration"
|
rlm@0
|
323 [#^dna-data d]
|
rlm@0
|
324 (assert (= (:glow-units d) "Raw Glow (Volts)"))
|
rlm@0
|
325 (assoc d
|
rlm@0
|
326 :name (str (:name d) " |glow-units| ")
|
rlm@0
|
327 :glow-units "DNA-Fracton"
|
rlm@0
|
328 :glow (normalize (:glow d))))
|
rlm@0
|
329
|
rlm@0
|
330
|
rlm@0
|
331 (defn thermal-efficiency
|
rlm@0
|
332 [#^dna-data d]
|
rlm@0
|
333 (assert (= (:temp-units d) "Temperature (Kelvin)"))
|
rlm@0
|
334 (assert (= (:glow-units d) "Raw Glow (Volts)"))
|
rlm@0
|
335 (assoc d
|
rlm@0
|
336 :name (str (:name d) " |thermal-efficiency| ")
|
rlm@0
|
337 :glow nil))
|
rlm@0
|
338
|
rlm@0
|
339
|
rlm@0
|
340
|
rlm@0
|
341 (defn remove-first-values
|
rlm@0
|
342 "remove the first minute or so of data so that the initial temperature of
|
rlm@0
|
343 the cuvette itself doesn't matter"
|
rlm@0
|
344 [#^dna-data d]
|
rlm@0
|
345 (assoc d
|
rlm@0
|
346 :name (str (:name d) " |remove-first| ")
|
rlm@0
|
347 :temp (drop 600 (:temp d))
|
rlm@0
|
348 :glow (drop 600 (:glow d))))
|
rlm@0
|
349
|
rlm@0
|
350 (defn diff [coll] (map - (rest coll) coll ))
|
rlm@0
|
351 (defn temp-transform
|
rlm@0
|
352 "my assumption will be that at the very end of the temperature data, the
|
rlm@0
|
353 cuvette is exactly the same temperature as the RTD. Then, I use newton's
|
rlm@0
|
354 heat transfer and numerical integration to derive the heat of the cuvette
|
rlm@0
|
355 for previous times"
|
rlm@0
|
356 [k-value temperature-data]
|
rlm@0
|
357 (let [T1 (vec (reverse temperature-data))
|
rlm@0
|
358 ;; want fast indexed access, and instead of doing the problem in reverse
|
rlm@0
|
359 ;; just reverse the vector at the start, treat it normally, and
|
rlm@0
|
360 ;; then reverse at the end
|
rlm@0
|
361 initial-cuvette-temp (first T1)
|
rlm@0
|
362 initial-RTD-temp (first T1)
|
rlm@0
|
363 step (fn [prev-T2 current-T1]
|
rlm@0
|
364 (let [new-T2 (+ prev-T2 (* k-value (- current-T1 prev-T2)))]
|
rlm@0
|
365 new-T2))]
|
rlm@0
|
366 (reverse (reductions step initial-cuvette-temp T1))))
|
rlm@0
|
367
|
rlm@0
|
368 (def k-value 0.1)
|
rlm@0
|
369
|
rlm@0
|
370 (defn temp-corr-2 [k-value T1*]
|
rlm@0
|
371 (let [T1 (reverse T1*)]
|
rlm@0
|
372 (reverse (cons (first T1) (map + (map (partial * k-value) (diff T1)) (rest T1))))))
|
rlm@0
|
373
|
rlm@0
|
374 (defn temp-corr-3 [initial-cuvette-temp k-value T1]
|
rlm@0
|
375 (let [step
|
rlm@0
|
376 (fn [prev-T2 current-T1]
|
rlm@0
|
377 (+ prev-T2 (* k-value (- current-T1 prev-T2))))]
|
rlm@0
|
378 (rest (reductions step initial-cuvette-temp T1))))
|
rlm@0
|
379
|
rlm@0
|
380 (defn newton-temp-correction* [k-value #^dna-data d]
|
rlm@0
|
381 (assoc d
|
rlm@0
|
382 :name (str (:name d) " |newton-temp| ")
|
rlm@0
|
383 :temp (temp-transform k-value (:temp d))))
|
rlm@0
|
384 (defn newton-temp-correction2* [k-value #^dna-data d]
|
rlm@0
|
385 (assoc d
|
rlm@0
|
386 :name (str (:name d) " |newton-temp-reverse| ")
|
rlm@0
|
387 :temp (temp-corr-2 k-value (:temp d))))
|
rlm@0
|
388 (defn newton-temp-correction3* [init-temp k-value #^dna-data d]
|
rlm@0
|
389 (assoc d
|
rlm@0
|
390 :name (str (:name d) " |newton-temp-reverse| ")
|
rlm@0
|
391 :temp (temp-corr-3 init-temp k-value (:temp d))))
|
rlm@0
|
392
|
rlm@0
|
393 (defn turn-e [n] (format "%5.2e" n))
|
rlm@0
|
394
|
rlm@0
|
395 (defn newton-temp-correction4*
|
rlm@0
|
396 "now with AIR!!!"
|
rlm@0
|
397 [k-block-sample k-sample-air initial-sample-temp room-temp #^dna-data d]
|
rlm@0
|
398 (let [old-temp (:temp d)
|
rlm@0
|
399 step (fn [prev-T2 current-T1]
|
rlm@0
|
400 (+ prev-T2
|
rlm@0
|
401 (* k-block-sample (- current-T1 prev-T2))
|
rlm@0
|
402 (* k-sample-air (- room-temp prev-T2))))
|
rlm@0
|
403 new-temp (rest (reductions step initial-sample-temp old-temp))]
|
rlm@0
|
404 (assoc d
|
rlm@0
|
405 :name (str (:name d)
|
rlm@0
|
406 " |temp-adjusted {"
|
rlm@0
|
407 "initial-sample:" (turn initial-sample-temp) " "
|
rlm@0
|
408 "room-temp:" (turn room-temp) " "
|
rlm@0
|
409 "k-block-sample:" (turn-e k-block-sample) " "
|
rlm@0
|
410 "k-sample-air:" (turn-e k-sample-air) "}")
|
rlm@0
|
411 :temp new-temp)))
|
rlm@0
|
412
|
rlm@0
|
413 (defn plot-adjusted-temp
|
rlm@0
|
414 "visualize the temperature correction"
|
rlm@0
|
415 [adjuster #^dna-data d]
|
rlm@0
|
416 (view
|
rlm@0
|
417 (add-lines
|
rlm@0
|
418 (visual (:temp d))
|
rlm@0
|
419 (vec (range (count (:temp d))))
|
rlm@0
|
420 (vec (:temp (adjuster d))))))
|
rlm@0
|
421
|
rlm@0
|
422
|
rlm@0
|
423
|
rlm@0
|
424
|
rlm@0
|
425
|
rlm@0
|
426 ;; process data
|
rlm@0
|
427
|
rlm@0
|
428 (defn process [#^dna-data d]
|
rlm@0
|
429 (cooling-phase
|
rlm@0
|
430 (unit-transform
|
rlm@0
|
431 (residual-fluorescence
|
rlm@0
|
432 (bleaching-correction
|
rlm@0
|
433 (low-pass-correction
|
rlm@0
|
434 ;(thermal-delay-correction
|
rlm@0
|
435 d))))))
|
rlm@0
|
436
|
rlm@0
|
437
|
rlm@0
|
438 (defn first-pass
|
rlm@0
|
439 "after this, the only thing stopping the melting curve from
|
rlm@0
|
440 being a straight line is thermal capacatance"
|
rlm@0
|
441 [#^dna-data data]
|
rlm@0
|
442 (low-pass-correction
|
rlm@0
|
443 ;(remove-first-values
|
rlm@0
|
444 (bleaching-correction
|
rlm@0
|
445 data)))
|
rlm@0
|
446
|
rlm@0
|
447
|
rlm@0
|
448
|
rlm@0
|
449 (defn explore [k-value]
|
rlm@0
|
450 (dna-melting (newton-temp-correction* k-value (first-pass A-2))))
|
rlm@0
|
451
|
rlm@0
|
452
|
rlm@0
|
453 (defn explore2 [k-value]
|
rlm@0
|
454 (dna-melting (newton-temp-correction2* k-value (first-pass A-2))))
|
rlm@0
|
455
|
rlm@0
|
456
|
rlm@0
|
457 (defn explore3 [init-temp k-value]
|
rlm@0
|
458 (dna-melting (newton-temp-correction3* init-temp k-value (unit-transform (first-pass A-2)))))
|
rlm@0
|
459
|
rlm@0
|
460
|
rlm@0
|
461
|
rlm@0
|
462 ;; get the three characteristic values we want.
|
rlm@0
|
463
|
rlm@0
|
464 (defn ideal-dna-fraction* [concentration [delta-S delta-H] temperature]
|
rlm@0
|
465 (let [R 8.31447215
|
rlm@0
|
466 t temperature
|
rlm@0
|
467 Ct concentration
|
rlm@0
|
468 K_eq (exp (- (/ delta-S R) (/ delta-H (* R t))))
|
rlm@0
|
469 f (/
|
rlm@0
|
470 (+ 1 (* Ct K_eq)
|
rlm@0
|
471 (- (sqrt (+ 1 (* 2 Ct K_eq)))))
|
rlm@0
|
472 (* Ct K_eq))]
|
rlm@0
|
473 f))
|
rlm@0
|
474
|
rlm@0
|
475 (def ideal-dna-fraction (partial ideal-dna-fraction* 30e-6))
|
rlm@0
|
476
|
rlm@0
|
477 (defvar starting-values
|
rlm@0
|
478 ;; [-184 -71e3]
|
rlm@0
|
479 [-904 -32e4]
|
rlm@0
|
480 "chosen from theoritical estimation")
|
rlm@0
|
481
|
rlm@0
|
482
|
rlm@0
|
483 (defn analyze* [iterations #^dna-data d]
|
rlm@0
|
484 (let [best-fit-results
|
rlm@0
|
485 (non-linear-model
|
rlm@0
|
486 ideal-dna-fraction
|
rlm@0
|
487 (:glow d)
|
rlm@0
|
488 (:temp d)
|
rlm@0
|
489 starting-values
|
rlm@0
|
490 :max-iter iterations)]
|
rlm@0
|
491 best-fit-results))
|
rlm@0
|
492
|
rlm@0
|
493
|
rlm@0
|
494 (def analyze (partial analyze* 200))
|
rlm@0
|
495
|
rlm@0
|
496 (defn plot-results [d]
|
rlm@0
|
497 (let [ans (analyze d)]
|
rlm@0
|
498 (view (doto (xy-plot (:x ans) (:y ans))
|
rlm@0
|
499 (set-title (:name d))
|
rlm@0
|
500 (add-lines (:x ans) (:fitted ans))))
|
rlm@0
|
501 ans))
|
rlm@0
|
502
|
rlm@0
|
503 (defn zzz [d [A B]]
|
rlm@0
|
504 (let [gar (vec (map (partial ideal-dna-fraction [A B]) (:temp d)))]
|
rlm@0
|
505 (view (doto (xy-plot (:temp d) (:glow d))
|
rlm@0
|
506 (set-title (:name d))
|
rlm@0
|
507 (add-lines (:temp d) gar)))))
|
rlm@0
|
508
|
rlm@0
|
509
|
rlm@0
|
510
|
rlm@0
|
511
|
rlm@0
|
512 (defn index-filter [pred coll]
|
rlm@0
|
513 (when pred
|
rlm@0
|
514 (for [[idx elt] (indexed coll) :when (pred elt)] idx)))
|
rlm@0
|
515
|
rlm@0
|
516
|
rlm@0
|
517
|
rlm@0
|
518
|
rlm@0
|
519 ;; we're going to estimate two correction factors,
|
rlm@0
|
520 ;; a constant capacative temperature delay,
|
rlm@0
|
521 ;; and a linear decrease in brightness. We'll do this
|
rlm@0
|
522 ;; before unit-transforms.
|
rlm@0
|
523
|
rlm@0
|
524 (defn zipmap-collect
|
rlm@0
|
525 "Returns a map with the keys mapped to the corresponding vals.
|
rlm@0
|
526 adds extra values to a list"
|
rlm@0
|
527 [keys vals]
|
rlm@0
|
528 (loop [map {}
|
rlm@0
|
529 ks (seq keys)
|
rlm@0
|
530 vs (seq vals)]
|
rlm@0
|
531 (if (and ks vs)
|
rlm@0
|
532 (recur (if (contains? map (first ks))
|
rlm@0
|
533 (assoc map (first ks) (cons (first vs) (map (first ks))))
|
rlm@0
|
534 (assoc map (first ks) (list (first vs))))
|
rlm@0
|
535 (next ks)
|
rlm@0
|
536 (next vs))
|
rlm@0
|
537 map)))
|
rlm@0
|
538
|
rlm@0
|
539
|
rlm@0
|
540 (defn little-process [thermal-cap bleaching #^dna-data d]
|
rlm@0
|
541 ((partial thermal-delay-correction thermal-cap)
|
rlm@0
|
542 ((partial bleaching-correction bleaching) d)))
|
rlm@0
|
543
|
rlm@0
|
544 (def little-process (memoize little-process))
|
rlm@0
|
545
|
rlm@0
|
546 ;(defn lookup-map [#^dna-data d]
|
rlm@0
|
547
|
rlm@0
|
548
|
rlm@0
|
549 (defn converge-function
|
rlm@0
|
550 "this is the vertical discrepancy after applying the two transforms"
|
rlm@0
|
551 [#^dna-data d [thermal-cap bleaching] x]
|
rlm@0
|
552 (let [new-d (little-process thermal-cap bleaching d)
|
rlm@0
|
553 spreads
|
rlm@0
|
554 (map (vec (:glow new-d))
|
rlm@0
|
555 (index-filter (partial = x) (vec (:temp new-d))))
|
rlm@0
|
556 spread (- (apply max spreads) (apply min spreads))]
|
rlm@0
|
557 spread))
|
rlm@0
|
558
|
rlm@0
|
559
|
rlm@0
|
560
|
rlm@0
|
561
|
rlm@0
|
562
|
rlm@0
|
563
|
rlm@0
|
564
|
rlm@0
|
565
|
rlm@0
|
566
|
rlm@0
|
567
|
rlm@0
|
568 (defn blah [k] (comp
|
rlm@0
|
569 ( partial newton-temp-correction3* 280 k ) temp-unit-transform first-pass))
|
rlm@0
|
570
|
rlm@0
|
571 ;(def ggg (comp
|
rlm@0
|
572 ; dna-melting ( partial thermal-delay-correction* 120 ) (blah 0.0018)))
|
rlm@0
|
573
|
rlm@0
|
574
|
rlm@0
|
575
|
rlm@0
|
576
|
rlm@0
|
577
|
rlm@0
|
578 ;; Ensure that the melting function is uniquely valued by
|
rlm@0
|
579 ;; combining samples with identical temperature values.
|
rlm@0
|
580 (defn cooling-function
|
rlm@0
|
581 "just put it in a map for this!"
|
rlm@0
|
582 [source]
|
rlm@0
|
583 (let [data (cooling-phase source)
|
rlm@0
|
584 lookup (apply sorted-map
|
rlm@0
|
585 (interleave
|
rlm@0
|
586 (:temperature data) (:dna-fraction data)))]
|
rlm@0
|
587 {:temperature (vec (keys lookup))
|
rlm@0
|
588 :dna-fraction (vec (map lookup (keys lookup)))}))
|
rlm@0
|
589
|
rlm@0
|
590
|
rlm@0
|
591 (defn turn [n]
|
rlm@0
|
592 (format "%5.2f" n))
|
rlm@0
|
593
|
rlm@0
|
594 (defn bleach-correction-2 [bleaching-constant #^dna-data d]
|
rlm@0
|
595 (assert (= (:glow-units d) "Raw Glow (Volts)"))
|
rlm@0
|
596 (let [old-glow (vec (:glow d))
|
rlm@0
|
597 integral-old-glow (vec (reductions + old-glow))
|
rlm@0
|
598 correction (fn [t] (+
|
rlm@0
|
599 (old-glow t)
|
rlm@0
|
600 (* bleaching-constant (integral-old-glow t))))
|
rlm@0
|
601 new-glow (map + (map correction (range 0 (count old-glow))) old-glow)]
|
rlm@0
|
602 ;;glow is related to the number of particles,
|
rlm@0
|
603 ;;and the number destroyed is some fraction of
|
rlm@0
|
604 ;;the number present and glowing.
|
rlm@0
|
605 (assoc d
|
rlm@0
|
606 :name (str (:name d) " |bleaching-corrected constant = " (turn bleaching-constant) "| ")
|
rlm@0
|
607 :glow new-glow)))
|
rlm@0
|
608
|
rlm@0
|
609
|
rlm@0
|
610
|
rlm@0
|
611 (defn ultimate-correction
|
rlm@0
|
612 ([k-block-sample k-sample-air initial-sample-temp room-temp bleach-constant
|
rlm@0
|
613 #^dna-data d
|
rlm@0
|
614 ]
|
rlm@0
|
615 ((comp
|
rlm@0
|
616 (partial newton-temp-correction4*
|
rlm@0
|
617 k-block-sample k-sample-air initial-sample-temp room-temp)
|
rlm@0
|
618
|
rlm@0
|
619 temp-unit-transform
|
rlm@0
|
620 (partial bleach-correction-2 bleach-constant)
|
rlm@0
|
621 low-pass-correction) d))
|
rlm@0
|
622
|
rlm@0
|
623 ([{:keys [test-k-sample-air
|
rlm@0
|
624 test-k-block-sample
|
rlm@0
|
625 test-sample-temp
|
rlm@0
|
626 test-room-temp
|
rlm@0
|
627 test-bleach]}
|
rlm@0
|
628 #^dna-data d]
|
rlm@0
|
629 (ultimate-correction
|
rlm@0
|
630 (mean test-k-block-sample) (mean test-k-sample-air)
|
rlm@0
|
631 (mean test-sample-temp) (mean test-room-temp) (mean test-bleach) d)))
|
rlm@0
|
632
|
rlm@0
|
633
|
rlm@0
|
634 ;;((uui 0.0020 0.0001 280 293) (bleach-correction-2 1.1e-5 ?-3))
|
rlm@0
|
635
|
rlm@0
|
636
|
rlm@0
|
637 ;;; this is for ?-3 after 35 iterations
|
rlm@0
|
638 (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
|
639
|
rlm@0
|
640
|
rlm@0
|
641 ;;; C-1 after 20 iterations
|
rlm@0
|
642 (def opt2
|
rlm@0
|
643 {: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
|
644
|
rlm@0
|
645
|
rlm@0
|
646 ;;(time (def result-?-3 (doall (take 10 (iterate (partial improve ?-3) start-values-overestimates)))))
|
rlm@0
|
647
|
rlm@0
|
648
|
rlm@0
|
649 ;;(doall (map (comp dna-melting #(ultimate-correction % ?-3)) result-?-3))
|
rlm@0
|
650
|
rlm@0
|
651
|
rlm@0
|
652 (defn subsample
|
rlm@0
|
653 "subsample v every n points"
|
rlm@0
|
654 [n v]
|
rlm@0
|
655
|
rlm@0
|
656 (vec (map (vec v) (vec (range 0 (count (vec v)) n)))))
|
rlm@0
|
657
|
rlm@0
|
658
|
rlm@0
|
659 (defn max-difference [v] (- (apply max v) (apply min v)))
|
rlm@0
|
660
|
rlm@0
|
661 (defn alignment [#^dna-data d]
|
rlm@0
|
662 (let [distribution (vec (vals
|
rlm@0
|
663 (apply sorted-map
|
rlm@0
|
664 (interleave (:temp d) (:glow d)))))]
|
rlm@0
|
665 (mean (map variance (partition 10 distribution)))))
|
rlm@0
|
666
|
rlm@0
|
667 (defn garigh [d]
|
rlm@0
|
668 (let [distribution (vec (vals
|
rlm@0
|
669 (apply sorted-map
|
rlm@0
|
670 (interleave (:temp d) (:glow d)))))]
|
rlm@0
|
671 (view (visual (vec (map max-difference (partition 100 distribution)))))))
|
rlm@0
|
672
|
rlm@0
|
673
|
rlm@0
|
674 (defn alignment-4 [#^dna-data d]
|
rlm@0
|
675 (let [distribution (vec (vals
|
rlm@0
|
676 (apply sorted-map
|
rlm@0
|
677 (interleave (:temp d) (:glow d)))))]
|
rlm@0
|
678 (reduce + (map sq (map max-difference (partition 100 distribution))))))
|
rlm@0
|
679
|
rlm@0
|
680
|
rlm@0
|
681 (defn alignment-2 [#^dna-data d]
|
rlm@0
|
682 (let [distribution
|
rlm@0
|
683 (vec
|
rlm@0
|
684 (vals
|
rlm@0
|
685 (apply sorted-map
|
rlm@0
|
686 (interleave (subsample 2 (:temp d)) (subsample 2(:glow d))))))]
|
rlm@0
|
687 (reduce + (map max-difference (partition 30 distribution)))))
|
rlm@0
|
688 (defn alignment-3 [#^dna-data d]
|
rlm@0
|
689 (let [distribution
|
rlm@0
|
690 (vec
|
rlm@0
|
691 (vals
|
rlm@0
|
692 (apply sorted-map
|
rlm@0
|
693 (interleave (subsample 2 (:temp d)) (subsample 2(:glow d))))))]
|
rlm@0
|
694 (reduce max (map max-difference (partition 30 distribution)))))
|
rlm@0
|
695
|
rlm@0
|
696
|
rlm@0
|
697 (comment
|
rlm@0
|
698 ;; this takes a while
|
rlm@0
|
699 (def opts (zipmap all-samples (map #(last (take 20 (iterate (partial improve %) start-values))) all-samples))))
|
rlm@0
|
700
|
rlm@0
|
701 ;(def opts (zipmap all-samples opts-vals))
|
rlm@0
|
702
|
rlm@0
|
703 (defn mean-map [v s]
|
rlm@0
|
704 (map mean (map (partial vector v) s)))
|
rlm@0
|
705
|
rlm@0
|
706
|
rlm@0
|
707 (def start-values
|
rlm@0
|
708 {:test-k-sample-air [0.000001 0.001]
|
rlm@0
|
709 :test-k-block-sample [0.001 0.03]
|
rlm@0
|
710 :test-room-temp [290 300]
|
rlm@0
|
711 :test-sample-temp [285 310]
|
rlm@0
|
712 :test-bleach [0.00001 0.0009]})
|
rlm@0
|
713
|
rlm@0
|
714 (def start-values-2
|
rlm@0
|
715 {:test-k-sample-air [1e-15 0.0005]
|
rlm@0
|
716 :test-k-block-sample [0.001 0.002]
|
rlm@0
|
717 :test-room-temp [290 300]
|
rlm@0
|
718 :test-sample-temp [290 300]
|
rlm@0
|
719 :test-bleach [1e-7 0.0005]})
|
rlm@0
|
720
|
rlm@0
|
721
|
rlm@0
|
722
|
rlm@0
|
723 (def start-values-3
|
rlm@0
|
724 {:test-k-block-sample [0.0005 0.00125 0.002]
|
rlm@0
|
725 :test-k-sample-air [0.0000001 0.000001 0.0000019]
|
rlm@0
|
726 :test-sample-temp [280 300]
|
rlm@0
|
727 :test-room-temp [300 306]
|
rlm@0
|
728 :test-bleach [0.000001 0.00001]})
|
rlm@0
|
729
|
rlm@0
|
730 (defn condense [[a b :as v]]
|
rlm@0
|
731 (if (< (/ (abs (- (apply max v) (apply min v))) (apply max v)) 1e-2) (vector (mean v)) v))
|
rlm@0
|
732
|
rlm@0
|
733 (defn mass-shift
|
rlm@0
|
734 [value v]
|
rlm@0
|
735 (let [frac 2/3]
|
rlm@0
|
736 (cond (>= (count v) 2)
|
rlm@0
|
737 (let [[a b] v]
|
rlm@0
|
738 (condense (vec (map #(float (+ (* frac %) (* (- 1 frac) value))) v))))
|
rlm@0
|
739 (= (count v) 1)
|
rlm@0
|
740 v)))
|
rlm@0
|
741
|
rlm@0
|
742
|
rlm@0
|
743 (defn improve* [#^dna-data d
|
rlm@0
|
744 {:keys [test-k-sample-air
|
rlm@0
|
745 test-k-block-sample
|
rlm@0
|
746 test-room-temp
|
rlm@0
|
747 test-sample-temp
|
rlm@0
|
748 test-bleach] :as domain}]
|
rlm@0
|
749 (do (print \newline)
|
rlm@0
|
750 (reduce (fn [m1 m2] (if (< (:align m1) (:align m2)) m1 m2))
|
rlm@0
|
751
|
rlm@0
|
752
|
rlm@0
|
753 (for [k-sample-air (condense test-k-sample-air)
|
rlm@0
|
754 k-block-sample (condense test-k-block-sample)
|
rlm@0
|
755 room-temp (condense test-room-temp)
|
rlm@0
|
756 initial-sample-temp (condense test-sample-temp)
|
rlm@0
|
757 bleach-constant (condense test-bleach)]
|
rlm@0
|
758
|
rlm@0
|
759 (do
|
rlm@0
|
760 (print ".")
|
rlm@0
|
761 (let [align
|
rlm@0
|
762 (alignment-2
|
rlm@0
|
763 (ultimate-correction
|
rlm@0
|
764 k-block-sample k-sample-air initial-sample-temp room-temp
|
rlm@0
|
765 bleach-constant d))
|
rlm@0
|
766 align (if (Double/isNaN align) Double/POSITIVE_INFINITY align)]
|
rlm@0
|
767 {:test-k-sample-air (mass-shift k-sample-air test-k-sample-air)
|
rlm@0
|
768 :test-k-block-sample(mass-shift k-block-sample test-k-block-sample)
|
rlm@0
|
769 :test-room-temp (mass-shift room-temp test-room-temp)
|
rlm@0
|
770 :test-sample-temp (mass-shift initial-sample-temp test-sample-temp)
|
rlm@0
|
771 :test-bleach (mass-shift bleach-constant test-bleach)
|
rlm@0
|
772 :align align}))))))
|
rlm@0
|
773
|
rlm@0
|
774
|
rlm@0
|
775
|
rlm@0
|
776 (def println-repl (bound-fn [& args] (apply println args)))
|
rlm@0
|
777 (def print-repl (bound-fn [& args] (apply print args)))
|
rlm@0
|
778
|
rlm@0
|
779
|
rlm@0
|
780
|
rlm@0
|
781 (defn improve [#^dna-data d
|
rlm@0
|
782 {:keys [test-k-block-sample
|
rlm@0
|
783 test-k-sample-air
|
rlm@0
|
784 test-sample-temp
|
rlm@0
|
785 test-room-temp
|
rlm@0
|
786 test-bleach] :as domain}]
|
rlm@0
|
787 (reduce
|
rlm@0
|
788 (fn [m1 m2] (if (< (:align m1) (:align m2)) m1 m2))
|
rlm@0
|
789 (let [return-hash
|
rlm@0
|
790 (fn
|
rlm@0
|
791 [[k-block-sample k-sample-air initial-sample-temp
|
rlm@0
|
792 room-temp bleach-constant]]
|
rlm@0
|
793 (do
|
rlm@0
|
794 (print-repl ".")
|
rlm@0
|
795 (let [align
|
rlm@0
|
796 (alignment-4
|
rlm@0
|
797 (ultimate-correction
|
rlm@0
|
798 k-block-sample k-sample-air initial-sample-temp room-temp
|
rlm@0
|
799 bleach-constant d))
|
rlm@0
|
800 align (if (Double/isNaN align) Double/POSITIVE_INFINITY align)]
|
rlm@0
|
801 {:test-k-sample-air (mass-shift k-sample-air test-k-sample-air)
|
rlm@0
|
802 :test-k-block-sample(mass-shift k-block-sample test-k-block-sample)
|
rlm@0
|
803 :test-room-temp (mass-shift room-temp test-room-temp)
|
rlm@0
|
804 :test-sample-temp (mass-shift initial-sample-temp test-sample-temp)
|
rlm@0
|
805 :test-bleach (mass-shift bleach-constant test-bleach)
|
rlm@0
|
806 :align align})))]
|
rlm@0
|
807 (pmap return-hash
|
rlm@0
|
808 (cartesian-product
|
rlm@0
|
809 test-k-block-sample
|
rlm@0
|
810 test-k-sample-air
|
rlm@0
|
811 test-sample-temp
|
rlm@0
|
812 test-room-temp
|
rlm@0
|
813 test-bleach)))))
|
rlm@0
|
814
|
rlm@0
|
815 (def improve (memoize improve))
|
rlm@0
|
816
|
rlm@0
|
817 (defn iterate-test [n #^dna-data d]
|
rlm@0
|
818 (last (take n (iterate (partial improve d) start-values-3))))
|
rlm@0
|
819
|
rlm@0
|
820
|
rlm@0
|
821 (defn nested-vals [key coll]
|
rlm@0
|
822 (for [x (tree-seq coll? seq coll) :when (contains? x key)]
|
rlm@0
|
823 (get x key)))
|
rlm@0
|
824
|
rlm@0
|
825
|
rlm@0
|
826 (defn quantum-efficiency* [q-constant #^dna-data d]
|
rlm@0
|
827 ;; 0.8 at room temperature
|
rlm@0
|
828 ;; quantum efficiency is somehow
|
rlm@0
|
829 ;; inversely proportional to temperature
|
rlm@0
|
830 (let [correction
|
rlm@0
|
831
|
rlm@0
|
832
|
rlm@0
|
833
|
rlm@0
|
834 (fn [t g] (/ g (* 0.8 (exp (* q-constant (- 292 t))))))]
|
rlm@0
|
835 (assoc d
|
rlm@0
|
836 :name (str (:name d) " |quantum-correction| ")
|
rlm@0
|
837 :glow (map correction (:temp d) (:glow d)))))
|
rlm@0
|
838
|
rlm@0
|
839 (def quantum-efficiency (partial quantum-efficiency* 1.2e-2))
|
rlm@0
|
840
|
rlm@0
|
841 (defn exp-decay [[A] x]
|
rlm@0
|
842 (* 0.8 (exp (* A (- 298 x)))))
|
rlm@0
|
843
|
rlm@0
|
844 ;(defn final-join [#^dna-data d]
|
rlm@0
|
845
|
rlm@0
|
846
|
rlm@0
|
847
|
rlm@0
|
848 (def join-opts
|
rlm@0
|
849 {
|
rlm@0
|
850 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
|
851 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
|
852 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
|
853
|
rlm@0
|
854 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
|
855 B-3
|
rlm@0
|
856 {: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
|
857
|
rlm@0
|
858 C-1
|
rlm@0
|
859 {: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
|
860 C-2
|
rlm@0
|
861 {: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
|
862 C-3
|
rlm@0
|
863 {: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
|
864
|
rlm@0
|
865 ?-1
|
rlm@0
|
866 {: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
|
867 ?-3
|
rlm@0
|
868 {: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
|
869
|
rlm@0
|
870 })
|
rlm@0
|
871
|
rlm@0
|
872
|
rlm@0
|
873 (defn final-join [#^dna-data d]
|
rlm@0
|
874
|
rlm@0
|
875 (let [
|
rlm@0
|
876 m (apply sorted-map (interleave (:temp d) (:glow d)))
|
rlm@0
|
877 distribution (vec (vals m))
|
rlm@0
|
878 ordinates (vec (keys m))
|
rlm@0
|
879 ]
|
rlm@0
|
880 (assoc d
|
rlm@0
|
881 :name (str (:name d) " |joined| ")
|
rlm@0
|
882 :temp (map (fn [region] (mean [(apply min region) (apply max region)]))
|
rlm@0
|
883 (partition 30 ordinates))
|
rlm@0
|
884
|
rlm@0
|
885 :glow (map (fn [region] (mean [(apply min region) (apply max region)]))
|
rlm@0
|
886 (partition 30 distribution)))))
|
rlm@0
|
887
|
rlm@0
|
888
|
rlm@0
|
889 (defn melting-point [#^dna-data d]
|
rlm@0
|
890 (let [t (drop 150 (:temp d))
|
rlm@0
|
891 g (drop 150 (:glow d))]
|
rlm@0
|
892 (nth t (find-min (map / (diff g) (diff t) )))))
|
rlm@0
|
893
|
rlm@0
|
894
|
rlm@0
|
895 (def target-dir (file-str "/home/r/MIT/6/6.122/DNA/output"))
|
rlm@0
|
896
|
rlm@0
|
897 (defn x-form-dna [#^dna-data d]
|
rlm@0
|
898 (final-join (ultimate-correction (join-opts d) d)))
|
rlm@0
|
899
|
rlm@0
|
900 (defn save-raw-dna-curve [#^dna-data d]
|
rlm@0
|
901 (save (dna-melting d) (str target-dir File/separator (:name d) ".png")))
|
rlm@0
|
902
|
rlm@0
|
903 (defn save-dna-curve-with-temp [#^dna-data d*]
|
rlm@0
|
904 (let [
|
rlm@0
|
905
|
rlm@0
|
906 melt (melting-point d*)
|
rlm@0
|
907 d (assoc d* :name (str (:name d*) " |M.P. = " (turn melt) " |"))
|
rlm@0
|
908 ]
|
rlm@0
|
909 (save (add-lines (dna-melting d) (repeat (count (:glow d)) melt) (:glow d))
|
rlm@0
|
910 (str target-dir File/separator (:name d) ".png"))))
|
rlm@0
|
911
|
rlm@0
|
912
|
rlm@0
|
913
|
rlm@0
|
914 (defn process-all []
|
rlm@0
|
915 (doall (pmap save-dna-curve all-samples))
|
rlm@0
|
916 (doall (pmap (comp save-dna-curve-with-temp x-form-dna) all-samples)))
|
rlm@0
|
917
|
rlm@0
|
918
|
rlm@0
|
919 (def ?-temps (vec (map (comp melting-point x-form-dna) [?-1 ?-3])))
|
rlm@0
|
920 (def A-temps (vec (map (comp melting-point x-form-dna) [A-1 A-2 A-3])))
|
rlm@0
|
921 (def B-temps (vec (map (comp melting-point x-form-dna) [B-2 B-3])))
|
rlm@0
|
922 (def C-temps (vec (map (comp melting-point x-form-dna) [C-1 C-2 C-3])))
|
rlm@0
|
923
|
rlm@0
|
924 (defn report [temps]
|
rlm@0
|
925
|
rlm@0
|
926 (sorted-map :data temps :std (sqrt (variance temps)) :mean (mean temps)))
|
rlm@0
|
927
|
rlm@0
|
928
|
rlm@0
|
929 ;;; now it's time to start data reporting
|
rlm@0
|
930
|
rlm@0
|
931 ;(defn process [#^dna-data d]
|
rlm@0
|
932 ; (ultimate-correction (opts d) d)
|
rlm@0
|
933
|
rlm@0
|
934
|
rlm@0
|
935
|
rlm@0
|
936 (report A-temps)
|
rlm@0
|
937 {:data [348.8771481779248 349.3989434763391 348.76079920968573],
|
rlm@0
|
938 :mean 349.0122969546499,
|
rlm@0
|
939 :std 0.33986161910548834}
|
rlm@0
|
940
|
rlm@0
|
941 (report B-temps)
|
rlm@0
|
942 {:data [349.0655539116902 349.92155530431506],
|
rlm@0
|
943 :mean 349.4935546080026,
|
rlm@0
|
944 :std 0.6052843894485812}
|
rlm@0
|
945
|
rlm@0
|
946 (report C-temps)
|
rlm@0
|
947 {:data [349.92341967511976 349.68832687337056 348.2872189461977],
|
rlm@0
|
948 :mean 349.299655164896,
|
rlm@0
|
949 :std 0.884639745344559}
|
rlm@0
|
950
|
rlm@0
|
951 (report ?-temps)
|
rlm@0
|
952 {:data [350.4525305192239 349.87377093318975],
|
rlm@0
|
953 :mean 350.1631507262068,
|
rlm@0
|
954 :std 0.4092448279511269}
|
rlm@0
|
955
|
rlm@0
|
956
|