diff --git a/common/LICENCE b/common/LICENCE
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/common/README.md b/common/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/common/common/L1L2_cells_ids.py b/common/common/L1L2_cells_ids.py
new file mode 100644
index 0000000000000000000000000000000000000000..0d33309d0b03b5c9862ff8b9298a049368562b4c
--- /dev/null
+++ b/common/common/L1L2_cells_ids.py
@@ -0,0 +1,10 @@
+L1_10h = [514, 515, 516, 517, 519, 520, 521, 522, 523, 527, 528, 529, 530, 532, 533, 534, 535, 536, 537, 538, 539, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 554, 555, 558, 122, 126, 129, 149, 158, 173, 183, 196, 197, 205, 211, 215, 232, 240, 247, 249, 265, 269, 274, 278, 280, 281, 285, 292, 293, 295, 307, 310, 322, 323, 332, 337, 342, 345, 352, 355, 362, 367, 372, 384, 385, 387, 391, 396, 401, 402, 405, 408, 409, 411, 417, 424, 430, 435, 439, 440, 447, 448, 451, 455, 456, 460, 461, 464, 465, 466, 470, 473, 474, 476, 479, 484, 485, 487, 489, 491, 494, 499, 501, 502, 504, 505, 506, 509, 510]
+L2_10h = [518, 525, 120, 130, 135, 142, 148, 154, 162, 164, 170, 181, 182, 186, 193, 202, 204, 209, 210, 221, 223, 228, 229, 233, 238, 239, 242, 254, 279, 282, 286, 298, 302, 303, 305, 309, 316, 317, 318, 333, 339, 340, 346, 349, 350, 356, 358, 366, 368, 369, 371, 374, 379, 383, 386, 390, 398, 399, 404, 406, 407, 410, 412, 413, 415, 419, 426, 427, 429, 431, 432, 433, 434, 438, 445, 449, 453, 454, 457, 458, 459, 463, 467, 471, 472, 475, 477, 481, 482, 488, 492, 493, 497, 498, 500, 508, 511]
+L1_40h = [112, 123, 126, 131, 134, 135, 140, 146, 149, 151, 166, 167, 185, 187, 195, 196, 205, 221, 222, 224, 227, 229, 232, 238, 239, 243, 246, 253, 255, 256, 260, 261, 272, 274, 275, 276, 282, 286, 288, 290, 303, 308, 317, 318, 319, 321, 322, 323, 326, 328, 333, 334, 343, 344, 350, 351, 354, 357, 358, 360, 361, 362, 363, 364, 366, 371, 372, 376, 378, 382, 386, 387, 388, 389, 390, 391, 392, 398, 400, 401, 402, 405, 406, 408, 409, 410, 412, 415, 417, 418, 420, 421, 422, 425, 428]
+L2_40h = [257, 258, 259, 179, 266, 267, 269, 270, 237, 277, 278, 153, 218, 281, 154, 156, 285, 240, 291, 292, 293, 294, 369, 168, 297, 170, 299, 300, 301, 302, 306, 307, 310, 311, 313, 249, 193, 327, 200, 330, 331, 335, 336, 210, 340, 341, 342, 215, 345, 346, 207, 349, 223, 352, 356, 103, 242, 171, 359, 365, 110, 368, 113, 370, 235, 245, 169, 121, 217, 252, 125]
+L1_96h = [28, 96, 101, 115, 116, 119, 122, 134, 149, 168, 172, 175, 185, 186, 192, 195, 209, 213, 215, 224, 229, 233, 234, 243, 260, 265, 266, 267, 274, 278, 279, 280, 281, 282, 283, 294, 299, 302, 303, 307, 310, 323, 333, 334, 342, 346, 347, 348, 349, 351, 352, 353, 360, 361, 364, 365, 369, 372, 373, 375, 377, 378, 379, 380, 381, 382, 385, 391, 392, 394, 397, 398, 399, 401, 402, 403, 406, 407, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430]
+L2_96h = [384, 311, 386, 388, 389, 390, 263, 264, 393, 395, 396, 270, 400, 146, 404, 277, 408, 153, 284, 248, 197, 240, 164, 167, 296, 405, 300, 173, 142, 200, 306, 253, 309, 137, 314, 318, 321, 322, 324, 325, 198, 328, 329, 330, 331, 332, 205, 206, 336, 337, 212, 341, 343, 344, 345, 225, 354, 357, 358, 359, 232, 235, 237, 366, 367, 112, 371, 244, 245, 118, 376, 121, 374, 127, 125, 383]
+L1_120h = [515, 699, 518, 769, 522, 526, 921, 531, 533, 534, 857, 536, 543, 544, 545, 546, 554, 555, 560, 561, 862, 566, 567, 569, 576, 577, 586, 781, 598, 599, 601, 783, 605, 607, 609, 99, 786, 112, 625, 627, 628, 633, 634, 637, 639, 640, 129, 875, 132, 135, 649, 139, 140, 657, 658, 659, 662, 663, 664, 155, 704, 672, 624, 676, 679, 683, 684, 685, 688, 177, 179, 692, 693, 695, 891, 697, 187, 188, 189, 192, 709, 710, 711, 201, 887, 915, 740, 212, 213, 726, 728, 890, 222, 739, 228, 742, 743, 233, 234, 747, 237, 751, 752, 753, 755, 756, 757, 246, 760, 766, 257, 261, 776, 265, 778, 779, 268, 269, 271, 785, 274, 787, 788, 279, 745, 644, 286, 901, 816, 845, 295, 808, 297, 811, 777, 818, 814, 815, 304, 817, 306, 819, 308, 311, 820, 827, 830, 319, 320, 834, 837, 838, 327, 328, 908, 330, 823, 333, 847, 848, 850, 853, 825, 856, 345, 911, 860, 350, 863, 352, 865, 866, 867, 868, 357, 358, 913, 872, 363, 876, 367, 880, 369, 882, 883, 372, 885, 375, 888, 378, 831, 892, 893, 894, 895, 896, 897, 386, 899, 900, 773, 902, 903, 904, 905, 906, 907, 396, 909, 398, 399, 912, 401, 914, 403, 916, 917, 918, 919, 920, 409, 922, 411, 412, 413, 926, 416, 798, 923, 924, 429, 430, 925, 432, 884, 448, 449, 454, 458, 461, 462, 463, 466, 473, 477, 910, 482, 749, 484, 489, 499, 502, 869, 506, 509, 510, 511]
+L2_120h = [512, 516, 520, 523, 525, 528, 529, 542, 347, 548, 551, 553, 559, 562, 564, 570, 575, 578, 68, 582, 585, 588, 589, 79, 855, 597, 804, 600, 186, 606, 608, 97, 610, 611, 613, 615, 616, 619, 620, 621, 110, 623, 626, 629, 638, 127, 642, 643, 645, 648, 137, 656, 154, 156, 157, 158, 671, 673, 678, 167, 680, 682, 689, 691, 457, 698, 701, 702, 705, 706, 197, 289, 712, 718, 208, 723, 724, 214, 727, 729, 730, 731, 733, 735, 225, 738, 229, 744, 748, 238, 242, 244, 758, 761, 762, 252, 765, 768, 770, 771, 772, 774, 264, 780, 782, 272, 789, 278, 791, 280, 281, 794, 284, 861, 801, 802, 803, 292, 805, 806, 807, 296, 813, 302, 305, 821, 310, 824, 313, 314, 315, 829, 832, 323, 324, 822, 839, 840, 841, 842, 843, 332, 846, 852, 854, 343, 344, 346, 859, 826, 864, 870, 873, 874, 877, 878, 879, 881, 886, 376, 889, 382, 391, 392, 397, 835, 750, 406, 407, 408, 836, 414, 420, 421, 422, 427, 428, 754, 946, 670, 441, 444, 451, 844, 468, 471, 472, 476, 480, 481, 485, 486, 764, 490, 500, 501, 505, 507, 508, 784]
+L1_132h = [71, 73, 89, 90, 92, 108, 109, 120, 121, 123, 127, 142, 158, 162, 164, 176, 177, 181, 182, 186, 197, 206, 208, 209, 215, 221, 222, 224, 234, 238, 248, 250, 252, 255, 258, 265, 267, 276, 277, 289, 291, 293, 300, 307, 308, 309, 315, 319, 322, 323, 324, 325, 328, 338, 352, 353, 356, 377, 378, 382, 383, 385, 387, 389, 394, 395, 396, 397, 400, 402, 403, 406, 407, 414, 420, 426, 427, 429, 436, 444, 446, 447, 454, 459, 466, 467, 473, 475, 480, 482, 483, 484, 492, 493, 494, 497, 498, 503, 507, 508, 513, 516, 517, 518, 520, 521, 527, 529, 530, 532, 536, 539, 541, 542, 543, 544, 545, 549, 553, 555, 557, 558, 559, 560, 562, 563, 577, 580, 582, 585, 587, 590, 591, 592, 601, 606, 607, 615, 616, 618, 619, 626, 628, 630, 635, 639, 640, 641, 643, 644, 646, 651, 661, 666, 677, 678, 682, 683, 688, 689, 697, 703, 707, 710, 711, 715, 716, 721, 723, 727, 730, 736, 739, 742, 744, 747, 748, 754, 758, 759, 761, 770, 771, 775, 788, 792, 793, 799, 800, 813, 821, 825, 827, 828, 829, 832, 835, 841, 842, 848, 849, 852, 853, 858, 859, 861, 868, 869, 879, 882, 886, 887, 890, 896, 899, 901, 905, 910, 911, 913, 921, 922, 923, 927, 936, 939, 943, 946, 947, 950, 957, 958, 960, 963, 968, 969, 970, 977, 978, 980, 985, 986, 996, 1000, 1004, 1008, 1013, 1014, 1016, 1019, 1024, 1029, 1031, 1033, 1039, 1040, 1045, 1047, 1048, 1049, 1052, 1053, 1067, 1071, 1079, 1086, 1087, 1088, 1092, 1093, 1094, 1097, 1100, 1102, 1104, 1105, 1106, 1117, 1118, 1123, 1126, 1127, 1128, 1129, 1132, 1134, 1139, 1142, 1143, 1144, 1145, 1151, 1155, 1163, 1164, 1166, 1169, 1177, 1179, 1183, 1184, 1191, 1192, 1194, 1195, 1196, 1200, 1201, 1203, 1205, 1208, 1210, 1211, 1215, 1219, 1221, 1223, 1224, 1228, 1230, 1231, 1234, 1235, 1237, 1239, 1242, 1247, 1254, 1255, 1258, 1259, 1260, 1263, 1266, 1268, 1269, 1270, 1272, 1273, 1274, 1275, 1276, 1279, 1281, 1282, 1283, 1284, 1285, 1286, 1288, 1292, 1294, 1300, 1301, 1302, 1303, 1304, 1306, 1308, 1309, 1310, 1311, 1313, 1314, 1315, 1316, 1317, 1318, 1319, 1322, 1323, 1324, 1325, 1326, 1327, 1328, 1329, 1330, 1331, 1332, 1333, 1334, 1335, 1336, 1337, 1338, 1340, 1344, 1345, 1346, 1348, 1350, 1351, 1354, 1355, 1356, 1357, 1360, 1363, 1369, 1370, 1371, 1374, 1375, 1376, 1377, 1379, 1382, 1383, 1384, 1386, 1387, 1389, 1392, 1393, 1394, 1395, 1396, 1400, 1401, 1402, 1403, 1404, 1405, 1406, 1407, 1408, 1409, 1410, 1411, 1412, 1413, 1414, 1417, 1418, 1452]
+L2_132h = [512, 1025, 392, 1028, 598, 519, 1034, 1035, 524, 525, 526, 1197, 1041, 988, 534, 535, 1050, 602, 1054, 1055, 1058, 1059, 1060, 380, 550, 1063, 1065, 1066, 556, 1214, 1070, 1073, 1074, 53, 566, 1081, 1082, 575, 576, 1089, 1091, 1366, 1121, 72, 586, 1207, 1378, 79, 81, 594, 1107, 596, 1110, 612, 1114, 1115, 1116, 361, 1120, 97, 98, 1124, 614, 1297, 617, 1130, 1298, 110, 111, 624, 113, 1140, 629, 105, 122, 1147, 1148, 1149, 126, 1154, 363, 1156, 1157, 1158, 705, 1160, 1388, 935, 489, 141, 654, 1168, 146, 660, 1173, 663, 152, 1220, 1178, 667, 156, 1182, 671, 1185, 674, 1187, 1188, 1189, 1193, 1341, 685, 1198, 1199, 1165, 1399, 1202, 1204, 1206, 183, 184, 185, 1213, 702, 373, 192, 193, 1218, 195, 708, 228, 1222, 502, 200, 1225, 1226, 203, 1229, 720, 211, 1236, 213, 1240, 1241, 1244, 1245, 1246, 1248, 1232, 1250, 1252, 1253, 232, 233, 1261, 1262, 1264, 1265, 242, 243, 638, 1271, 760, 249, 762, 251, 1277, 1278, 725, 1280, 959, 773, 774, 1287, 1289, 778, 779, 1293, 785, 274, 1299, 1103, 790, 791, 282, 1307, 796, 797, 798, 437, 801, 290, 390, 487, 1320, 297, 298, 1135, 304, 305, 820, 734, 312, 826, 317, 830, 1343, 198, 321, 1347, 326, 1249, 1352, 329, 330, 844, 846, 653, 1361, 335, 339, 1364, 1365, 854, 569, 1368, 740, 1372, 349, 741, 367, 866, 867, 256, 1381, 871, 360, 1385, 362, 875, 364, 877, 1390, 1391, 1342, 372, 1397, 1398, 745, 888, 933, 374, 1257, 892, 381, 384, 878, 902, 833, 904, 548, 909, 237, 912, 920, 245, 916, 348, 408, 409, 855, 412, 413, 926, 419, 345, 241, 623, 327, 428, 430, 925, 1096, 755, 949, 951, 440, 500, 955, 445, 816, 1321, 964, 965, 753, 457, 971, 1380, 462, 976, 931, 1358, 984, 676, 987, 476, 478, 479, 1170, 763, 485, 998, 999, 1001, 490, 1003, 1367, 1005, 1009, 1011, 1217, 766, 1227, 681, 668, 1176, 1020, 1021, 945]
diff --git a/common/common/__init__.py b/common/common/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/common/common/__pycache__/L1L2_cells_ids.cpython-37.pyc b/common/common/__pycache__/L1L2_cells_ids.cpython-37.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..b625b0d491786ec38d5b49865e9914b726b6cf61
Binary files /dev/null and b/common/common/__pycache__/L1L2_cells_ids.cpython-37.pyc differ
diff --git a/common/common/__pycache__/__init__.cpython-37.pyc b/common/common/__pycache__/__init__.cpython-37.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..19fec233ee28f7a3f8a9731001549ede26956b8a
Binary files /dev/null and b/common/common/__pycache__/__init__.cpython-37.pyc differ
diff --git a/common/common/__pycache__/lin.cpython-37.pyc b/common/common/__pycache__/lin.cpython-37.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..eac3d3c486ca32d2f01e700fa68d5f251a9f26f8
Binary files /dev/null and b/common/common/__pycache__/lin.cpython-37.pyc differ
diff --git a/common/common/__pycache__/seg.cpython-37.pyc b/common/common/__pycache__/seg.cpython-37.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..51c47c8b7428722ffce143168c611ed13dbca604
Binary files /dev/null and b/common/common/__pycache__/seg.cpython-37.pyc differ
diff --git a/common/common/edict.py b/common/common/edict.py
new file mode 100644
index 0000000000000000000000000000000000000000..8df594153ea1b32e24bda9ddc3ea77d5cebcdda5
--- /dev/null
+++ b/common/common/edict.py
@@ -0,0 +1,58 @@
+def unionWith(m1, m2, f):
+ combined = dict()
+
+ combKeys = set.union(set(m1.keys()), set(m2.keys()))
+
+ for k in combKeys:
+ if k in m1 and k in m2:
+ combined[k] = f(m1[k], m2[k])
+ elif k in m1 and k not in m2:
+ combined[k] = m1[k]
+ elif k not in m1 and k in m2:
+ combined[k] = m2[k]
+
+ return combined
+
+def unionsWith(ms, f, df):
+ combined = dict()
+
+ combKeys = set.union(*[set(m.keys()) for m in ms])
+
+ for k in combKeys:
+ combined[k] = f([m.get(k, df) for m in ms])
+
+ return combined
+
+def union(m1, m2):
+ combined = dict()
+
+ if set(m1.keys()).intersection(m2.keys()):
+ raise ValueError("non distinct keys")
+
+ combKeys = set.union(set(m1.keys()), set(m2.keys()))
+
+ for k in combKeys:
+ if k in m1:
+ combined[k] = m1[k]
+ elif k in m2:
+ combined[k] = m2[k]
+
+ return combined
+
+
+def map(m, f):
+ return dict([(k, f(v)) for k, v in m.items()])
+
+def mapWithKey(m, f):
+ return dict([(k, f(k, v)) for k, v in m.items()])
+
+def relabel(m, kms):
+ m1 = dict()
+ for k, v in m:
+ for km in kms.get(k, []):
+ m1[km] = v
+
+ return m1
+
+def filter(m, f):
+ return dict([(k, v) for k, v in m.items() if f(v)])
diff --git a/common/common/lin.py b/common/common/lin.py
new file mode 100644
index 0000000000000000000000000000000000000000..b607cca5875c180d994ea8af8847ecc8d36082d0
--- /dev/null
+++ b/common/common/lin.py
@@ -0,0 +1,498 @@
+from collections import defaultdict
+from os.path import join
+import os
+import _pickle as cPickle
+import common.seg as seg
+from itertools import repeat
+import common.L1L2_cells_ids as layers
+
+def getSurvivingBetween(tss, linss, t1, t2):
+ linsL = filterLinssBetween(linss, (t1, t2))
+
+ cids = list()
+
+ for c in tss[t1]:
+ if tSearchN(linsL, c.cid, len(linsL), 0):
+ cids.append(c.cid)
+
+ return cids
+
+def filterSurvivingUntil(tss, linss, t2):
+ for t, ts in tss.items():
+ cids = getSurvivingBetween(tss, linss, t, t2)
+ filterT(ts, cids)
+
+ return
+
+def fst(x): return x[0]
+
+def snd(x): return x[1]
+
+def filterT(ts, cids):
+ scids = set(cids)
+ ts.cells = dict([(c.cid, c) for c in ts if c.cid in scids])
+
+ neighsL1 = dict()
+ for cid, ds in ts.neighs.items():
+ if cid in scids:
+ neighsL1[cid] = dict([(n, d) for n, d in ds.items()
+ if n in scids])
+
+ ts.neighs = neighsL1
+
+def filterL1(tss):
+ L1_cids = {10: layers.L1_10h,
+ 40: layers.L1_40h,
+ 96: layers.L1_96h,
+ 120: layers.L1_120h,
+ 132: layers.L1_132h}
+
+ for t, ts in tss.items():
+ filterT(ts, L1_cids[t])
+
+ return
+
+
+def filterL1L2(tss):
+ L1_cids = {10: layers.L1_10h,
+ 40: layers.L1_40h,
+ 96: layers.L1_96h,
+ 120: layers.L1_120h,
+ 132: layers.L1_132h}
+
+ L2_cids = {10: layers.L2_10h,
+ 40: layers.L2_40h,
+ 96: layers.L2_96h,
+ 120: layers.L2_120h,
+ 132: layers.L2_132h}
+
+ for t, ts in tss.items():
+ filterT(ts, L1_cids[t] + L2_cids[t])
+
+ return
+
+
+def getTs(linFn):
+ fn, _ = os.path.splitext(linFn)
+ (_, t1, _, t2) = fn.split("_")
+
+ return (int(t1[:-1]), int(t2[:-1]))
+
+
+def readLin(linFn):
+ with open(linFn, 'rb') as linF:
+ lin = cPickle.load(linF, encoding='latin1')
+
+ return lin, getTs(linFn)
+
+
+def readLins(linDataLoc):
+ lins = [readLin(join(linDataLoc, f)) for f in os.listdir(linDataLoc)
+ if os.path.isfile(join(linDataLoc, f))]
+
+ return lins
+
+
+def getLinMap(lin):
+ cMap = defaultdict(list)
+ for c1, c2 in lin[0]:
+ cMap[c2].append(c1)
+
+ return cMap
+
+
+def mkLinsSeries(lins):
+ return [(ts, getLinMap(lin))
+ for lin, ts in sorted(lins, key=lambda lt: lt[1][0])]
+
+
+def after(t, lt):
+ ts, _ = lt
+ (t1, t2) = ts
+
+ return t >= t1
+
+
+def before(t, lt):
+ ts, _ = lt
+ (t1, t2) = ts
+
+ return t <= t2
+
+
+def between(linSeries, t1, t2):
+ afters = [i for i, lt in enumerate(linSeries) if after(t1, lt)]
+ befores = [i for i, lt in enumerate(linSeries) if before(t2, lt)]
+
+ s = afters[-1]
+ e = befores[0]
+
+ return [(ts, lin) for ts, lin in linSeries[s:(e+1)]]
+
+
+def mergeLins(lin1, lin2):
+ mLin = dict()
+
+ for m, ds in lin1.items():
+ mLin[m] = sum([lin2[d] for d in ds], [])
+
+ return mLin
+
+
+def mergeLinss1(lins, acc):
+ if not lins: return acc
+
+ return mergeLinss1(lins[1:], mergeLins(acc, lins[0]))
+
+
+def mergeLinss(lins):
+ if len(lins) > 2: return mergeLinss1(lins[2:], mergeLins(lins[0], lins[1]))
+ elif len(lins) == 2: return mergeLins(lins[0], lins[1])
+ elif len(lins) == 1: return lins[0]
+ else: return []
+
+
+######################################################
+def concat(xss):
+ return sum(xss, [])
+
+
+def tSearchN(mss, r, n, c):
+ if c == n: return [r]
+ else:
+ return concat([tSearchN(mss, d, n, c+1) for d in mss[c].get(r, [])])
+
+
+def tSearch1(mss, r):
+ try:
+ return mss[0][r]
+ except KeyError:
+ return []
+
+
+def tSearch2(mss, r):
+ try:
+ return sum([mss[1][d] for d in mss[0][r]], [])
+ except KeyError:
+ print(r)
+ return []
+
+
+def cvol(c):
+ if c: return c.vol
+ else: return 0.0
+
+
+def sumVol(cs, f=cvol):
+ return sum([cvol(c) for c in cs])
+
+
+def mergeCellVolsRefN(tss, mss):
+ volss = []
+
+ for i, ts in enumerate(tss):
+ volss.append(dict())
+ for c in tss[0]:
+ volss[i][c.cid] = sum(cvol(ts.cells.get(ci, None))
+ for ci in tSearchN(mss, c.cid, i, 0))
+
+ return volss
+
+
+def const1(c):
+ if c: return 1
+ else: return 0
+
+
+def mergeCellNsRefN(tss, linss):
+ attrss = []
+
+ for i, ts in enumerate(tss):
+ attrss.append(dict())
+ for c in tss[0]:
+ attrss[i][c.cid] = sum(const1(ts.cells.get(ci, None))
+ for ci in tSearchN(linss, c.cid, i, 0))
+
+ return attrss
+
+
+def mergeCellVolsRef3(tss, mss):
+ (ts1, ts2, ts3) = tss
+
+ vols1 = dict()
+ vols2 = dict()
+ vols3 = dict()
+
+ for c in ts1:
+ vols1[c.cid] = c.vol
+ vols2[c.cid] = sumVol([ts2.cells[i] for i in tSearch1(mss, c.cid)])
+ vols3[c.cid] = sumVol([ts3.cells[i] for i in tSearch2(mss, c.cid)])
+
+ return vols1, vols2, vols3
+
+
+def mergeCellVolsRef2(tss, mss):
+ (ts1, ts2) = tss
+
+ vols1 = dict()
+ vols2 = dict()
+
+ for c in ts1:
+ print(c.cid)
+ vols1[c.cid] = c.vol
+ vols2[c.cid] = sumVol([ts2.cells[i] for i in tSearch1(mss, c.cid)])
+
+ return vols1, vols2
+
+
+def divOr0(n1, n2):
+ if n2 == 0: return 0
+ else: return n1 / n2
+
+
+def growthRates3(volss, m12):
+ vols1, vols2, vols3 = volss
+ gr2 = dict()
+ for i, vol in vols1.items():
+ gr12 = divOr0((vols2[i] - vols1[i]), vols1[i])
+ gr23 = divOr0((vols3[i] - vols2[i]), vols2[i])
+ gr2[i] = 0.5 * (gr12 + gr23)
+
+ gr22 = [zip(m12.get(i, []), repeat(gr, len(m12.get(i, []))))
+ for i, gr in gr2.items()]
+
+ return dict(sum(gr22, []))
+
+
+def growthRatesL2(volss):
+ vols1, vols2 = volss
+ gr1 = dict()
+ for i, vol in vols1.items():
+ gr1[i] = (vols2[i] - vols1[i]) / vols1[i]
+
+ return gr1
+
+
+def growthRatesR2(volss, m12):
+ vols1, vols2 = volss
+ gr2 = dict()
+ for i, vol in vols1.items():
+ gr2[i] = (vols2[i] - vols1[i]) / vols1[i]
+
+ gr22 = [zip(m12.get(i, []), repeat(gr, len(m12.get(i, []))))
+ for i, gr in gr2.items()]
+
+ return dict(sum(gr22, []))
+
+
+def mkSeries1(dataDir):
+ d = "organism"
+ dExprs = "geneExpression/new_patterns"
+ linDataLoc = "tracking"
+
+
+ timepoints = [10, 40, 96, 120, 132]
+ lins = mkLinsSeries(readLins(join(dataDir, linDataLoc)))
+
+ tss = dict()
+ linss = dict()
+
+ for t1, t2 in zip(timepoints, timepoints[1:]):
+ tss[t1] = readDataT1(join(dataDir, d), join(dataDir, dExprs), t1)
+ tss[t2] = readDataT1(join(dataDir, d), join(dataDir, dExprs), t2)
+
+ linss[(t1, t2)] = mergeLinss(list(map(lambda x: x[1],
+ between(lins, t1, t2))))
+
+ return tss, linss
+
+def mkSeriesIm0(dataDir):
+ d = "yr01"
+ linDataLoc = "yr01/tracking"
+
+ timepoints = [0, 10, 40, 96, 120, 132]
+ lins = mkLinsSeries(readLins(join(dataDir, linDataLoc)))
+
+ tss = dict()
+ linss = dict()
+
+ for t in timepoints:
+ tss[t] = readDataT1Im(join(dataDir, d), t)
+
+ for t1, t2 in zip(timepoints, timepoints[1:]):
+ linss[(t1, t2)] = mergeLinss(map(lambda x: x[1],
+ between(lins, t1, t2)))
+
+ return tss, linss
+
+
+def mkGrowthRates():
+ tss, linss = mkSeries()
+ grs = dict()
+
+ grs[10] = growthRatesL2(mergeCellVolsRef2((tss[10], tss[40]),
+ [linss[(10, 40)]]))
+ grs[40] = growthRates3(mergeCellVolsRef3((tss[10], tss[40], tss[96]),
+ [linss[(10, 40)], linss[(40, 96)]]),
+ linss[(10, 40)])
+ grs[96] = growthRates3(mergeCellVolsRef3((tss[40], tss[96], tss[120]),
+ [linss[(40, 96)], linss[(96, 120)]]),
+ linss[(40, 96)])
+ grs[120] = growthRates3(mergeCellVolsRef3((tss[96], tss[120], tss[132]),
+ [linss[(96, 120)], linss[(120, 132)]]),
+ linss[(96, 120)])
+ grs[132] = growthRatesR2(mergeCellVolsRef2((tss[120], tss[132]),
+ [linss[(120, 132)]]),
+ linss[(120, 132)])
+
+ return grs
+
+def readGeomImT(dataDir, fid, t):
+ ressFn = "yr" + fid + "/resolutions.txt"
+ segFn = ("yr" +
+ fid +
+ "/segmentations/YR_" +
+ fid + "_" + str(t) +
+ "h_segmented.tif")
+
+ ress = readRess(join(dataDir, ressFn))
+ img = readImage(join(dataDir, segFn))
+ res = ress[t]
+
+ ts = Tissue.fromImage(data=img, res=res)
+ for c in ts:
+ c.swapZX()
+
+ return ts
+
+def mkSeriesGeomIm(dataDir, fid):
+ ressFn = "yr" + fid + "/resolutions.txt"
+ linDataLoc = "yr" + fid + "/trackingFiles/"
+ ress = readRess(join(dataDir, ressFn))
+ tpoints = sorted(ress.keys())
+
+ lins = mkLinsSeries(readLins(join(dataDir, linDataLoc)))
+ tss = dict()
+ linss = dict()
+
+ for t in tpoints:
+ tss[t] = readGeomImT(dataDir, fid, t)
+
+ for t1, t2 in zip(tpoints, tpoints[1:]):
+ linss[(t1, t2)] = mergeLinss(map(lambda x:
+ x[1], between(lins, t1, t2)))
+
+ return tss, linss
+
+
+###################################
+def invertMap(ln):
+ iLin = [list(zip(ds, repeat(m, len(ds)))) for m, ds in ln.items()]
+
+ return dict(sum(iLin, []))
+
+def updateCells(ts, ts1, iLn):
+ for c1 in ts1:
+ try:
+ c = ts.cells[iLn[c1.cid]]
+ except KeyError:
+ print(c1.cid)
+ c = Cell(0, Vec3(0, 0, 0), 0, dict([(g, False) for g in geneNms]))
+ c1.exprs = dict([(g, c.exprs[g]) for g, v in c1.exprs.items()])
+
+ return
+
+def prev(t):
+ if t == 40: return 10
+ elif t == 96: return 40
+ elif t == 120: return 96
+ elif t == 132: return 120
+ else: return None
+
+def succ(t):
+ if t==10: return 40
+ elif t==40: return 96
+ elif t==96: return 120
+ elif t==120: return 132
+ else: return None
+
+def updateTs(tss, linss, tpoints):
+ #evolve patterns based on lineage
+ tssLin = {}
+ for t in tpoints[1:]:
+ tssLin[(prev(t), t)] = tss[t].updateExprsTs(tss[prev(t)],
+ invertMap(linss[(prev(t), t)]))
+
+ return tssLin
+
+def updateTsBF(tss, fs):
+ return dict((t, tss[t].updateExprsBF(fs[t])) for t in tss.keys())
+
+def toNewmanFn(fn):
+ (fn, ext) = os.path.splitext(fn)
+ return fn + "000" + ".tif"
+
+
+def tFn(t, ref=""):
+ return "ts_t" + str(t) + ref + ".txt"
+
+
+def snd(x):
+ return x[1]
+
+def betweenTs(t, bounds):
+ ts, te = bounds
+ return t >= ts and t <= te
+
+def filterLinssBetween(linss, bounds):
+ t1, t2 = bounds
+ linssL = sorted([(ts, ln) for ts, ln in linss.items()
+ if betweenTs(ts[0], (t1, t2)) and
+ betweenTs(ts[1], (t1, t2))], key=lambda p: p[0][0])
+
+ return list(map(snd, linssL))
+
+def filterTssBetween(tss, tbounds):
+ t1, t2 = tbounds
+
+ tssL = sorted([(t, ts) for t, ts in tss.items() if betweenTs(t, (t1, t2))], key=fst)
+
+ return list(map(snd, tssL))
+
+def readDataT1(d, dExprs, t):
+ from os.path import join
+
+ geomF = "_segmented_tvformat_0_volume_position.txt"
+ neighF = "_segmented_tvformat_0_neighbors_py.txt"
+ exprF = "h.txt"
+
+ fnGeom = join(d, "t" + str(t) + geomF)
+ fnNeigh = join(d, "t" + str(t) + neighF)
+ fnExpr = join(dExprs, "t_" + str(t) + exprF)
+
+ gexprs = seg.readExprs(fnExpr)
+ ts = seg.STissue.fromTV(fnGeom, fnNeigh, gexprs, seg.geneNms)
+
+ return ts
+
+def readDataT1Im(d, dExprs, t):
+ from os.path import join
+
+ exprF = "h.txt"
+ fnRes = join(d, "resolutions.txt")
+ fnSeg = join(d, "segmentations/t"+str(t)+"_segmented.tif")
+ fnExpr = join(dExprs, "t_" + str(t) + exprF)
+
+
+ data, _ = seg.readImage(fnSeg)
+ ress = seg.readRess(fnRes)
+ res = ress[t]
+
+ gexprs = seg.readExprs(fnExpr)
+ ts = seg.STissue.fromImage(data, res, gexprs, seg.geneNms)
+
+ for c in ts:
+ c.swapZX()
+
+ return ts
+
diff --git a/common/common/lin.py~ b/common/common/lin.py~
new file mode 100644
index 0000000000000000000000000000000000000000..3cd910b4b85dba6e884fc51bbdd0c67b86643f41
--- /dev/null
+++ b/common/common/lin.py~
@@ -0,0 +1,546 @@
+from __future__ import print_function
+from __future__ import division
+from collections import defaultdict
+from os.path import join
+import os
+import _pickle as cPickle
+from seg import *
+import analyseExprs as a
+from itertools import repeat
+import boolFs as b
+import L1L2_cells_ids as layers
+
+linDataLoc = "../data/tracking/"
+
+def getSurvivingBetween(tss, linss, t1, t2):
+ linsL = filterLinssBetween(linss, (t1, t2))
+
+ cids = list()
+
+
+ for c in tss[t1]:
+ if tSearchN(linsL, c.cid, len(linsL), 0):
+ cids.append(c.cid)
+
+ return cids
+
+def filterSurvivingUntil(tss, linss, t2):
+ for t, ts in tss.items():
+ cids = getSurvivingBetween(tss, linss, t, t2)
+ filterT(ts, cids)
+
+ return
+
+def fst(x): return x[0]
+
+def snd(x): return x[1]
+
+def filterT(ts, cids):
+ scids = set(cids)
+ ts.cells = dict([(c.cid, c) for c in ts if c.cid in scids])
+
+ neighsL1 = dict()
+ for cid, ds in ts.neighs.items():
+ if cid in scids:
+ neighsL1[cid] = dict([(n, d) for n, d in ds.items() if n in scids])
+
+ ts.neighs = neighsL1
+
+def filterL1(tss):
+ L1_cids = {10: layers.L1_10h,
+ 40: layers.L1_40h,
+ 96: layers.L1_96h,
+ 120: layers.L1_120h,
+ 132: layers.L1_132h}
+
+ for t, ts in tss.items():
+ filterT(ts, L1_cids[t])
+
+ return
+
+
+def filterL1L2(tss):
+ L1_cids = {10: layers.L1_10h,
+ 40: layers.L1_40h,
+ 96: layers.L1_96h,
+ 120: layers.L1_120h,
+ 132: layers.L1_132h}
+
+ L2_cids = {10: layers.L2_10h,
+ 40: layers.L2_40h,
+ 96: layers.L2_96h,
+ 120: layers.L2_120h,
+ 132: layers.L2_132h}
+
+ for t, ts in tss.items():
+ #scids1 = set(L1_cids[t])
+ #scids2 = set(L2_cids[t])
+
+ filterT(ts, L1_cids[t] + L2_cids[t])
+
+# ts.cells = dict([(c.cid, c) for c in ts
+# if c.cid in scids1 or c.cid in scids2])
+
+ return
+
+
+def getTs(linFn):
+ fn, _ = os.path.splitext(linFn)
+ (_, t1, _, t2) = fn.split("_")
+
+ return (int(t1[:-1]), int(t2[:-1]))
+
+
+def readLin(linFn):
+ with open(linFn, 'rb') as linF:
+ lin = cPickle.load(linF, encoding='latin1')
+
+ return lin, getTs(linFn)
+
+
+def readLins(linDataLoc):
+ lins = [readLin(join(linDataLoc, f)) for f in os.listdir(linDataLoc)
+ if os.path.isfile(join(linDataLoc, f))]
+
+ return lins
+
+
+def getLinMap(lin):
+ cMap = defaultdict(list)
+ for c1, c2 in lin[0]:
+ cMap[c2].append(c1)
+
+ return cMap
+
+
+def mkLinsSeries(lins):
+ return [(ts, getLinMap(lin))
+ for lin, ts in sorted(lins, key=lambda lt: lt[1][0])]
+
+
+def after(t, lt):
+ ts, _ = lt
+ (t1, t2) = ts
+
+ return t >= t1
+
+
+def before(t, lt):
+ ts, _ = lt
+ (t1, t2) = ts
+
+ return t <= t2
+
+
+def between(linSeries, t1, t2):
+ afters = [i for i, lt in enumerate(linSeries) if after(t1, lt)]
+ befores = [i for i, lt in enumerate(linSeries) if before(t2, lt)]
+
+ s = afters[-1]
+ e = befores[0]
+
+ return [(ts, lin) for ts, lin in linSeries[s:(e+1)]]
+
+
+def mergeLins(lin1, lin2):
+ mLin = dict()
+
+ for m, ds in lin1.items():
+ mLin[m] = sum([lin2[d] for d in ds], [])
+
+ return mLin
+
+
+def mergeLinss1(lins, acc):
+ if not lins: return acc
+
+ return mergeLinss1(lins[1:], mergeLins(acc, lins[0]))
+
+
+def mergeLinss(lins):
+ if len(lins) > 2: return mergeLinss1(lins[2:], mergeLins(lins[0], lins[1]))
+ elif len(lins) == 2: return mergeLins(lins[0], lins[1])
+ elif len(lins) == 1: return lins[0]
+ else: return []
+
+
+######################################################
+def concat(xss):
+ return sum(xss, [])
+
+
+def tSearchN(mss, r, n, c):
+ if c == n: return [r]
+ else:
+ return concat([tSearchN(mss, d, n, c+1) for d in mss[c].get(r, [])])
+
+
+def tSearch1(mss, r):
+ try:
+ return mss[0][r]
+ except KeyError:
+ return []
+
+
+def tSearch2(mss, r):
+ try:
+ return sum([mss[1][d] for d in mss[0][r]], [])
+ except KeyError:
+ print(r)
+ return []
+
+
+def cvol(c):
+ if c: return c.vol
+ else: return 0.0
+
+
+def sumVol(cs, f=cvol):
+ return sum([cvol(c) for c in cs])
+
+
+def mergeCellVolsRefN(tss, mss):
+ volss = []
+
+ for i, ts in enumerate(tss):
+ volss.append(dict())
+ for c in tss[0]:
+ volss[i][c.cid] = sum(cvol(ts.cells.get(ci, None)) for ci in tSearchN(mss, c.cid, i, 0))
+
+ return volss
+
+
+def const1(c):
+ if c: return 1
+ else: return 0
+
+
+def mergeCellNsRefN(tss, linss):
+ attrss = []
+
+ for i, ts in enumerate(tss):
+ attrss.append(dict())
+ for c in tss[0]:
+ attrss[i][c.cid] = sum(const1(ts.cells.get(ci, None)) for ci in tSearchN(linss, c.cid, i, 0))
+
+ return attrss
+
+
+def mergeCellVolsRef3(tss, mss):
+ (ts1, ts2, ts3) = tss
+
+ vols1 = dict()
+ vols2 = dict()
+ vols3 = dict()
+
+ for c in ts1:
+ vols1[c.cid] = c.vol
+ vols2[c.cid] = sumVol([ts2.cells[i] for i in tSearch1(mss, c.cid)])
+ vols3[c.cid] = sumVol([ts3.cells[i] for i in tSearch2(mss, c.cid)])
+
+ return vols1, vols2, vols3
+
+
+def mergeCellVolsRef2(tss, mss):
+ (ts1, ts2) = tss
+
+ vols1 = dict()
+ vols2 = dict()
+
+ for c in ts1:
+ print(c.cid)
+ vols1[c.cid] = c.vol
+ vols2[c.cid] = sumVol([ts2.cells[i] for i in tSearch1(mss, c.cid)])
+
+ return vols1, vols2
+
+
+def divOr0(n1, n2):
+ if n2 == 0: return 0
+ else: return n1 / n2
+
+
+def growthRates3(volss, m12):
+ vols1, vols2, vols3 = volss
+ gr2 = dict()
+ for i, vol in vols1.items():
+ gr12 = divOr0((vols2[i] - vols1[i]), vols1[i])
+ gr23 = divOr0((vols3[i] - vols2[i]), vols2[i])
+ gr2[i] = 0.5 * (gr12 + gr23)
+
+ gr22 = [zip(m12.get(i, []), repeat(gr, len(m12.get(i, []))))
+ for i, gr in gr2.items()]
+
+ return dict(sum(gr22, []))
+
+
+def growthRatesL2(volss):
+ vols1, vols2 = volss
+ gr1 = dict()
+ for i, vol in vols1.items():
+ gr1[i] = (vols2[i] - vols1[i]) / vols1[i]
+
+ return gr1
+
+
+def growthRatesR2(volss, m12):
+ vols1, vols2 = volss
+ gr2 = dict()
+ for i, vol in vols1.items():
+ gr2[i] = (vols2[i] - vols1[i]) / vols1[i]
+
+ gr22 = [zip(m12.get(i, []), repeat(gr, len(m12.get(i, []))))
+ for i, gr in gr2.items()]
+
+ return dict(sum(gr22, []))
+
+
+def mkSeries():
+ d = "../data/organism"
+ dExprs = "../data/geneExpression/new_patterns"
+
+ timepoints = [10, 40, 96, 120, 132]
+ lins = mkLinsSeries(readLins(linDataLoc))
+
+ tss = dict()
+ linss = dict()
+
+ for t1, t2 in zip(timepoints, timepoints[1:]):
+ tss[t1] = a.readDataT1(d, dExprs, t1)
+ tss[t2] = a.readDataT1(d, dExprs, t2)
+
+ linss[(t1, t2)] = mergeLinss(list(map(lambda x: x[1], between(lins, t1, t2))))
+
+ return tss, linss
+
+def mkSeriesIm0():
+ d = "../data/yr01"
+ linDataLoc = "../data/yr01/tracking"
+
+ timepoints = [0, 10, 40, 96, 120, 132]
+ lins = mkLinsSeries(readLins(linDataLoc))
+
+ tss = dict()
+ linss = dict()
+
+ for t in timepoints:
+ tss[t] = a.readDataT1Im(d, t)
+
+ for t1, t2 in zip(timepoints, timepoints[1:]):
+ linss[(t1, t2)] = mergeLinss(map(lambda x: x[1],
+ between(lins, t1, t2)))
+
+ return tss, linss
+
+
+def mkGrowthRates():
+ tss, linss = mkSeries()
+ grs = dict()
+
+ grs[10] = growthRatesL2(mergeCellVolsRef2((tss[10], tss[40]),
+ [linss[(10, 40)]]))
+ grs[40] = growthRates3(mergeCellVolsRef3((tss[10], tss[40], tss[96]),
+ [linss[(10, 40)], linss[(40, 96)]]),
+ linss[(10, 40)])
+ grs[96] = growthRates3(mergeCellVolsRef3((tss[40], tss[96], tss[120]),
+ [linss[(40, 96)], linss[(96, 120)]]),
+ linss[(40, 96)])
+ grs[120] = growthRates3(mergeCellVolsRef3((tss[96], tss[120], tss[132]),
+ [linss[(96, 120)], linss[(120, 132)]]),
+ linss[(96, 120)])
+ grs[132] = growthRatesR2(mergeCellVolsRef2((tss[120], tss[132]),
+ [linss[(120, 132)]]),
+ linss[(120, 132)])
+
+ return grs
+
+def readGeomImT(fid, t):
+ ressFn = "../data/yr" + fid + "/resolutions.txt"
+ segFn = ("../data/yr" +
+ fid +
+ "/segmentations/YR_" +
+ fid + "_" + str(t) +
+ "h_segmented.tif")
+
+ ress = readRess(ressFn)
+ img = readImage(segFn)
+ res = ress[t]
+
+ ts = Tissue.fromImage(data=img, res=res)
+ for c in ts:
+ c.swapZX()
+
+ return ts
+
+def mkSeriesGeomIm(fid):
+ ressFn = "../data/yr" + fid + "/resolutions.txt"
+ linDataLoc = "../data/yr" + fid + "/trackingFiles/"
+ ress = readRess(ressFn)
+ tpoints = sorted(ress.keys())
+
+ lins = mkLinsSeries(readLins(linDataLoc))
+ tss = dict()
+ linss = dict()
+
+ for t in tpoints:
+ tss[t] = readGeomImT(fid, t)
+
+ for t1, t2 in zip(tpoints, tpoints[1:]):
+ linss[(t1, t2)] = mergeLinss(map(lambda x:
+ x[1], between(lins, t1, t2)))
+
+ return tss, linss
+
+
+###################################
+def invertMap(ln):
+ iLin = [list(zip(ds, repeat(m, len(ds)))) for m, ds in ln.items()]
+
+ return dict(sum(iLin, []))
+
+def updateCells(ts, ts1, iLn):
+ for c1 in ts1:
+ try:
+ c = ts.cells[iLn[c1.cid]]
+ except KeyError:
+ print(c1.cid)
+ c = Cell(0, Vec3(0, 0, 0), 0, dict([(g, False) for g in geneNms]))
+ c1.exprs = dict([(g, c.exprs[g]) for g, v in c1.exprs.items()])
+
+ return
+
+def prev(t):
+ if t == 40: return 10
+ elif t == 96: return 40
+ elif t == 120: return 96
+ elif t == 132: return 120
+ else: return None
+
+def succ(t):
+ if t==10: return 40
+ elif t==40: return 96
+ elif t==96: return 120
+ elif t==120: return 132
+ else: return None
+
+def updateTs(tss, linss, tpoints):
+ #evolve patterns based on lineage
+ tssLin = {}
+ for t in tpoints[1:]:
+ tssLin[(prev(t), t)] = tss[t].updateExprsTs(tss[prev(t)],
+ invertMap(linss[(prev(t), t)]))
+
+ return tssLin
+
+def updateTsBF(tss, fs):
+ return dict((t, tss[t].updateExprsBF(fs[t])) for t in tss.keys())
+
+def toNewmanFn(fn):
+ (fn, ext) = os.path.splitext(fn)
+ return fn + "000" + ".tif"
+
+def plotNewman1(tss, i, geneNm, outDir="."):
+ from subprocess import call
+ import os
+
+ confFile = "newman.conf"
+ visCmd = "/home2/Organism/tools/plot/bin/newman"
+ mergeCmd = "convert"
+
+ b.cellsToDat("ts1.dat", list(tss))
+
+ call([visCmd,
+ "-shape", "sphereVolume",
+ "-d", "3",
+ "ts1.dat",
+ "-column", str(i+1),
+ "-output", "tiff",
+ geneNm,
+ "-camera", confFile,
+ "-min", str(0),
+ "-max", str(1)])
+
+ os.remove("ts1.dat")
+
+ return
+
+def plotNewman2(tss, tss1, i, geneNm, outDir="."):
+ from os.path import join
+ from subprocess import call
+ import os
+
+ confFile = "newman.conf"
+ visCmd = "/home2/Organism/tools/plot/bin/newman"
+ mergeCmd = "convert"
+
+ b.cellsToDat("ts1.dat", list(tss))
+ b.cellsToDat("ts2.dat", list(tss1))
+
+ call([visCmd,
+ "-shape", "sphereVolume",
+ "-d", "3",
+ "ts1.dat",
+ "-column", str(i+1),
+ "-output", "tiff",
+ "ts1",
+ "-camera", confFile,
+ "-min", str(0),
+ "-max", str(1)])
+
+ call([visCmd,
+ "-shape", "sphereVolume",
+ "-d", "3",
+ "ts2.dat",
+ "-column", str(i+1),
+ "-output", "tiff",
+ "ts2",
+ "-camera", confFile,
+ "-min", str(0),
+ "-max", str(1)])
+
+ call([mergeCmd,
+ toNewmanFn("ts1.tif"),
+ toNewmanFn("ts2.tif"),
+ "+append",
+ "-pointsize", "20",
+ "-fill", "white",
+ "-annotate", "+10+40", geneNm,
+ join(outDir, geneNm+".tif")])
+
+ os.remove("ts1.dat")
+ os.remove("ts2.dat")
+ os.remove(toNewmanFn("ts1.tif"))
+ os.remove(toNewmanFn("ts2.tif"))
+
+
+def tFn(t, ref=""):
+ return "ts_t" + str(t) + ref + ".txt"
+
+
+def toOrgs(tss, ref="", outDir="."):
+ for t, ts in tss.items():
+ b.cellsToOrg(join(outDir, tFn(t, ref)), list(ts))
+
+ return
+
+def snd(x):
+ return x[1]
+
+def betweenTs(t, bounds):
+ ts, te = bounds
+ return t >= ts and t <= te
+
+def filterLinssBetween(linss, bounds):
+ t1, t2 = bounds
+ linssL = sorted([(ts, ln) for ts, ln in linss.items()
+ if betweenTs(ts[0], (t1, t2)) and betweenTs(ts[1], (t1, t2))], key=lambda p: p[0][0])
+
+ return list(map(snd, linssL))
+
+def filterTssBetween(tss, tbounds):
+ t1, t2 = tbounds
+
+ tssL = sorted([(t, ts) for t, ts in tss.items() if betweenTs(t, (t1, t2))], key=fst)
+
+ return list(map(snd, tssL))
+
diff --git a/common/common/seg.py b/common/common/seg.py
new file mode 100644
index 0000000000000000000000000000000000000000..b4b9e6533052295ed5d0076984385e4d01132405
--- /dev/null
+++ b/common/common/seg.py
@@ -0,0 +1,594 @@
+import numpy as np
+import ast
+from tifffile import TiffFile
+import matplotlib.pyplot as plt
+from copy import deepcopy
+from os.path import join
+
+geneNms = ['AG',
+ 'AHP6',
+ 'ANT',
+ 'AP1',
+ 'AP2',
+ 'AP3',
+ 'AS1',
+ 'ATML1',
+ 'CLV3',
+ 'CUC1_2_3',
+ 'ETTIN',
+ 'FIL',
+ 'LFY',
+ 'MP',
+ 'PHB_PHV',
+ 'PI',
+ 'PUCHI',
+ 'REV',
+ 'SEP1',
+ 'SEP2',
+ 'SEP3',
+ 'STM',
+ 'SUP',
+ 'SVP',
+ 'WUS']
+
+class Vec3(object):
+ def __init__(self, x, y, z):
+ self.x = x
+ self.y = y
+ self.z = z
+
+ def dist(self, v1):
+ return np.sqrt((self.x - v1.x)**2 +
+ (self.y - v1.y)**2 +
+ (self.z - v1.z)**2)
+
+ def toOrganism(self, sep=' '):
+ return sep.join([str(self.x), str(self.y), str(self.z)])
+
+ def swapZX(self):
+ xt = self.x
+ self.x = self.z
+ self.z = xt
+
+ def toArray(self):
+ return np.array([self.x, self.y, self.z])
+
+ def __repr__(self):
+ return self.toOrganism()
+
+class Cell(object):
+ def __init__(self, cid, pos, vol, exprs=None):
+ self.cid = cid #Int
+ self.pos = pos #Vec3
+ self.vol = vol #Real
+ self.exprs = exprs #Map String Bool
+
+ def dist(self, c1):
+ return self.pos.dist(c1.pos)
+
+ def addExprs(self, exprs, geneNms):
+ for gn in geneNms:
+ if not gn in exprs.keys():
+ exprs[gn] = False
+
+ self.exprs = exprs
+ self.geneNms = geneNms
+
+ return
+
+ def updateExprs(self, gfs):
+ gexprsUpd = dict()
+
+ for g, _ in self.exprs.items():
+ if g in gfs:
+ gf = gfs[g]
+ (acts, reprs, out, fsAct, fsRepr) = gf
+ inputs = acts + reprs
+ vals = [self.exprs[i] for i in inputs]
+ gexprsUpd[g] = evalI((acts, reprs, out), vals, fsAct, fsRepr)
+ else:
+ gexprsUpd[g] = self.exprs[g]
+
+ ncell = Cell(self.cid, self.pos, self.vol)
+ ncell.addExprs(gexprsUpd, deepcopy(self.geneNms))
+
+ return ncell
+
+ def updateExprsC(self, c):
+ exprsUpd = dict([(g, c.exprs[g]) for g, v in self.exprs.items()])
+ nc = Cell(self.cid, self.pos, self.vol, exprsUpd)
+ nc.geneNms = self.geneNms
+ return nc
+
+ def _boolToInt(self, b):
+ if b: return 1
+
+ return 0
+
+ def _exprsToOrg(self, sep=' '):
+ exprs = [str(self._boolToInt(self.exprs[gn])) for gn in self.geneNms]
+
+ return sep.join(exprs)
+
+ def swapZX(self):
+ self.pos.swapZX()
+
+ @classmethod
+ def fromTV(cls, nelems):
+ def parsePos(spos):
+ pos = list(np.fromstring(spos.replace('[', '').replace(']', ''),
+ sep=' '))
+ return pos
+
+ cid = int(nelems['cid'])
+ vol = float(nelems['volume'])
+ [xc, yc, zc] = parsePos(nelems['center'])
+ pos = Vec3(xc, yc, zc)
+
+ return cls(cid, pos, vol)
+
+ @classmethod
+ def fromImage(cls, cd, res):
+ volVoxel = np.prod(res)
+ xc, yc, zc = cd.centroid
+ nVoxels = cd.area
+
+ xc, yc, zc = np.multiply(np.array([xc, yc, zc]), res)
+
+ vol = nVoxels * volVoxel
+ pos = Vec3(xc, yc, zc)
+ cid = cd.label
+
+ return cls(cid, pos, vol)
+
+ def toOrganism(self):
+ if self.exprs:
+ return (' '.join([self.pos.toOrganism(), str(self.vol),
+ str(self.cid), self._exprsToOrg()]))
+ else:
+ return (' '.join([self.pos.toOrganism(), str(self.vol),
+ str(self.cid)]))
+
+ def toText(self, sep=" "):
+ if self.exprs:
+ return (sep.join([str(self.cid),
+ self.pos.toOrganism(sep=sep),
+ str(self.vol),
+ self._exprsToOrg(sep=sep)]))
+ else:
+ return sep.join([str(self.cid), self.pos.toOrganism(sep=sep)])
+
+ def centreToArray(self):
+ return self.pos.toArray()
+
+ def getOnGenes(self):
+ return frozenset([g for g, st in self.exprs.items() if st]) - frozenset(['WUS'])
+
+ def __repr__(self):
+ return self.toOrganism()
+
+class Tissue(object):
+ def __init__(self, cells, neighs, ecellIds = None):
+ self.cells = cells #Map Int Cell
+ self.neighs = neighs #Map Int (Map Int Dist)
+ self.L1_cells = ecellIds #Set Int
+
+ def __iter__(self):
+ #iterate over a list instead of map
+ return iter([v for _, v in self.cells.items()])
+
+ @classmethod
+ def _extractCells(cls, data, res):
+ from skimage.measure import regionprops
+ cellsSeg = regionprops(data)
+
+ cells = dict()
+ for cd in cellsSeg:
+ cells[int(cd.label)] = Cell.fromImage(cd, res)
+
+ return cells
+
+ @classmethod
+ def _getConn(cls, cells, data):
+ from skimage.future.graph import RAG
+
+ connG = RAG(data)
+ neigh = dict()
+
+ for n in connG.nodes:
+ nbrs = connG.adj[n].keys()
+ ds = [cells[n].dist(cells[nb]) for nb in nbrs]
+
+ neigh[n] = dict(zip(nbrs, ds))
+
+ return neigh
+
+ @classmethod
+ def fromTV(cls, fnGeom, fnNeigh):
+ def parseLn(ln):
+ ln = str.strip(ln)
+ nelems = dict()
+ elems = str.split(ln, ',')
+ for elm in elems:
+ [name, val] = str.split(elm, ':')
+ name = str.strip(name)
+ val = str.strip(val)
+ nelems[name] = val
+
+ return nelems
+
+ cells = dict()
+ neighs = dict()
+ ecellIds = []
+
+ with open(fnGeom, 'r') as f:
+ for ln in f:
+ nelems = parseLn(ln)
+ c = Cell.fromTV(nelems)
+ cells[c.cid] = c
+
+ with open(fnNeigh, 'r') as f:
+ for ln in f:
+ nelems = ast.literal_eval(ln)
+ cid = nelems['cid']
+ ns = nelems['neighbors ids']
+
+ if 1 in ns:
+ ecellIds.append(cid)
+ ns.remove(1)
+
+ neighs[cid] = dict((n, cells[cid].dist(cells[n])) for n in ns)
+
+ return cls(cells, neighs, set(ecellIds))
+
+ @classmethod
+ def fromImage(cls, data, res):
+ cells = cls._extractCells(data, res)
+ conn = cls._getConn(cells, data)
+
+ return cls(cells, conn)
+
+ def adjView(self):
+ return dict([(cid, ds.keys()) for cid, ds in self.neighs.items()])
+
+ def relabelFromMap(self, relabelM):
+ neighsR = []
+
+ for c in self:
+ c.cid = relabelM[c.cid]
+
+ for c, dists in self.neighs.items():
+ reNs = dict([(relabelM[n], d) for n, d in dists.items()])
+ neighsR.append((relabelM[c], reNs))
+
+ self.neighs = dict(neighsR)
+
+ return
+
+ def relabel(self):
+ relabelM = dict([(c.cid, i) for i, c in enumerate(self)])
+
+ self.relabelFromMap(relabelM)
+
+ return
+
+ def filterL1(self):
+ if not self.neighs:
+ return
+
+ L1_cell_ids = self.neighs[1]
+
+ L1_cells = dict([(c.cid, c) for c in self if c.cid in L1_cell_ids])
+ L1_neighs = dict([(cid, nbrs) for cid, nbrs in deepcopy(self.neighs.items())
+ if cid in L1_cells_ids])
+
+ return Tissue(L1_cells, L1_neighs)
+
+ def toNeigh(self):
+
+ def cellToNeigh(cid, nbs):
+ nbsIds = nbs.keys()
+ nNbs = len(nbsIds)
+ nbsDists = [nbs[nid] for nid in nbsIds]
+
+ return " ".join([str(cid), str(nNbs)]
+ + map(str, nbsIds)
+ + map(str, nbsDists))
+
+ ncells = sum(1 for _ in self)
+ header = str(ncells) + " " + "1"
+
+ cellReprs = ([cellToNeigh(cid, nbs)
+ for cid, nbs in self.neighs.items()])
+
+ return "\n".join([header] + cellReprs)
+
+ def toOrganism(self):
+ headerComment = "#x y z vol id"
+ ns = str(len(self.cells)) + " 5"
+ content = '\n'.join([cell.toOrganism() for cell in self])
+
+ return "\n".join([headerComment, ns, content])
+
+ def centresToArray(self):
+ return np.array([c.centreToArray() for c in self])
+
+
+class STissue(Tissue):
+ def __init__(self, cells, neighs, gexprs=None, geneNms=None):
+ super(STissue, self).__init__(cells, neighs)
+ if gexprs and geneNms:
+ self.addExprs(gexprs, geneNms)
+
+ def _defGeneExprs(self):
+ return dict([(gn, False) for gn in self.geneNms])
+
+ def addExprs(self, gexprs, geneNms):
+ self.geneNms = geneNms
+ for cell in self:
+ try:
+ cell.addExprs(gexprs[cell.cid], geneNms)
+ except KeyError:
+ cell.addExprs(self._defGeneExprs(), geneNms)
+
+ return
+
+ def relabel(self):
+ super(STissue, self).relabel()
+
+ @classmethod
+ def fromImage(cls, data, res, gexprs, geneNms):
+ ts = Tissue.fromImage(data, res)
+ return cls(ts.cells, ts.neighs, gexprs, geneNms)
+
+ @classmethod
+ def fromTV(cls, fnGeom, fnNeigh, gexprs, geneNms):
+ ts = Tissue.fromTV(fnGeom, fnNeigh)
+ return cls(ts.cells, ts.neighs, gexprs, geneNms)
+
+ def toOrganism(self):
+ headerComment = "#x y z vol id " + ' '.join(self.geneNms)
+ ns = str(len(self.cells)) + " " + str(len(self.geneNms)+5)
+ content = '\n'.join([cell.toOrganism() for cell in self])
+
+ return "\n".join([headerComment, ns, content])
+
+ def updateExprsBF(self, gfs):
+ ncells = dict([(c.cid, c.updateExprs(gfs)) for c in self])
+ return STissue(ncells, deepcopy(self.neighs))
+
+ def updateExprsTs(self, ts, iLn):
+ ncells = {}
+ for c1 in self:
+ try:
+ c = ts.cells[iLn[c1.cid]]
+ except KeyError:
+ c = Cell(0, Vec3(0, 0, 0), 0, dict([(g, False) for g in geneNms]))
+
+ ncells[c1.cid] = c1.updateExprsC(c)
+
+ return STissue(ncells, deepcopy(self.neighs))
+
+ def _boolToInt(b):
+ if b: return 1
+ else: return 0
+
+ def toDict(self):
+ tsDict = {}
+
+ for c in self:
+ tsDict[c.cid] = dict(((g, self._boolToInt(v)) for g, v in c.exprs))
+
+ return tsDict
+
+ def toCSV(self):
+ headerComment = "id,x,y,z,vol," + ','.join(self.geneNms)
+
+ content = '\n'.join([cell.toText(sep=",") for cell in self])
+
+ return "\n".join([headerComment, content])
+
+
+############### reading
+def getResolution(f):
+ try:
+ x1, x2 = f.pages.pages[0].tags['XResolution']
+ y1, y2 = f.pages.pages[0].tags['YResolution']
+ xres = x2 / x1
+ yres = y2 / y1
+ zres = 1.0
+ except:
+ xres = 1.0
+ yres = 1.0
+ zres = 1.0
+
+ return np.array([xres, yres, zres])
+
+def readImage(fseg):
+ f = TiffFile(fseg)
+ seg = f.asarray().astype(int)
+
+ return seg
+
+def readImageRegionPoints(fseg, res=[1, 1, 1], d=2):
+ from skimage.measure import regionprops
+
+ cell_regions = dict()
+
+ seg = readImage(fseg)
+ cellsSeg = regionprops(seg)
+
+ for cregion in cellsSeg:
+ cell_regions[cregion.label] = np.multiply(cregion.coords[:, :d], res[:d])
+
+ return cell_regions
+
+def readRegionVolume(fseg, res=[1, 1, 1]):
+ from skimage.measure import regionprops
+
+ cell_regions = dict()
+
+ seg = readImage(fseg)
+ cellsSeg = regionprops(seg)
+
+ for cregion in cellsSeg:
+ cell_regions[cregion.label] = cregion.area*np.product(res)
+
+ return cell_regions
+
+
+def mergePoints(cell_regions, cids):
+ return np.vstack(tuple([cell_regions[cid] for cid in cids]))
+
+def getCellRegions(dataDir, fid="01", d=2):
+ segFn = "yr{fid}/segmentations/YR_{fid}_{t}h_segmented.tif"
+ ressFn = "yr" + fid + "/resolutions.txt"
+
+ ress = readRess(join(dataDir, ressFn))
+ tpoints = sorted(ress.keys())
+
+ fsegs = dict([(t, join(dataDir, segFn.format(fid=fid, t=str(t)))) for t in tpoints])
+
+ return dict([(t, readImageRegionPoints(fn, res=ress[t], d=d))
+ for t, fn in fsegs.items()])
+
+def getCellRegionVolumesF1(dataDir, ress):
+ fsegs = {10: "yr01/segmentations/YR_01_10h_segmented.tif",
+ 40: "yr01/segmentations/YR_01_40h_segmented.tif",
+ 96: "yr01/segmentations/YR_01_96h_segmented.tif",
+ 120: "yr01/segmentations/YR_01_120h_segmented.tif",
+ 132: "yr01/segmentations/YR_01_132h_segmented.tif"}
+
+ return dict([(t, readRegionVolume(join(dataDir, fn), res=ress[t]))
+ for t, fn in fsegs.items()])
+
+def readExprs(fn):
+ #returns: Map Int (Map String Bool) Cid -> (GeneName -> {On/Off})
+ gexprs = dict()
+
+ with open(fn) as f:
+ geneNms = f.readline().rstrip().split(' ')[1:]
+
+ for line in f:
+ els = line.rstrip().split(' ')
+ cid = els[0]
+ exprs = els[1:]
+ gexprs[int(cid)] = dict(zip(geneNms, [bool(int(exp))
+ for exp in exprs]))
+
+ return gexprs
+
+def writeOrg(fnGeom, fnNeigh, ts):
+ with open(fnGeom, 'w+') as f:
+ f.write(ts.toOrganism())
+
+ with open(fnNeigh, 'w+') as f:
+ f.write(ts.toNeigh())
+
+ return
+
+def readRess(fn):
+ import ast
+
+ with open(fn, 'r') as f:
+ sress = ("{" + f.read().replace("=", ":").replace("h", "").
+ replace("\n", "").replace("t", "") + "}")
+ ress = ast.literal_eval(sress)
+
+ return ress
+
+#functions for common tasks
+def mkOrgs(d, dExprs, dOut):
+ from os.path import join
+
+ timepoints = [10, 40, 96, 120, 132]
+ geomF = "_segmented_tvformat_0_volume_position.txt"
+ neighF = "_segmented_tvformat_0_neighbors_py.txt"
+ exprF = "h.txt"
+
+ initF = ".init"
+ orgNeighF = ".neigh"
+
+ for t in timepoints:
+ print(t)
+ fnGeom = join(d, "t" + str(t) + geomF)
+ fnNeigh = join(d, "t" + str(t) + neighF)
+ fnExpr = join(dExprs, "t_" + str(t) + exprF)
+ fnGeomOut = join(dOut, "t" + str(t) + initF)
+ fnNeighOut = join(dOut, "t" + str(t) + orgNeighF)
+
+ gexprs = readExprs(fnExpr)
+ ts = STissue.fromTV(fnGeom, fnNeigh, gexprs, geneNms)
+
+ writeOrg(fnGeomOut, fnNeighOut, ts)
+
+ return
+
+def mkOrgGeomState(fim, fseg, fexpr, fout):
+ (data, res) = readImage(fim, fseg)
+ gexprs = readExprs(fexpr)
+
+ ts = STissue.fromImage(data, res, gexprs, geneNms)
+
+ writeOrg(fout, ts)
+
+def mkOrgGeom(fim, fseg, fout):
+ (data, res) = readImage(fim, fseg)
+ ts = Tissue.fromImage(data, res)
+
+ writeOrg(fout, ts)
+
+
+def plotTissue(ts, hl):
+ xs = [c.pos.x for c in ts]
+ ys = [c.pos.y for c in ts]
+ zs = [c.pos.z for c in ts]
+
+ cs = ['red' if c.cid in hl else 'blue' for c in ts]
+
+ fig = plt.figure()
+ ax = fig.add_subplot(211, projection='3d')
+ ax.view_init(elev=95, azim=270)
+ ax.scatter(xs, ys, zs, c=cs, alpha=0.9, s=200,
+ edgecolors='black', linewidths=0.5)
+
+ plt.show()
+
+def filterL1(ts):
+ if not ts.neighs:
+ return
+
+ L1_cell_ids = ts.neighs[1]
+
+ L1_cells = dict([(c.cid, c) for c in ts if c.cid in L1_cell_ids])
+ L1_neighs = dict([(cid, nbrs) for cid, nbrs in deepcopy(list(ts.neighs.items()))
+ if cid in L1_cell_ids])
+
+ return Tissue(L1_cells, L1_neighs)
+
+def const2(a, b): return b
+
+def folds1(vals, fs, iVal):
+ acc = iVal
+ for i, v in enumerate(vals):
+ acc = fs[i](acc, v)
+
+ return acc
+
+def folds(vals, fs):
+ fs = (const2, ) + fs
+
+ return folds1(vals, fs, False)
+
+def evalI(intr, vals, fsAct, fsRepr):
+ acts = intr[0]
+ reprs = intr[1]
+
+ nActs = len(acts)
+ nReprs = len(reprs)
+
+ if vals[nActs:]:
+ return (folds(vals[:nActs], fsAct) and
+ not folds(vals[nActs:], fsRepr))
+ else:
+ return folds(vals[:nActs], fsAct)
diff --git a/common/common/seg.py~ b/common/common/seg.py~
new file mode 100644
index 0000000000000000000000000000000000000000..8fe48542254b939da38243fc05e7587c0dc0170e
--- /dev/null
+++ b/common/common/seg.py~
@@ -0,0 +1,601 @@
+from __future__ import division
+from __future__ import print_function
+import numpy as np
+import ast
+from tifffile import TiffFile
+import matplotlib.pyplot as plt
+import boolFs as b
+from copy import deepcopy
+from mpl_toolkits.mplot3d import axes3d
+
+
+geneNms = ['AG',
+ 'AHP6',
+ 'ANT',
+ 'AP1',
+ 'AP2',
+ 'AP3',
+ 'AS1',
+ 'ATML1',
+ 'CLV3',
+ 'CUC1_2_3',
+ 'ETTIN',
+ 'FIL',
+ 'LFY',
+ 'MP',
+ 'PHB_PHV',
+ 'PI',
+ 'PUCHI',
+ 'REV',
+ 'SEP1',
+ 'SEP2',
+ 'SEP3',
+ 'STM',
+ 'SUP',
+ 'SVP',
+ 'WUS']
+
+class Vec3(object):
+ def __init__(self, x, y, z):
+ self.x = x
+ self.y = y
+ self.z = z
+
+ def dist(self, v1):
+ return np.sqrt((self.x - v1.x)**2 +
+ (self.y - v1.y)**2 +
+ (self.z - v1.z)**2)
+
+ def toOrganism(self, sep=' '):
+ return sep.join([str(self.x), str(self.y), str(self.z)])
+
+ def swapZX(self):
+ xt = self.x
+ self.x = self.z
+ self.z = xt
+
+ def toArray(self):
+ return np.array([self.x, self.y, self.z])
+
+ def __repr__(self):
+ return self.toOrganism()
+
+class Cell(object):
+ def __init__(self, cid, pos, vol, exprs=None):
+ self.cid = cid #Int
+ self.pos = pos #Vec3
+ self.vol = vol #Real
+ self.exprs = exprs #Map String Bool
+
+ def dist(self, c1):
+ return self.pos.dist(c1.pos)
+
+ def addExprs(self, exprs, geneNms):
+ for gn in geneNms:
+ if not gn in exprs.keys():
+ exprs[gn] = False
+
+ self.exprs = exprs
+ self.geneNms = geneNms
+
+ return
+
+ def updateExprs(self, gfs):
+ gexprsUpd = dict()
+
+ for g, _ in self.exprs.items():
+ if g in gfs:
+ gf = gfs[g]
+ (acts, reprs, out, fsAct, fsRepr) = gf
+ inputs = acts + reprs
+ vals = [self.exprs[i] for i in inputs]
+ gexprsUpd[g] = b.evalI((acts, reprs, out), vals, fsAct, fsRepr)
+ else:
+ gexprsUpd[g] = self.exprs[g]
+
+ ncell = Cell(self.cid, self.pos, self.vol)
+ ncell.addExprs(gexprsUpd, deepcopy(self.geneNms))
+
+ return ncell
+
+ def updateExprsC(self, c):
+ exprsUpd = dict([(g, c.exprs[g]) for g, v in self.exprs.items()])
+ nc = Cell(self.cid, self.pos, self.vol, exprsUpd)
+ nc.geneNms = self.geneNms
+ return nc
+
+ def _boolToInt(self, b):
+ if b: return 1
+
+ return 0
+
+ def _exprsToOrg(self, sep=' '):
+ exprs = [str(self._boolToInt(self.exprs[gn])) for gn in self.geneNms]
+
+ return sep.join(exprs)
+
+ def swapZX(self):
+ self.pos.swapZX()
+
+ @classmethod
+ def fromTV(cls, nelems):
+ def parsePos(spos):
+ pos = list(np.fromstring(spos.replace('[', '').replace(']', ''),
+ sep=' '))
+ return pos
+
+ cid = int(nelems['cid'])
+ vol = float(nelems['volume'])
+ [xc, yc, zc] = parsePos(nelems['center'])
+ pos = Vec3(xc, yc, zc)
+
+ return cls(cid, pos, vol)
+
+ @classmethod
+ def fromImage(cls, cd, res):
+ volVoxel = np.prod(res)
+ xc, yc, zc = cd.centroid
+ nVoxels = cd.area
+
+ xc, yc, zc = np.multiply(np.array([xc, yc, zc]), res)
+
+ vol = nVoxels * volVoxel
+ pos = Vec3(xc, yc, zc)
+ cid = cd.label
+
+ return cls(cid, pos, vol)
+
+ def toOrganism(self):
+ if self.exprs:
+ return (' '.join([self.pos.toOrganism(), str(self.vol),
+ str(self.cid), self._exprsToOrg()]))
+ else:
+ return (' '.join([self.pos.toOrganism(), str(self.vol),
+ str(self.cid)]))
+
+ def toText(self, sep=" "):
+ if self.exprs:
+ return (sep.join([str(self.cid),
+ self.pos.toOrganism(sep=sep),
+ str(self.vol),
+ self._exprsToOrg(sep=sep)]))
+ else:
+ return sep.join([str(self.cid), self.pos.toOrganism(sep=sep)])
+
+ def centreToArray(self):
+ return self.pos.toArray()
+
+ def getOnGenes(self):
+ return frozenset([g for g, st in self.exprs.items() if st]) - frozenset(['WUS'])
+
+ def __repr__(self):
+ return self.toOrganism()
+
+class Tissue(object):
+ def __init__(self, cells, neighs, ecellIds = None):
+ self.cells = cells #Map Int Cell
+ self.neighs = neighs #Map Int (Map Int Dist)
+ self.L1_cells = ecellIds #Set Int
+
+ def __iter__(self):
+ #iterate over a list instead of map
+ return iter([v for _, v in self.cells.items()])
+
+ @classmethod
+ def _extractCells(cls, data, res):
+ from skimage.measure import regionprops
+ cellsSeg = regionprops(data)
+
+ cells = dict()
+ for cd in cellsSeg:
+ cells[int(cd.label)] = Cell.fromImage(cd, res)
+
+ return cells
+
+ @classmethod
+ def _getConn(cls, cells, data):
+ from skimage.future.graph import RAG
+
+ connG = RAG(data)
+ neigh = dict()
+
+ for n in connG.nodes:
+ nbrs = connG.adj[n].keys()
+ ds = [cells[n].dist(cells[nb]) for nb in nbrs]
+
+ neigh[n] = dict(zip(nbrs, ds))
+
+ return neigh
+
+ @classmethod
+ def fromTV(cls, fnGeom, fnNeigh):
+ def parseLn(ln):
+ ln = str.strip(ln)
+ nelems = dict()
+ elems = str.split(ln, ',')
+ for elm in elems:
+ [name, val] = str.split(elm, ':')
+ name = str.strip(name)
+ val = str.strip(val)
+ nelems[name] = val
+
+ return nelems
+
+ cells = dict()
+ neighs = dict()
+ ecellIds = []
+
+ with open(fnGeom, 'r') as f:
+ for ln in f:
+ nelems = parseLn(ln)
+ c = Cell.fromTV(nelems)
+ cells[c.cid] = c
+
+ with open(fnNeigh, 'r') as f:
+ for ln in f:
+ nelems = ast.literal_eval(ln)
+ cid = nelems['cid']
+ ns = nelems['neighbors ids']
+
+ if 1 in ns:
+ ecellIds.append(cid)
+ ns.remove(1)
+
+ neighs[cid] = dict((n, cells[cid].dist(cells[n])) for n in ns)
+
+ return cls(cells, neighs, set(ecellIds))
+
+ @classmethod
+ def fromImage(cls, data, res):
+ cells = cls._extractCells(data, res)
+ conn = cls._getConn(cells, data)
+
+ return cls(cells, conn)
+
+ def adjView(self):
+ return dict([(cid, ds.keys()) for cid, ds in self.neighs.items()])
+
+ def relabelFromMap(self, relabelM):
+ neighsR = []
+
+ for c in self:
+ c.cid = relabelM[c.cid]
+
+ for c, dists in self.neighs.items():
+ reNs = dict([(relabelM[n], d) for n, d in dists.items()])
+ neighsR.append((relabelM[c], reNs))
+
+ self.neighs = dict(neighsR)
+
+ return
+
+ def relabel(self):
+ relabelM = dict([(c.cid, i) for i, c in enumerate(self)])
+
+ self.relabelFromMap(relabelM)
+
+ return
+
+ def filterL1(self):
+ if not self.neighs:
+ return
+
+ L1_cell_ids = self.neighs[1]
+
+ L1_cells = dict([(c.cid, c) for c in self if c.cid in L1_cell_ids])
+ L1_neighs = dict([(cid, nbrs) for cid, nbrs in deepcopy(self.neighs.items())
+ if cid in L1_cells_ids])
+
+ return Tissue(L1_cells, L1_neighs)
+
+ def toNeigh(self):
+
+ def cellToNeigh(cid, nbs):
+ nbsIds = nbs.keys()
+ nNbs = len(nbsIds)
+ nbsDists = [nbs[nid] for nid in nbsIds]
+
+ return " ".join([str(cid), str(nNbs)]
+ + map(str, nbsIds)
+ + map(str, nbsDists))
+
+ ncells = sum(1 for _ in self)
+ header = str(ncells) + " " + "1"
+
+ cellReprs = ([cellToNeigh(cid, nbs)
+ for cid, nbs in self.neighs.items()])
+
+ return "\n".join([header] + cellReprs)
+
+ def toOrganism(self):
+ headerComment = "#x y z vol id"
+ ns = str(len(self.cells)) + " 5"
+ content = '\n'.join([cell.toOrganism() for cell in self])
+
+ return "\n".join([headerComment, ns, content])
+
+ def centresToArray(self):
+ return np.array([c.centreToArray() for c in self])
+
+
+class STissue(Tissue):
+ def __init__(self, cells, neighs, gexprs=None, geneNms=None):
+ super(STissue, self).__init__(cells, neighs)
+ if gexprs and geneNms:
+ self.addExprs(gexprs, geneNms)
+
+ def _defGeneExprs(self):
+ return dict([(gn, False) for gn in self.geneNms])
+
+ def addExprs(self, gexprs, geneNms):
+ self.geneNms = geneNms
+ for cell in self:
+ try:
+ cell.addExprs(gexprs[cell.cid], geneNms)
+ except KeyError:
+ cell.addExprs(self._defGeneExprs(), geneNms)
+
+ return
+
+ def relabel(self):
+ super(STissue, self).relabel()
+
+ @classmethod
+ def fromImage(cls, data, res, gexprs, geneNms):
+ ts = Tissue.fromImage(data, res)
+ return cls(ts.cells, ts.neighs, gexprs, geneNms)
+
+ @classmethod
+ def fromTV(cls, fnGeom, fnNeigh, gexprs, geneNms):
+ ts = Tissue.fromTV(fnGeom, fnNeigh)
+ return cls(ts.cells, ts.neighs, gexprs, geneNms)
+
+ def toOrganism(self):
+ headerComment = "#x y z vol id " + ' '.join(self.geneNms)
+ ns = str(len(self.cells)) + " " + str(len(self.geneNms)+5)
+ content = '\n'.join([cell.toOrganism() for cell in self])
+
+ return "\n".join([headerComment, ns, content])
+
+ def updateExprsBF(self, gfs):
+ ncells = dict([(c.cid, c.updateExprs(gfs)) for c in self])
+ return STissue(ncells, deepcopy(self.neighs))
+
+ def updateExprsTs(self, ts, iLn):
+ ncells = {}
+ for c1 in self:
+ try:
+ c = ts.cells[iLn[c1.cid]]
+ except KeyError:
+ c = Cell(0, Vec3(0, 0, 0), 0, dict([(g, False) for g in geneNms]))
+
+ ncells[c1.cid] = c1.updateExprsC(c)
+
+ return STissue(ncells, deepcopy(self.neighs))
+
+ def toDict(self):
+ tsDict = {}
+
+ for c in self:
+ tsDict[c.cid] = dict(((g, b.boolToInt(v)) for g, v in c.exprs))
+
+ return tsDict
+
+ def toCSV(self):
+ headerComment = "id,x,y,z,vol," + ','.join(self.geneNms)
+
+ content = '\n'.join([cell.toText(sep=",") for cell in self])
+
+ return "\n".join([headerComment, content])
+
+
+
+
+#plotting
+class IndexTracker(object):
+ def __init__(self, ax, X):
+ self.ax = ax
+ ax.set_title('use scroll wheel to navigate images')
+
+ self.X = X
+ self.slices, rows, cols = X.shape
+ self.ind = self.slices//2
+
+ self.im = ax.imshow(self.X[self.ind, :, :])
+ self.update()
+
+ def onscroll(self, event):
+ if event.button == 'up':
+ self.ind = (self.ind + 1) % self.slices
+ else:
+ self.ind = (self.ind - 1) % self.slices
+ self.update()
+
+ def update(self):
+ self.im.set_data(self.X[self.ind, :, :])
+ self.ax.set_ylabel('slice %s' % self.ind)
+ self.im.axes.figure.canvas.draw()
+
+def plot3d(data):
+ fig, ax = plt.subplots(1, 1)
+ tracker = IndexTracker(ax, data)
+ fig.canvas.mpl_connect('scroll_event', tracker.onscroll)
+ plt.show()
+
+
+############### reading
+def getResolution(f):
+ try:
+ x1, x2 = f.pages.pages[0].tags['XResolution']
+ y1, y2 = f.pages.pages[0].tags['YResolution']
+ xres = x2 / x1
+ yres = y2 / y1
+ zres = 1.0
+ except:
+ xres = 1.0
+ yres = 1.0
+ zres = 1.0
+
+ return np.array([xres, yres, zres])
+
+def readImage(fseg):
+ f = TiffFile(fseg)
+ seg = f.asarray().astype(int)
+
+ return seg
+
+def readImageRegionPoints(fseg, res=[1, 1, 1], d=2):
+ from skimage.measure import regionprops
+
+ cell_regions = dict()
+
+ seg = readImage(fseg)
+ cellsSeg = regionprops(seg)
+
+ for cregion in cellsSeg:
+ cell_regions[cregion.label] = np.multiply(cregion.coords[:, :d], res[:d])
+
+ return cell_regions
+
+def readRegionVolume(fseg, res=[1, 1, 1]):
+ from skimage.measure import regionprops
+
+ cell_regions = dict()
+
+ seg = readImage(fseg)
+ cellsSeg = regionprops(seg)
+
+ for cregion in cellsSeg:
+ cell_regions[cregion.label] = cregion.area*np.product(res)
+
+ return cell_regions
+
+
+
+def mergePoints(cell_regions, cids):
+ return np.vstack(tuple([cell_regions[cid] for cid in cids]))
+
+def getCellRegions(fid="01", d=2):
+ segFn = "../data/yr{fid}/segmentations/YR_{fid}_{t}h_segmented.tif"
+ ressFn = "../data/yr" + fid + "/resolutions.txt"
+ ress = readRess(ressFn)
+ tpoints = sorted(ress.keys())
+
+ fsegs = dict([(t, segFn.format(fid=fid, t=str(t))) for t in tpoints])
+
+ return dict([(t, readImageRegionPoints(fn, res=ress[t], d=d))
+ for t, fn in fsegs.items()])
+
+def getCellRegionVolumes(fid="01"):
+ fsegs = {10: "../data/yr01/segmentations/YR_01_10h_segmented.tif",
+ 40: "../data/yr01/segmentations/YR_01_40h_segmented.tif",
+ 96: "../data/yr01/segmentations/YR_01_96h_segmented.tif",
+ 120: "../data/yr01/segmentations/YR_01_120h_segmented.tif",
+ 132: "../data/yr01/segmentations/YR_01_132h_segmented.tif"}
+
+ return dict([(t, readRegionVolume(fn, res=ress[t]))
+ for t, fn in fsegs.items()])
+
+def readExprs(fn):
+ #returns: Map Int (Map String Bool) Cid -> (GeneName -> {On/Off})
+ gexprs = dict()
+
+ with open(fn) as f:
+ geneNms = f.readline().rstrip().split(' ')[1:]
+
+ for line in f:
+ els = line.rstrip().split(' ')
+ cid = els[0]
+ exprs = els[1:]
+ gexprs[int(cid)] = dict(zip(geneNms, [bool(int(exp))
+ for exp in exprs]))
+
+ return gexprs
+
+def writeOrg(fnGeom, fnNeigh, ts):
+ with open(fnGeom, 'w+') as f:
+ f.write(ts.toOrganism())
+
+ with open(fnNeigh, 'w+') as f:
+ f.write(ts.toNeigh())
+
+ return
+
+def readRess(fn):
+ import ast
+
+ with open(fn, 'r') as f:
+ sress = ("{" + f.read().replace("=", ":").replace("h", "").
+ replace("\n", "").replace("t", "") + "}")
+ ress = ast.literal_eval(sress)
+
+ return ress
+
+#functions for common tasks
+def mkOrgs(d, dExprs, dOut):
+ from os.path import join
+
+ timepoints = [10, 40, 96, 120, 132]
+ geomF = "_segmented_tvformat_0_volume_position.txt"
+ neighF = "_segmented_tvformat_0_neighbors_py.txt"
+ exprF = "h.txt"
+
+ initF = ".init"
+ orgNeighF = ".neigh"
+
+ for t in timepoints:
+ print(t)
+ fnGeom = join(d, "t" + str(t) + geomF)
+ fnNeigh = join(d, "t" + str(t) + neighF)
+ fnExpr = join(dExprs, "t_" + str(t) + exprF)
+ fnGeomOut = join(dOut, "t" + str(t) + initF)
+ fnNeighOut = join(dOut, "t" + str(t) + orgNeighF)
+
+ gexprs = readExprs(fnExpr)
+ ts = STissue.fromTV(fnGeom, fnNeigh, gexprs, geneNms)
+
+ writeOrg(fnGeomOut, fnNeighOut, ts)
+
+ return
+
+def mkOrgGeomState(fim, fseg, fexpr, fout):
+ (data, res) = readImage(fim, fseg)
+ gexprs = readExprs(fexpr)
+
+ ts = STissue.fromImage(data, res, gexprs, geneNms)
+
+ writeOrg(fout, ts)
+
+def mkOrgGeom(fim, fseg, fout):
+ (data, res) = readImage(fim, fseg)
+ ts = Tissue.fromImage(data, res)
+
+ writeOrg(fout, ts)
+
+
+def plotTissue(ts, hl):
+ xs = [c.pos.x for c in ts]
+ ys = [c.pos.y for c in ts]
+ zs = [c.pos.z for c in ts]
+
+ cs = ['red' if c.cid in hl else 'blue' for c in ts]
+
+ fig = plt.figure()
+ ax = fig.add_subplot(211, projection='3d')
+ ax.view_init(elev=95, azim=270)
+ ax.scatter(xs, ys, zs, c=cs, alpha=0.9, s=200,
+ edgecolors='black', linewidths=0.5)
+
+ plt.show()
+
+def filterL1(ts):
+ if not ts.neighs:
+ return
+
+ L1_cell_ids = ts.neighs[1]
+
+ L1_cells = dict([(c.cid, c) for c in ts if c.cid in L1_cell_ids])
+ L1_neighs = dict([(cid, nbrs) for cid, nbrs in deepcopy(list(ts.neighs.items()))
+ if cid in L1_cell_ids])
+
+ return Tissue(L1_cells, L1_neighs)
diff --git a/common/setup.py b/common/setup.py
new file mode 100644
index 0000000000000000000000000000000000000000..f75eef42d31f00e094818ba9d0148f4b3ae46e2c
--- /dev/null
+++ b/common/setup.py
@@ -0,0 +1,18 @@
+import setuptools
+
+with open("README.md", "r") as fh:
+ long_description = fh.read()
+
+setuptools.setup(
+ name="common", # Replace with your own username
+ version="0.0.1",
+ author="Argyris Zardilis",
+ author_email="argyris.zardilis@slcu.cam.ac.uk",
+ packages=setuptools.find_packages(),
+ classifiers=[
+ "Programming Language :: Python :: 3",
+ "License :: OSI Approved :: MIT License",
+ "Operating System :: OS Independent",
+ ],
+ python_requires='>=3.6',
+)
diff --git a/common/setup.py~ b/common/setup.py~
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391