| blob | c0843efb26b6b03b51f17b4a30b616483379610f |
1 (ns lineage.identify
2 (:import
3 (lineage LineageClassifier)
4 (java.io InputStreamReader FileInputStream PushbackReader)
5 (javax.vecmath Point3d)
6 (javax.swing JTable JFrame JScrollPane JPanel JLabel BoxLayout JPopupMenu JMenuItem)
7 (javax.swing.table AbstractTableModel DefaultTableCellRenderer)
8 (java.awt.event MouseAdapter ActionListener)
9 (java.awt Color Dimension Component)
10 (mpicbg.models AffineModel3D)
11 (ij.measure Calibration)
12 (ini.trakem2.utils Utils)
13 (ini.trakem2.display Display Display3D LayerSet)
14 (ini.trakem2.analysis Compare)
15 (ini.trakem2.vector VectorString3D Editions))
16 (:use [clojure.set :only (intersection)]))
18 (defmacro report
19 [& args]
20 `(ij.IJ/log (str ~@args)))
22 (defn- as-VectorString3D
23 "Convert a map of {\"name\" {:x (...) :y (...) :z (...)}} into a map of {\"name\" (VectorString3D. ...)}"
24 [SATs]
25 (map
26 (fn [e]
27 {(key e)
28 (let [data (val e)]
29 (VectorString3D. (into-array Double/TYPE (data :x))
30 (into-array Double/TYPE (data :y))
31 (into-array Double/TYPE (data :z))
32 false))})
33 SATs))
35 (defn- load-SAT-lib-as-vs
36 "Returns the SATs library with VectorString3D instances"
37 [filepath]
38 (reduce
39 (fn [m e]
40 (let [label (key e)
41 data (val e)]
42 (assoc m
43 label
44 (assoc data :SATs (as-VectorString3D (data :SATs))))))
45 {}
46 (read
47 (PushbackReader.
48 (InputStreamReader. (FileInputStream. filepath))
49 4096))))
51 (defn- fids-as-Point3d
52 "Convert a map like {\"blva_joint\" [70.62114008906023 140.07819545137772 -78.72010436407604] ...}
53 into a sorted map of {\"blva_joint\" (Point3d ...)}"
54 [fids]
55 (reduce
56 (fn [m e]
57 (let [xyz (val e)]
58 (assoc m (key e) (Point3d. (double (get xyz 0))
59 (double (get xyz 1))
60 (double (get xyz 2))))))
61 {}
62 fids))
64 (defn- register-SATs
65 "Returns a map of SAT name vs VectorString3D,
66 where the VectorString3D is registered to the target fids."
67 [brain-label brain-data target-fids]
68 (let [source-fids (fids-as-Point3d (brain-data :fiducials))
69 SATs (into (sorted-map) (brain-data :SATs))
70 common-fid-keys (intersection (into #{} (keys source-fids)) (into #{} (keys target-fids)))
71 vs (Compare/transferVectorStrings (vals SATs)
72 (vals (into (sorted-map) (select-keys source-fids common-fid-keys)))
73 (vals (into (sorted-map) (select-keys target-fids common-fid-keys)))
74 AffineModel3D)]
75 (zipmap
76 (map #(str (.replaceAll % "\\[.*\\] " " [") \space brain-label \]) (keys SATs))
77 vs)))
79 (defn- load-mb
80 [SAT-lib reference-brain]
81 (let [r (re-pattern (str "peduncle . (dorsal|medial) lobe.*" reference-brain ".*"))]
82 (loop [sl SAT-lib
83 mb {}]
84 (if (= 2 (count mb))
85 mb
86 (if-let [[k v] (first sl)]
87 (recur (next sl)
88 (if (re-matches r k)
89 (into mb {k v})
90 mb)))))))
92 (def *libs* (ref {}))
94 (defn- load-SAT-lib
95 "Load the named SAT-lib from filepath into *libs* and returns it.
96 Will setup the reference set of fiducial points and the mushroom body lobes from the entry named by reference-brain,
97 as we as register all SATs to the reference brain."
98 [lib-name filepath reference-brain]
99 (let [SAT-lib (load-SAT-lib-as-vs filepath)
100 target-fids (fids-as-Point3d ((SAT-lib reference-brain) :fiducials))
101 SATs (reduce
102 (fn [m e]
103 (conj m (register-SATs (key e) (val e) target-fids)))
104 {}
105 SAT-lib)]
106 (dosync
107 (commute *libs* (fn [libs]
108 (assoc libs lib-name {:filepath filepath
109 :SATs SATs
110 :fids target-fids
111 :mb (load-mb SATs reference-brain)})))))
112 (@*libs* lib-name))
116 (defn- register-vs
117 "Register a singe VectorString3D from source-fids to target-fids."
118 [vs source-fids target-fids]
119 (let [common-fid-keys (intersection (into #{} (keys source-fids)) (into #{} (keys target-fids)))]
120 (if (empty? common-fid-keys)
121 nil
122 (first (Compare/transferVectorStrings [vs]
123 (vals (into (sorted-map) (select-keys source-fids common-fid-keys)))
124 (vals (into (sorted-map) (select-keys target-fids common-fid-keys)))
125 AffineModel3D)))))
127 (defn resample
128 "Returns a resampled copy of vs."
129 [vs delta]
130 (let [copy (.clone vs)]
131 (.resample copy delta)
132 copy))
134 (defn- match-all
135 "Returns a vector of two elements; the first element is the list of matches
136 of the query-vs against all SAT vs in the library sorted by mean euclidean distance,
137 and labeled as correct of incorrect matches according to the Random Forest classifier.
138 The second element is a list of corresponding SAT names for each match.
139 Assumes the query-vs is already registered to the corresponding SAT-lib reference fiducials."
140 [SATs query-vs delta direct substring]
141 (let [vs1 (resample query-vs delta) ; query-vs is already registered into FRT42-fids
142 matches (sort
143 #(int (- (%1 :med) (%2 :med)))
144 (doall (pmap
145 (fn [e]
146 (let [vs2 (let [copy (.clone (val e))]
147 (.resample copy delta)
148 copy)
149 c (Compare/findBestMatch vs1 vs2
150 (double delta) false (int 5) (float 0.5) Compare/AVG_PHYS_DIST
151 direct substring
152 (double 1.1) (double 1.1) (double 1))
153 stats (.getStatistics (get c 0) false (int 0) (float 0) false)]
154 {:SAT-name (key e)
155 :stats stats
156 :med (get stats 0)
157 :correct (LineageClassifier/classify stats)}))
158 SATs)))]
159 ; Cleanup thread table
160 (LineageClassifier/flush)
162 [matches
163 (map (fn [match]
164 (match :SAT-name))
165 matches)]))
167 (defn text-width
168 "Measure the width, in pixels, of the String text, by the Font and FontMetrics of component."
169 [#^String text #^Component component]
170 ; Get the int[] of widths of the first 256 chars
171 (let [#^ints ws (.getWidths (.getFontMetrics component (.getFont component)))]
172 (reduce (fn [sum c]
173 (if (< (int c) (int 256))
174 (+ sum (aget ws (int c)))
175 sum))
176 0
177 text)))
179 (defn- identify-SAT
180 "Takes a calibrated VectorString3D and a list of fiducial points, and checks against the library for identity.
181 For consistency in the usage of the Random Forest classifier, the registration is done into the FRT42D-BP106 brain."
182 [pipe SAT-lib query-vs fids delta direct substring]
183 (let [SATs (SAT-lib :SATs)
184 vs1 (register-vs query-vs fids (SAT-lib :fids))
185 [matches names] (match-all SATs vs1 delta direct substring)
186 SAT-names (vec names)
187 indexed (vec matches)
188 column-names ["SAT" "Match" "Seq sim %" "Lev Dist" "Med Dist" "Avg Dist" "Cum Dist" "Std Dev" "Prop Mut" "Prop Lengths" "Proximity" "Prox Mut" "Tortuosity"]
189 worker (agent nil) ; for event dispatch
190 table (JTable. (proxy [AbstractTableModel] []
191 (getColumnName [col]
192 (get column-names col))
193 (getRowCount []
194 (count matches))
195 (getColumnCount []
196 (count column-names))
197 (getValueAt [row col]
198 (let [match (get indexed row)
199 stats (match :stats)]
200 (cond
201 (= col 0) (get SAT-names row)
202 (= col 1) (str (match :correct)) ; Whether the classifier considered it correct or not
203 true (Utils/cutNumber
204 (cond
205 (= col 2) (* 100 (get stats 6)) ; Similarity
206 (= col 3) (get stats 5) ; Levenshtein
207 (= col 4) (get stats 3) ; Median Physical Distance
208 (= col 5) (get stats 0) ; Average Physical Distance
209 (= col 6) (get stats 1) ; Cummulative Physical Distance
210 (= col 7) (get stats 2) ; Std Dev
211 (= col 8) (get stats 4) ; Prop Mut
212 (= col 9) (get stats 9) ; Prop Lengths
213 (= col 10) (get stats 7) ; Proximity
214 (= col 11) (get stats 8) ; Prox Mut
215 (= col 12) (get stats 10)) ; Tortuosity
216 2))))
217 (isCellEditable [row col]
218 false)
219 (setValueAt [ob row col] nil)))
220 frame (JFrame. (str "Matches for " pipe))
221 dummy-ls (LayerSet. (.. Display getFront getProject) (long -1) "Dummy" (double 0) (double 0) (double 0) (double 0) (double 0) (double 512) (double 512) false (int 0) (java.awt.geom.AffineTransform.))]
222 (.setCellRenderer (.getColumn table "Match")
223 (proxy [DefaultTableCellRenderer] []
224 (getTableCellRendererComponent [t v sel foc row col]
225 (proxy-super setText (str v))
226 (proxy-super setBackground
227 (if (Boolean/parseBoolean v)
228 (Color. 166 255 166)
229 (if sel
230 (Color. 184 207 229)
231 Color/white)))
233 this)))
234 (.add frame (JScrollPane. table))
235 (.setSize frame (int 950) (int 550))
236 (.addMouseListener table
237 (proxy [MouseAdapter] []
238 (mousePressed [ev]
239 (try
240 (clear-agent-errors worker) ; I don't care what it was
241 (catch Exception e))
242 (send-off worker
243 (fn [_]
244 (let [match (indexed (.rowAtPoint table (.getPoint ev)))
245 show-match (fn []
246 (Display3D/addMesh dummy-ls
247 (resample (.clone vs1) delta)
248 "Query"
249 Color/yellow)
250 (Display3D/addMesh dummy-ls
251 (resample (SATs (match :SAT-name)) delta)
252 (match :SAT-name)
253 (if (match :correct)
254 Color/red
255 Color/blue)))]
256 (cond
257 ; On double-click, show 3D view of the match:
258 (= 2 (.getClickCount ev))
259 (show-match)
260 ; On right-click, show menu
261 (Utils/isPopupTrigger ev)
262 (let [popup (JPopupMenu.)
263 new-command (fn [title action]
264 (let [item (JMenuItem. title)]
265 (.addActionListener item (proxy [ActionListener] []
266 (actionPerformed [evt]
267 (send-off worker (fn [_] (action))))))
268 item))]
269 (doto popup
270 (.add (new-command "Show match in 3D"
271 show-match))
272 (.add (new-command "Show Mushroom body"
273 #(doseq [[k v] (SAT-lib :mb)]
274 (Display3D/addMesh dummy-ls v k Color/gray))))
275 (.add (new-command "Show interpolated"
276 #(Display3D/addMesh dummy-ls
277 (VectorString3D/createInterpolatedPoints
278 (Editions. (SATs (match :SAT-name)) (.clone vs1) delta false (double 1.1) (double 1.1) (double 1)) (float 0.5))
279 (str "Interpolated with " (match :SAT-name)) Color/magenta)))
280 (.add (new-command "Show stdDev plot"
281 #(let [cp (ini.trakem2.analysis.Compare$CATAParameters.)]
282 (if (.setup cp false nil true true)
283 (let [condensed (VectorString3D/createInterpolatedPoints
284 (let [v1 (.clone vs1)
285 v2 (SATs (match :SAT-name))]
286 (.resample v1 delta true)
287 (.resample v2 delta true)
288 (Editions. v1 v2 delta false (double 1.1) (double 1.1) (double 1)))
289 (float 0.5))
290 _ (let [c1 (Calibration.); Dummy default calibration but with same units as the pipe.
291 ; VectorString3D instances are already calibrated.
292 c2 (.. pipe getLayerSet getCalibrationCopy)]
293 (.setUnit c1 (.getUnit c2))
294 (.calibrate condensed c1))
295 _ (set! (. cp plot_max_x) (.length condensed))
296 plot (Compare/makePlot cp (str "Query versus " (match :SAT-name)) condensed)]
297 ; There are lazy evaluations ... that need to be promoted: lets print some data:
298 (println "plot is" plot " and condensed has " (.getStdDevAtEachPoint condensed))
299 (if plot
300 (.show plot)
301 (report "ERROR: could not make a plot!")))))))
302 (.show table (.getX ev) (.getY ev)))))))))))
304 ; Enlarge the cell width of the first column
305 (.setMinWidth (.. table getColumnModel (getColumn 0)) (int 250))
306 (doto frame
307 ;(.pack)
308 (.setVisible true))))
311 (def known-libs (ref {"Drosophila-3rd-instar" {:filepath "plugins/SAT-lib-Drosophila-3rd-instar.clj"
312 :reference-brain "FRT42 new"}}))
314 (defn fetch-lib
315 "Returns the SAT lib for the given name, loading it if not there yet. Nil otherwise."
316 [lib-name]
317 (if-let [cached (@*libs* lib-name)]
318 cached
319 (try
320 (if-let [m (@known-libs lib-name)]
321 (load-SAT-lib lib-name
322 (m :filepath)
323 (m :reference-brain))
324 (do
325 (report "Unknown lib " lib-name)
326 nil))
327 (catch Exception e
328 (do
329 (report "An error ocurred while loading the SAT library: " e "\nCheck the terminal output.")
330 (.printStackTrace e))))))
333 (defn extract-fiducial-points
334 "Find a set of fiducial points in the project.
335 Will find the first node of type 'fiducial_points' and recursively inspect
336 its children nodes for nodes of type 'ball'. To each found 'ball', it will
337 assign the name of the ball node itself or the first superior node that
338 has a name different that this type.
339 Returns nil if none found."
340 [project]
341 (if-let [fp-node (.. project getRootProjectThing (findChildrenOfTypeR "fiducial_points"))]
342 (if-let [ball-nodes (.findChildrenOfTypeR (first fp-node) "ball")]
343 (let [fids (reduce
344 (fn [fids pt-ball]
345 (let [title (.. pt-ball getObject getTitle) ; the title of the Ball ob itself, if any
346 b (.. pt-ball getObject getWorldBalls)]
347 (if (= 0 (count b))
348 (do
349 (println "WARNING: empty ball object" (.getObject pt-ball))
350 fids)
351 (do
352 (if (> (count b) 1)
353 (println "WARNING: ball object with more than one ball:" (.getObject pt-ball)))
354 (assoc fids
355 (-> (if (not (= title (.getType pt-ball)))
356 title
357 (.. pt-ball getParent getTitle))
358 (.toLowerCase)
359 (.replace \space \_))
360 (let [b0 (get b 0)]
361 (Point3d. (get b0 0) (get b0 1) (get b0 2))))))))
362 {}
363 ball-nodes)]
364 (println "Found" (count fids) "fiducial points:\n" fids)
365 fids)
366 (println "No ball objects found for" fp-node))
367 (println "No fiducial points found")))
369 (defn identify
370 "Identify a Pipe or Polyline (which implement Line3D) that represent a SAT."
371 ([p lib-name]
372 (identify p lib-name 1.0 true false))
373 ([p lib-name delta direct substring]
374 (if p
375 (if-let [SAT-lib (fetch-lib lib-name)]
376 (if-let [fids (extract-fiducial-points (.getProject p))]
377 (identify-SAT
378 p
379 SAT-lib
380 (let [vs (.asVectorString3D p)]
381 (.calibrate vs (.. p getLayerSet getCalibrationCopy))
382 vs)
383 fids
384 delta
385 direct
386 substring)
387 (report "Cannot find fiducial points for project" (.getProject p)))
388 (report "Cannot find a SAT library for" lib-name))
389 (report "Cannot identify a null pipe or polyline!"))))
392 (defn identify-without-gui
393 "Identify a Pipe or Polyline (which implement Line3D) that represent a SAT.
394 No GUI is shown. Returns the vector containing the list of matches and the list of names."
395 ([p lib-name]
396 (identify-without-gui p lib-name 1.0 true false))
397 ([p lib-name delta direct substring]
398 (if p
399 (if-let [SAT-lib (fetch-lib lib-name)]
400 (let [SATs (SAT-lib :SATs)
401 query-vs (let [vs (.asVectorString3D p)]
402 (.calibrate vs (.. p getLayerSet getCalibrationCopy))
403 vs)
404 fids (extract-fiducial-points (.getProject p)) ; (first (.. p getProject getRootProjectThing (findChildrenOfTypeR "fiducial_points"))))
405 vs1 (register-vs query-vs fids (SAT-lib :fids))]
406 (match-all SATs vs1 delta direct substring)))
407 (report "Cannot identify a null pipe or polyline!"))))
410 ;(defn lib-stats
411 ; "Returns a map with the number of SATs in the library, a list of names vs. sorted sequence lengths, the median length, the average length."
412 ; [lib-name delta]
413 ; (if-let [lib (fetch-lib lib-name)]
414 ; (let [SATs (SAT-lib :SATs)
415 ; lengths (reduce
416 ; (sorted-map-by
417 ; #(<
418 ;
419 ;
420 ; {:n (count SATs)
421 ; :
422 ; (report (str "Unknown library " lib-name))))
424 (defn- ready-vs
425 "Return a calibrate and registered VectorString3D from a chain."
426 [chain source-fids target-fids]
427 (let [vs (.clone (.vs chain))]
428 (.calibrate vs (.. chain getRoot getLayerSet getCalibrationCopy))
429 (register-vs vs source-fids target-fids)))
431 (defn quantify-all
432 "Take all pipes in project and score/classify them.
433 Returns a sorted map of name vs. a vector with:
434 - if the top 1,2,3,4,5 have a homonymous
435 - the number of positives: 'true' and homonymous
436 - the number of false positives: 'true' and not homonymous
437 - the number of false negatives: 'false' and homonymous
438 - the number of true negatives: 'false' and not homonymous
439 - the FPR: false positive rate: false positives / ( false positives + true negatives )
440 - the FNR: false negative rate: false negatives / ( false negatives + true positives )
441 - the TPR: true positive rate: 1 - FNR
442 - the length of the sequence queried"
443 [project lib-name regex-exclude delta direct substring]
444 (let [fids (extract-fiducial-points (first (.. project getRootProjectThing (findChildrenOfTypeR "fiducial_points"))))
445 tops (int 5)
446 SAT-lib (fetch-lib lib-name)]
447 (if SAT-lib
448 (reduce
449 (fn [m chain]
450 (let [vs (resample (ready-vs chain fids (SAT-lib :fids)) delta)
451 [matches names] (match-all (SAT-lib :SATs) vs delta direct substring)
452 names (vec (take tops names))
453 #^String SAT-name (let [t (.getCellTitle chain)]
454 (.substring t (int 0) (.indexOf t (int \space))))
455 ;has-top-match (fn [n] (some #(.startsWith % SAT-name) (take names n)))
456 dummy (println "###\nSAT-name: " SAT-name "\ntop 5 names: " names "\ntop 5 meds: " (map #(% :med) (take 5 matches)))
457 top-matches (loop [i (int 0)
458 r []]
459 (if (>= i tops)
460 r
461 (if (.startsWith (names i) SAT-name)
462 (into r (repeat (- tops i) true)) ; the rest are all true
463 (recur (inc i) (into [false] r)))))
464 true-positives (filter #(and
465 (% :correct)
466 (.startsWith (% :SAT-name) SAT-name))
467 matches)
468 true-negatives (filter #(and
469 (not (% :correct))
470 (not (.startsWith (% :SAT-name) SAT-name)))
471 matches)
472 false-positives (filter #(and
473 (% :correct)
474 (not (.startsWith (% :SAT-name) SAT-name)))
475 matches)
476 false-negatives (filter #(and
477 (not (% :correct))
478 (.startsWith (% :SAT-name) SAT-name))
479 matches)]
480 (println "top matches for " SAT-name " : " (count top-matches) top-matches)
481 (assoc m
482 SAT-name
483 [(top-matches 0)
484 (top-matches 1)
485 (top-matches 2)
486 (top-matches 3)
487 (top-matches 4)
488 (count true-positives)
489 (count false-positives)
490 (count true-negatives)
491 (count false-negatives)
492 (/ (count false-positives) (+ (count false-positives) (count true-negatives))) ; False positive rate
493 (let [divisor (+ (count false-negatives) (count true-positives))]
494 (if (= 0 divisor)
495 Double/MAX_VALUE
496 (/ (count false-negatives) divisor))) ; False negative rate
497 (let [divisor (+ (count false-negatives) (count true-positives))]
498 (if (= 0 divisor)
499 Double/MAX_VALUE
500 (- 1 (/ (count false-negatives) divisor)))) ; True positive rate
501 (.length vs)])))
502 (sorted-map)
503 (Compare/createPipeChains (.getRootProjectThing project) (.getRootLayerSet project) regex-exclude)))))
505 (defn summarize-as-confusion-matrix
506 "Takes the results of quantify-all and returns the confusion matrix as a map."
507 [qa]
508 (reduce
509 (fn [m results]
510 (merge-with + m {:true-positives (results 5)
511 :false-positives (results 6)
512 :true-negatives (results 7)
513 :false-negatives (results 8)}))
514 {:true-positives 0
515 :false-positives 0
516 :true-negatives 0
517 :false-negatives 0}
518 (vals qa)))
520 (defn summarize-as-distributions
521 "Takes the results of quantify-all and returns four vectors representing four histograms,
522 one for each distribution of true-positives, false-positives, true-negatives and false-negatives."
523 [qa]
524 (reduce
525 (fn [m results]
526 (merge-with (fn [m1 m2]
527 (merge-with + m1 m2))
528 m
529 {:true-positives {(results 5) 1}
530 :false-positives {(results 6) 1}
531 :true-negatives {(results 7) 1}
532 :false-negatives {(results 8) 1}}))
533 {:true-positives (sorted-map)
534 :false-positives (sorted-map)
535 :true-negatives (sorted-map)
536 :false-negatives (sorted-map)}
537 (vals qa)))
539 (defn summarize-quantify-all
540 "Takes the results of quantify-all and returns a map of two maps,
541 one with the confusion matrix and one with the distribution of true/false-positive/negative matches."
542 [project lib-name regex-exclude delta direct substring]
543 (let [qa (quantify-all project lib-name regex-exclude delta direct substring)]
544 {:confusion-matrix (summarize-as-confusion-matrix qa)
545 :distributions (summarize-as-distributions)}))
548 (defn print-quantify-all [t]
549 (doseq [[k v] t]
550 (print k \tab)
551 (doseq [x (take 5 v)] (print x \tab))
552 (doseq [x (nthnext v 5)] (print (float x) \tab))
553 (print \newline)))
556 (defn find-rev-match
557 "Take a set of non-homonymous VectorString3D.
558 Do pairwise comparisons, in this fashion:
559 1. reverse vs2
560 2. substring length SL = 0.5 * max(len(vs1), len(vs2))
561 3. for all substrings of len SL of vs1, match against all of vs2.
562 4. Record best result in a set sorted by mean euclidean distance of the best match.
563 Assumes chains are registered and resampled.
564 Will ignore cases where the shorter chain is shorter than half the longer chain."
565 [chains delta]
566 (reduce
567 (fn [s chain1]
568 (let [rvs1 (let [v (.clone (.vs chain1))]
569 (.reverse v)
570 (.resample v delta)
571 v)
572 len1 (.length rvs1)]
573 (conj s (reduce
574 ; Find best reverse submatch for chain1 in chains
575 (fn [best chain2]
576 (println "Processing " (.getShortCellTitle chain1) " versus " (.getShortCellTitle chain2))
577 (let [vs2 (.vs chain2)
578 len2 (.length vs2)
579 SL (int (/ (max len1 len2) 2))]
580 (if (or (= chain1 chain2) (< len1 SL) (< len2 SL))
581 best ; A chain has to be longer than half the length of the other to be worth looking at
582 (let [b (reduce
583 ; Find best reverse submatch between chain1 and chain2
584 (fn [best-rsub start]
585 (let [bm (Compare/findBestMatch (.substring rvs1 start (+ start SL))
586 vs2
587 delta
588 false 0 0
589 Compare/AVG_PHYS_DIST ; Also known as mean euclidean distance
590 true true
591 1.1 1.1 1.0)]
592 (if (and
593 best-rsub
594 (< (best-rsub :med) (get bm 1)))
595 best-rsub
596 {:chain1 chain1 :chain2 chain2 :med (get bm 1)})))
597 {:med Double/MAX_VALUE}
598 (range (- len1 SL)))]
599 (if (and
600 best
601 (< (best :med) (b :med)))
602 best
603 b)))))
604 {:med Double/MAX_VALUE} ; bogus initial element
605 chains))))
606 (sorted-set-by
607 #(int (- (%1 :med) (%2 :med)))
608 {:med Double/MAX_VALUE}) ; adding a bogus initial element
609 chains))
612 (defn find-reverse-match
613 "Takes all the chains of a project and compares them all-to-all pairwise with one of the two reversed.
614 Will calibrate the chains and resample them to delta."
615 [project delta regex-exclude]
616 (let [cal (.. project getRootLayerSet getCalibrationCopy)]
617 (find-rev-match (map
618 (fn [chain]
619 (.calibrate (.vs chain) cal)
620 (.resample (.vs chain) delta)
621 chain)
622 (Compare/createPipeChains (.getRootProjectThing project) (.getRootLayerSet project) regex-exclude))
623 delta)))
625 (defn print-reversed-match
626 "Prints the sorted set returned by find-reversed-similar"
627 [s]
628 (doseq [m s]
629 (let [c1 (m :chain1)
630 c2 (m :chain2)]
631 (if c1
632 (println (.getShortCellTitle c1) (.getShortCellTitle c2) (m :med))))))
635 (defn print-stats
636 "Prints the number of unique SATs and the number of unique lineages in a library of SATs.
637 Returns the two sets in a map."
638 [lib-name]
639 (let [lib (fetch-lib lib-name)
640 unique-sats (reduce
641 (fn [unique clone]
642 (if (> (.indexOf clone (int \space)) -1)
643 (conj unique (.substring clone 0 (.indexOf clone (int \space))))
644 unique))
645 (sorted-set)
646 (keys (lib :SATs)))
647 unique-lineages (reduce
648 (fn [unique clone]
649 (if (> (.indexOf clone (int \:)) -1)
650 (conj unique (.substring clone 0 (.indexOf clone (int \:))))
651 unique))
652 (sorted-set)
653 unique-sats)]
654 (println "Unique SATs: " (count unique-sats) \newline "Unique lineages: " (count unique-lineages))
655 {:unique-sats unique-sats
656 :unique-lineages unique-lineages}))