hypothesisB.py (7015B)
1 import pymultinest 2 import numpy as np 3 from scipy import stats 4 from scipy.special import erf 5 import os 6 import sys 7 from datetime import datetime 8 9 startTime = datetime.now() 10 11 ### 12 # Directory of current script 13 # Required by Condor jobs to find relevant samples etc. 14 topDirectory = os.path.dirname(os.path.realpath(__file__)) + '/' 15 16 # Load BNS mass samples 17 # (eg. 17x10000 array) 18 pulsarSamples = np.loadtxt(topDirectory + 'Samples/mr_samples.txt') 19 companionSamples = np.loadtxt(topDirectory + 'Samples/ms_samples.txt') 20 21 # For each BNS, group pulsar & companion samples into pairs. 22 # Creates array of shape like (17x10000x2) 23 bothMassSamples = np.stack((pulsarSamples, companionSamples), axis=-1) 24 25 # Only use first 1000 samples to increase speed 26 massSamples = bothMassSamples[:,:1000,:] 27 28 # Define the nMeasurements and nSamples per mass measurement 29 nSamples = len(massSamples[0]) 30 nMeasurements = len(massSamples) 31 32 33 34 ### Probability Distribution Functions 35 # Single Gaussian Functions 36 def evalSingleGaussian(theta, x): 37 mu, sig = theta[0], theta[1] 38 normalisingTerm = (0.5*(1+erf((2-mu)/(sig*2**0.5)) - (1+erf((0.8-mu)/(sig*2**0.5))))) 39 return stats.norm(mu, sig).pdf(x) * 1.0/normalisingTerm 40 41 # Two Gaussian Functions 42 def evalTwoGaussian(theta, x): 43 mu1, mu2, sig1, sig2, alpha = theta 44 normalisingTerm1 = (0.5*(1+erf((2-mu1)/(sig1*2**0.5)) - (1+erf((0.8-mu1)/(sig1*2**0.5))))) 45 normalisingTerm2 = (0.5*(1+erf((2-mu2)/(sig2*2**0.5)) - (1+erf((0.8-mu2)/(sig2*2**0.5))))) 46 return alpha * stats.norm(mu1, sig1).pdf(x) * 1.0/normalisingTerm1 + (1-alpha) * stats.norm(mu2, sig2).pdf(x) * 1.0/normalisingTerm2 47 48 # Uniform Functions 49 def evalUniform(theta, x): 50 mMin, mMax = theta[0], theta[1] 51 return stats.uniform(mMin, mMax-mMin).pdf(x) 52 53 54 55 ### Model Information Dictionaries 56 singleGaussianModel = { 57 "name": "singleGaussian", 58 "pdf": evalSingleGaussian, 59 "ndim": 2, 60 "params": ['mu', 'sigma'] 61 } 62 63 twoGaussianModel = { 64 "name": "twoGaussian", 65 "pdf": evalTwoGaussian, 66 "ndim": 5, 67 "params": ['mu1', 'mu2', 'sigma1', 'sigma2', 'alpha'] 68 } 69 70 uniformModel = { 71 "name": "uniform", 72 "pdf": evalUniform, 73 "ndim": 2, 74 "params": ['mMin', 'mMax'] 75 } 76 77 78 79 ### Prior for each model. 80 def prior(cube, ndim, nparams): 81 # cube is initially a unit hypercube which is to be mapped onto the relevant prior space. 82 83 84 # Hyperparameter Index. We map the priors beginning with the parameters of model1, 85 # and then incriment j by the number of parameters in that model1. Then mapping 86 # the parameters of the next model2. 87 j = 0 88 89 # Loop over both models. 90 for modelBeingMapped in [modelName1, modelName2]: 91 if modelBeingMapped == 'singleGaussian': 92 cube[j] = 0.8 + cube[j] * (2 - 0.8) 93 cube[j+1] = 0.005 + cube[j+1] * (0.5 - 0.005) 94 j += 2 95 if modelBeingMapped == 'twoGaussian': 96 cube[j] = 0.8 + cube[j] * (2 - 0.8) 97 cube[j+1] = cube[j] + cube[j+1] * (2 - cube[j]) 98 cube[j+2] = 0.005 + cube[j+2] * (0.5 - 0.005) 99 cube[j+3] = 0.005 + cube[j+3] * (0.5 - 0.005) 100 cube[j+4] = cube[j+4] * 1 101 j += 5 102 if modelBeingMapped == 'uniform': 103 cube[j] = 0.8 + cube[j] * (2 - 0.8) 104 cube[j+1] = cube[j] + cube[j+1] * (2 - cube[j]) 105 j += 2 106 107 # After this process the number of parameters mapped, j, should be equal to the number of dimensions for the problem (hyperparameters of model1 + model2). 108 if j != ndim: 109 print("SOME PARAMETERS WERE NOT MAPPED TO!!!") 110 print(ndim) 111 print(j) 112 113 return 114 115 116 117 ### Likelihood Function (Same as Farr et al.) 118 def likelihood(cube, ndim, nparams): 119 120 # Create lists of the parameters for each model. Model1 has parameters in cube from 0 to ndim1-1, Model2 has parameters in cube from ndim1 to ndim-1. 121 paramList1 = [cube[i] for i in range(ndim1)] 122 paramList2 = [cube[i] for i in range(ndim1, ndim)] 123 124 # Initial list to contain the sum of the products of the probability for each m_r and m_s sample in their respective models. 125 pdfProductSumList = [] 126 127 # For the m_r and m_s pairs in each BNS system. (eg. 1000x2) 128 for massSample in massSamples: 129 130 # Evaluate the PDF function down the m_r and m_s samples of the BNS 131 mrProbabilities = modelEval1(paramList1, massSample[:,0]) 132 msProbabilities = modelEval2(paramList2, massSample[:,1]) 133 134 # Evaluate the product of the m_r and m_s probability for each pair. 135 probabilityProduct = mrProbabilities*msProbabilities 136 137 # Append the sum over all the probability products of each pair. 138 pdfProductSumList.append(np.sum(probabilityProduct)) 139 140 # If either all the m_r or all the m_s samples are completely outside their model then return a log-likelihood of -inf. 141 if 0 in pdfProductSumList: 142 print("Zero probability value - Parameters: {}, {}".format(paramList1,paramList2)) 143 return -np.inf 144 145 # The log-likelihood is the log of the normalised sum over the log of each pdfProductSum 146 loglikelihood = nMeasurements * np.log(1.0/nSamples) + np.sum(np.log(pdfProductSumList)) 147 return loglikelihood 148 149 150 151 ### Models 152 # This list contains the set of all model dictionary combinations 153 # Model names, pdf functions, nDimensions, ParameterNames 154 modelSet = [[singleGaussianModel, singleGaussianModel], 155 [singleGaussianModel, twoGaussianModel], 156 [singleGaussianModel, uniformModel], 157 [twoGaussianModel, singleGaussianModel], 158 [twoGaussianModel, twoGaussianModel], 159 [twoGaussianModel, uniformModel], 160 [uniformModel, singleGaussianModel], 161 [uniformModel, twoGaussianModel], 162 [uniformModel, uniformModel]] 163 164 165 ### System input (models choice) 166 # Takes the system argument given. 167 # eg hypothesisB.py 2 will select the model pair with index 1, singleGuassian and twoGaussian. 168 modelSetIndex = int(sys.argv[1]) - 1 169 170 # Set relevant variables which are to be used in this sampling 171 modelDictionary1, modelDictionary2 = modelSet[modelSetIndex] 172 modelName1, modelEval1, ndim1, paramNames1 = modelDictionary1['name'], modelDictionary1['pdf'], modelDictionary1['ndim'], modelDictionary1['params'] 173 modelName2, modelEval2, ndim2, paramNames2 = modelDictionary2['name'], modelDictionary2['pdf'], modelDictionary2['ndim'], modelDictionary2['params'] 174 175 # Combined parameter list 176 paramNames = paramNames1 + paramNames2 177 178 # Define the total number of parameters 179 ndim = ndim1 + ndim2 180 181 182 183 ### Inference 184 # Directory to send output to. Create it if it does not exist. 185 directoryName = topDirectory + 'hypB_out/' + modelName1[:4] + "/" + modelName2[:4] 186 if not os.path.exists(directoryName): 187 os.makedirs(directoryName) 188 189 # Run pymultinest sampling 190 pymultinest.run(likelihood, prior, ndim, n_live_points=1000, sampling_efficiency=0.3, importance_nested_sampling=False, outputfiles_basename=directoryName + '/', resume=False, verbose=True) 191 192 print("Code took {} seconds to run.".format(datetime.now() - startTime)) 193 print("Exiting!") 194