From e6d6bc2b2e08744664d284849c5f971b24cd78fe Mon Sep 17 00:00:00 2001 From: Miguel Zumalacarregui Date: Tue, 22 Mar 2016 13:31:30 -0700 Subject: [PATCH 01/35] fixed typo in analyze.py --- montepython/analyze.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/montepython/analyze.py b/montepython/analyze.py index c302bbdc..abe5806c 100644 --- a/montepython/analyze.py +++ b/montepython/analyze.py @@ -683,7 +683,7 @@ def compute_posterior(information_instances): # store the coordinates of the points for further # plotting. store_contour_coordinates( - conf, standard_name, second_standard_name, cotours) + conf, standard_name, second_standard_name, contours) for info in information_instances: if not info.ignore_param and info.has_second_param: From 7d72cfc8ccb48436df8930883bf7681e9e1b75dd Mon Sep 17 00:00:00 2001 From: Miguel Zumalacarregui Date: Wed, 12 Oct 2016 21:38:35 -0700 Subject: [PATCH 02/35] writes error_log for failed models --- montepython/data.py | 5 +++++ montepython/io_mp.py | 40 ++++++++++++++++++++++++++++++++++++++++ montepython/mcmc.py | 2 +- montepython/sampler.py | 6 +++++- 4 files changed, 51 insertions(+), 2 deletions(-) diff --git a/montepython/data.py b/montepython/data.py index 05382341..714e62c7 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: diff --git a/montepython/io_mp.py b/montepython/io_mp.py index 7eecba6c..82bce9fa 100644 --- a/montepython/io_mp.py +++ b/montepython/io_mp.py @@ -211,6 +211,36 @@ 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:: + + It is the `current` point that is printed, which produced the error + No derived parameters are written (maybe not produced) + """ + out = data.err + for elem in data.get_mcmc_parameters(['varying']): + out.write('%.6e\t' % + data.mcmc_parameters[elem]['current']) + out.write(error_msg) + out.write('\n') + + def refresh_file(data): """ Closes and reopen the output file to write any buffered quantities @@ -218,6 +248,9 @@ 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): @@ -273,6 +306,13 @@ def create_output_files(command_line, data): if command_line.restart is not None: for line in open(command_line.restart, 'r'): data.out.write(line) + + #open/create error file + data.err_name = os.path.join(command_line.folder,"error_log.txt") + data.err = open(data.err_name,'a') + print 'Creating error file %s\n' % data.err_name + + def get_tex_name(name, number=1): diff --git a/montepython/mcmc.py b/montepython/mcmc.py index b1ef1af7..84227d72 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/sampler.py b/montepython/sampler.py index e29c8aba..9318c117 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( From 7ed9391b5c2ca8908a38296ab45a0ff78964e6a4 Mon Sep 17 00:00:00 2001 From: Miguel Zumalacarregui Date: Thu, 13 Oct 2016 18:13:52 -0700 Subject: [PATCH 03/35] added error codes to record error_log --- montepython/io_mp.py | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/montepython/io_mp.py b/montepython/io_mp.py index 82bce9fa..65705b42 100644 --- a/montepython/io_mp.py +++ b/montepython/io_mp.py @@ -229,15 +229,43 @@ def print_error(error_msg, data): 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 'Gradient instability' and 'scalar' in error_msg: + err_code = 1 + if 'Gradient instability' and 'tensor' in error_msg: + err_code = 2 + if 'Ghost instability' and 'scalar' in error_msg: + err_code = 3 + if 'Ghost instability' and 'tensor' in error_msg: + err_code = 4 + + if "Isnan v_X" in error_msg: + err_code = 10 + + out = data.err for elem in data.get_mcmc_parameters(['varying']): out.write('%.6e\t' % data.mcmc_parameters[elem]['current']) - out.write(error_msg) + # if error unidentified write full message + if err_code == -1: + out.write(str(err_code) + "\t" + error_msg) + else: + out.write(str(err_code)) out.write('\n') From 5c493d780a514530ae3dceee053905c98683e934 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Tue, 6 Dec 2016 12:35:08 +0100 Subject: [PATCH 04/35] fixed conditionals for err_code selection --- montepython/io_mp.py | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/montepython/io_mp.py b/montepython/io_mp.py index 65705b42..a2fd7692 100644 --- a/montepython/io_mp.py +++ b/montepython/io_mp.py @@ -223,51 +223,51 @@ def print_error(error_msg, data): 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, + 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 + # rewrite with different possibilities # errors in background are 1-10, perturbations 10-99... if 'Shooting failed' in error_msg: err_code = 0 - if 'Gradient instability' and 'scalar' in error_msg: + if all(x in error_msg for x in ['Gradient instability', 'scalar']): err_code = 1 - if 'Gradient instability' and 'tensor' in error_msg: + if all(x in error_msg for x in ['Gradient instability', 'tensor']): err_code = 2 - if 'Ghost instability' and 'scalar' in error_msg: + if all(x in error_msg for x in ['Ghost instability', 'scalar']): err_code = 3 - if 'Ghost instability' and 'tensor' in error_msg: + if all(x in error_msg for x in ['Ghost instability', 'tensor']): err_code = 4 - + if "Isnan v_X" in error_msg: err_code = 10 - - + + 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 error unidentified write full message if err_code == -1: out.write(str(err_code) + "\t" + error_msg) else: out.write(str(err_code)) out.write('\n') - + def refresh_file(data): """ @@ -276,7 +276,7 @@ def refresh_file(data): """ data.out.close() data.out = open(data.out_name, 'a') - + data.err.close() data.err = open(data.err_name, 'a') @@ -334,13 +334,13 @@ def create_output_files(command_line, data): if command_line.restart is not None: for line in open(command_line.restart, 'r'): data.out.write(line) - - #open/create error file + + #open/create error file data.err_name = os.path.join(command_line.folder,"error_log.txt") data.err = open(data.err_name,'a') print 'Creating error file %s\n' % data.err_name - - + + def get_tex_name(name, number=1): From 7211aa3ceab66871b550605e603d76e1bc1cfb96 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Tue, 6 Dec 2016 12:35:36 +0100 Subject: [PATCH 05/35] fixed missing __init__.py file --- montepython/likelihoods/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 montepython/likelihoods/__init__.py diff --git a/montepython/likelihoods/__init__.py b/montepython/likelihoods/__init__.py new file mode 100644 index 00000000..e69de29b From 30286564521c40362371a1d1f38244a68d9229c8 Mon Sep 17 00:00:00 2001 From: Miguel Zumalacarregui Date: Tue, 6 Dec 2016 10:11:12 -0800 Subject: [PATCH 06/35] changes in print_error --- montepython/io_mp.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/montepython/io_mp.py b/montepython/io_mp.py index 65705b42..fab600f4 100644 --- a/montepython/io_mp.py +++ b/montepython/io_mp.py @@ -244,17 +244,19 @@ def print_error(error_msg, data): # errors in background are 1-10, perturbations 10-99... if 'Shooting failed' in error_msg: err_code = 0 - if 'Gradient instability' and 'scalar' in error_msg: + if 'Gradient' and 'scalar' in error_msg: err_code = 1 - if 'Gradient instability' and 'tensor' in error_msg: + if 'Gradient' and 'tensor' in error_msg: err_code = 2 - if 'Ghost instability' and 'scalar' in error_msg: + if 'Ghost' and 'scalar' in error_msg: err_code = 3 - if 'Ghost instability' and 'tensor' in error_msg: + if 'Ghost' and 'tensor' in error_msg: err_code = 4 - if "Isnan v_X" in error_msg: + if 'Isnan v_X' in error_msg: err_code = 10 + if 'early_smg' and 'adiabatic' in error_msg: + err_code = 11 out = data.err @@ -265,6 +267,7 @@ def print_error(error_msg, data): 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') From e5498ec9eb7200348fe1eab48f803428eadeb22f Mon Sep 17 00:00:00 2001 From: Miguel Zumalacarregui Date: Tue, 6 Dec 2016 13:07:09 -0800 Subject: [PATCH 07/35] commented out CMASS data (use bao_boss_aniso instead) --- data/bao_2014.txt | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/data/bao_2014.txt b/data/bao_2014.txt index bd879e5a..5a2072e2 100644 --- a/data/bao_2014.txt +++ b/data/bao_2014.txt @@ -19,13 +19,14 @@ LOWZ 0.32 8.47 0.17 3 # # BOSS CMASS DR10&11 (D_V in Mpc divided by r_s^fid in Mpc) Anderson et al. 1312.4877 # found in Section 8.2 (directly using the exact sound horizon) -CMASS 0.57 13.77 0.13 3 +#CMASS 0.57 13.77 0.13 3 +#NOTE: commented out because of using anisotropic BOSS-CMASS # # 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 From 15652943a9e6c7ba3054bf2731624a7cad7c46fa Mon Sep 17 00:00:00 2001 From: Miguel Zumalacarregui Date: Tue, 6 Dec 2016 18:07:49 -0800 Subject: [PATCH 08/35] uncommented data: exclusion was made automatic --- data/bao_2014.txt | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/data/bao_2014.txt b/data/bao_2014.txt index 5a2072e2..0dec83bc 100644 --- a/data/bao_2014.txt +++ b/data/bao_2014.txt @@ -19,8 +19,7 @@ LOWZ 0.32 8.47 0.17 3 # # BOSS CMASS DR10&11 (D_V in Mpc divided by r_s^fid in Mpc) Anderson et al. 1312.4877 # found in Section 8.2 (directly using the exact sound horizon) -#CMASS 0.57 13.77 0.13 3 -#NOTE: commented out because of using anisotropic BOSS-CMASS +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 From 17cf1ec710180f2079f73bd61c4d5a0386b86e45 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Wed, 14 Dec 2016 18:52:05 +0100 Subject: [PATCH 09/35] improved input param handling: error with param scale 0 --- montepython/data.py | 7 ++++++- montepython/io_mp.py | 5 +++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/montepython/data.py b/montepython/data.py index 714e62c7..158ae434 100644 --- a/montepython/data.py +++ b/montepython/data.py @@ -1000,7 +1000,12 @@ def __init__(self, array, key): dict.__init__(self) self['initial'] = array[0:4] - self['scale'] = array[4] + 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[-1] self['tex_name'] = io_mp.get_tex_name(key) if array[3] == 0: diff --git a/montepython/io_mp.py b/montepython/io_mp.py index a2fd7692..998a5ef5 100644 --- a/montepython/io_mp.py +++ b/montepython/io_mp.py @@ -622,3 +622,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 From a69f5bc94cdb237c57ef4789b76b6f375368e40f Mon Sep 17 00:00:00 2001 From: ardok-m Date: Mon, 17 Apr 2017 17:12:44 +0200 Subject: [PATCH 10/35] fixed TypeError np.max for numpy 1.12 --- montepython/analyze.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/montepython/analyze.py b/montepython/analyze.py index a7c6df0a..2402875c 100644 --- a/montepython/analyze.py +++ b/montepython/analyze.py @@ -788,8 +788,8 @@ def minimum_credible_intervals(info): bounds = np.zeros((len(levels), 2)) j = 0 delta = bincenters[1]-bincenters[0] - left_edge = np.max(histogram[0] - 0.5*(histogram[1]-histogram[0]), 0.) - right_edge = np.max(histogram[-1] + 0.5*(histogram[-1]-histogram[-2]), 0.) + left_edge = np.max(histogram[0] - 0.5*(histogram[1]-histogram[0]), 0) + right_edge = np.max(histogram[-1] + 0.5*(histogram[-1]-histogram[-2]), 0) failed = False for level in levels: norm = float( From 9674a286acade1cd936d6b2b70b6e2f2cbafff85 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Tue, 18 Apr 2017 10:57:04 +0200 Subject: [PATCH 11/35] fixed contours level error handling. Matplotlib must have promoted it from a Warning to a Value Error --- montepython/analyze.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/montepython/analyze.py b/montepython/analyze.py index 2402875c..26ce89f0 100644 --- a/montepython/analyze.py +++ b/montepython/analyze.py @@ -672,6 +672,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]) if index == len(plotted_parameters)-1: From aa0498499cad2e5653d09a09845690970d7089e6 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Mon, 22 May 2017 11:24:10 +0200 Subject: [PATCH 12/35] added Planck_compressed likelihood arXiv:1502.01590 --- .../Planck_compressed/Planck_compressed.data | 24 +++++++++++++ .../likelihoods/Planck_compressed/__init__.py | 34 +++++++++++++++++++ 2 files changed, 58 insertions(+) create mode 100644 montepython/likelihoods/Planck_compressed/Planck_compressed.data create mode 100644 montepython/likelihoods/Planck_compressed/__init__.py diff --git a/montepython/likelihoods/Planck_compressed/Planck_compressed.data b/montepython/likelihoods/Planck_compressed/Planck_compressed.data new file mode 100644 index 00000000..4229176d --- /dev/null +++ b/montepython/likelihoods/Planck_compressed/Planck_compressed.data @@ -0,0 +1,24 @@ +# 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 = np.linalg.inv([[1., 0.54, -0.63, -0.86], [0.54, 1., -0.43, -0.48], [-0.63, -0.43, 1., 0.58], [-0.86, -0.48, 0.58, 1.]]) + +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 = np.array([[1.7488, 0.0074], [301.76, 0.14], [0.02228, 0.00023], [0.9660, 0.0061]]) +Planck_compressed.planck_values = [[1.7488, 0.0074], [301.76, 0.14], [0.02228, 0.00023], [0.9660, 0.0061]] + + diff --git a/montepython/likelihoods/Planck_compressed/__init__.py b/montepython/likelihoods/Planck_compressed/__init__.py new file mode 100644 index 00000000..7b0bf127 --- /dev/null +++ b/montepython/likelihoods/Planck_compressed/__init__.py @@ -0,0 +1,34 @@ +from montepython.likelihood_class import Likelihood_prior +import numpy as np + + +class Planck_compressed(Likelihood_prior): + + # 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): + + z_star = cosmo.z_optical_depth_unity() # z at which the optical depth is unity + D_A_star = cosmo.angular_distance(z_star) + rs_star = cosmo.rs_at_z(z_star) + + 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() + + values = np.array([R, l_a, omega_b, n_s]) + planck_values = zip(*self.planck_values) + + # 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 = ((values[:][0] - planck_values[0]) / + planck_values[1]) + + loglkl = -0.5 * np.dot(differences_over_sigma, + np.dot(self.inverse_correlation_matrix, + differences_over_sigma)) + return loglkl From a6cb5935ab6871e391f4873a310596653176d84b Mon Sep 17 00:00:00 2001 From: ardok-m Date: Mon, 22 May 2017 13:15:31 +0200 Subject: [PATCH 13/35] fix angular diamter distance. Needed comoving --- montepython/likelihoods/Planck_compressed/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/montepython/likelihoods/Planck_compressed/__init__.py b/montepython/likelihoods/Planck_compressed/__init__.py index 7b0bf127..b3f0c4de 100644 --- a/montepython/likelihoods/Planck_compressed/__init__.py +++ b/montepython/likelihoods/Planck_compressed/__init__.py @@ -11,7 +11,8 @@ class Planck_compressed(Likelihood_prior): def loglkl(self, cosmo, data): z_star = cosmo.z_optical_depth_unity() # z at which the optical depth is unity - D_A_star = cosmo.angular_distance(z_star) + # comoving angular diameter distance + D_A_star = cosmo.angular_distance(z_star) * (z_star + 1) rs_star = cosmo.rs_at_z(z_star) R = np.sqrt(cosmo.Omega_m()) * cosmo.Hubble(0) * D_A_star From a1ce2739afe61c1296b214b99bba85fe0e5f3b71 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Tue, 23 May 2017 10:04:43 +0200 Subject: [PATCH 14/35] separed mean and sigma values --- .../likelihoods/Planck_compressed/Planck_compressed.data | 7 +++---- montepython/likelihoods/Planck_compressed/__init__.py | 5 ++--- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/montepython/likelihoods/Planck_compressed/Planck_compressed.data b/montepython/likelihoods/Planck_compressed/Planck_compressed.data index 4229176d..9d6bb8b2 100644 --- a/montepython/likelihoods/Planck_compressed/Planck_compressed.data +++ b/montepython/likelihoods/Planck_compressed/Planck_compressed.data @@ -8,8 +8,6 @@ # Covariance matrix_ij = sigma_i sigma_j Correlation_matrix_ij -#Planck_compressed.inverse_correlation_matrix = np.linalg.inv([[1., 0.54, -0.63, -0.86], [0.54, 1., -0.43, -0.48], [-0.63, -0.43, 1., 0.58], [-0.86, -0.48, 0.58, 1.]]) - 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) @@ -18,7 +16,8 @@ Planck_compressed.inverse_correlation_matrix = [[ 4.51266537, -0.59697118, 0.755 #Planck_compressed.omega_b = [0.02228, 0.00023] #Planck_compressed.n_s = [0.9660, 0.0061] -#Planck_compressed.planck_values = np.array([[1.7488, 0.0074], [301.76, 0.14], [0.02228, 0.00023], [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 = [[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 index b3f0c4de..37b2e09c 100644 --- a/montepython/likelihoods/Planck_compressed/__init__.py +++ b/montepython/likelihoods/Planck_compressed/__init__.py @@ -21,15 +21,14 @@ def loglkl(self, cosmo, data): n_s = cosmo.n_s() values = np.array([R, l_a, omega_b, n_s]) - planck_values = zip(*self.planck_values) # 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 = ((values[:][0] - planck_values[0]) / - planck_values[1]) + differences_over_sigma = (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 From a5c6cf2d23b677caacbd8d057c1b24372477cf94 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Tue, 23 May 2017 13:59:23 +0200 Subject: [PATCH 15/35] added triangular prior and fixed some bugs --- montepython/data.py | 2 +- montepython/prior.py | 23 ++++++++++++++++++++--- 2 files changed, 21 insertions(+), 4 deletions(-) diff --git a/montepython/data.py b/montepython/data.py index 158ae434..9c6975c8 100644 --- a/montepython/data.py +++ b/montepython/data.py @@ -1006,7 +1006,7 @@ def __init__(self, array, key): raise io_mp.ParameterError("You chose scale 0 to param {}. Correct\ it.".format(key)) - self['role'] = array[-1] + self['role'] = array[5] self['tex_name'] = io_mp.get_tex_name(key) if array[3] == 0: self['status'] = 'fixed' 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])) From 3954be28584f00060643ad66d79fb9f36eda3cf0 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Thu, 25 May 2017 10:27:42 +0200 Subject: [PATCH 16/35] changed name of parent likelihood, behavior is the same --- montepython/likelihoods/Planck_compressed/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/montepython/likelihoods/Planck_compressed/__init__.py b/montepython/likelihoods/Planck_compressed/__init__.py index 37b2e09c..ca9be4ae 100644 --- a/montepython/likelihoods/Planck_compressed/__init__.py +++ b/montepython/likelihoods/Planck_compressed/__init__.py @@ -1,8 +1,8 @@ -from montepython.likelihood_class import Likelihood_prior +from montepython.likelihood_class import Likelihood import numpy as np -class Planck_compressed(Likelihood_prior): +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 From 20cba93886c9933f5e298861bb321fbf66740019 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Wed, 31 May 2017 10:46:23 +0200 Subject: [PATCH 17/35] store values as object to be accesible by wic data_mp_likelihoods.py --- montepython/likelihoods/Planck_compressed/__init__.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/montepython/likelihoods/Planck_compressed/__init__.py b/montepython/likelihoods/Planck_compressed/__init__.py index ca9be4ae..f378df47 100644 --- a/montepython/likelihoods/Planck_compressed/__init__.py +++ b/montepython/likelihoods/Planck_compressed/__init__.py @@ -20,12 +20,16 @@ def loglkl(self, cosmo, data): omega_b = cosmo.omega_b() n_s = cosmo.n_s() - values = np.array([R, l_a, omega_b, 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 = (values - np.array(self.planck_values_mean)) / self.planck_values_sigma + 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, From 814138565ef564ef9682182dd46aa09eabc7144f Mon Sep 17 00:00:00 2001 From: ardok-m Date: Thu, 27 Jul 2017 13:35:59 +0200 Subject: [PATCH 18/35] added only_error optiont to create error file when using NS --- montepython/initialise.py | 1 + montepython/io_mp.py | 18 +++++++++++------- 2 files changed, 12 insertions(+), 7 deletions(-) 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 index cd85713f..fca41ab5 100644 --- a/montepython/io_mp.py +++ b/montepython/io_mp.py @@ -255,12 +255,12 @@ def print_error(error_msg, data): 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' % @@ -286,7 +286,7 @@ def refresh_file(data): 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. @@ -302,6 +302,14 @@ def create_output_files(command_line, data): changing things here. """ + + if only_error: + #open/create error file + data.err_name = os.path.join(command_line.folder,"error_log.txt") + data.err = open(data.err_name,'a') + print 'Creating error file %s\n' % data.err_name + return + if command_line.restart is None: number = command_line.N else: @@ -340,10 +348,6 @@ def create_output_files(command_line, data): for line in open(command_line.restart, 'r'): data.out.write(line) - #open/create error file - data.err_name = os.path.join(command_line.folder,"error_log.txt") - data.err = open(data.err_name,'a') - print 'Creating error file %s\n' % data.err_name From fa8d9746990231375944d196f8a3aaa34610430d Mon Sep 17 00:00:00 2001 From: ardok-m Date: Thu, 27 Jul 2017 15:49:40 +0200 Subject: [PATCH 19/35] error_log file per chain --- montepython/io_mp.py | 99 ++++++++++++++++++++++++-------------------- 1 file changed, 55 insertions(+), 44 deletions(-) diff --git a/montepython/io_mp.py b/montepython/io_mp.py index fca41ab5..a88e201c 100644 --- a/montepython/io_mp.py +++ b/montepython/io_mp.py @@ -303,53 +303,64 @@ def create_output_files(command_line, data, only_error=False): """ - if only_error: - #open/create error file - data.err_name = os.path.join(command_line.folder,"error_log.txt") - data.err = open(data.err_name,'a') - print 'Creating error file %s\n' % data.err_name - return + 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 - 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: + 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): From 343b5b0f060815e7bef40e909cee0c8390d2bdcc Mon Sep 17 00:00:00 2001 From: ardok-m Date: Tue, 8 Aug 2017 15:07:59 +0200 Subject: [PATCH 20/35] exclude from analysis error_log files --- montepython/analyze.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/montepython/analyze.py b/montepython/analyze.py index 26ce89f0..5c8e7ce3 100644 --- a/montepython/analyze.py +++ b/montepython/analyze.py @@ -1144,6 +1144,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]): @@ -1151,6 +1155,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: @@ -1164,6 +1169,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 From b86ddd62599c9d285840e9349758af3032a56486 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Thu, 14 Sep 2017 15:36:02 +0200 Subject: [PATCH 21/35] fixed missing '' --- montepython/io_mp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) mode change 100644 => 100755 montepython/io_mp.py diff --git a/montepython/io_mp.py b/montepython/io_mp.py old mode 100644 new mode 100755 index a88e201c..26babacb --- a/montepython/io_mp.py +++ b/montepython/io_mp.py @@ -352,7 +352,7 @@ def create_output_files(command_line, data, only_error=False): 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'): + for line in open(command_line.restart[:-4]+suffix_err+'txt', 'r'): open_file.write(line) if f == 'out': From b81e79126ff6684f2670591dca3031c4a2187177 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Thu, 14 Sep 2017 15:38:02 +0200 Subject: [PATCH 22/35] fixed missing . in suffix --- montepython/io_mp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/montepython/io_mp.py b/montepython/io_mp.py index 26babacb..f7f3f8c7 100755 --- a/montepython/io_mp.py +++ b/montepython/io_mp.py @@ -352,7 +352,7 @@ def create_output_files(command_line, data, only_error=False): 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'): + for line in open(command_line.restart[:-4]+suffix_err+'.txt', 'r'): open_file.write(line) if f == 'out': From f9c326915116e86ee7e670d6ebdad200dbe89bd3 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Wed, 11 Oct 2017 10:42:22 +0200 Subject: [PATCH 23/35] modified bao_boss to interact with work_in_class --- montepython/likelihoods/bao_boss/__init__.py | 5 +++++ 1 file changed, 5 insertions(+) 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 From a55a5517c0026da07e90327bc57d1d904af92cab Mon Sep 17 00:00:00 2001 From: ardok-m Date: Mon, 23 Oct 2017 12:39:43 +0200 Subject: [PATCH 24/35] updated to match stable CosmoHammer version --- montepython/cosmo_hammer.py | 10 ++++++---- montepython/run.py | 2 ++ 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/montepython/cosmo_hammer.py b/montepython/cosmo_hammer.py index efec34d4..50339f76 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 @@ -127,7 +129,7 @@ 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, 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) + From afeb21d06f16e8fa4c97977f1343f288bcfc3135 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Tue, 24 Oct 2017 10:14:52 +0200 Subject: [PATCH 25/35] fixed option input cosmo_hammer --- montepython/cosmo_hammer.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/montepython/cosmo_hammer.py b/montepython/cosmo_hammer.py index 50339f76..b0f8e311 100644 --- a/montepython/cosmo_hammer.py +++ b/montepython/cosmo_hammer.py @@ -103,7 +103,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: @@ -133,9 +135,6 @@ def run(cosmo, data, command_line): params=params, likelihoodComputationChain=chain, filePrefix=file_prefix, - walkersRatio=50, - burninIterations=10, - sampleIterations=30, storageUtil=derived_util, threadCount=num_threads, **data.CH_arguments) From a738c0c5374650d3add8fbdab6c8d3d0caded67e Mon Sep 17 00:00:00 2001 From: ardok-m Date: Wed, 25 Oct 2017 13:19:30 +0200 Subject: [PATCH 26/35] fixed derived CH when not computed --- montepython/cosmo_hammer.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/montepython/cosmo_hammer.py b/montepython/cosmo_hammer.py index b0f8e311..16f45a1f 100644 --- a/montepython/cosmo_hammer.py +++ b/montepython/cosmo_hammer.py @@ -50,6 +50,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): """ @@ -59,6 +63,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( @@ -196,8 +203,11 @@ 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]) final = np.concatenate((pos, derived), axis=1) posFile.write("\n".join( From b70913f0dad68c10502332f0d3898637dc0b78ed Mon Sep 17 00:00:00 2001 From: ardok-m Date: Wed, 25 Oct 2017 15:04:49 +0200 Subject: [PATCH 27/35] save errors in different file --- montepython/cosmo_hammer.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/montepython/cosmo_hammer.py b/montepython/cosmo_hammer.py index 16f45a1f..de05bab8 100644 --- a/montepython/cosmo_hammer.py +++ b/montepython/cosmo_hammer.py @@ -35,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_' @@ -166,6 +167,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) @@ -186,11 +189,20 @@ def from_CH_output_to_chains(folder): ## Create the array of ones ones = np.array([[1] for _ in range(len(lkl))]) + for x, p in [chains, lkl]: + if p is not -np.inf: + chainsClean.append(x) + lklClean.append(p) + else: + error.append(x) + ## Concatenate everything and save to file - final = np.concatenate((ones, lkl, chains), axis=1) + final = np.concatenate((ones, lkl, 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): From 50ad1506630032faae6c7352afe1a65d98e3f05a Mon Sep 17 00:00:00 2001 From: ardok-m Date: Wed, 25 Oct 2017 15:54:53 +0200 Subject: [PATCH 28/35] fixed problems saving errors diff files --- montepython/cosmo_hammer.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/montepython/cosmo_hammer.py b/montepython/cosmo_hammer.py index de05bab8..67f12a1c 100644 --- a/montepython/cosmo_hammer.py +++ b/montepython/cosmo_hammer.py @@ -183,21 +183,22 @@ 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. + for x, p in zip(chains, lkl): + if p != -np.inf: + chainsClean.append(x) + lklClean.append([p]) + else: + error.append(x) + ## First, reshape the lkl array - lkl = np.array([[elem] for elem in lkl]) + #lkl = np.array([[elem] for elem in lkl]) ## Create the array of ones - ones = np.array([[1] for _ in range(len(lkl))]) + ones = np.array([[1] for _ in range(len(lklClean))]) - for x, p in [chains, lkl]: - if p is not -np.inf: - chainsClean.append(x) - lklClean.append(p) - else: - error.append(x) ## Concatenate everything and save to file - final = np.concatenate((ones, lkl, chainsClean), 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) From 13f0ab793d044adddb800260ab77f161778481b8 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Fri, 27 Oct 2017 11:59:35 +0200 Subject: [PATCH 29/35] added function to store derived paramters from cosmo --- montepython/cosmo_hammer.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/montepython/cosmo_hammer.py b/montepython/cosmo_hammer.py index 67f12a1c..9595b612 100644 --- a/montepython/cosmo_hammer.py +++ b/montepython/cosmo_hammer.py @@ -101,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(): @@ -231,3 +232,30 @@ 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. Copied from 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'])) + for name, value in derived.iteritems(): + data.mcmc_parameters[name]['current'] = value + 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") + for elem in data.get_mcmc_parameters(['derived']): + data.mcmc_parameters[elem]['current'] /= \ + data.mcmc_parameters[elem]['scale'] From 47b2fa44c467479c074b47bbedd5ab301c9a0adf Mon Sep 17 00:00:00 2001 From: ardok-m Date: Fri, 27 Oct 2017 13:55:54 +0200 Subject: [PATCH 30/35] remove nan's from error-log file --- montepython/cosmo_hammer.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/montepython/cosmo_hammer.py b/montepython/cosmo_hammer.py index 9595b612..0d5bb184 100644 --- a/montepython/cosmo_hammer.py +++ b/montepython/cosmo_hammer.py @@ -184,15 +184,13 @@ 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, separe errors from accepted points: for x, p in zip(chains, lkl): if p != -np.inf: chainsClean.append(x) lklClean.append([p]) else: - error.append(x) - - ## First, reshape the lkl array - #lkl = np.array([[elem] for elem in lkl]) + 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(lklClean))]) From 16085a36b8f51a4ec5a93ed4a511936f8d7f474d Mon Sep 17 00:00:00 2001 From: ardok-m Date: Tue, 14 Nov 2017 14:52:36 +0100 Subject: [PATCH 31/35] convert to -loglkl the CH loglkl when analyzing --- montepython/cosmo_hammer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/montepython/cosmo_hammer.py b/montepython/cosmo_hammer.py index 0d5bb184..f73de175 100644 --- a/montepython/cosmo_hammer.py +++ b/montepython/cosmo_hammer.py @@ -188,7 +188,7 @@ def from_CH_output_to_chains(folder): for x, p in zip(chains, lkl): if p != -np.inf: chainsClean.append(x) - lklClean.append([p]) + 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 From f8b41424040b691acb22975a3939737b56b1a3eb Mon Sep 17 00:00:00 2001 From: ardok-m Date: Mon, 14 May 2018 13:07:36 +0200 Subject: [PATCH 32/35] sort by name the derived parameters (placed at the end of the param list) --- montepython/cosmo_hammer.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/montepython/cosmo_hammer.py b/montepython/cosmo_hammer.py index f73de175..28be6c8f 100644 --- a/montepython/cosmo_hammer.py +++ b/montepython/cosmo_hammer.py @@ -219,7 +219,8 @@ def persistValues(self, posFile, probFile, pos, prob, data): derived_not_computed = [np.NaN] * len(CH_derived_names) derived = np.array( - [[a for a in elem.itervalues()] if elem else derived_not_computed 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( From 62315b912c18bbd83787de99d6a392a2398bca1b Mon Sep 17 00:00:00 2001 From: ardok-m Date: Mon, 28 May 2018 10:49:07 +0200 Subject: [PATCH 33/35] fixed Planck_compressed. Use default CLASS vars --- montepython/likelihoods/Planck_compressed/__init__.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/montepython/likelihoods/Planck_compressed/__init__.py b/montepython/likelihoods/Planck_compressed/__init__.py index f378df47..0390461b 100644 --- a/montepython/likelihoods/Planck_compressed/__init__.py +++ b/montepython/likelihoods/Planck_compressed/__init__.py @@ -10,10 +10,14 @@ class Planck_compressed(Likelihood): def loglkl(self, cosmo, data): - z_star = cosmo.z_optical_depth_unity() # z at which the optical depth is unity + 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 = 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 From 81788c56e89749f7bff0831c5b2dffe5c8d7da75 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Wed, 20 Jun 2018 10:42:10 +0200 Subject: [PATCH 34/35] fixed error storing derived parameters --- montepython/cosmo_hammer.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/montepython/cosmo_hammer.py b/montepython/cosmo_hammer.py index 28be6c8f..93d3a108 100644 --- a/montepython/cosmo_hammer.py +++ b/montepython/cosmo_hammer.py @@ -235,7 +235,7 @@ def persistValues(self, posFile, probFile, pos, prob, data): def store_cosmo_derived(ctx): """ - Store derived parameters. Copied from sampler.compute_lkl + Store derived parameters. Based on sampler.compute_lkl """ from classy import CosmoSevereError @@ -246,8 +246,6 @@ def store_cosmo_derived(ctx): try: derived = cosmo.get_current_derived_parameters( data.get_mcmc_parameters(['derived'])) - for name, value in derived.iteritems(): - data.mcmc_parameters[name]['current'] = value except AttributeError: # This happens if the classy wrapper is still using the old # convention, expecting data as the input parameter @@ -255,6 +253,9 @@ def store_cosmo_derived(ctx): except CosmoSevereError: raise io_mp.CosmologicalModuleError( "Could not write the current derived parameters") - for elem in data.get_mcmc_parameters(['derived']): - data.mcmc_parameters[elem]['current'] /= \ - data.mcmc_parameters[elem]['scale'] + + # 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) + + ctx.add('key_data', derived) From bbcedf87904ecfc7033e4fbbf866ac55109dbdd6 Mon Sep 17 00:00:00 2001 From: ardok-m Date: Wed, 20 Jun 2018 11:25:50 +0200 Subject: [PATCH 35/35] fixed identation in store_cosmo_derived --- montepython/cosmo_hammer.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/montepython/cosmo_hammer.py b/montepython/cosmo_hammer.py index 93d3a108..a384f71e 100644 --- a/montepython/cosmo_hammer.py +++ b/montepython/cosmo_hammer.py @@ -254,8 +254,8 @@ def store_cosmo_derived(ctx): raise io_mp.CosmologicalModuleError( "Could not write the current derived parameters") - # 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) + ctx.add('key_data', derived) - 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)