compareRegPerfAA1.py 9.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233
  1. import os
  2. import pandas as pd
  3. import seaborn as sns
  4. from matplotlib import pyplot as plt
  5. from regmaxsn.core.matplotlibRCParams import mplPars
  6. from regmaxsn.core.occupancyBasedMeasure import occupancyEMD
  7. from regmaxsn.core.farthestPointStats import maxDistStats
  8. import numpy as np
  9. import sys
  10. plt.ion()
  11. sns.set(rc=mplPars)
  12. def standardizedExpNameLambda(x):
  13. if x.endswith("_Standardized"):
  14. return x[:-len("_Standardized")]
  15. else:
  16. return x
  17. homeFolder = "/home/aj/DataAndResults/morphology"
  18. expNames = [
  19. "VGlut-F-300181_Standardized",
  20. "VGlut-F-300181.CNG",
  21. "VGlut-F-400545.CNG",
  22. "VGlut-F-500778.CNG",
  23. "VGlut-F-300196.CNG",
  24. "VGlut-F-300288.CNG",
  25. "VGlut-F-600290.CNG",
  26. "VGlut-F-600499.CNG",
  27. "VGlut-F-400665.CNG",
  28. "VGlut-F-300142.CNG",
  29. "VGlut-F-500147.CNG",
  30. "VGlut-F-600181.CNG",
  31. "VGlut-F-700190.CNG"
  32. ]
  33. case1 = {'resDirs': {
  34. "PCA": os.path.join(homeFolder, "PCA-Based", "chiangAA1"),
  35. "BlastNeuron": os.path.join(homeFolder, "BlastNeuron", "chiangAA1"),
  36. "PCA + RobartsICP": os.path.join(homeFolder, "RobartsICP", "chiangAA1"),
  37. "Reg-MaxS": os.path.join(homeFolder, "Reg-MaxS", "chiangAA1"),
  38. "Reg-MaxS-N": os.path.join(homeFolder, "Reg-MaxS-N", "chiangAA1"),
  39. "Standardized": os.path.join(homeFolder, "Registered", "chiangAA1",),
  40. },
  41. 'initRef': "VGlut-F-300181\n-Standardized",
  42. 'expNameLambdas': {
  43. "PCA": lambda x: x,
  44. "BlastNeuron": lambda x: x,
  45. "PCA + RobartsICP": lambda x: x,
  46. "Reg-MaxS": lambda x: x,
  47. "Reg-MaxS-N": lambda x: x,
  48. "Standardized": standardizedExpNameLambda}
  49. }
  50. case2 = {'resDirs': {
  51. "PCA": os.path.join(homeFolder, "PCA-Based", "chiangAA1_VGlut-F-300181.CNG"),
  52. "BlastNeuron": os.path.join(homeFolder, "BlastNeuron", "chiangAA1_VGlut-F-300181.CNG"),
  53. "PCA + RobartsICP": os.path.join(homeFolder, "RobartsICP", "chiangAA1_VGlut-F-300181.CNG"),
  54. "Reg-MaxS": os.path.join(homeFolder, "Reg-MaxS", "chiangAA1_VGlut-F-300181.CNG"),
  55. "Reg-MaxS-N": os.path.join(homeFolder, "Reg-MaxS-N", "chiangAA1_VGlut-F-300181.CNG"),
  56. "Standardized": os.path.join(homeFolder, "Registered", "chiangAA1"),
  57. },
  58. 'initRef': "VGlut-F-300181",
  59. 'expNameLambdas': {
  60. "PCA": lambda x: x,
  61. "BlastNeuron": lambda x: x,
  62. "PCA + RobartsICP": lambda x: x,
  63. "Reg-MaxS": lambda x: x,
  64. "Reg-MaxS-N": lambda x: x,
  65. "Standardized": standardizedExpNameLambda}
  66. }
  67. case3 = {'resDirs': {
  68. "PCA": os.path.join(homeFolder, "PCA-Based", "chiangAA1_VGlut-F-600290.CNG"),
  69. "BlastNeuron": os.path.join(homeFolder, "BlastNeuron", "chiangAA1_VGlut-F-600290.CNG"),
  70. "PCA + RobartsICP": os.path.join(homeFolder, "RobartsICP", "chiangAA1_VGlut-F-600290.CNG"),
  71. "Reg-MaxS": os.path.join(homeFolder, "Reg-MaxS", "chiangAA1_VGlut-F-600290.CNG"),
  72. "Reg-MaxS-N": os.path.join(homeFolder, "Reg-MaxS-N", "chiangAA1_VGlut-F-600290.CNG"),
  73. "Standardized": os.path.join(homeFolder, "Registered", "chiangAA1"),
  74. },
  75. 'initRef': "VGlut-F-600290",
  76. 'expNameLambdas': {
  77. "PCA": lambda x: x,
  78. "BlastNeuron": lambda x: x,
  79. "PCA + RobartsICP": lambda x: x,
  80. "Reg-MaxS": lambda x: x,
  81. "Reg-MaxS-N": lambda x: x,
  82. "Standardized": standardizedExpNameLambda}
  83. }
  84. case4 = {'resDirs': {
  85. "PCA": os.path.join(homeFolder, "PCA-Based", "chiangAA1_VGlut-F-500147.CNG"),
  86. "BlastNeuron": os.path.join(homeFolder, "BlastNeuron", "chiangAA1_VGlut-F-500147.CNG"),
  87. "PCA + RobartsICP": os.path.join(homeFolder, "RobartsICP", "chiangAA1_VGlut-F-500147.CNG"),
  88. "Reg-MaxS": os.path.join(homeFolder, "Reg-MaxS", "chiangAA1_VGlut-F-500147.CNG"),
  89. "Reg-MaxS-N": os.path.join(homeFolder, "Reg-MaxS-N", "chiangAA1_VGlut-F-500147.CNG"),
  90. "Standardized": os.path.join(homeFolder, "Registered", "chiangAA1"),
  91. },
  92. 'initRef': "VGlut-F-500147",
  93. 'expNameLambdas': {
  94. "PCA": lambda x: x,
  95. "BlastNeuron": lambda x: x,
  96. "PCA + RobartsICP": lambda x: x,
  97. "Reg-MaxS": lambda x: x,
  98. "Reg-MaxS-N": lambda x: x,
  99. "Standardized": standardizedExpNameLambda}
  100. }
  101. cases = [case1, case2, case3, case4]
  102. voxelSize = 10
  103. def saveData(outXLFile):
  104. metricsDF = pd.DataFrame()
  105. maxDistStatsDF = pd.DataFrame()
  106. for case in cases:
  107. resDirs = case["resDirs"]
  108. initRef = case["initRef"]
  109. expNameLambdas = case['expNameLambdas']
  110. for (resDirLabel, resDir) in resDirs.iteritems():
  111. outFiles = []
  112. expNameLambda = expNameLambdas[resDirLabel]
  113. for expName in expNames:
  114. outFile = os.path.join(resDir, "{}.swc".format(expNameLambda(expName)))
  115. if os.path.isfile(outFile):
  116. outFiles.append(outFile)
  117. else:
  118. print("{} not found. Ignoring it.".format(outFile))
  119. if outFiles:
  120. print("Collecting data for resDirLabel={}, initRef={}".format(resDirLabel, initRef))
  121. metric = occupancyEMD(outFiles, voxelSize)
  122. if resDirLabel == "Reg-MaxS-N":
  123. finalRef = os.path.join(resDir, "finalRef.swc")
  124. initialRef = os.path.join(resDir, "ref-1.swc")
  125. runtime = os.stat(finalRef).st_mtime - os.stat(initialRef).st_mtime
  126. else:
  127. outFileModTimes = [os.stat(outFile).st_mtime for outFile in outFiles]
  128. outFileModTimesSorted = sorted(outFileModTimes)
  129. nFiles = float(len(outFiles))
  130. runtime = (outFileModTimesSorted[-1] - outFileModTimesSorted[0]) * nFiles / (nFiles - 1)
  131. tempDict = {"Initial Reference": initRef,
  132. "Occupancy Based Dissimilarity Measure": metric,
  133. "Total runtime (s)": runtime,
  134. "Method": resDirLabel}
  135. metricsDF = metricsDF.append(tempDict, ignore_index=True)
  136. # crdMaxDistsStatsDFFull = maxDistStats(outFiles)
  137. # aggDictRename = {"mean": "mean of \nmaximum distances",
  138. # "std": "standard deviation of \nmaximum distances"}
  139. # crdMaxDistsStatsDF = crdMaxDistsStatsDFFull.loc[:, ("maximum distance", "source point")]
  140. # crdMaxDistsStatsDF.set_index("source point", inplace=True)
  141. # crdMaxDistMeanStd = crdMaxDistsStatsDF.groupby("source point")["maximum distance"]\
  142. # .aggregate([np.mean, np.std])
  143. # crdMaxDistMeanStd = crdMaxDistMeanStd.rename(columns=aggDictRename)
  144. # tempDF = crdMaxDistMeanStd.reset_index()
  145. # tempDF['Initial Reference'] = initRef
  146. # tempDF['Method'] = resDirLabel
  147. # maxDistStatsDF = maxDistStatsDF.append(tempDF, ignore_index=True)
  148. else:
  149. print("No usable SWCs found in {}".format(resDir))
  150. metricsDF.to_excel(outXLFile)
  151. def plotData(inFile):
  152. [darkblue, green, red, violet, yellow, lightblue] = sns.color_palette()
  153. metricsDF = pd.read_excel(inFile)
  154. fig1, ax1 = plt.subplots(figsize=(14, 11.2))
  155. sns.barplot(data=metricsDF, x="Initial Reference",
  156. y="Occupancy Based Dissimilarity Measure", hue="Method",
  157. ax=ax1, hue_order=["PCA", "PCA + RobartsICP", "BlastNeuron",
  158. "Reg-MaxS", "Reg-MaxS-N", "Standardized"],
  159. palette=[red, violet, yellow, lightblue, darkblue, green])
  160. ax1.legend(loc='best', ncol=3)
  161. ax1.set_ylabel("Occupancy Based Dissimilarity Measure")
  162. temp = ax1.get_ylim()
  163. ax1.set_ylim(temp[0], temp[1] + 2)
  164. # fig2, ax2 = plt.subplots(figsize=(14, 11.2))
  165. # sns.boxplot(data=maxDistStatsDF, x="Initial Reference", y="mean of \nmaximum distances",
  166. # hue="Method", whis=np.inf,
  167. # ax=ax2, hue_order=["PCA", "BlastNeuron","PCA + RobartsICP",
  168. # "Reg-MaxS", "Reg-MaxS-N", "Standardized"])
  169. # ax2.legend(loc="best", ncol=3)
  170. #
  171. # fig3, ax3 = plt.subplots(figsize=(14, 11.2))
  172. # sns.boxplot(data=maxDistStatsDF, x="Initial Reference", y="standard deviation of \nmaximum distances",
  173. # hue="Method", whis=np.inf,
  174. # ax=ax3, hue_order=["PCA", "BlastNeuron","PCA + RobartsICP",
  175. # "Reg-MaxS", "Reg-MaxS-N", "Standardized"])
  176. # ax3.legend(loc="best", ncol=3)
  177. # for fig in [fig1, fig2, fig3]:
  178. for fig in [fig1]:
  179. fig.tight_layout()
  180. return fig1
  181. if __name__ == "__main__":
  182. assert len(sys.argv) == 3, "Improper Usage! Please use as:\n" \
  183. "python {fName} save <outFile> or python {fName} plot <inFile>".format(fName=sys.argv[0])
  184. if sys.argv[1] == "save":
  185. saveData(sys.argv[2])
  186. elif sys.argv[1] == "plot":
  187. fig = plotData(sys.argv[2])
  188. else:
  189. raise(ValueError("Improper Usage! Please use as:\n"
  190. "python {fName} save <outFile> or python {fName} plot <inFile>".format(fName=sys.argv[0])))