diff --git a/data/bao_2014.txt b/data/bao_2014.txt index bd879e5a..0dec83bc 100644 --- a/data/bao_2014.txt +++ b/data/bao_2014.txt @@ -22,10 +22,10 @@ LOWZ 0.32 8.47 0.17 3 CMASS 0.57 13.77 0.13 3 # # BOSS LyaQSO (DA/rs) Font-Ribera et al. 1311.1767 -# LyaQSO 2.36 10.8 0.4 5 + LyaQSO 2.36 10.8 0.4 5 # # BOSS LyaQSO c/(H rs) Font-Ribera et al. 1311.1767 -# LyaQSO 2.36 9.0 0.3 6 + LyaQSO 2.36 9.0 0.3 6 # # SDSS DR7 MGS, Ross et al. 1409.3242v1 MGS 0.15 4.47 0.16 3 diff --git a/montepython/analyze.py b/montepython/analyze.py index 8da81f05..f710b05c 100644 --- a/montepython/analyze.py +++ b/montepython/analyze.py @@ -832,6 +832,23 @@ def compute_posterior(information_instances): "'%s-%s' 2d-plot" % ( info.plotted_parameters[info.native_index], info.plotted_parameters[info.native_second_index])) + except ValueError as e: + if str(e) == "Contour levels must be increasing": + warnings.warn( + "The routine could not find the contour of the " + + "'%s-%s' 2d-plot. \n " % ( + info.plotted_parameters[info.native_index], + info.plotted_parameters[info.native_second_index]) + + 'The error is: "Contour levels must be increasing"' + + " but " + str(ctr_level(info.n, info.levels[:2])) + + " were found. This may happen when most" + + " points fall in the same bin.") + else: + warnings.warn( + "The routine could not find the contour of the " + + "'%s-%s' 2d-plot" % ( + info.plotted_parameters[info.native_index], + info.plotted_parameters[info.native_second_index])) ax2dsub.set_xticks(info.ticks[info.native_second_index]) ax2dsub.set_yticks(info.ticks[info.native_index]) @@ -1382,6 +1399,10 @@ def recover_folder_and_files(files): # The following list defines the substring that a chain should contain for # the code to recognise it as a proper chain. substrings = ['.txt', '__'] + # The following variable defines the substring that identify error_log + # files and therefore there must not be taken into account in the analysis. + substring_err = 'error_log' + limit = 10 # If the first element is a folder, grab all chain files inside if os.path.isdir(files[0]): @@ -1389,6 +1410,7 @@ def recover_folder_and_files(files): files = [os.path.join(folder, elem) for elem in os.listdir(folder) if not os.path.isdir(os.path.join(folder, elem)) and not os.path.getsize(os.path.join(folder, elem)) < limit + and (substring_err not in elem) and all([x in elem for x in substrings])] # Otherwise, extract the folder from the chain file-name. else: @@ -1402,6 +1424,7 @@ def recover_folder_and_files(files): if os.path.join(folder, elem) in np.copy(files) and not os.path.isdir(os.path.join(folder, elem)) and not os.path.getsize(os.path.join(folder, elem)) < limit + and (substring_err not in elem) and all([x in elem for x in substrings])] basename = os.path.basename(folder) return folder, files, basename diff --git a/montepython/cosmo_hammer.py b/montepython/cosmo_hammer.py index efec34d4..a384f71e 100644 --- a/montepython/cosmo_hammer.py +++ b/montepython/cosmo_hammer.py @@ -22,9 +22,11 @@ import io_mp import sampler -from cosmoHammer.likelihood.chain.LikelihoodComputationChain import ( - LikelihoodComputationChain) -from cosmoHammer.sampler.CosmoHammerSampler import CosmoHammerSampler +from cosmoHammer import MpiCosmoHammerSampler +from cosmoHammer import LikelihoodComputationChain +#from cosmoHammer.likelihood.chain.LikelihoodComputationChain import ( +# LikelihoodComputationChain) +#from cosmoHammer.sampler.CosmoHammerSampler import CosmoHammerSampler from cosmoHammer.util.SampleFileUtil import SampleFileUtil # Cosmo Hammer subfolder and name separator @@ -33,6 +35,7 @@ # Cosmo Hammer file names ending, after the defined 'base_name' name_arguments = '.arguments' name_chain = 'chain_CH__sampling.txt' +name_error_chain = 'chain_CH__sampling-error_log.txt' # Cosmo Hammer option prefix CH_prefix = 'CH_' @@ -48,6 +51,10 @@ {'help': 'Number of sample iterations', 'type': int}} +# List of strings that will contain the derived parameters names. +# Filled in run using the MontePython class data +# Used in DerivedUtil.persistValues +CH_derived_names = [] def run(cosmo, data, command_line): """ @@ -57,6 +64,9 @@ def run(cosmo, data, command_line): # Store the parameters inside the format expected by CosmoHammer # TODO: about the derived params? parameter_names = data.get_mcmc_parameters(["varying"]) + # Fill now CH_derived_names to be used in DerivedUtil.persistValues + global CH_derived_names + CH_derived_names= data.get_mcmc_parameters(["derived"]) # Ensure that their prior is bound and flat is_flat, is_bound = sampler.check_flat_bound_priors( @@ -91,6 +101,7 @@ def run(cosmo, data, command_line): # here, since data must be called before cosmo. chain.addCoreModule(data) chain.addCoreModule(cosmo) + chain.addCoreModule(store_cosmo_derived) # Add each likelihood class as a LikelihoodModule for likelihood in data.lkl.itervalues(): @@ -101,7 +112,9 @@ def run(cosmo, data, command_line): file_prefix = os.path.join(command_line.folder, CH_subfolder, chain_name) # Recover the User options - data.CH_arguments = {} + data.CH_arguments = {'walkersRatio': 50, + 'burninIterations': 10, + 'sampleIterations': 30} for arg in CH_user_arguments: value = getattr(command_line, CH_prefix+arg) if value != -1: @@ -127,13 +140,10 @@ def run(cosmo, data, command_line): num_threads = 1 # Create the Sampler object - sampler_hammer = CosmoHammerSampler( + sampler_hammer = MpiCosmoHammerSampler( params=params, likelihoodComputationChain=chain, filePrefix=file_prefix, - walkersRatio=50, - burninIterations=10, - sampleIterations=30, storageUtil=derived_util, threadCount=num_threads, **data.CH_arguments) @@ -158,6 +168,8 @@ def from_CH_output_to_chains(folder): This function will be called by the module :mod:`analyze`. """ + chainsClean, lklClean = [], [] # Arrays to store the valid steps. + error = [] # Array to store the step did not could not be computed. chain_name = [a for a in folder.split(os.path.sep) if a][-2] base_name = os.path.join(folder, chain_name) @@ -172,17 +184,25 @@ def from_CH_output_to_chains(folder): # multiplicity. This does not mean that the acceptance rate is one, but # points where the code stayed are duplicated in the file. - ## First, reshape the lkl array - lkl = np.array([[elem] for elem in lkl]) + # First, separe errors from accepted points: + for x, p in zip(chains, lkl): + if p != -np.inf: + chainsClean.append(x) + lklClean.append([-p]) # MP stores -loglkl in its files. + else: + error.append([val for val in x if not np.isnan(val)]) # Remove nan's ## Create the array of ones - ones = np.array([[1] for _ in range(len(lkl))]) + ones = np.array([[1] for _ in range(len(lklClean))]) + ## Concatenate everything and save to file - final = np.concatenate((ones, lkl, chains), axis=1) + final = np.concatenate((ones, lklClean, chainsClean), axis=1) output_folder = os.path.join(folder, '..') output_chain_path = os.path.join(output_folder, name_chain) + output_error_chain_path = os.path.join(output_folder, name_error_chain) np.savetxt(output_chain_path, final) + np.savetxt(output_error_chain_path, error) class DerivedUtil(SampleFileUtil): @@ -195,8 +215,12 @@ def persistValues(self, posFile, probFile, pos, prob, data): """ # extend the pos array to also contain the value of the derived # parameters + + derived_not_computed = [np.NaN] * len(CH_derived_names) + derived = np.array( - [[a for a in elem.itervalues()] for elem in data]) + #[[a for a in elem.itervalues()] if elem else derived_not_computed for elem in data]) + [[elem[key] for key in sorted(elem.iterkeys())] if elem else derived_not_computed for elem in data]) final = np.concatenate((pos, derived), axis=1) posFile.write("\n".join( @@ -207,3 +231,31 @@ def persistValues(self, posFile, probFile, pos, prob, data): probFile.write("\n".join([str(p) for p in prob])) probFile.write("\n") probFile.flush() + + +def store_cosmo_derived(ctx): + """ + Store derived parameters. Based on sampler.compute_lkl + """ + from classy import CosmoSevereError + + cosmo = ctx.get("cosmo") + data = ctx.get("data") + + if data.get_mcmc_parameters(['derived']) != []: + try: + derived = cosmo.get_current_derived_parameters( + data.get_mcmc_parameters(['derived'])) + except AttributeError: + # This happens if the classy wrapper is still using the old + # convention, expecting data as the input parameter + cosmo.get_current_derived_parameters(data) + except CosmoSevereError: + raise io_mp.CosmologicalModuleError( + "Could not write the current derived parameters") + + ctx.add('key_data', derived) + + # If we want the derived parameters stored rescaled + # derived_scaled = {elem: value/data.mcmc_parameters[elem]['scale'] for elem, value in derived.iteritems()} + # ctx.add('key_data', derived_scaled) diff --git a/montepython/data.py b/montepython/data.py index 05382341..9c6975c8 100644 --- a/montepython/data.py +++ b/montepython/data.py @@ -189,6 +189,11 @@ def __init__(self, command_line, path): self.out = None self.out_name = '' + # Create the variable err, and err_name to keep the errors + # will be initialised later by the :mod:`io_mp` module + self.err = None + self.err_name = '' + # If the parameter file is not a log.param, the path will be read # before reading the parameter file. if self.param.find('log.param') == -1: @@ -995,8 +1000,13 @@ def __init__(self, array, key): dict.__init__(self) self['initial'] = array[0:4] - self['scale'] = array[4] - self['role'] = array[-1] + if array[4]: + self['scale'] = array[4] + else: + raise io_mp.ParameterError("You chose scale 0 to param {}. Correct\ + it.".format(key)) + + self['role'] = array[5] self['tex_name'] = io_mp.get_tex_name(key) if array[3] == 0: self['status'] = 'fixed' diff --git a/montepython/initialise.py b/montepython/initialise.py index ca09f782..009fe172 100644 --- a/montepython/initialise.py +++ b/montepython/initialise.py @@ -84,6 +84,7 @@ def initialise(custom_command=''): io_mp.create_output_files(command_line, data) # NS: Creating the NS subfolder and the MultiNest arguments elif command_line.method == 'NS': + io_mp.create_output_files(command_line, data, only_error=True) from nested_sampling import initialise as initialise_ns initialise_ns(cosmo, data, command_line) diff --git a/montepython/io_mp.py b/montepython/io_mp.py old mode 100644 new mode 100755 index 7eecba6c..f7f3f8c7 --- a/montepython/io_mp.py +++ b/montepython/io_mp.py @@ -211,6 +211,69 @@ def print_vector(out, N, loglkl, data): out[j].write('\n') +def print_error(error_msg, data): + """ + Register parameters leading to class_error to :code:`out` + + Parameters + ---------- + out : list + Array containing both standard output and the output file. + + This way, if you run in interactive mode, you will be able to monitor + the progress of the chain. + error_msg : string + Error message: the idea is not to keep the whole error message, + but to read part of it and asign it an error tag + 0 => error not identified + n>0 => different error codes + + .. note:: + + We use error codes instead of recording the full error messages + unless the error is not identified + + It is the `current` point that is printed, which produced the error + No derived parameters are written (maybe not produced) + + """ + + # -1 is for unmarked error + err_code = -1 + # rewrite with different possibilities + # errors in background are 1-10, perturbations 10-99... + if 'Shooting failed' in error_msg: + err_code = 0 + if all(x in error_msg for x in ['Gradient', 'scalar']): + err_code = 1 + if all(x in error_msg for x in ['Gradient', 'tensor']): + err_code = 2 + if all(x in error_msg for x in ['Ghost', 'scalar']): + err_code = 3 + if all(x in error_msg for x in ['Ghost', 'tensor']): + err_code = 4 + + if "Isnan v_X" in error_msg: + err_code = 10 + + if 'Isnan v_X' in error_msg: + err_code = 10 + if 'early_smg' and 'adiabatic' in error_msg: + err_code = 11 + + out = data.err + for elem in data.get_mcmc_parameters(['varying']): + out.write('%.6e\t' % + data.mcmc_parameters[elem]['current']) + # if error unidentified write full message + if err_code == -1: + out.write(str(err_code) + "\t" + error_msg) + else: + #print "error",err_code + out.write(str(err_code)) + out.write('\n') + + def refresh_file(data): """ Closes and reopen the output file to write any buffered quantities @@ -219,8 +282,11 @@ def refresh_file(data): data.out.close() data.out = open(data.out_name, 'a') + data.err.close() + data.err = open(data.err_name, 'a') -def create_output_files(command_line, data): + +def create_output_files(command_line, data, only_error=False): """ Automatically create a new name for the chain. @@ -236,43 +302,65 @@ def create_output_files(command_line, data): changing things here. """ - if command_line.restart is None: - number = command_line.N - else: - number = int( - command_line.restart.split(os.path.sep)[-1].split('__')[0]. - split('_')[1]) + command_line.N - - # output file - outname_base = '{0}_{1}__'.format(date.today(), number) - suffix = 0 - trying = True - if command_line.chain_number is None: - for files in os.listdir(command_line.folder): - if files.find(outname_base) != -1: - if int(files.split('__')[-1].split('.')[0]) > suffix: - suffix = int(files.split('__')[-1].split('.')[0]) - suffix += 1 - while trying: - data.out = open(os.path.join( - command_line.folder, outname_base)+str(suffix)+'.txt', 'w') - try: - lock(data.out, fcntl.LOCK_EX | fcntl.LOCK_NB) - trying = False - except LockError: + + for f in ('out', 'err'): + if (f == 'out') and only_error: + continue + + if command_line.restart is None: + number = command_line.N + else: + number = int( + command_line.restart.split(os.path.sep)[-1].split('__')[0]. + split('_')[1]) + command_line.N + + # output file + outname_base = '{0}_{1}__'.format(date.today(), number) + suffix = 0 + suffix_err = "" + if f == 'err': + suffix_err = "-error_log" + + trying = True + if command_line.chain_number is None: + for files in os.listdir(command_line.folder): + if files.find(outname_base) != -1: + matched_number = re.search('\d*', files.split('__')[-1]) + number = int(matched_number.group(0)) + if number > suffix: + suffix = number + + if (f == 'out') or only_error: suffix += 1 - data.out_name = os.path.join( - command_line.folder, outname_base)+str(suffix)+'.txt' - print 'Creating %s\n' % data.out_name - else: - data.out_name = os.path.join( - command_line.folder, outname_base)+command_line.chain_number+'.txt' - data.out = open(data.out_name, 'w') - print 'Creating %s\n' % data.out_name - # in case of a restart, copying the whole thing in the new file - if command_line.restart is not None: - for line in open(command_line.restart, 'r'): - data.out.write(line) + + while trying: + open_file = open(os.path.join( + command_line.folder, outname_base)+str(suffix)+suffix_err+'.txt', 'w') + + try: + lock(open_file, fcntl.LOCK_EX | fcntl.LOCK_NB) + trying = False + except LockError: + suffix += 1 + name = os.path.join( + command_line.folder, outname_base)+str(suffix)+ suffix_err+'.txt' + print 'Creating %s\n' % name + else: + name = os.path.join( + command_line.folder, outname_base)+command_line.chain_number+suffix_err+'.txt' + open_file = open(name, 'w') + print 'Creating %s\n' % name + # in case of a restart, copying the whole thing in the new file + if command_line.restart is not None: + for line in open(command_line.restart[:-4]+suffix_err+'.txt', 'r'): + open_file.write(line) + + if f == 'out': + data.out = open_file + data.out_name = name + else: + data.err = open_file + data.err_name = name def get_tex_name(name, number=1): @@ -554,3 +642,8 @@ class FiducialModelWritten(MyError): class AnalyzeError(MyError): """Used when encountering a fatal mistake in analyzing chains""" pass + + +class ParameterError(MyError): + """Used when some MontePython parameter is not correctly set""" + pass diff --git a/montepython/likelihoods/Planck_compressed/Planck_compressed.data b/montepython/likelihoods/Planck_compressed/Planck_compressed.data new file mode 100644 index 00000000..9d6bb8b2 --- /dev/null +++ b/montepython/likelihoods/Planck_compressed/Planck_compressed.data @@ -0,0 +1,23 @@ +# Covriance matrix of the compressed likelihood for Planck TT + lowP +# (arXiv:1502.01590v2, Table 4) + # R, l_A, omega_b, n_s +#Planck_compressed.correlation_matrix = [[1., 0.54, -0.63, -0.86],\ # R +# [0.54, 1., -0.43, -0.48],\ # l_A +# [-0.63, -0.43, 1., 0.58],\ # omega_b +# [-0.86, -0.48, 0.58, 1.]] # n_s + +# Covariance matrix_ij = sigma_i sigma_j Correlation_matrix_ij + +Planck_compressed.inverse_correlation_matrix = [[ 4.51266537, -0.59697118, 0.75581806, 3.15597158], [-0.59697118, 1.43957767, 0.21084593, 0.05531143], [ 0.75581806, 0.21084593, 1.70453221, -0.2374191 ], [ 3.15597158, 0.05531143, -0.2374191 , 3.87838813]] + + # mean, std. dev. (sigma) +#Planck_compressed.R = [1.7488, 0.0074] +#Planck_compressed.l_a = [301.76, 0.14] +#Planck_compressed.omega_b = [0.02228, 0.00023] +#Planck_compressed.n_s = [0.9660, 0.0061] + +#Planck_compressed.planck_values = [[1.7488, 0.0074], [301.76, 0.14], [0.02228, 0.00023], [0.9660, 0.0061]] + +Planck_compressed.planck_values_mean = [1.7488, 301.76, 0.02228, 0.9660] +Planck_compressed.planck_values_sigma = [0.0074, 0.14, 0.00023, 0.0061] + diff --git a/montepython/likelihoods/Planck_compressed/__init__.py b/montepython/likelihoods/Planck_compressed/__init__.py new file mode 100644 index 00000000..0390461b --- /dev/null +++ b/montepython/likelihoods/Planck_compressed/__init__.py @@ -0,0 +1,42 @@ +from montepython.likelihood_class import Likelihood +import numpy as np + + +class Planck_compressed(Likelihood): + + # initialisation of the class is done within the parent Likelihood_prior. For + # this case, it does not differ, actually, from the __init__ method in + # Likelihood class. + + def loglkl(self, cosmo, data): + + derived_params = cosmo.get_current_derived_parameters(['z_rec', 'rs_rec']) + + #z_star = cosmo.z_optical_depth_unity() # z at which the optical depth is unity + z_star = derived_params['z_rec'] # z at which the optical depth is unity + # comoving angular diameter distance + D_A_star = cosmo.angular_distance(z_star) * (z_star + 1) + #rs_star = cosmo.rs_at_z(z_star) + rs_star = derived_params['rs_rec'] + + R = np.sqrt(cosmo.Omega_m()) * cosmo.Hubble(0) * D_A_star + l_a = np.pi * D_A_star / rs_star + omega_b = cosmo.omega_b() + n_s = cosmo.n_s() + + # Using self in order to be accesible by data_mp_likelihoods.py + # work-in-class module. + # https://github.com/miguelzuma/work-in-class && + # https://github.com/ardok-m/work-in-class && + self.values = np.array([R, l_a, omega_b, n_s]) + + # loglkl = -1/2 * (x_i - x_mean_i) * 1/sigma_i * corr_mat_ij^-1 * + # 1/sigma_j * (x_j - x_mean_j) + + differences_over_sigma = (self.values - np.array(self.planck_values_mean)) / self.planck_values_sigma + + loglkl = -0.5 * np.dot(differences_over_sigma, + np.dot(self.inverse_correlation_matrix, + differences_over_sigma)) + + return loglkl diff --git a/montepython/likelihoods/bao_boss/__init__.py b/montepython/likelihoods/bao_boss/__init__.py index faf5dd81..76028d79 100644 --- a/montepython/likelihoods/bao_boss/__init__.py +++ b/montepython/likelihoods/bao_boss/__init__.py @@ -60,6 +60,9 @@ def loglkl(self, cosmo, data): chi2 = 0. + # This object is meant to be used with data_mp_likelihoods.py + self.theo = [] + # for each point, compute angular distance da, radial distance dr, # volume distance dv, sound horizon at baryon drag rs_d, # theoretical prediction and chi2 contribution @@ -92,6 +95,8 @@ def loglkl(self, cosmo, data): chi2 += ((theo - self.data[i]) / self.error[i]) ** 2 + self.theo.append(theo) + # return ln(L) lkl = - 0.5 * chi2 diff --git a/montepython/mcmc.py b/montepython/mcmc.py index 9df8b1e8..6644c54b 100644 --- a/montepython/mcmc.py +++ b/montepython/mcmc.py @@ -452,7 +452,7 @@ def chain(cosmo, data, command_line): # just computed ('current' flag), but really the previous one.) # with its proper multiplicity (number of times the system stayed # there). - io_mp.print_vector(outputs, N, loglike, data) + io_mp.print_vector(outputs, N, loglike, data) # Report the 'current' point to the 'last_accepted' sampler.accept_step(data) diff --git a/montepython/prior.py b/montepython/prior.py index f86db5a6..12e15a19 100644 --- a/montepython/prior.py +++ b/montepython/prior.py @@ -50,11 +50,25 @@ def __init__(self, array): "You asked for a gaussian prior, but provided no " + "mean nor sigma. Please add them in the parameter " + "file.") + # in case of a triangular prior, one expects one more entry, the + # position of the mode + elif self.prior_type == 'triangular': + try: + self.mode = array[7] + # TODO: Consider unifying prior variable names. + except IndexError: + raise io_mp.ConfigurationError( + "You asked for a triangular prior, but provided no " + + "position for the mode. Please add it in the parameter " + + "file.") # Store boundaries for convenient access later # Put all fields that are -1 to None to avoid confusion later on. self.prior_range = [a if not((a is -1) or (a is None)) else None for a in deepcopy(array[1:3])] + # TODO: As implemented above, you wouldn't be able to set -1 as + # boundary for one of the parameters. Consider leaving the None part + # only def draw_from_prior(self): """ @@ -76,10 +90,13 @@ def draw_from_prior(self): while not within_bounds: value = rd.gauss(self.mu, self.sigma) # Check for boundaries problem - within_bounds = calue_within_prior_range(value) + within_bounds = self.value_within_prior_range(value) return value - + + elif self.prior_type == 'triangular': + return rd.triangular(self.prior_range[0], self.prior_range[1], self.mode) + def value_within_prior_range(self, value): """ Check for a value being in or outside the prior range @@ -112,5 +129,5 @@ def map_from_unit_interval(self, value): which should have been previously checked with :func:`is_bound` """ - return (self.prior_range[0] + + return (self.prior_range[0] + value * (self.prior_range[1] - self.prior_range[0])) diff --git a/montepython/run.py b/montepython/run.py index c74a8068..3787f9de 100644 --- a/montepython/run.py +++ b/montepython/run.py @@ -115,6 +115,8 @@ def mpi_run(custom_command=""): status = suffix elif command_line.method == "NS": status = 1 + elif command_line.method == "CH": + status = 1 else: warnings.warn( "The method '%s' is not supported"%(command_line.method) + diff --git a/montepython/sampler.py b/montepython/sampler.py index 5caeca99..c9276b3e 100644 --- a/montepython/sampler.py +++ b/montepython/sampler.py @@ -489,10 +489,14 @@ def compute_lkl(cosmo, data): # point, but with minimum likelihood, so will be rejected, resulting in # the choice of a new point. try: - cosmo.compute(["lensing"]) + cosmo.compute(["lensing"]) except CosmoComputationError as failure_message: sys.stderr.write(str(failure_message)+'\n') sys.stderr.flush() + + #write the error file + io_mp.print_error(str(failure_message), data) + return data.boundary_loglike except CosmoSevereError as critical_message: raise io_mp.CosmologicalModuleError(