hypothesisAnalyse.py (7641B)
1 import sys 2 import os 3 import matplotlib 4 import datetime 5 matplotlib.use('Agg') 6 7 8 import numpy as np 9 from numpy import exp, log 10 import matplotlib.pyplot as plt 11 import sys, os 12 import json 13 import pymultinest 14 15 16 17 18 modelSetB = [['singleGaussian', 'singleGaussian', 2, 2, ['mu', 'sigma', 'mu', 'sigma']], 19 ['singleGaussian', 'twoGaussian', 2, 5, ['mu', 'sigma', 'mu1', 'mu2', 'sigma1', 'sigma2', 'alpha']], 20 ['singleGaussian', 'uniform', 2, 2, ['mu', 'sigma', 'mMin', 'mMax']], 21 ['twoGaussian', 'singleGaussian', 5, 2, ['mu1', 'mu2', 'sigma1', 'sigma2', 'alpha', 'mu', 'sigma']], 22 ['twoGaussian', 'twoGaussian', 5, 5, ['mu1', 'mu2', 'sigma1', 'sigma2', 'alpha', 'mu1', 'mu2', 'sigma1', 'sigma2', 'alpha']], 23 ['twoGaussian', 'uniform', 5, 2, ['mu1', 'mu2', 'sigma1', 'sigma2', 'alpha', 'mMin', 'mMax']], 24 ['uniform', 'singleGaussian', 2, 2, ['mMin', 'mMax', 'mu', 'sigma']], 25 ['uniform', 'twoGaussian', 2, 5, ['mMin', 'mMax', 'mu1', 'mu2', 'sigma1', 'sigma2', 'alpha']], 26 ['uniform', 'uniform', 2, 2, ['mMin', 'mMax', 'mMin', 'mMax']]] 27 28 29 modelSetA = [['singleGaussian', 2, ['mu', 'sigma']], 30 ['twoGaussian', 5, ['mu1', 'mu2', 'sigma1', 'sigma2', 'alpha']], 31 ['uniform', 2, ['mMin', 'mMax']]] 32 33 def roundNumber(number): 34 return round(number, 4) 35 36 37 def cleanRound(number, dec=3): 38 newNum = str(round(number, 3)) 39 print(newNum) 40 missing = dec - len(newNum.split('.')[1]) 41 return newNum + '0' * missing 42 43 44 45 def readStats(statsFile): 46 globalEvidenceLine = statsFile.readline() 47 nestedGlobalEvidenceLine = statsFile.readline() 48 49 globalEvidenceLine 50 51 print(globalEvidenceLine.split()) 52 globalEvidence = float(globalEvidenceLine.split()[-3]) 53 print(globalEvidenceLine.split()[-1]) 54 evidenceStd = float(globalEvidenceLine.split()[-1]) 55 print(globalEvidence, evidenceStd) 56 57 return globalEvidence, evidenceStd 58 59 60 def fractionalEvidence(resultList, resultDirectory, saveName="fracResults.txt"): 61 with open(resultDirectory + saveName,"w+") as z: 62 63 totalEv = 0 64 for result in resultList: 65 totalEv += result[-2] 66 67 68 ### Hypo B 69 if len(resultList[0]) == 4: 70 for result1 in resultList: 71 tb = ' & ' 72 tableString = ' ' + result1[0][:4] + '-' + result1[1][:4] + tb + str(result1[2]/totalEv) + "\\" * 2 73 z.write(tableString + '\n') 74 75 else: 76 for result1 in resultList: 77 tb = ' & ' 78 tableString = ' ' + result1[0][:4] + tb + str(result1[1]/totalEv) + "\\" * 2 79 z.write(tableString + '\n') 80 return 81 82 83 84 def makeTableB(resultList, resultDirectory, saveName="tableResults.txt"): 85 with open(resultDirectory + saveName,"w+") as z: 86 87 reverseNameTitle = ' ' * 4 + ' & '.join([result[0][:4] + '-' + result[1][:4] for result in resultList[::-1]]) + "\\" * 2 88 z.write(reverseNameTitle + '\n') 89 90 for index, result1 in enumerate(resultList): 91 print(index, result1) 92 evidence1 = result1[2] 93 94 lineBFs = [] 95 #for result2 in resultList[:-index - 1]: 96 for result2 in resultList[:index:-1]: 97 evidence2 = result2[2] 98 bayesFactor = evidence1/evidence2 99 lineBFs.append(cleanRound(bayesFactor)) 100 101 tb = ' & ' 102 tableString = ' ' + result1[0][:4] + '-' + result1[1][:4] + tb + tb.join(lineBFs) + "\\" * 2 103 z.write(tableString + '\n') 104 return 105 106 107 108 def makeTableA(resultList, resultDirectory, saveName="tableResults.txt"): 109 with open(resultDirectory + saveName,"w+") as z: 110 111 reverseNameTitle = ' ' * 3 + ' & '.join([result[0][:4] for result in resultList[::-1]]) + "\\" * 2 112 z.write(reverseNameTitle + '\n') 113 114 for index, result1 in enumerate(resultList): 115 print(index, result1) 116 evidence1 = result1[1] 117 118 lineBFs = [] 119 #for result2 in resultList[:-index - 1]: 120 for result2 in resultList[:index:-1]: 121 evidence2 = result2[1] 122 bayesFactor = evidence1/evidence2 123 lineBFs.append(cleanRound(bayesFactor)) 124 125 tb = ' & ' 126 tableString = ' ' + result1[0][:4] + tb + tb.join(lineBFs) + "\\" * 2 127 z.write(tableString + '\n') 128 return 129 130 131 132 133 134 def analyseMainB(resultDirectory): 135 with open(resultDirectory + "collectedResults.txt","w+") as g: 136 g.write("Model1 ------ Model2 ------ Evidence ------ LogEvidence std\n") 137 resultsList = [] 138 for modelName1, modelName2, ndim1, ndim2, paramNames in modelSetB: 139 ndim = ndim1 + ndim2 140 print(modelName1, modelName2) 141 142 prefix = resultDirectory + modelName1[:4] + "/" + modelName2[:4] + "/" 143 144 with open(prefix + "stats.dat","r") as statsFile: 145 logevidence, evidenceStd = readStats(statsFile) 146 evidence = np.exp(logevidence) 147 print("{}, {}: {} +- {}\n".format(modelName1, modelName2, evidence, evidenceStd)) 148 149 150 g.write("{}, {}: {} +- {}\n".format(modelName1, modelName2, evidence, evidenceStd)) 151 152 resultsList.append([modelName1, modelName2, evidence, evidenceStd]) 153 154 155 evidenceSum = 0 156 for i in range(len(resultsList)): 157 evidenceSum += float(resultsList[i][2]) 158 g.write("Total evidence: {}".format(evidenceSum)) 159 160 #makeTableB(resultsList, resultDirectory) 161 162 163 sortedResults = sorted(resultsList, key = lambda x: x[2], reverse=True) 164 makeTableB(sortedResults, resultDirectory, "sortedResultsB.txt") 165 166 fractionalEvidence(sortedResults, resultDirectory, "fracResultsB.txt") 167 168 169 return 170 171 def analyseMainA(resultDirectory): 172 with open(resultDirectory + "collectedResults.txt","w+") as g: 173 g.write("Model------ Evidence ------ LogEvidence std\n") 174 175 resultsList = [] 176 for modelName, ndim, paramNames in modelSetA: 177 print(modelName) 178 179 prefix = resultDirectory + modelName[:4] + "/" 180 181 with open(prefix + "stats.dat","r") as statsFile: 182 logevidence, evidenceStd = readStats(statsFile) 183 evidence = np.exp(logevidence) 184 print("{}: {} +- {}\n".format(modelName, evidence, evidenceStd)) 185 g.write("{}: {} +- {}\n".format(modelName, evidence, evidenceStd)) 186 resultsList.append([modelName, evidence, evidenceStd]) 187 188 189 evidenceSum = 0 190 for i in range(len(resultsList)): 191 evidenceSum += float(resultsList[i][1]) 192 g.write("Total evidence: {}".format(evidenceSum)) 193 194 sortedResults = sorted(resultsList, key = lambda x: x[1], reverse=True) 195 makeTableA(sortedResults, resultDirectory, "sortedResultsA.txt") 196 197 fractionalEvidence(sortedResults, resultDirectory, "fracResultsA.txt") 198 return 199 200 201 ### Get results directory 202 203 if len(sys.argv) != 3: 204 print("ERROR") 205 print("Usage: hypothesisAnalyse.py [-A, -B] /out_directory/") 206 else: 207 outputfileName = sys.argv[2] 208 if sys.argv[1] == 'A': 209 analyseMainA(outputfileName) 210 elif sys.argv[1] == 'B': 211 analyseMainB(outputfileName) 212 else: 213 print("Argument Error") 214 215