From 3bb57afeee6de73cf637fb6f02d8df3e59dac5d2 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 19 Sep 2019 18:19:58 -0300 Subject: [PATCH 01/27] mark branch as develop --- packages/_version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/_version.py b/packages/_version.py index d17e4edd..abf19198 100644 --- a/packages/_version.py +++ b/packages/_version.py @@ -1 +1 @@ -__version__ = "v0.2.6" \ No newline at end of file +__version__ = "v0.2.6-dev" \ No newline at end of file From 1a2005beec0d9e909005ce79fbacbcb740e6168b Mon Sep 17 00:00:00 2001 From: Gabriel Date: Mon, 23 Sep 2019 11:06:56 -0300 Subject: [PATCH 02/27] small fix when parallel runs are used --- packages/out/inparams_out.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/packages/out/inparams_out.py b/packages/out/inparams_out.py index 9aeea7f5..8595697a 100644 --- a/packages/out/inparams_out.py +++ b/packages/out/inparams_out.py @@ -1,4 +1,5 @@ +from os.path import isfile from time import strftime @@ -6,8 +7,12 @@ def main(npd, file_end, **kwargs): """ Save input 'params_input.dat' file as output file. """ - # Read - with open('./params_input' + file_end + '.dat', 'r') as f: + + fname = './params_input.dat' + if isfile('./params_input' + file_end + '.dat'): + fname = './params_input' + file_end + '.dat' + + with open(fname, 'r') as f: # Read file into data var. data = f.readlines() From ce5bd9fe529c2a74ca4db68a06469ae660f1b2d0 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Tue, 24 Sep 2019 12:00:12 -0300 Subject: [PATCH 03/27] intermediate stage, testing phase --- packages/best_fit/zaInterp.py | 10 +- packages/inp/input_params.py | 2 +- packages/synth_clust/imf.py | 106 +++++++++++++++++----- packages/synth_clust/mass_distribution.py | 5 +- 4 files changed, 91 insertions(+), 32 deletions(-) diff --git a/packages/best_fit/zaInterp.py b/packages/best_fit/zaInterp.py index 979489d0..39b8c018 100644 --- a/packages/best_fit/zaInterp.py +++ b/packages/best_fit/zaInterp.py @@ -46,11 +46,11 @@ def main(theor_tracks, fundam_params, varIdxs, model): al = ah - 1 model_proper.append(par[ah]) a_model = model[i - j] - # If it is the parameter mass. - elif i == 4: - # Select the closest value in the array of allowed values. - model_proper.append(min( - par, key=lambda x: abs(x - model[i - j]))) + # # If it is the parameter mass. + # elif i == 4: + # # Select the closest value in the array of allowed values. + # model_proper.append(min( + # par, key=lambda x: abs(x - model[i - j]))) else: model_proper.append(model[i - j]) else: diff --git a/packages/inp/input_params.py b/packages/inp/input_params.py index 3c19aeaf..9aa5d552 100644 --- a/packages/inp/input_params.py +++ b/packages/inp/input_params.py @@ -256,7 +256,7 @@ def main(mypath, pars_f_path): 'tolstoy', 'duong', 'dolphin', 'mighell', 'kdeKL', 'dolphin_kde') # Accepted IMF functions. imf_funcs = ('chabrier_2001_exp', 'chabrier_2001_log', 'kroupa_1993', - 'kroupa_2002') + 'kroupa_2002', 'salpeter_1955') # Optimizing algorithm # TODO 'brute', 'emcee', 'abc' optimz_algors = ('ptemcee', 'boot+GA', 'n') diff --git a/packages/synth_clust/imf.py b/packages/synth_clust/imf.py index c7c53907..4b3dbc43 100644 --- a/packages/synth_clust/imf.py +++ b/packages/synth_clust/imf.py @@ -1,6 +1,7 @@ from ..core_imp import np from scipy.integrate import quad +from scipy.interpolate import interp1d def main(IMF_name, m_high, masses): @@ -27,19 +28,19 @@ def main(IMF_name, m_high, masses): """ print("Sampling selected IMF ({})".format(IMF_name)) + # Low mass limits are defined for each IMF to avoid numerical # issues when integrating. - imfs_dict = {'chabrier_2001_exp': (0.01), 'chabrier_2001_log': (0.01), - 'kroupa_1993': (0.081), 'kroupa_2002': (0.011)} - - # Set IMF low mass limit. + imfs_dict = { + 'chabrier_2001_exp': 0.01, 'chabrier_2001_log': 0.01, + 'kroupa_1993': 0.081, 'kroupa_2002': 0.011, 'salpeter_1955': 0.3} + # IMF low mass limit. m_low = imfs_dict[IMF_name] - # Set IMF mass interpolation step. - # The step (mass_step) should not be too small as it will have an impact - # on the performance of the 'mass_distribution' function later on. - mass_step = 0.1 + # IMF mass interpolation step. + mass_step = 0.05 + ################## # Obtain normalization constant. This is equivalent to 'k/M_total' in # Eq. (7) of Popescu & Hanson 2009 (138:1724-1740; PH09). # For m_high > 100 Mo the differences in the resulting normalization @@ -62,40 +63,93 @@ def main(IMF_name, m_high, masses): # Normalize number of stars by constant. N_dist = np.asarray(N_dist) * norm_const - st_dist_mass = massDist(m_low, mass_up, N_dist, masses) + # st_dist_mass = massDist(m_low, mass_up, N_dist, masses) + ################## + + ################## + # Obtain normalization constant (k = \int_{m_low}^{m_up} \xi(m) dm). This + # makes the IMF behave like a PDF. + norm_const = quad(integral_IMF_N, m_low, m_high, args=(IMF_name))[0] + + # The CDF is defined as: $F(m)= \int_{m_low}^{m} PDF(m) dm$ + # Sample the CDF + mass_values = np.arange(m_low, m_high, mass_step) + CDF_samples = [] + for m in mass_values: + CDF_samples.append(quad(integral_IMF_N, m_low, m, args=(IMF_name))[0]) + CDF_samples = np.array(CDF_samples) / norm_const + + inv_cdf = interp1d(CDF_samples, mass_values) + CDF_min, CDF_max = CDF_samples.min(), CDF_samples.max() + + def sampled_inv_cdf(N): + mr = np.random.rand(N) + mr = mr[(mr >= CDF_min) & (mr <= CDF_max)] + return inv_cdf(mr) + + # sampled_IMF = sampled_inv_cdf(10000) + + # st_dist_mass2 = (sampled_IMF, np.cumsum(sampled_IMF)) + ################## + + def return_intersection(hist_1, hist_2): + minima = np.minimum(hist_1, hist_2) + intersection = np.true_divide(np.sum(minima), np.sum(hist_2)) + return intersection + + import matplotlib.pyplot as plt + dd = [] + for mass in np.random.uniform(100., 1000., 1000): + data_1 = massDist(m_low, mass_up, N_dist, [mass]) + data_1 = data_1[mass] + sampled_IMF = sampled_inv_cdf(10000) + data_2 = sampled_IMF[:np.searchsorted(np.cumsum(sampled_IMF), mass)] + hist_1, _ = np.histogram(data_1, bins=100, range=[0., 150]) + hist_2, _ = np.histogram(data_2, bins=100, range=[0., 150]) + rr = return_intersection(hist_1, hist_2) + dd.append([mass, rr]) + if rr < .7: + print(mass, rr, len(data_1), len(data_2)) + plt.hist(data_1, 100, alpha=.5, density=True, label='old') + plt.hist(data_2, 100, alpha=.5, density=True, label='new') + plt.xscale('log');plt.yscale('log') + plt.legend();plt.show() + + dd = np.array(dd).T + plt.scatter(*dd) + plt.show() return st_dist_mass def integral_IMF_M(m_star, IMF_name): ''' - Return the properly normalized function to perform the integration of the - selected IMF. Returns mass values. + Returns mass values: $$ ''' - imf_val = m_star * imfs(IMF_name, m_star) - return imf_val + return m_star * imfs(IMF_name, m_star) def integral_IMF_N(m_star, IMF_name): ''' - Return the properly normalized function to perform the integration of the - selected IMF. Returns number of stars. + Returns number of stars: $$ ''' - imf_val = imfs(IMF_name, m_star) - return imf_val + return imfs(IMF_name, m_star) def imfs(IMF_name, m_star): - ''' + """ Define any number of IMFs. - ''' + + The package https://github.com/keflavich/imf has some more (I think, + 24-09-2019). + """ + if IMF_name == 'kroupa_1993': # Kroupa, Tout & Gilmore. (1993) piecewise IMF. # http://adsabs.harvard.edu/abs/1993MNRAS.262..545K # Eq. (13), p. 572 (28) alpha = [-1.3, -2.2, -2.7] - m_i = [0.08, 0.5, 1.] - m0, m1, m2 = m_i + m0, m1, m2 = [0.08, 0.5, 1.] factor = [0.035, 0.019, 0.019] if m0 < m_star <= m1: i = 0 @@ -106,9 +160,8 @@ def imfs(IMF_name, m_star): imf_val = factor[i] * (m_star ** alpha[i]) elif IMF_name == 'kroupa_2002': - # Kroupa (2002) piecewise IMF (taken from MASSCLEAN article). - # http://adsabs.harvard.edu/abs/2002Sci...295...82K - # Eq. (2) & (3), p. 1725 + # Kroupa (2002) Salpeter (1995) piecewise IMF taken from MASSCLEAN + # article, Eq. (2) & (3), p. 1725 alpha = [-0.3, -1.3, -2.3] m0, m1, m2 = [0.01, 0.08, 0.5] factor = [(1. / m1) ** alpha[0], (1. / m1) ** alpha[1], @@ -135,6 +188,11 @@ def imfs(IMF_name, m_star): # Eq (8) imf_val = 3. * m_star ** (-3.3) * np.exp(-(716.4 / m_star) ** 0.25) + elif IMF_name == 'salpeter_1955': + # Salpeter (1955) IMF. + # https://ui.adsabs.harvard.edu/abs/1955ApJ...121..161S/ + imf_val = m_star ** -2.35 + return imf_val diff --git a/packages/synth_clust/mass_distribution.py b/packages/synth_clust/mass_distribution.py index db1f18c7..f5d1e340 100644 --- a/packages/synth_clust/mass_distribution.py +++ b/packages/synth_clust/mass_distribution.py @@ -1,5 +1,5 @@ -# import numpy as np +import numpy as np def main(st_dist_mass, M_total): @@ -23,6 +23,7 @@ def main(st_dist_mass, M_total): # base, scale, N_stars = st_dist_mass[M_total][:-1] # mass_dist = np.random.random(N_stars) * scale + base # else: - mass_dist = st_dist_mass[M_total] + # mass_dist = st_dist_mass[M_total] + mass_dist = st_dist_mass[0][:np.searchsorted(st_dist_mass[1], M_total)] return mass_dist From f151cbc64720b38b52974d385a4169205083c3d2 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 26 Sep 2019 12:27:35 -0300 Subject: [PATCH 04/27] minor, in place for #437 --- packages/out/cluster_members_file.py | 1 + 1 file changed, 1 insertion(+) diff --git a/packages/out/cluster_members_file.py b/packages/out/cluster_members_file.py index d5dec219..7ca09c1c 100644 --- a/packages/out/cluster_members_file.py +++ b/packages/out/cluster_members_file.py @@ -53,6 +53,7 @@ def main( st[:3] + mag + emag + cols + ecols + kinem + ekinem + [st[9], idx[i]]) + # TODO: this block gets really slow for large clusters # Add "incomplete" data in cluster region to file. ids_data = list(zip(*data))[0] for i, st in enumerate(clp['cl_region_i']): From 63f717daf2f5435451493f0b012c1d6417256204 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 26 Sep 2019 12:28:41 -0300 Subject: [PATCH 05/27] much faster algor, see stackoverflow.com/a/58118907/1391441 --- packages/decont_algors/decont_algors.py | 8 +++++--- packages/out/prep_plots.py | 13 ++++++------- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/packages/decont_algors/decont_algors.py b/packages/decont_algors/decont_algors.py index 44aa81d2..9b853218 100644 --- a/packages/decont_algors/decont_algors.py +++ b/packages/decont_algors/decont_algors.py @@ -58,10 +58,12 @@ def filterMPs(memb_probs_cl_region, cl_region_i, cl_region_c): ids_c, ids_i = list(zip(*cl_region_c))[0], list(zip(*cl_region_i))[0] memb_probs_cl_region_c = np.zeros(len(ids_c)) + dict_ids_c = {_: i for i, _ in enumerate(ids_c)} for i, id_i in enumerate(ids_i): - if id_i in ids_c: - j = ids_c.index(id_i) - memb_probs_cl_region_c[j] = memb_probs_cl_region[i] + try: + memb_probs_cl_region_c[dict_ids_c[id_i]] = memb_probs_cl_region[i] + except KeyError: + pass return memb_probs_cl_region_c diff --git a/packages/out/prep_plots.py b/packages/out/prep_plots.py index 67399257..7489c9aa 100644 --- a/packages/out/prep_plots.py +++ b/packages/out/prep_plots.py @@ -554,14 +554,13 @@ def complSeparate(cl_region_i, membs_i, memb_prob_avrg_sort): """ mags_c, mags_i, cols_c, cols_i, colors_c, colors_i = [[] for _ in range(6)] ids_c = list(zip(*memb_prob_avrg_sort))[0] + dict_ids_c = {_: i for i, _ in enumerate(ids_c)} for i, star in enumerate(cl_region_i): - id_i = star[0] - if id_i in ids_c: - j = ids_c.index(id_i) - mags_c.append(memb_prob_avrg_sort[j][3][0]) - cols_c.append(memb_prob_avrg_sort[j][5]) - colors_c.append(memb_prob_avrg_sort[j][9]) - else: + try: + mags_c.append(memb_prob_avrg_sort[dict_ids_c[star[0]]][3][0]) + cols_c.append(memb_prob_avrg_sort[dict_ids_c[star[0]]][5]) + colors_c.append(memb_prob_avrg_sort[dict_ids_c[star[0]]][9]) + except KeyError: mags_i.append(star[3][0]) cols_i.append(star[5]) colors_i.append(membs_i[i]) From 74e1583bc4316aa7b2b2f3a3e555216be16ed4ff Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 26 Sep 2019 12:29:02 -0300 Subject: [PATCH 06/27] useful when only IDs are passed in the membership file --- packages/decont_algors/read_da.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/packages/decont_algors/read_da.py b/packages/decont_algors/read_da.py index d5ad3ee5..efad0b02 100644 --- a/packages/decont_algors/read_da.py +++ b/packages/decont_algors/read_da.py @@ -1,4 +1,5 @@ +import numpy as np from astropy.io import ascii @@ -13,8 +14,12 @@ def main(cl_region, memb_file, readda_idcol=0, readda_mpcol=-2): # Read IDs and MPs from file. data = ascii.read(memb_file) # Read IDs as strings since that is how they are stored in 'cl_region' - id_list, memb_probs = [str(_) for _ in data.columns[readda_idcol]],\ - data.columns[readda_mpcol] + id_list = [str(_) for _ in data.columns[readda_idcol]] + try: + memb_probs = data.columns[readda_mpcol] + except IndexError: + print(" WARNING: MPs column not found. Assigned MP=1. to all stars") + memb_probs = np.ones(len(data)) N_not, memb_probs_cl_region = 0, [] # Assign probabilities read from file according to the star's IDs. From 15be4b1052b0879b7ffbe567f0cc0bb8700a931c Mon Sep 17 00:00:00 2001 From: Gabriel Date: Sun, 29 Sep 2019 12:11:51 -0300 Subject: [PATCH 07/27] remove not used files from emcee folder --- packages/best_fit/emcee3rc2/mpi_pool.py | 18 -- packages/best_fit/emcee3rc2/ptsampler.py | 18 -- packages/best_fit/emcee3rc2/tests/__init__.py | 0 .../emcee3rc2/tests/integration/__init__.py | 0 .../emcee3rc2/tests/integration/test_de.py | 20 -- .../tests/integration/test_de_snooker.py | 17 -- .../tests/integration/test_gaussian.py | 62 ----- .../emcee3rc2/tests/integration/test_kde.py | 20 -- .../tests/integration/test_proposal.py | 85 ------- .../tests/integration/test_stretch.py | 28 --- .../emcee3rc2/tests/integration/test_walk.py | 16 -- .../best_fit/emcee3rc2/tests/unit/__init__.py | 0 .../emcee3rc2/tests/unit/test_autocorr.py | 49 ---- .../emcee3rc2/tests/unit/test_backends.py | 215 ------------------ .../emcee3rc2/tests/unit/test_blobs.py | 51 ----- .../emcee3rc2/tests/unit/test_sampler.py | 186 --------------- .../emcee3rc2/tests/unit/test_state.py | 43 ---- .../emcee3rc2/tests/unit/test_stretch.py | 39 ---- 18 files changed, 867 deletions(-) delete mode 100644 packages/best_fit/emcee3rc2/mpi_pool.py delete mode 100644 packages/best_fit/emcee3rc2/ptsampler.py delete mode 100644 packages/best_fit/emcee3rc2/tests/__init__.py delete mode 100644 packages/best_fit/emcee3rc2/tests/integration/__init__.py delete mode 100644 packages/best_fit/emcee3rc2/tests/integration/test_de.py delete mode 100644 packages/best_fit/emcee3rc2/tests/integration/test_de_snooker.py delete mode 100644 packages/best_fit/emcee3rc2/tests/integration/test_gaussian.py delete mode 100644 packages/best_fit/emcee3rc2/tests/integration/test_kde.py delete mode 100644 packages/best_fit/emcee3rc2/tests/integration/test_proposal.py delete mode 100644 packages/best_fit/emcee3rc2/tests/integration/test_stretch.py delete mode 100644 packages/best_fit/emcee3rc2/tests/integration/test_walk.py delete mode 100644 packages/best_fit/emcee3rc2/tests/unit/__init__.py delete mode 100644 packages/best_fit/emcee3rc2/tests/unit/test_autocorr.py delete mode 100644 packages/best_fit/emcee3rc2/tests/unit/test_backends.py delete mode 100644 packages/best_fit/emcee3rc2/tests/unit/test_blobs.py delete mode 100644 packages/best_fit/emcee3rc2/tests/unit/test_sampler.py delete mode 100644 packages/best_fit/emcee3rc2/tests/unit/test_state.py delete mode 100644 packages/best_fit/emcee3rc2/tests/unit/test_stretch.py diff --git a/packages/best_fit/emcee3rc2/mpi_pool.py b/packages/best_fit/emcee3rc2/mpi_pool.py deleted file mode 100644 index a8c7fea5..00000000 --- a/packages/best_fit/emcee3rc2/mpi_pool.py +++ /dev/null @@ -1,18 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import print_function, absolute_import - -try: - from schwimmbad import MPIPool -except ImportError: - - class MPIPool(object): - - def __init__(self, *args, **kwargs): - raise ImportError( - "The MPIPool from emcee has been forked to " - "https://github.com/adrn/schwimmbad, " - "please install that package to continue using the MPIPool" - ) - -__all__ = ["MPIPool"] diff --git a/packages/best_fit/emcee3rc2/ptsampler.py b/packages/best_fit/emcee3rc2/ptsampler.py deleted file mode 100644 index ac4638f2..00000000 --- a/packages/best_fit/emcee3rc2/ptsampler.py +++ /dev/null @@ -1,18 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import print_function, absolute_import - -try: - from ptemcee import Sampler as PTSampler -except ImportError: - - class PTSampler(object): - - def __init__(self, *args, **kwargs): - raise ImportError( - "The PTSampler from emcee has been forked to " - "https://github.com/willvousden/ptemcee, " - "please install that package to continue using the PTSampler" - ) - -__all__ = ["PTSampler"] diff --git a/packages/best_fit/emcee3rc2/tests/__init__.py b/packages/best_fit/emcee3rc2/tests/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/packages/best_fit/emcee3rc2/tests/integration/__init__.py b/packages/best_fit/emcee3rc2/tests/integration/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/packages/best_fit/emcee3rc2/tests/integration/test_de.py b/packages/best_fit/emcee3rc2/tests/integration/test_de.py deleted file mode 100644 index 6ce2f0b1..00000000 --- a/packages/best_fit/emcee3rc2/tests/integration/test_de.py +++ /dev/null @@ -1,20 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -from emcee import moves -from .test_proposal import _test_normal, _test_uniform - -__all__ = ["test_normal_de", "test_normal_de_no_gamma", "test_uniform_de"] - - -def test_normal_de(**kwargs): - _test_normal(moves.DEMove(), **kwargs) - - -def test_normal_de_no_gamma(**kwargs): - _test_normal(moves.DEMove(gamma0=1.0), **kwargs) - - -def test_uniform_de(**kwargs): - _test_uniform(moves.DEMove(), **kwargs) diff --git a/packages/best_fit/emcee3rc2/tests/integration/test_de_snooker.py b/packages/best_fit/emcee3rc2/tests/integration/test_de_snooker.py deleted file mode 100644 index d4afc87f..00000000 --- a/packages/best_fit/emcee3rc2/tests/integration/test_de_snooker.py +++ /dev/null @@ -1,17 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -from emcee import moves -from .test_proposal import _test_normal, _test_uniform - -__all__ = ["test_normal_de_snooker", "test_uniform_de_snooker"] - - -def test_normal_de_snooker(**kwargs): - kwargs["nsteps"] = 4000 - _test_normal(moves.DESnookerMove(), **kwargs) - - -def test_uniform_de_snooker(**kwargs): - _test_uniform(moves.DESnookerMove(), **kwargs) diff --git a/packages/best_fit/emcee3rc2/tests/integration/test_gaussian.py b/packages/best_fit/emcee3rc2/tests/integration/test_gaussian.py deleted file mode 100644 index 6a11b0d1..00000000 --- a/packages/best_fit/emcee3rc2/tests/integration/test_gaussian.py +++ /dev/null @@ -1,62 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -import pytest -import numpy as np -from itertools import product - -from emcee import moves - -from .test_proposal import _test_normal, _test_uniform - -__all__ = ["test_normal_gaussian", "test_uniform_gaussian", - "test_normal_gaussian_nd"] - - -@pytest.mark.parametrize("mode,factor", product( - ["vector"], [None, 2.0, 5.0] -)) -def test_normal_gaussian(mode, factor, **kwargs): - _test_normal(moves.GaussianMove(0.5, mode=mode, factor=factor), **kwargs) - - -@pytest.mark.parametrize("mode,factor", product( - ["vector", "random", "sequential"], [None, 2.0] -)) -def test_normal_gaussian_nd(mode, factor, **kwargs): - ndim = 3 - kwargs["nsteps"] = 8000 - - # Isotropic. - _test_normal(moves.GaussianMove(0.5, factor=factor, mode=mode), ndim=ndim, - **kwargs) - - # Axis-aligned. - _test_normal(moves.GaussianMove(0.5 * np.ones(ndim), factor=factor, - mode=mode), ndim=ndim, **kwargs) - with pytest.raises(ValueError): - _test_normal(moves.GaussianMove(0.5 * np.ones(ndim-1), factor=factor, - mode=mode), ndim=ndim, - **kwargs) - - # Full matrix. - if mode == "vector": - _test_normal(moves.GaussianMove(np.diag(0.5 * np.ones(ndim)), - factor=factor, mode=mode), ndim=ndim, - **kwargs) - with pytest.raises(ValueError): - _test_normal(moves.GaussianMove(np.diag(0.5 * np.ones(ndim-1))), - ndim=ndim, **kwargs) - else: - with pytest.raises(ValueError): - _test_normal(moves.GaussianMove(np.diag(0.5 * np.ones(ndim)), - factor=factor, mode=mode), - ndim=ndim, **kwargs) - - -@pytest.mark.parametrize("mode,factor", product( - ["vector"], [None, 2.0, 5.0] -)) -def test_uniform_gaussian(mode, factor, **kwargs): - _test_uniform(moves.GaussianMove(0.5, factor=factor, mode=mode), **kwargs) diff --git a/packages/best_fit/emcee3rc2/tests/integration/test_kde.py b/packages/best_fit/emcee3rc2/tests/integration/test_kde.py deleted file mode 100644 index 647f8aa5..00000000 --- a/packages/best_fit/emcee3rc2/tests/integration/test_kde.py +++ /dev/null @@ -1,20 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -from emcee import moves -from .test_proposal import _test_normal, _test_uniform - -__all__ = ["test_normal_kde", "test_uniform_kde", "test_nsplits_kde"] - - -def test_normal_kde(**kwargs): - _test_normal(moves.KDEMove(), **kwargs) - - -def test_uniform_kde(**kwargs): - _test_uniform(moves.KDEMove(), **kwargs) - - -def test_nsplits_kde(**kwargs): - _test_normal(moves.KDEMove(nsplits=5), **kwargs) diff --git a/packages/best_fit/emcee3rc2/tests/integration/test_proposal.py b/packages/best_fit/emcee3rc2/tests/integration/test_proposal.py deleted file mode 100644 index 2c759a71..00000000 --- a/packages/best_fit/emcee3rc2/tests/integration/test_proposal.py +++ /dev/null @@ -1,85 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -import numpy as np -from scipy import stats - -import emcee - -__all__ = ["_test_normal", "_test_uniform"] - - -def normal_log_prob_blobs(params): - return -0.5 * np.sum(params**2), params - - -def normal_log_prob(params): - return -0.5 * np.sum(params**2) - - -def uniform_log_prob(params): - if np.any(params > 1) or np.any(params < 0): - return -np.inf - return 0.0 - - -def _test_normal(proposal, ndim=1, nwalkers=32, nsteps=2000, seed=1234, - check_acceptance=True, pool=None, blobs=False): - # Set up the random number generator. - np.random.seed(seed) - - # Initialize the ensemble and proposal. - coords = np.random.randn(nwalkers, ndim) - - if blobs: - lp = normal_log_prob_blobs - else: - lp = normal_log_prob - - sampler = emcee.EnsembleSampler(nwalkers, ndim, lp, - moves=proposal, pool=pool) - if hasattr(proposal, "ntune") and proposal.ntune > 0: - coords = sampler.run_mcmc(coords, proposal.ntune, tune=True) - sampler.reset() - sampler.run_mcmc(coords, nsteps) - - # Check the acceptance fraction. - if check_acceptance: - acc = sampler.acceptance_fraction - assert np.all((acc < 0.9) * (acc > 0.1)), \ - "Invalid acceptance fraction\n{0}".format(acc) - - # Check the resulting chain using a K-S test and compare to the mean and - # standard deviation. - samps = sampler.get_chain(flat=True) - mu, sig = np.mean(samps, axis=0), np.std(samps, axis=0) - assert np.all(np.abs(mu) < 0.08), "Incorrect mean" - assert np.all(np.abs(sig - 1) < 0.05), "Incorrect standard deviation" - - if ndim == 1: - ks, _ = stats.kstest(samps[:, 0], "norm") - assert ks < 0.05, "The K-S test failed" - - -def _test_uniform(proposal, nwalkers=32, nsteps=2000, seed=1234): - # Set up the random number generator. - np.random.seed(seed) - - # Initialize the ensemble and proposal. - coords = np.random.rand(nwalkers, 1) - - sampler = emcee.EnsembleSampler(nwalkers, 1, normal_log_prob, - moves=proposal) - sampler.run_mcmc(coords, nsteps) - - # Check the acceptance fraction. - acc = sampler.acceptance_fraction - assert np.all((acc < 0.9) * (acc > 0.1)), \ - "Invalid acceptance fraction\n{0}".format(acc) - - # Check that the resulting chain "fails" the K-S test. - samps = sampler.get_chain(flat=True) - np.random.shuffle(samps) - ks, _ = stats.kstest(samps[::100], "uniform") - assert ks > 0.1, "The K-S test failed" diff --git a/packages/best_fit/emcee3rc2/tests/integration/test_stretch.py b/packages/best_fit/emcee3rc2/tests/integration/test_stretch.py deleted file mode 100644 index d09289f5..00000000 --- a/packages/best_fit/emcee3rc2/tests/integration/test_stretch.py +++ /dev/null @@ -1,28 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -import pytest -from emcee import moves -from .test_proposal import _test_normal, _test_uniform - -__all__ = ["test_normal_stretch", "test_uniform_stretch", - "test_nsplits_stretch"] - - -@pytest.mark.parametrize("blobs", [True, False]) -def test_normal_stretch(blobs, **kwargs): - kwargs["blobs"] = blobs - _test_normal(moves.StretchMove(), **kwargs) - - -def test_uniform_stretch(**kwargs): - _test_uniform(moves.StretchMove(), **kwargs) - - -def test_nsplits_stretch(**kwargs): - _test_normal(moves.StretchMove(nsplits=5), **kwargs) - - -def test_randomize_stretch(**kwargs): - _test_normal(moves.StretchMove(randomize_split=True), **kwargs) diff --git a/packages/best_fit/emcee3rc2/tests/integration/test_walk.py b/packages/best_fit/emcee3rc2/tests/integration/test_walk.py deleted file mode 100644 index d09888ec..00000000 --- a/packages/best_fit/emcee3rc2/tests/integration/test_walk.py +++ /dev/null @@ -1,16 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -from emcee import moves -from .test_proposal import _test_normal, _test_uniform - -__all__ = ["test_normal_walk", "test_uniform_walk"] - - -def test_normal_walk(**kwargs): - _test_normal(moves.WalkMove(s=3), **kwargs) - - -def test_uniform_walk(**kwargs): - _test_uniform(moves.WalkMove(s=3), **kwargs) diff --git a/packages/best_fit/emcee3rc2/tests/unit/__init__.py b/packages/best_fit/emcee3rc2/tests/unit/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/packages/best_fit/emcee3rc2/tests/unit/test_autocorr.py b/packages/best_fit/emcee3rc2/tests/unit/test_autocorr.py deleted file mode 100644 index f78dc2ce..00000000 --- a/packages/best_fit/emcee3rc2/tests/unit/test_autocorr.py +++ /dev/null @@ -1,49 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -import pytest -import numpy as np - -from emcee.autocorr import integrated_time, AutocorrError - - -def get_chain(seed=1234, ndim=3, N=100000): - np.random.seed(seed) - a = 0.9 - x = np.empty((N, ndim)) - x[0] = np.zeros(ndim) - for i in range(1, N): - x[i] = x[i-1] * a + np.random.rand(ndim) - return x - - -def test_1d(seed=1234, ndim=1, N=250000): - x = get_chain(seed=seed, ndim=ndim, N=N) - tau = integrated_time(x) - assert np.all(np.abs(tau - 19.0) / 19. < 0.2) - - -def test_nd(seed=1234, ndim=3, N=150000): - x = get_chain(seed=seed, ndim=ndim, N=N) - tau = integrated_time(x) - assert np.all(np.abs(tau - 19.0) / 19. < 0.2) - - -def test_too_short(seed=1234, ndim=3, N=100): - x = get_chain(seed=seed, ndim=ndim, N=N) - with pytest.raises(AutocorrError): - integrated_time(x) - tau = integrated_time(x, quiet=True) # NOQA - - -def test_autocorr_multi_works(): - np.random.seed(42) - xs = np.random.randn(16384, 2) - - # This throws exception unconditionally in buggy impl's - acls_multi = integrated_time(xs) - acls_single = np.array([integrated_time(xs[:, i]) - for i in range(xs.shape[1])]) - - assert np.all(np.abs(acls_multi - acls_single) < 2) diff --git a/packages/best_fit/emcee3rc2/tests/unit/test_backends.py b/packages/best_fit/emcee3rc2/tests/unit/test_backends.py deleted file mode 100644 index 392aee6e..00000000 --- a/packages/best_fit/emcee3rc2/tests/unit/test_backends.py +++ /dev/null @@ -1,215 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -import os -from itertools import product - -import numpy as np - -import pytest - -from emcee import backends, EnsembleSampler - -import h5py - -__all__ = ["test_backend", "test_reload"] - -all_backends = backends.get_test_backends() -other_backends = all_backends[1:] -dtypes = [ - None, - [("log_prior", float), ("mean", int)] -] - - -def normal_log_prob(params): - return -0.5 * np.sum(params**2) - - -def normal_log_prob_blobs(params): - return normal_log_prob(params), 0.1, int(5) - - -def run_sampler(backend, nwalkers=32, ndim=3, nsteps=25, seed=1234, thin_by=1, - dtype=None, blobs=True, lp=None): - if lp is None: - lp = normal_log_prob_blobs if blobs else normal_log_prob - if seed is not None: - np.random.seed(seed) - coords = np.random.randn(nwalkers, ndim) - sampler = EnsembleSampler(nwalkers, ndim, lp, - backend=backend, blobs_dtype=dtype) - sampler.run_mcmc(coords, nsteps, thin_by=thin_by) - return sampler - - -def _custom_allclose(a, b): - if a.dtype.fields is None: - assert np.allclose(a, b) - else: - for n in a.dtype.names: - assert np.allclose(a[n], b[n]) - - -def test_uninit(tmpdir): - fn = str(tmpdir.join("EMCEE_TEST_FILE_DO_NOT_USE.h5")) - if os.path.exists(fn): - os.remove(fn) - - with backends.HDFBackend(fn) as be: - run_sampler(be) - - assert os.path.exists(fn) - os.remove(fn) - - -@pytest.mark.parametrize("backend", all_backends) -def test_uninit_errors(backend): - with backend() as be: - with pytest.raises(AttributeError): - be.get_last_sample() - - for k in ["chain", "log_prob", "blobs"]: - with pytest.raises(AttributeError): - getattr(be, "get_" + k)() - - -@pytest.mark.parametrize("backend", all_backends) -def test_blob_usage_errors(backend): - with backend() as be: - run_sampler(be, blobs=True) - with pytest.raises(ValueError): - run_sampler(be, blobs=False) - - with backend() as be: - run_sampler(be, blobs=False) - with pytest.raises(ValueError): - run_sampler(be, blobs=True) - - -@pytest.mark.parametrize("backend,dtype,blobs", - product(other_backends, dtypes, [True, False])) -def test_backend(backend, dtype, blobs): - # Run a sampler with the default backend. - sampler1 = run_sampler(backends.Backend(), dtype=dtype, blobs=blobs) - - with backend() as be: - sampler2 = run_sampler(be, dtype=dtype, blobs=blobs) - - values = ["chain", "log_prob"] - if blobs: - values += ["blobs"] - else: - assert sampler1.get_blobs() is None - assert sampler2.get_blobs() is None - - # Check all of the components. - for k in values: - a = getattr(sampler1, "get_" + k)() - b = getattr(sampler2, "get_" + k)() - _custom_allclose(a, b) - - last1 = sampler1.get_last_sample() - last2 = sampler2.get_last_sample() - assert np.allclose(last1.coords, last2.coords) - assert np.allclose(last1.log_prob, last2.log_prob) - assert all(np.allclose(l1, l2) - for l1, l2 in zip(last1.random_state[1:], - last2.random_state[1:])) - if blobs: - _custom_allclose(last1.blobs, last2.blobs) - else: - assert last1.blobs is None and last2.blobs is None - - a = sampler1.acceptance_fraction - b = sampler2.acceptance_fraction - assert np.allclose(a, b), "inconsistent acceptance fraction" - - -@pytest.mark.parametrize("backend,dtype", product(other_backends, dtypes)) -def test_reload(backend, dtype): - with backend() as backend1: - run_sampler(backend1, dtype=dtype) - - # Test the state - state = backend1.random_state - np.random.set_state(state) - - # Load the file using a new backend object. - backend2 = backends.HDFBackend(backend1.filename, backend1.name, - read_only=True) - - with pytest.raises(RuntimeError): - backend2.reset(32, 3) - - assert state[0] == backend2.random_state[0] - assert all(np.allclose(a, b) - for a, b in zip(state[1:], backend2.random_state[1:])) - - # Check all of the components. - for k in ["chain", "log_prob", "blobs"]: - a = backend1.get_value(k) - b = backend2.get_value(k) - _custom_allclose(a, b) - - last1 = backend1.get_last_sample() - last2 = backend2.get_last_sample() - assert np.allclose(last1.coords, last2.coords) - assert np.allclose(last1.log_prob, last2.log_prob) - assert all(np.allclose(l1, l2) - for l1, l2 in zip(last1.random_state[1:], - last2.random_state[1:])) - _custom_allclose(last1.blobs, last2.blobs) - - a = backend1.accepted - b = backend2.accepted - assert np.allclose(a, b), "inconsistent accepted" - - -@pytest.mark.parametrize("backend,dtype", product(other_backends, dtypes)) -def test_restart(backend, dtype): - # Run a sampler with the default backend. - b = backends.Backend() - run_sampler(b, dtype=dtype) - sampler1 = run_sampler(b, seed=None, dtype=dtype) - - with backend() as be: - run_sampler(be, dtype=dtype) - sampler2 = run_sampler(be, seed=None, dtype=dtype) - - # Check all of the components. - for k in ["chain", "log_prob", "blobs"]: - a = getattr(sampler1, "get_" + k)() - b = getattr(sampler2, "get_" + k)() - _custom_allclose(a, b) - - last1 = sampler1.get_last_sample() - last2 = sampler2.get_last_sample() - assert np.allclose(last1.coords, last2.coords) - assert np.allclose(last1.log_prob, last2.log_prob) - assert all(np.allclose(l1, l2) - for l1, l2 in zip(last1.random_state[1:], - last2.random_state[1:])) - _custom_allclose(last1.blobs, last2.blobs) - - a = sampler1.acceptance_fraction - b = sampler2.acceptance_fraction - assert np.allclose(a, b), "inconsistent acceptance fraction" - - -def test_multi_hdf5(): - with backends.TempHDFBackend() as backend1: - run_sampler(backend1) - - backend2 = backends.HDFBackend(backend1.filename, name='mcmc2') - run_sampler(backend2) - chain2 = backend2.get_chain() - - with h5py.File(backend1.filename, 'r') as f: - assert set(f.keys()) == {backend1.name, 'mcmc2'} - - backend1.reset(10, 2) - assert np.allclose(backend2.get_chain(), chain2) - with pytest.raises(AttributeError): - backend1.get_chain() diff --git a/packages/best_fit/emcee3rc2/tests/unit/test_blobs.py b/packages/best_fit/emcee3rc2/tests/unit/test_blobs.py deleted file mode 100644 index 1f3695f3..00000000 --- a/packages/best_fit/emcee3rc2/tests/unit/test_blobs.py +++ /dev/null @@ -1,51 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -import pytest -import numpy as np - -from emcee import backends, EnsembleSampler - -__all__ = ["test_blob_shape"] - - -class BlobLogProb(object): - - def __init__(self, blob_function): - self.blob_function = blob_function - - def __call__(self, params): - return -0.5 * np.sum(params**2), self.blob_function(params) - - -@pytest.mark.parametrize("backend", backends.get_test_backends()) -@pytest.mark.parametrize("blob_spec", [ - (True, 5, lambda x: np.random.randn(5)), - (True, 0, lambda x: np.random.randn()), - (False, 2, lambda x: (1.0, np.random.randn(3))), - (False, 0, lambda x: "face"), - (False, 2, lambda x: (np.random.randn(5), "face")), -]) -def test_blob_shape(backend, blob_spec): - # HDF backends don't support the object type - if backend in (backends.TempHDFBackend, ) and not blob_spec[0]: - return - - with backend() as be: - np.random.seed(42) - - model = BlobLogProb(blob_spec[2]) - coords = np.random.randn(32, 3) - nwalkers, ndim = coords.shape - - sampler = EnsembleSampler(nwalkers, ndim, model, backend=be) - nsteps = 10 - - sampler.run_mcmc(coords, nsteps) - - shape = [nsteps, nwalkers] - if blob_spec[1] > 0: - shape += [blob_spec[1]] - - assert sampler.get_blobs().shape == tuple(shape) diff --git a/packages/best_fit/emcee3rc2/tests/unit/test_sampler.py b/packages/best_fit/emcee3rc2/tests/unit/test_sampler.py deleted file mode 100644 index 38775f4d..00000000 --- a/packages/best_fit/emcee3rc2/tests/unit/test_sampler.py +++ /dev/null @@ -1,186 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -import pickle -from itertools import product - -import pytest -import numpy as np - -from emcee import moves, backends, EnsembleSampler - -__all__ = ["test_shapes", "test_errors", "test_thin", "test_vectorize"] - -all_backends = backends.get_test_backends() - - -def normal_log_prob(params): - return -0.5 * np.sum(params**2) - - -@pytest.mark.parametrize("backend, moves", product( - all_backends, - [ - None, - moves.GaussianMove(0.5), - [moves.StretchMove(), moves.GaussianMove(0.5)], - [(moves.StretchMove(), 0.3), (moves.GaussianMove(0.5), 0.1)], - ]) -) -def test_shapes(backend, moves, nwalkers=32, ndim=3, nsteps=10, seed=1234): - # Set up the random number generator. - np.random.seed(seed) - - with backend() as be: - # Initialize the ensemble, moves and sampler. - coords = np.random.randn(nwalkers, ndim) - sampler = EnsembleSampler(nwalkers, ndim, normal_log_prob, - moves=moves, backend=be) - - # Run the sampler. - sampler.run_mcmc(coords, nsteps) - chain = sampler.get_chain() - assert len(chain) == nsteps, "wrong number of steps" - - tau = sampler.get_autocorr_time(quiet=True) - assert tau.shape == (ndim,) - - # Check the shapes. - assert sampler.chain.shape == (nwalkers, nsteps, ndim), \ - "incorrect coordinate dimensions" - assert sampler.get_chain().shape == (nsteps, nwalkers, ndim), \ - "incorrect coordinate dimensions" - assert sampler.lnprobability.shape == (nsteps, nwalkers), \ - "incorrect probability dimensions" - - assert sampler.acceptance_fraction.shape == (nwalkers,), \ - "incorrect acceptance fraction dimensions" - - # Check the shape of the flattened coords. - assert sampler.get_chain(flat=True).shape == \ - (nsteps * nwalkers, ndim), "incorrect coordinate dimensions" - assert sampler.get_log_prob(flat=True).shape == \ - (nsteps*nwalkers,), "incorrect probability dimensions" - - -@pytest.mark.parametrize("backend", all_backends) -def test_errors(backend, nwalkers=32, ndim=3, nsteps=5, seed=1234): - # Set up the random number generator. - np.random.seed(seed) - - with backend() as be: - # Initialize the ensemble, proposal, and sampler. - coords = np.random.randn(nwalkers, ndim) - sampler = EnsembleSampler(nwalkers, ndim, normal_log_prob, - backend=be) - - # Test for not running. - with pytest.raises(AttributeError): - sampler.chain - with pytest.raises(AttributeError): - sampler.lnprobability - - # What about not storing the chain. - sampler.run_mcmc(coords, nsteps, store=False) - with pytest.raises(AttributeError): - sampler.chain - - # Now what about if we try to continue using the sampler with an - # ensemble of a different shape. - sampler.run_mcmc(coords, nsteps, store=False) - - coords2 = np.random.randn(nwalkers, ndim+1) - with pytest.raises(ValueError): - list(sampler.run_mcmc(coords2, nsteps)) - - -def run_sampler(backend, nwalkers=32, ndim=3, nsteps=25, seed=1234, - thin=None, thin_by=1, progress=False, store=True): - np.random.seed(seed) - coords = np.random.randn(nwalkers, ndim) - sampler = EnsembleSampler(nwalkers, ndim, normal_log_prob, - backend=backend) - sampler.run_mcmc(coords, nsteps, thin=thin, thin_by=thin_by, - progress=progress, store=store) - return sampler - - -@pytest.mark.parametrize("backend", all_backends) -def test_thin(backend): - with backend() as be: - with pytest.raises(ValueError): - run_sampler(be, thin=-1) - with pytest.raises(ValueError): - run_sampler(be, thin=0.1) - thinby = 3 - sampler1 = run_sampler(None) - sampler2 = run_sampler(be, thin=thinby) - for k in ["get_chain", "get_log_prob"]: - a = getattr(sampler1, k)()[thinby-1::thinby] - b = getattr(sampler2, k)() - c = getattr(sampler1, k)(thin=thinby) - assert np.allclose(a, b), "inconsistent {0}".format(k) - assert np.allclose(a, c), "inconsistent {0}".format(k) - - -@pytest.mark.parametrize("backend,progress", - product(all_backends, [True, False])) -def test_thin_by(backend, progress): - with backend() as be: - with pytest.raises(ValueError): - run_sampler(be, thin_by=-1) - with pytest.raises(ValueError): - run_sampler(be, thin_by=0.1) - nsteps = 25 - thinby = 3 - sampler1 = run_sampler(None, nsteps=nsteps*thinby, progress=progress) - sampler2 = run_sampler(be, thin_by=thinby, progress=progress, - nsteps=nsteps) - for k in ["get_chain", "get_log_prob"]: - a = getattr(sampler1, k)()[thinby-1::thinby] - b = getattr(sampler2, k)() - c = getattr(sampler1, k)(thin=thinby) - assert np.allclose(a, b), "inconsistent {0}".format(k) - assert np.allclose(a, c), "inconsistent {0}".format(k) - assert sampler1.iteration == sampler2.iteration*thinby - - -@pytest.mark.parametrize("backend", all_backends) -def test_restart(backend): - with backend() as be: - sampler = run_sampler(be, nsteps=0) - with pytest.raises(ValueError): - sampler.run_mcmc(None, 10) - - sampler = run_sampler(be) - sampler.run_mcmc(None, 10) - - with backend() as be: - sampler = run_sampler(be, store=False) - sampler.run_mcmc(None, 10) - - -def test_vectorize(): - def lp_vec(p): - return -0.5 * np.sum(p**2, axis=1) - - np.random.seed(42) - nwalkers, ndim = 32, 3 - coords = np.random.randn(nwalkers, ndim) - sampler = EnsembleSampler(nwalkers, ndim, lp_vec, vectorize=True) - sampler.run_mcmc(coords, 10) - - assert sampler.get_chain().shape == (10, nwalkers, ndim) - - -@pytest.mark.parametrize("backend", all_backends) -def test_pickle(backend): - with backend() as be: - sampler1 = run_sampler(be) - s = pickle.dumps(sampler1, -1) - sampler2 = pickle.loads(s) - for k in ["get_chain", "get_log_prob"]: - a = getattr(sampler1, k)() - b = getattr(sampler2, k)() - assert np.allclose(a, b), "inconsistent {0}".format(k) diff --git a/packages/best_fit/emcee3rc2/tests/unit/test_state.py b/packages/best_fit/emcee3rc2/tests/unit/test_state.py deleted file mode 100644 index 22d678aa..00000000 --- a/packages/best_fit/emcee3rc2/tests/unit/test_state.py +++ /dev/null @@ -1,43 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -import numpy as np -from emcee.state import State -from emcee import EnsembleSampler - - -def test_back_compat(seed=1234): - np.random.seed(seed) - coords = np.random.randn(16, 3) - log_prob = np.random.randn(len(coords)) - blobs = np.random.randn(len(coords)) - rstate = np.random.get_state() - - state = State(coords, log_prob, blobs, rstate) - c, l, r, b = state - assert np.allclose(coords, c) - assert np.allclose(log_prob, l) - assert np.allclose(blobs, b) - assert all(np.allclose(a, b) for a, b in zip(rstate[1:], r[1:])) - - state = State(coords, log_prob, None, rstate) - c, l, r = state - assert np.allclose(coords, c) - assert np.allclose(log_prob, l) - assert all(np.allclose(a, b) for a, b in zip(rstate[1:], r[1:])) - - -def test_overwrite(seed=1234): - np.random.seed(seed) - - def ll(x): - return -0.5 * np.sum(x**2) - - nwalkers = 64 - p0 = np.random.normal(size=(nwalkers, 1)) - init = np.copy(p0) - - sampler = EnsembleSampler(nwalkers, 1, ll) - sampler.run_mcmc(p0, 10) - assert np.allclose(init, p0) diff --git a/packages/best_fit/emcee3rc2/tests/unit/test_stretch.py b/packages/best_fit/emcee3rc2/tests/unit/test_stretch.py deleted file mode 100644 index 5dd64c85..00000000 --- a/packages/best_fit/emcee3rc2/tests/unit/test_stretch.py +++ /dev/null @@ -1,39 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import division, print_function - -import warnings - -import pytest -import numpy as np - -from emcee import moves -from emcee.state import State -from emcee.model import Model - -__all__ = ["test_live_dangerously"] - - -def test_live_dangerously(nwalkers=32, nsteps=3000, seed=1234): - warnings.filterwarnings("error") - - # Set up the random number generator. - np.random.seed(seed) - state = State(np.random.randn(nwalkers, 2 * nwalkers), - log_prob=np.random.randn(nwalkers)) - model = Model( - None, - lambda x: (np.zeros(len(x)), None), - map, - np.random - ) - proposal = moves.StretchMove() - - # Test to make sure that the error is thrown if there aren't enough - # walkers. - with pytest.raises(RuntimeError): - proposal.propose(model, state) - - # Living dangerously... - proposal.live_dangerously = True - proposal.propose(model, state) From 7d4b721a87a2300a5b133b0c8c37dc3225889a36 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Sun, 29 Sep 2019 12:12:28 -0300 Subject: [PATCH 08/27] use sigma_clipped_stats func to remove outliers --- packages/aux_funcs.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/packages/aux_funcs.py b/packages/aux_funcs.py index 16d88ab2..189747f0 100644 --- a/packages/aux_funcs.py +++ b/packages/aux_funcs.py @@ -1,17 +1,16 @@ import numpy as np from scipy import stats +from astropy.stats import sigma_clipped_stats def reject_outliers(data, m=4.): """ Outlier rejection. - Source: https://stackoverflow.com/a/16562028/1391441 """ - d = np.abs(data - np.median(data)) - mdev = np.median(d) - s = d / (mdev if mdev else 1.) - return data[s < m] + mean, median, std = sigma_clipped_stats(data) + msk = (data > mean - m * std) & (data < mean + m * std) + return data[msk] def kde1D(data): From 2899a9c3c875f1eb4440b68d84e22393e3e24b57 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Sun, 29 Sep 2019 12:19:30 -0300 Subject: [PATCH 09/27] dont use 'auto' in np.histogram, causes MemoryError if outliers are present --- packages/out/mp_kinem_plx.py | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/packages/out/mp_kinem_plx.py b/packages/out/mp_kinem_plx.py index 568f1430..8541302e 100644 --- a/packages/out/mp_kinem_plx.py +++ b/packages/out/mp_kinem_plx.py @@ -122,7 +122,8 @@ def plx_vs_mag( ax.grid(b=True, which='major', color='gray', linestyle='--', lw=.5, zorder=1) - ax.set_title("Plx clip " + r"$(med\pm2\sigma,\;N={})$".format( + # HARDCODED in plx_analysis: 3 sigma outlier rejection + ax.set_title("Plx clip " + r"$(med\pm3\sigma,\;N={})$".format( len(plx_clp)), fontsize=9) plt.xlabel('Plx [mas]', fontsize=12) plt.ylabel(y_ax, fontsize=12) @@ -197,11 +198,9 @@ def pms_bys_params( ax1 = plt.subplot(gs[2:4, 0:2]) plt.xlabel(r"$Plx_{Bay}$") - plx_outlr = reject_outliers(plx_samples.flatten()) - # Obtain the bin values and edges using numpy - hist, bin_edges = np.histogram(1. / plx_outlr, bins='auto') - if len(bin_edges) > 25: - hist, bin_edges = np.histogram(1. / plx_outlr, bins=20) + plx_no_outlr = reject_outliers(plx_samples.flatten()) + # Obtain the bin edges using numpy + hist, bin_edges = np.histogram(1. / plx_no_outlr, bins=25) # Plot bars with the proper positioning, height, and width. plt.bar( (bin_edges[1:] + bin_edges[:-1]) * .5, hist / float(hist.max()), @@ -230,20 +229,29 @@ def pms_bys_params( # Traceplot ax = plt.subplot(gs[2:3, 2:6]) N_tot = plx_samples.shape[0] - plt.plot(1. / plx_samples, c='k', lw=.8, ls='-', alpha=0.5) - # HARDCODED: 25% of trace is burned - plt.axvline(x=.25 * N_tot, linestyle=':', color='r', zorder=4) + trace = 1. / plx_samples + plt.plot(trace, c='k', lw=.8, ls='-', alpha=0.5) + # HARDCODED in plx_analysis: 25% of trace is burned + Nburn = .25 * N_tot + plt.axvline(x=Nburn, linestyle=':', color='r', zorder=4) # 16th and 84th percentiles + median. plt.axhline(y=1. / plx_Bys[0], linestyle=':', color='orange', zorder=4) plt.axhline(y=1. / plx_Bys[1], linestyle=':', color='blue', zorder=4) plt.axhline(y=1. / plx_Bys[2], linestyle=':', color='orange', zorder=4) plt.xlabel("Steps") plt.ylabel(r"$Plx_{Bay}$") + + # Use the last 10% of the chains. + N = int(trace.shape[0] * .1) + std = np.std(trace[-N:]) + pmin, pmax = np.min(trace[-N:]), np.max(trace[-N:]) + ymin, ymax = pmin - std, pmax + std + ax.set_ylim(ymin, ymax) ax.set_xlim(0, N_tot) # Tau plt.subplot(gs[3:4, 2:4]) - # HARDCODED: store samples every 10 steps + # HARDCODED in plx_analysis: store samples every 10 steps plt.plot( 10 * np.arange(plx_tau_autocorr.size), plx_tau_autocorr, label=r"$N_{{ESS}}\approx${:.0f}".format(plx_ess)) From 2af637df500f40af03f915a171e21bcbc12eaad4 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Sun, 29 Sep 2019 14:13:08 -0300 Subject: [PATCH 10/27] minor minor --- packages/data_analysis/luminosity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/data_analysis/luminosity.py b/packages/data_analysis/luminosity.py index 8b9ae9ec..9d91c481 100644 --- a/packages/data_analysis/luminosity.py +++ b/packages/data_analysis/luminosity.py @@ -52,7 +52,7 @@ def main(clp, **kwargs): np.array([0.]))) else: print(" WARNING: no field regions defined. Luminosity function\n" - " is not cleaned from field star contamination.") + " is not cleaned from field star contamination") # Pass dummy lists. x_fl, y_fl = [], [] From a107b92366316ade697a2573b61930e0ab9bc43b Mon Sep 17 00:00:00 2001 From: Gabriel Date: Sun, 29 Sep 2019 14:13:55 -0300 Subject: [PATCH 11/27] limit N for faster PMs KDE estimation --- packages/data_analysis/pms_analysis.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/packages/data_analysis/pms_analysis.py b/packages/data_analysis/pms_analysis.py index c28a51a0..54d5b551 100644 --- a/packages/data_analysis/pms_analysis.py +++ b/packages/data_analysis/pms_analysis.py @@ -241,6 +241,18 @@ def kde_2d(xarr, xsigma, yarr, ysigma, Nstd=3, grid_dens=50): Take an array of x,y data with their errors, create a grid of points in x,y and return the 2D KDE density map. ''' + # For better performance, use max 10000 random stars when estimating + # the KDE. + N_max = 10000 + if xarr.size > N_max: + print((" WARNING: used {} stars to estimate the PMs KDE\n" + " instead of the {} total").format(N_max, xarr.size)) + xarr, xsigma, yarr, ysigma =\ + np.random.choice(xarr, N_max, replace=False),\ + np.random.choice(xsigma, N_max, replace=False),\ + np.random.choice(yarr, N_max, replace=False),\ + np.random.choice(ysigma, N_max, replace=False) + _, x_median, x_std = sigma_clipped_stats(xarr) _, y_median, y_std = sigma_clipped_stats(yarr) From 5110f8ab4d6d6c3f67085ae9a0f4b6c94ee65e25 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Sun, 29 Sep 2019 15:05:36 -0300 Subject: [PATCH 12/27] fixed a few things, particularly the 'pos0' generation that caused issues --- packages/data_analysis/plx_analysis.py | 80 +++++++++++++++----------- 1 file changed, 46 insertions(+), 34 deletions(-) diff --git a/packages/data_analysis/plx_analysis.py b/packages/data_analysis/plx_analysis.py index 00e9492c..2f3402c6 100644 --- a/packages/data_analysis/plx_analysis.py +++ b/packages/data_analysis/plx_analysis.py @@ -9,14 +9,14 @@ def main( clp, plx_bayes_flag, plx_offset, plx_chains, plx_runs, flag_plx_mp, - flag_make_plot, **kwargs): + flag_make_plot, outlr_std=3., **kwargs): """ Bayesian parallax distance using the Bailer-Jones (2015) model with the 'shape parameter' marginalized. Hardcoded choices: - * 2 sigma outliers are rejected + * outlr_std sigma outliers are rejected * MPs are used * Bayesian prior is a Gaussian with a fixed standard deviation """ @@ -38,9 +38,9 @@ def main( if plx_flag_clp: print("Processing parallaxes") - # Reject 2\sigma outliers. - max_plx, min_plx = np.nanmedian(plx) + 2. * np.nanstd(plx),\ - np.nanmedian(plx) - 2. * np.nanstd(plx) + # Reject outlr_std*\sigma outliers. + max_plx, min_plx = np.nanmedian(plx) + outlr_std * np.nanstd(plx),\ + np.nanmedian(plx) - outlr_std * np.nanstd(plx) # Suppress Runtimewarning issued when 'plx' contains 'nan' values. with np.warnings.catch_warnings(): @@ -111,43 +111,24 @@ def plxBayes( ndim, nwalkers, nruns = 1, plx_chains, plx_runs print(" Bayesian Plx model ({} runs)".format(nruns)) - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - - # Define the 'r_i' values used to evaluate the integral. - int_max = 20. - N = int(int_max / 0.01) - x = np.linspace(.1, int_max, N).reshape(-1, 1) - B1 = ((plx_clp - (1. / x)) / e_plx_clp)**2 - B2 = (np.exp(-.5 * B1) / e_plx_clp) - - # Use DE to estimate the ML - def DEdist(model): - return -lnlike(model, x, B2, mp_clp) - bounds = [[0., 20.]] - result = DE(DEdist, bounds, popsize=20, maxiter=100) - # print(result) + # DE initial mu position + mu_p = DE_mu_sol(plx_clp, e_plx_clp, mp_clp) - # Prior parameters. - mu_p = result.x # Define the 'r_i' values used to evaluate the integral. int_max = mu_p + 5. - N = int(int_max / 0.01) - x = np.linspace(.1, int_max, N).reshape(-1, 1) - B1 = ((plx_clp - (1. / x)) / e_plx_clp)**2 - B2 = (np.exp(-.5 * B1) / e_plx_clp) + x, B2 = r_iVals(int_max, plx_clp, e_plx_clp) # emcee sampler sampler = ensemble.EnsembleSampler( nwalkers, ndim, lnprob, args=(x, B2, mp_clp, mu_p)) # Random initial guesses. - # pos0 = [np.random.uniform(0., 1., ndim) for i in range(nwalkers)] # Ball of initial guesses around 'mu_p' - pos0 = [mu_p + .5 * np.random.normal() for i in range(nwalkers)] + pos0 = np.clip( + np.array([[mu_p + .05 * np.random.normal()] for i in range(nwalkers)]), + a_min=0., a_max=None) tau_index, autocorr_vals, afs = 0, np.empty(nruns), np.empty(nruns) - try: with warnings.catch_warnings(): warnings.simplefilter("ignore") @@ -185,7 +166,7 @@ def DEdist(model): tau_mean = np.mean(sampler.get_autocorr_time(tol=0)) plx_ess = samples.size / tau_mean - # For plotting + # For plotting, (nsteps, nchains, ndim) plx_samples = sampler.get_chain()[:, :, 0] print("Bayesian plx estimated: " + @@ -194,17 +175,48 @@ def DEdist(model): except Exception as e: print(e) print("\n ERROR: could not process Plx data with emcee") - plx_samples, plx_Bys, plx_bayes_flag_clp, plx_ess = [], np.array([]),\ - False, np.nan + plx_samples, plx_Bys, plx_bayes_flag_clp, plx_ess, tau_autocorr,\ + mean_afs = [], np.array([]), False, np.nan, np.nan, np.nan return plx_samples, plx_Bys, plx_bayes_flag_clp, tau_autocorr,\ mean_afs, plx_ess +def DE_mu_sol(plx_clp, e_plx_clp, mp_clp, int_max=20., psize=20, maxi=100): + """ + Use the Differential Evolution algorithm to approximate the best solution + used as the mean of the prior. + """ + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + + # Define the 'r_i' values used to evaluate the integral. + x, B2 = r_iVals(int_max, plx_clp, e_plx_clp) + + # Use DE to estimate the ML + def DEdist(model): + return -lnlike(model, x, B2, mp_clp) + bounds = [[0., 20.]] + result = DE(DEdist, bounds, popsize=psize, maxiter=maxi) + + return result.x[0] + + +def r_iVals(int_max, plx, e_plx): + """ + The 'r_i' values used to evaluate the integral. + """ + N = int(int_max / 0.01) + x = np.linspace(.1, int_max, N).reshape(-1, 1) + B1 = ((plx - (1. / x)) / e_plx)**2 + B2 = (np.exp(-.5 * B1) / e_plx) + return x, B2 + + def lnprob(mu, x, B2, MPs, mu_p): lp = lnprior(mu, mu_p) if np.isinf(lp): - return lp + return -np.inf return lp + lnlike(mu, x, B2, MPs) From 3b4d1fb061c1b614aa6bd05900e8261274e045c6 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Sun, 29 Sep 2019 21:42:40 -0300 Subject: [PATCH 13/27] minor, same parameter name for consistency --- packages/synth_clust/binarity.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/packages/synth_clust/binarity.py b/packages/synth_clust/binarity.py index 6eec62c7..eb83ca4e 100644 --- a/packages/synth_clust/binarity.py +++ b/packages/synth_clust/binarity.py @@ -24,7 +24,7 @@ def mag_combine(m1, m2): def binarGen( - binar_fracs, N_interp, mags_theor, cols_theor, mags_cols_theor, + binar_fracs, N_mass_interp, mags_theor, cols_theor, mags_cols_theor, extra_pars, bin_mass_ratio): ''' Called by isoch_params(). @@ -51,7 +51,7 @@ def binarGen( # All theoretical isochrones are interpolated with the same length, # assign unique binarity probabilities to each star randomly. - unq_b_probs = np.arange(N_interp) / float(N_interp) + unq_b_probs = np.arange(N_mass_interp) / float(N_mass_interp) mags_binar, cols_binar, probs_binar, mass_binar = [], [], [], [] # For each metallicity defined. From 399c0faf7197a3f909d3c74006154aae46b2a2fc Mon Sep 17 00:00:00 2001 From: Gabriel Date: Mon, 30 Sep 2019 12:34:23 -0300 Subject: [PATCH 14/27] in place for #439 --- packages/best_fit/zaInterp.py | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/packages/best_fit/zaInterp.py b/packages/best_fit/zaInterp.py index 39b8c018..06e8fa69 100644 --- a/packages/best_fit/zaInterp.py +++ b/packages/best_fit/zaInterp.py @@ -63,9 +63,10 @@ def main(theor_tracks, fundam_params, varIdxs, model): model_proper.append(par[0]) j += 1 - # # Minimum and maximum initial mass for each of the four isochrones. - # mmin = np.min(isochs[:, -6, :], axis=1) - # mmax = np.max(isochs[:, -6, :], axis=1) + # TODO: IN PLACE FOR #439 + # # If (z, a) are both fixed, just return the single processed isochrone + # if ml == al == mh == ah == 0: + # return theor_tracks[ml][al], model_proper # Values of the four points in the (z, age) grid that contain the model # value (z_model, a_model) @@ -77,10 +78,24 @@ def main(theor_tracks, fundam_params, varIdxs, model): # to the grid (z, age) points. with warnings.catch_warnings(): warnings.simplefilter("ignore") + # Inverse distances between the (z, a) points in the 'model', and the # four points in the (z, a) grid that contain the model point. # Fast euclidean distance: https://stackoverflow.com/a/47775357/1391441 a_min_b = np.array([(z_model, a_model)]) - pts + + # TODO: IN PLACE FOR #439 + # # If the model has a 0. distance in (z,a) to the closest isochrone, + # # then just return that isochrone. + # try: + # # Which one is faster? + # idx = np.nonzero(a_min_b == 0.)[0][0] + # idx = np.where(a_min_b == 0.)[0][0] + # return theor_tracks[pts[idx][0]][pts[idx][1]], model_proper + # except IndexError: + # pass + + # Inverse distance. inv_d = 1. / np.sqrt(np.einsum('ij,ij->i', a_min_b, a_min_b)) # Order: (z1, a1), (z1, a2), (z2, a1), (z2, a2) From 8f782e712bb42f3db0fe7cd770dd38ca47b9b0f0 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Mon, 30 Sep 2019 15:55:20 -0300 Subject: [PATCH 15/27] close #434 --- packages/best_fit/best_fit_synth_cl.py | 3 +- packages/best_fit/bf_common.py | 51 +++---- packages/best_fit/bootstrap.py | 9 +- packages/best_fit/ptemcee_algor.py | 10 +- packages/best_fit/zaInterp.py | 5 - packages/check/params_match.py | 16 -- packages/check/read_met_files.py | 2 +- packages/defvals/params_input.dat | 32 ++-- packages/inp/input_params.py | 12 +- packages/inp/isoch_params.py | 5 +- packages/inp/met_ages_values.py | 25 +--- packages/synth_clust/imf.py | 172 +++------------------- packages/synth_clust/mass_distribution.py | 2 +- packages/synth_clust/synth_cl_plot.py | 4 +- 14 files changed, 86 insertions(+), 262 deletions(-) diff --git a/packages/best_fit/best_fit_synth_cl.py b/packages/best_fit/best_fit_synth_cl.py index cf1cba47..efec11ab 100644 --- a/packages/best_fit/best_fit_synth_cl.py +++ b/packages/best_fit/best_fit_synth_cl.py @@ -34,8 +34,7 @@ def main(clp, pd): # Obtain mass distribution using the selected IMF. We run it once # because the array only depends on the IMF selected. - st_dist_mass = imf.main( - pd['IMF_name'], pd['m_high'], pd['fundam_params'][4]) + st_dist_mass = imf.main(pd['IMF_name'], pd['fundam_params'][4]) # Store the number of defined filters and colors. N_fc = [len(pd['filters']), len(pd['colors'])] diff --git a/packages/best_fit/bf_common.py b/packages/best_fit/bf_common.py index 51ba6696..4d457d10 100644 --- a/packages/best_fit/bf_common.py +++ b/packages/best_fit/bf_common.py @@ -136,31 +136,32 @@ def synthClust(fundam_params, varIdxs, model, synthcl_args): R_V, ext_coefs, N_fc, cmpl_rnd, err_rnd, model_proper) -def discreteParams(fundam_params, varIdxs, chains_nruns, pushidxs): - """ - Push values in each chain for each discrete parameter in the 'pushidxs' - list to the closest grid value. - - chains_nruns.shape: (runs, nwalkers, ndim) - """ - params, j = [], 0 - for i, par in enumerate(fundam_params): - p = np.array(par) - # If this parameter is one of the 'free' parameters. - if i in varIdxs: - # If it is the parameter metallicity, age or mass. - if i in pushidxs: - pc = chains_nruns.T[j] - chains = [] - for c in pc: - chains.append( - p[abs(c[None, :] - p[:, None]).argmin(axis=0)]) - params.append(chains) - else: - params.append(chains_nruns.T[j]) - j += 1 - - return np.array(params).T +# DEPRECATED 24-09-2019 +# def discreteParams(fundam_params, varIdxs, chains_nruns, pushidxs): +# """ +# Push values in each chain for each discrete parameter in the 'pushidxs' +# list to the closest grid value. + +# chains_nruns.shape: (runs, nwalkers, ndim) +# """ +# params, j = [], 0 +# for i, par in enumerate(fundam_params): +# p = np.array(par) +# # If this parameter is one of the 'free' parameters. +# if i in varIdxs: +# # If it is the parameter metallicity, age or mass. +# if i in pushidxs: +# pc = chains_nruns.T[j] +# chains = [] +# for c in pc: +# chains.append( +# p[abs(c[None, :] - p[:, None]).argmin(axis=0)]) +# params.append(chains) +# else: +# params.append(chains_nruns.T[j]) +# j += 1 + +# return np.array(params).T def rangeCheck(model, ranges, varIdxs): diff --git a/packages/best_fit/bootstrap.py b/packages/best_fit/bootstrap.py index e697b4d7..df1b17d7 100644 --- a/packages/best_fit/bootstrap.py +++ b/packages/best_fit/bootstrap.py @@ -6,7 +6,7 @@ from . import obs_clust_prepare from .. import update_progress from .bf_common import random_population, varPars, modeKDE, fillParams,\ - r2Dist, closeSol + r2Dist def main( @@ -120,8 +120,6 @@ def main( isoch_fit_params['mode_sol'] = [np.nan] * 6 isoch_fit_params['param_r2'] = [np.nan] * 6 - isoch_fit_params['map_sol'] = closeSol( - pd['fundam_params'], isoch_fit_params['map_sol'], [4]) isoch_fit_params['btstrp_t'] = t.time() - btstrp_start_t isoch_fit_params['bf_elapsed'] = t.time() - start_t isoch_fit_params['N_total'] = float( @@ -213,11 +211,6 @@ def getSols(pd, isoch_fit_params): mode_boot_sol.append(par[0]) j += 1 - # Push Mass to grid values. - mean_boot_sol = closeSol(pd['fundam_params'], mean_boot_sol, [4]) - median_boot_sol = closeSol(pd['fundam_params'], median_boot_sol, [4]) - mode_boot_sol = closeSol(pd['fundam_params'], mode_boot_sol, [4]) - isoch_fit_params['mean_sol'] = mean_boot_sol isoch_fit_params['median_sol'] = median_boot_sol isoch_fit_params['mode_sol'] = mode_boot_sol diff --git a/packages/best_fit/ptemcee_algor.py b/packages/best_fit/ptemcee_algor.py index d4e1bb3a..84b16d60 100644 --- a/packages/best_fit/ptemcee_algor.py +++ b/packages/best_fit/ptemcee_algor.py @@ -6,7 +6,7 @@ # from .emcee3rc2 import autocorr from .mcmc_convergence import convergenceVals from .bf_common import initPop, varPars, synthClust, rangeCheck, fillParams,\ - closeSol, discreteParams, r2Dist, modeKDE # , thinChain + r2Dist, modeKDE # , thinChain from .ptemcee import sampler @@ -183,14 +183,11 @@ def main( bi_steps = int(nburn_pt * cold_chain.shape[0]) # chains_nruns.shape: (bi_steps, nchains, ndim) chains_nruns = cold_chain[:bi_steps] - # The Mass parameter is not interpolated, use its grid values. - chains_nruns = discreteParams(fundam_params, varIdxs, chains_nruns, [4]) # pars_chains_bi.shape: (ndim, nchains, bi_steps) pars_chains_bi = chains_nruns.T # After burn-in chains_nruns = cold_chain[bi_steps:] - chains_nruns = discreteParams(fundam_params, varIdxs, chains_nruns, [4]) pars_chains = chains_nruns.T # Convergence parameters. @@ -217,11 +214,6 @@ def main( mode_sol = fillParams(fundam_params, varIdxs, mode_sol) median_sol = fillParams(fundam_params, varIdxs, median_sol) - # Push Mass value to grid value for mean, map, and mode solutions. - mean_sol = closeSol(fundam_params, mean_sol, [4]) - mode_sol = closeSol(fundam_params, mode_sol, [4]) - map_sol = closeSol(fundam_params, map_sol, [4]) - # Total number of values used to estimate the parameter's distributions. N_total = mcmc_trace.shape[-1] diff --git a/packages/best_fit/zaInterp.py b/packages/best_fit/zaInterp.py index 06e8fa69..7867fcb4 100644 --- a/packages/best_fit/zaInterp.py +++ b/packages/best_fit/zaInterp.py @@ -46,11 +46,6 @@ def main(theor_tracks, fundam_params, varIdxs, model): al = ah - 1 model_proper.append(par[ah]) a_model = model[i - j] - # # If it is the parameter mass. - # elif i == 4: - # # Select the closest value in the array of allowed values. - # model_proper.append(min( - # par, key=lambda x: abs(x - model[i - j]))) else: model_proper.append(model[i - j]) else: diff --git a/packages/check/params_match.py b/packages/check/params_match.py index 58710a83..7ca68cc3 100644 --- a/packages/check/params_match.py +++ b/packages/check/params_match.py @@ -4,10 +4,6 @@ def check(mypath, pd): - # par_ranges, lkl_method, lkl_methods, lkl_binning, - # evol_track, max_mag, IMF_name, R_V, N_mass, bin_mr, bin_methods, - # lkl_weight, bin_weights, cmd_evol_tracks, iso_paths, imf_funcs, hmax, - # inst_packgs_lst, ntemps, nwalkers_ptm, nburn_ptm """ Check all parameters related to the search for the best synthetic cluster match. @@ -198,18 +194,6 @@ def checkSynthClustParams(pd): sys.exit("ERROR: Name of IMF ({}) is incorrect.".format( pd['IMF_name'])) - m_range = pd['par_ranges'][4] - # Check mass range. - try: - if pd['N_mass'] > 100 and m_range[1] > 1e5: - print(" WARNING: the number of masses defined is > 100 and\n" - " the max mass is large ({:.0f}). This could cause\n" - " memory issues when sampling the IMF.\n".format( - m_range[-1])) - except IndexError: - # Single value - pass - if not 0. <= pd['bin_mr'] <= 1.: sys.exit("ERROR: Binary mass ratio set ('{}') is out of\n" "boundaries. Please select a value in the range [0., 1.]". diff --git a/packages/check/read_met_files.py b/packages/check/read_met_files.py index 107cd088..cee124c3 100644 --- a/packages/check/read_met_files.py +++ b/packages/check/read_met_files.py @@ -22,7 +22,7 @@ def check_get(pd): params_values, met_f_filter, met_values, age_values, met_vals_all,\ age_vals_all = met_ages_values.main( pd['iso_paths'], pd['best_fit_algor'], pd['par_ranges'], - pd['za_steps'], pd['N_mass']) + pd['za_steps']) except Exception: print(traceback.format_exc()) sys.exit("\nERROR: error storing metallicity files.") diff --git a/packages/defvals/params_input.dat b/packages/defvals/params_input.dat index 42dc7838..8b73d107 100644 --- a/packages/defvals/params_input.dat +++ b/packages/defvals/params_input.dat @@ -475,31 +475,27 @@ LK dolphin knuth mean # evol_track colibri ET PAR12 No # -# Input [z, log(age)] steps to read the stored isochrones. +# * z, log(age): [floats] +# Input [z, log(age)] steps to read the stored isochrones. # -# z log(age) -PS 0.005 0.5 +# * N_interp: [int] +# Number of extra values to interpolate into the isochrones. +# +# z log(age) N_interp +PS 0.005 0.5 1000 # IMF used to populate the synthetic clusters. # # * IMF: [selection] # -# - kroupa_2002: Kroupa (2002) 295.82K, Eq. (2) & (3) +# - kroupa_2002 : Kroupa (2002) 295.82K, Eq. (2) & (3) # - chabrier_2001_exp: Chabrier (2001) 554.1274C, Eq (8) # - chabrier_2001_log: Chabrier (2001) 554.1274C, Eq (7) -# - kroupa_1993: Kroupa et al. (1993) 262.545K, Eq. (13) -# -# * m_high: [float] -# Upper mass value for the IMF (in solar masses). -# -# * N_interp: [int] -# Number of extra values to interpolate into the isochrones. -# -# * N_mass: [int] -# Number of masses defined within the mass range. +# - kroupa_1993 : Kroupa et al. (1993) 262.545K, Eq. (13) +# - salpeter_1955 : Salpeter (1955) # -# IMF m_high N_interp N_mass -MF kroupa_2002 150 2000 100 +# IMF +MF kroupa_2002 # Binaries mass ratio. # @@ -508,11 +504,11 @@ MF kroupa_2002 150 2000 100 # 'Binary fraction' parameter is fixed to '0.' # # min_mass_ratio -BI_m .7 +BR .7 # Ratio of total to selective absorption (R_V). # -RV 3.1 +RV 3.1 # Maximum magnitude cut. # diff --git a/packages/inp/input_params.py b/packages/inp/input_params.py index 9aa5d552..b0a3a053 100644 --- a/packages/inp/input_params.py +++ b/packages/inp/input_params.py @@ -216,13 +216,11 @@ def main(mypath, pars_f_path): evol_track = str(reader[1]) colibri = str(reader[2]) elif reader[0] == 'PS': - za_steps = list(map(float, reader[1:])) + za_steps = list(map(float, reader[1:3])) + N_mass_interp = int(reader[3]) elif reader[0] == 'MF': IMF_name = str(reader[1]) - m_high = float(reader[2]) - N_mass_interp = int(reader[3]) - N_mass = float(reader[4]) - elif reader[0] == 'BI_m': + elif reader[0] == 'BR': bin_mr = float(reader[1]) elif reader[0] == 'RV': R_V = float(reader[1]) @@ -323,8 +321,8 @@ def main(mypath, pars_f_path): 'lkl_weight': lkl_weight, # Synthetic cluster parameters 'evol_track': evol_track, 'za_steps': za_steps, - 'max_mag': max_mag, 'IMF_name': IMF_name, 'm_high': m_high, - 'N_mass_interp': N_mass_interp, 'N_mass': N_mass, 'R_V': R_V, + 'max_mag': max_mag, 'IMF_name': IMF_name, + 'N_mass_interp': N_mass_interp, 'R_V': R_V, 'bin_mr': bin_mr, # parameters ranges 'par_ranges': par_ranges, diff --git a/packages/inp/isoch_params.py b/packages/inp/isoch_params.py index 7c40d8b1..a3ba6d52 100644 --- a/packages/inp/isoch_params.py +++ b/packages/inp/isoch_params.py @@ -118,9 +118,8 @@ def main(met_f_filter, age_values, cmd_evol_tracks, evol_track, bin_mr, # the final isochrones. plot_isoch_data = np.concatenate((mags_theor, cols_theor), axis=2) - lens = [len(_) for _ in fundam_params] - print("\nGrid values: {} [z], {} [log(age)], {} [Mass]".format( - lens[0], lens[1], lens[4])) + print("\nGrid values: {} [z], {} [log(age)]".format( + len(fundam_params[0]), len(fundam_params[1]))) # # In place for #415 # filename = 'temp.memmap' diff --git a/packages/inp/met_ages_values.py b/packages/inp/met_ages_values.py index 1e3a92d3..b230fc3a 100644 --- a/packages/inp/met_ages_values.py +++ b/packages/inp/met_ages_values.py @@ -7,7 +7,7 @@ from . import isochs_format -def main(iso_paths, best_fit_algor, par_ranges, za_steps, N_mass): +def main(iso_paths, best_fit_algor, par_ranges, za_steps): ''' Obtain the correct metallicities and ages used by the code. ''' @@ -24,7 +24,7 @@ def main(iso_paths, best_fit_algor, par_ranges, za_steps, N_mass): age_vals_all = CMDAges(met_files[0][0]) # Get parameters ranges. - params_values = getParamVals(best_fit_algor, par_ranges, za_steps, N_mass) + params_values = getParamVals(best_fit_algor, par_ranges, za_steps) # Match values in metallicity and age ranges given by the user, with # those available in the theoretical isochrones. @@ -88,7 +88,7 @@ def CMDAges(met_file): return isoch_a -def getParamVals(best_fit_algor, par_ranges, za_steps, N_mass): +def getParamVals(best_fit_algor, par_ranges, za_steps): ''' Obtain parameter ranges to be used by the selected best fit method. ''' @@ -109,20 +109,11 @@ def getParamVals(best_fit_algor, par_ranges, za_steps, N_mass): if i in (0, 1): p_rang = np.arange(param[0], param[1], za_steps[i]) - # If processing either (E_BV, dm, b_fr) parameter, just store the - # limits since these are *continuous* variables. - elif i in (2, 3, 5): + # If processing either (E_BV, dm, Mass,b_fr) parameter, just store + # the limits since these are *continuous* variables. + elif i in (2, 3, 4, 5): p_rang = np.asarray(param) - # If processing the Mass, use `np.linspace()` instead of - # `np.arange()` since this *discrete* variable does not use steps - # but a total number of allowed values. - elif i == 4: - p_rang = np.linspace(param[0], param[1], N_mass) - # Round to integer and remove possible duplicated elements - p_rang = np.array(list(set(np.round(p_rang, 0)))) - p_rang.sort() - # Skip if array is empty. Checker will catch this. if p_rang.size: # Add max value if not present. Check this way to avoid @@ -162,8 +153,8 @@ def match_ranges(met_vals_all, met_files, age_vals_all, z_range, a_range): age_values.append(age) if not age_values: - txt = "No age value read:\n\n{}\n\ncould be matched to the ages" +\ - " given as input:\n\n{}" + txt = "None of the log(age) values read:\n\n{}\n\ncould be matched" +\ + " to the ages given as input:\n\n{}" sys.exit(txt.format(age_vals_all, a_range)) return met_f_filter, met_values, age_values diff --git a/packages/synth_clust/imf.py b/packages/synth_clust/imf.py index 4b3dbc43..68e87030 100644 --- a/packages/synth_clust/imf.py +++ b/packages/synth_clust/imf.py @@ -4,7 +4,7 @@ from scipy.interpolate import interp1d -def main(IMF_name, m_high, masses): +def main(IMF_name, masses, m_high=150.): """ Returns the number of stars per interval of mass for the selected IMF. @@ -12,11 +12,10 @@ def main(IMF_name, m_high, masses): ---------- IMF_name : str Name of the IMF to be used. + masses: array + Array of floats containing the range of masses defined. m_high : float Maximum mass value to be sampled. - masses: array - Array of floats containing the masses defined within the limiting - values and with the given mass step. Returns ------- @@ -29,110 +28,52 @@ def main(IMF_name, m_high, masses): print("Sampling selected IMF ({})".format(IMF_name)) - # Low mass limits are defined for each IMF to avoid numerical - # issues when integrating. + # Low mass limits for each IMF. Defined slightly larger to avoid sampling + # issues. imfs_dict = { - 'chabrier_2001_exp': 0.01, 'chabrier_2001_log': 0.01, - 'kroupa_1993': 0.081, 'kroupa_2002': 0.011, 'salpeter_1955': 0.3} + 'chabrier_2001_exp': 0.011, 'chabrier_2001_log': 0.011, + 'kroupa_1993': 0.081, 'kroupa_2002': 0.011, 'salpeter_1955': 0.31} # IMF low mass limit. m_low = imfs_dict[IMF_name] - # IMF mass interpolation step. - mass_step = 0.05 - - ################## - # Obtain normalization constant. This is equivalent to 'k/M_total' in - # Eq. (7) of Popescu & Hanson 2009 (138:1724-1740; PH09). - # For m_high > 100 Mo the differences in the resulting normalization - # constant are negligible. This is because th IMF drops very rapidly for - # high masses. - norm_const = 1. / quad(integral_IMF_M, m_low, m_high, args=(IMF_name))[0] - - # Obtain number of stars in each mass interval. Equivalent to the upper - # fraction of Eq. (8) in PH09, without the total mass. - mass_up, N_dist = [], [] - m_upper = m_low - while m_upper < m_high: - m_lower = m_upper - m_upper += mass_step - # Number of stars in the (m_lower, m_upper) interval. - N_stars = quad(integral_IMF_N, m_lower, m_upper, args=(IMF_name))[0] - mass_up.append(m_upper) - N_dist.append(N_stars) - - # Normalize number of stars by constant. - N_dist = np.asarray(N_dist) * norm_const - - # st_dist_mass = massDist(m_low, mass_up, N_dist, masses) - ################## - - ################## # Obtain normalization constant (k = \int_{m_low}^{m_up} \xi(m) dm). This # makes the IMF behave like a PDF. - norm_const = quad(integral_IMF_N, m_low, m_high, args=(IMF_name))[0] + norm_const = quad(IMF_func, m_low, m_high, args=(IMF_name))[0] + + # IMF mass interpolation step and grid values. + mass_step = 0.05 + mass_values = np.arange(m_low, m_high, mass_step) # The CDF is defined as: $F(m)= \int_{m_low}^{m} PDF(m) dm$ # Sample the CDF - mass_values = np.arange(m_low, m_high, mass_step) CDF_samples = [] for m in mass_values: - CDF_samples.append(quad(integral_IMF_N, m_low, m, args=(IMF_name))[0]) + CDF_samples.append(quad(IMF_func, m_low, m, args=(IMF_name))[0]) + + # Normalize values CDF_samples = np.array(CDF_samples) / norm_const + CDF_min, CDF_max = CDF_samples.min(), CDF_samples.max() + # Inverse CDF inv_cdf = interp1d(CDF_samples, mass_values) - CDF_min, CDF_max = CDF_samples.min(), CDF_samples.max() def sampled_inv_cdf(N): mr = np.random.rand(N) mr = mr[(mr >= CDF_min) & (mr <= CDF_max)] return inv_cdf(mr) - # sampled_IMF = sampled_inv_cdf(10000) - - # st_dist_mass2 = (sampled_IMF, np.cumsum(sampled_IMF)) - ################## - - def return_intersection(hist_1, hist_2): - minima = np.minimum(hist_1, hist_2) - intersection = np.true_divide(np.sum(minima), np.sum(hist_2)) - return intersection + # Sample in chunks of 100 stars until the maximum defined mass is reached. + mass_samples = [] + while np.sum(mass_samples) < masses[-1]: + mass_samples += sampled_inv_cdf(100).tolist() + sampled_IMF = np.array(mass_samples) - import matplotlib.pyplot as plt - dd = [] - for mass in np.random.uniform(100., 1000., 1000): - data_1 = massDist(m_low, mass_up, N_dist, [mass]) - data_1 = data_1[mass] - sampled_IMF = sampled_inv_cdf(10000) - data_2 = sampled_IMF[:np.searchsorted(np.cumsum(sampled_IMF), mass)] - hist_1, _ = np.histogram(data_1, bins=100, range=[0., 150]) - hist_2, _ = np.histogram(data_2, bins=100, range=[0., 150]) - rr = return_intersection(hist_1, hist_2) - dd.append([mass, rr]) - if rr < .7: - print(mass, rr, len(data_1), len(data_2)) - plt.hist(data_1, 100, alpha=.5, density=True, label='old') - plt.hist(data_2, 100, alpha=.5, density=True, label='new') - plt.xscale('log');plt.yscale('log') - plt.legend();plt.show() - - dd = np.array(dd).T - plt.scatter(*dd) - plt.show() + st_dist_mass = (sampled_IMF, np.cumsum(sampled_IMF)) return st_dist_mass -def integral_IMF_M(m_star, IMF_name): - ''' - Returns mass values: $$ - ''' - return m_star * imfs(IMF_name, m_star) - - -def integral_IMF_N(m_star, IMF_name): - ''' - Returns number of stars: $$ - ''' +def IMF_func(m_star, IMF_name): return imfs(IMF_name, m_star) @@ -194,68 +135,3 @@ def imfs(IMF_name, m_star): imf_val = m_star ** -2.35 return imf_val - - -def massDist(m_low, mass_up, N_dist, masses): - """ - Given the mass distribution from the sampled IMF, obtain for each mass - defined: the total number of stars, and the scale and base factors - that will distribute them properly. - """ - st_dist_mass = {} - - for M_total in masses: - # Normalize number of stars per interval of mass according to total - # mass. This is equivalent to Eq. (8) in PH09. - N_stars = N_dist * M_total - # Distribute the N_stars for this total mass. - m_low_i, m_up_i, N_stars_i = starsInterv(m_low, mass_up, N_stars) - - # Store mass distribution parameters. - N_stars_total = np.sum(N_stars_i) - base = np.repeat(m_low_i, N_stars_i) - scale = np.repeat(m_up_i, N_stars_i) - base - - # DEPRECATED May 2019, all models need to be identical for the same - # input parameters. This mode could be useful for the 'synth_mode' - # mode to generate synthetic clusters (#239). - # - # # Either pass the values so that each synthetic cluster can generate - # # its own mass distribution, or generate it once here. - # if m_sample_flag is True: - # # Add flag to the list to avoid having to pass 'm_sample_flag' - # # across the code. - # st_dist_mass[M_total] = [base, scale, N_stars_total, True] - # else: - st_dist_mass[M_total] =\ - np.random.random(N_stars_total) * scale + base - - return st_dist_mass - - -def starsInterv(m_low, mass_up, N_stars): - """ - Distribute the N_stars in each mass interval, so that each interval - contains at least one star. - """ - m_low_i, m_up_i, N_stars_i = [], [], [] - N_st_add = 0. - for m_up, N_st in zip(*[mass_up, N_stars]): - N_st += N_st_add - # If the number of stars in the interval is less than 1, combine as - # many adjacent intervals as necessary until it reaches at least one - # star and then generate that star(s) with a random mass in the - # m_low, m_up interval. - if N_st < 1.: - # Store this fraction to be added with the next one. - N_st_add = N_st - else: - # The N_st stars will be randomly distributed between - # m_low and m_up limits. - m_low_i.append(m_low) - m_up_i.append(m_up) - N_stars_i.append(int(round(N_st))) - # Reset parameters and move on to the next interval. - N_st_add, m_low = 0., m_up - - return m_low_i, m_up_i, N_stars_i diff --git a/packages/synth_clust/mass_distribution.py b/packages/synth_clust/mass_distribution.py index f5d1e340..087383b4 100644 --- a/packages/synth_clust/mass_distribution.py +++ b/packages/synth_clust/mass_distribution.py @@ -23,7 +23,7 @@ def main(st_dist_mass, M_total): # base, scale, N_stars = st_dist_mass[M_total][:-1] # mass_dist = np.random.random(N_stars) * scale + base # else: - # mass_dist = st_dist_mass[M_total] + mass_dist = st_dist_mass[0][:np.searchsorted(st_dist_mass[1], M_total)] return mass_dist diff --git a/packages/synth_clust/synth_cl_plot.py b/packages/synth_clust/synth_cl_plot.py index 2e178225..facce82e 100644 --- a/packages/synth_clust/synth_cl_plot.py +++ b/packages/synth_clust/synth_cl_plot.py @@ -17,8 +17,8 @@ def main( elif best_fit_algor == 'ptemcee': # Use mean fit values for all parameters. synth_cl_params = isoch_fit_params['mean_sol'] - # Grid values for (z, a, M) - synth_cl_params_grid = closeSol(fundam_params, synth_cl_params, [0, 1, 4]) + # Grid values for (z, a) + synth_cl_params_grid = closeSol(fundam_params, synth_cl_params, [0, 1]) # Find indexes for metallicity and age. If indexes are not found due # to some difference in the significant figures, use the indexes From e2cb9a6880a312f65ea48051e55414d05790a3a3 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Mon, 30 Sep 2019 17:25:46 -0300 Subject: [PATCH 16/27] this issue was closed long ago --- packages/synth_clust/binarity.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/packages/synth_clust/binarity.py b/packages/synth_clust/binarity.py index eb83ca4e..019d9229 100644 --- a/packages/synth_clust/binarity.py +++ b/packages/synth_clust/binarity.py @@ -13,11 +13,6 @@ def mag_combine(m1, m2): """ c = 10 ** -.4 - # TODO - # This catches an overflow warning issued because some Marigo isochrones - # contain huge values in the U filter. Not sure if this happens with - # other systems/filters. See issue #375 - np.warnings.filterwarnings('ignore') mbin = -2.5 * (-.4 * m1 + np.log10(1. + c ** (m2 - m1))) return mbin From 3e84db79bfd99f8c3555a8d88857206504141b39 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Tue, 1 Oct 2019 10:02:00 -0300 Subject: [PATCH 17/27] close #439 and #440 --- packages/best_fit/best_fit_synth_cl.py | 8 +- packages/best_fit/bf_common.py | 6 +- packages/best_fit/bootstrap.py | 29 ++-- packages/best_fit/genetic_algorithm.py | 16 +- packages/best_fit/ptemcee_algor.py | 4 +- packages/best_fit/zaInterp.py | 218 +++++++++++++------------ packages/synth_clust/synth_cluster.py | 8 +- 7 files changed, 148 insertions(+), 141 deletions(-) diff --git a/packages/best_fit/best_fit_synth_cl.py b/packages/best_fit/best_fit_synth_cl.py index efec11ab..7a160fdd 100644 --- a/packages/best_fit/best_fit_synth_cl.py +++ b/packages/best_fit/best_fit_synth_cl.py @@ -39,6 +39,10 @@ def main(clp, pd): # Store the number of defined filters and colors. N_fc = [len(pd['filters']), len(pd['colors'])] + # Index of m_ini (theoretical initial mass), stored in the theoretical + # isochrones. + m_ini = 2 * N_fc[0] + 2 * N_fc[1] + 2 + # HARDCODED: generate random floats to use in the synthetic cluster # completeness removal and error adding. cmpl_rnd = np.random.uniform(0., 1., 1000000) @@ -54,7 +58,7 @@ def main(clp, pd): pd['lkl_method'] == 'dolphin' else pd['lkl_method'])) isoch_fit_params = bootstrap.main( pd, clp, cl_max_mag, max_mag_syn, obs_clust, ext_coefs, - st_dist_mass, N_fc, cmpl_rnd, err_rnd) + st_dist_mass, N_fc, m_ini, cmpl_rnd, err_rnd) # Assign uncertainties. isoch_fit_errors = params_errors(pd, isoch_fit_params) @@ -66,7 +70,7 @@ def main(clp, pd): isoch_fit_params = ptemcee_algor.main( clp['err_lst'], clp['completeness'], clp['em_float'], max_mag_syn, obs_clust, ext_coefs, st_dist_mass, N_fc, - cmpl_rnd, err_rnd, **pd) + m_ini, cmpl_rnd, err_rnd, **pd) # Assign uncertainties. isoch_fit_errors = params_errors(pd, isoch_fit_params) diff --git a/packages/best_fit/bf_common.py b/packages/best_fit/bf_common.py index 4d457d10..2d050ae0 100644 --- a/packages/best_fit/bf_common.py +++ b/packages/best_fit/bf_common.py @@ -124,16 +124,16 @@ def synthClust(fundam_params, varIdxs, model, synthcl_args): Generate a synthetic cluster. """ theor_tracks, e_max, err_lst, completeness, max_mag_syn, st_dist_mass,\ - R_V, ext_coefs, N_fc, cmpl_rnd, err_rnd = synthcl_args + R_V, ext_coefs, N_fc, m_ini, cmpl_rnd, err_rnd = synthcl_args # Interpolate (z, a) data isochrone, model_proper = zaInterp.main( - theor_tracks, fundam_params, varIdxs, model) + theor_tracks, fundam_params, varIdxs, model, m_ini) # Generate synthetic cluster. return synth_cluster.main( e_max, err_lst, completeness, max_mag_syn, st_dist_mass, isochrone, - R_V, ext_coefs, N_fc, cmpl_rnd, err_rnd, model_proper) + R_V, ext_coefs, N_fc, m_ini, cmpl_rnd, err_rnd, model_proper) # DEPRECATED 24-09-2019 diff --git a/packages/best_fit/bootstrap.py b/packages/best_fit/bootstrap.py index df1b17d7..65c3f7ca 100644 --- a/packages/best_fit/bootstrap.py +++ b/packages/best_fit/bootstrap.py @@ -11,7 +11,7 @@ def main( pd, clp, cl_max_mag, max_mag_syn, obs_clust, ext_coefs, st_dist_mass, - N_fc, cmpl_rnd, err_rnd): + N_fc, m_ini, cmpl_rnd, err_rnd): ''' Non-parametric bootstrap process, runs the selected algorithm a number of times each time generating a new observed cluster representation through @@ -40,15 +40,16 @@ def main( flag_print_perc = True argsOF = [ pd, clp, max_mag_syn, obs_clust, ext_coefs, st_dist_mass, N_fc, - cmpl_rnd, err_rnd, available_secs, init_pop, N_pop_init, + m_ini, cmpl_rnd, err_rnd, available_secs, init_pop, N_pop_init, N_gen, flag_print_perc] elif pd['best_fit_algor'] == 'boot+DE': popsize, maxiter = pd['N_pop'], pd['N_gen'] argsOF = [ - pd, clp, max_mag_syn, st_dist_mass, ext_coefs, N_fc, cmpl_rnd, - err_rnd, obs_clust, popsize, maxiter, available_secs, True] + pd, clp, max_mag_syn, st_dist_mass, ext_coefs, N_fc, m_ini, + cmpl_rnd, err_rnd, obs_clust, popsize, maxiter, available_secs, + True] # First run of the numerical optimizer function with the observed data. isoch_fit_params = optimizerFunc(pd['best_fit_algor'], argsOF) @@ -77,14 +78,14 @@ def main( if pd['best_fit_algor'] == 'boot+GA': argsOF = [ pd, clp, max_mag_syn, obs_cl, ext_coefs, st_dist_mass, - N_fc, cmpl_rnd, err_rnd, np.inf, isoch_fit_params[ + N_fc, m_ini, cmpl_rnd, err_rnd, np.inf, isoch_fit_params[ 'OF_final_generation'][:pd['N_pop_btstrp']], pd['N_pop_btstrp'], pd['N_step_btstrp'], False] elif pd['best_fit_algor'] == 'boot+DE': argsOF = [ pd, clp, max_mag_syn, st_dist_mass, ext_coefs, N_fc, - cmpl_rnd, err_rnd, obs_cl, pd['N_pop_btstrp'], + m_ini, cmpl_rnd, err_rnd, obs_cl, pd['N_pop_btstrp'], pd['N_step_btstrp'], np.nan, False] params_boot.append(optimizerFunc(pd['best_fit_algor'], argsOF)) @@ -134,25 +135,25 @@ def optimizerFunc(best_fit_algor, args): """ if best_fit_algor == 'boot+GA': pd, clp, max_mag_syn, obs_clust, ext_coefs, st_dist_mass, N_fc,\ - cmpl_rnd, err_rnd, available_secs, init_pop, N_pop, N_gen,\ + m_ini, cmpl_rnd, err_rnd, available_secs, init_pop, N_pop, N_gen,\ flag_print_perc = args if flag_print_perc: # Print advances. isoch_fit_params = genetic_algorithm.main( available_secs, init_pop, flag_print_perc, N_pop, N_gen, - max_mag_syn, obs_clust, ext_coefs, st_dist_mass, N_fc, + max_mag_syn, obs_clust, ext_coefs, st_dist_mass, N_fc, m_ini, cmpl_rnd, err_rnd, clp['em_float'], clp['err_lst'], clp['completeness'], **pd) else: isoch_fit_params = genetic_algorithm.main( available_secs, init_pop, flag_print_perc, N_pop, N_gen, - max_mag_syn, obs_clust, ext_coefs, st_dist_mass, N_fc, + max_mag_syn, obs_clust, ext_coefs, st_dist_mass, N_fc, m_ini, cmpl_rnd, err_rnd, clp['em_float'], clp['err_lst'], clp['completeness'], **pd) elif best_fit_algor == 'boot+DE': - pd, clp, max_mag_syn, st_dist_mass, ext_coefs, N_fc, cmpl_rnd,\ + pd, clp, max_mag_syn, st_dist_mass, ext_coefs, N_fc, m_ini, cmpl_rnd,\ err_rnd, obs_clust, popsize, maxiter, available_secs,\ flag_print_perc = args @@ -160,14 +161,14 @@ def optimizerFunc(best_fit_algor, args): # Print advances. isoch_fit_params = de_algorithm.main( clp['em_float'], clp['err_lst'], clp['completeness'], - max_mag_syn, st_dist_mass, ext_coefs, N_fc, cmpl_rnd, err_rnd, - obs_clust, popsize, maxiter, available_secs, + max_mag_syn, st_dist_mass, ext_coefs, N_fc, m_ini, cmpl_rnd, + err_rnd, obs_clust, popsize, maxiter, available_secs, flag_print_perc, **pd) else: isoch_fit_params = de_algorithm.main( clp['em_float'], clp['err_lst'], clp['completeness'], - max_mag_syn, st_dist_mass, ext_coefs, N_fc, cmpl_rnd, err_rnd, - obs_clust, popsize, maxiter, available_secs, + max_mag_syn, st_dist_mass, ext_coefs, N_fc, m_ini, cmpl_rnd, + err_rnd, obs_clust, popsize, maxiter, available_secs, flag_print_perc, **pd) return isoch_fit_params diff --git a/packages/best_fit/genetic_algorithm.py b/packages/best_fit/genetic_algorithm.py index 0d5a91b9..4f44fa3c 100644 --- a/packages/best_fit/genetic_algorithm.py +++ b/packages/best_fit/genetic_algorithm.py @@ -26,7 +26,7 @@ def main( available_secs, ran_pop, flag_print_perc, N_popl, N_gener, max_mag_syn, - obs_clust, ext_coefs, st_dist_mass, N_fc, cmpl_rnd, err_rnd, e_max, + obs_clust, ext_coefs, st_dist_mass, N_fc, m_ini, cmpl_rnd, err_rnd, e_max, err_lst, completeness, lkl_method, fundam_params, theor_tracks, R_V, fit_diff, cross_prob, cross_sel, mut_prob, N_el, N_ei, N_es, **kwargs): ''' @@ -61,8 +61,8 @@ def main( # Evaluate initial random solutions in the objective function. generation, lkl = evaluation( lkl_method, e_max, err_lst, completeness, max_mag_syn, fundam_params, - obs_clust, theor_tracks, R_V, ext_coefs, st_dist_mass, N_fc, cmpl_rnd, - err_rnd, varIdxs, ranges, ran_pop) + obs_clust, theor_tracks, R_V, ext_coefs, st_dist_mass, N_fc, m_ini, + cmpl_rnd, err_rnd, varIdxs, ranges, ran_pop) # Store best solution(s) for passing along in the 'Elitism' block. best_sol = generation[:N_el] @@ -120,7 +120,8 @@ def main( generation, lkl = evaluation( lkl_method, e_max, err_lst, completeness, max_mag_syn, fundam_params, obs_clust, theor_tracks, R_V, ext_coefs, - st_dist_mass, N_fc, cmpl_rnd, err_rnd, varIdxs, ranges, p_lst_e) + st_dist_mass, N_fc, m_ini, cmpl_rnd, err_rnd, varIdxs, ranges, + p_lst_e) if flag_print_perc: @@ -307,8 +308,8 @@ def decode(fundam_params, n_bin, p_delta, p_mins, mut_chrom): def evaluation( lkl_method, e_max, err_lst, completeness, max_mag_syn, fundam_params, - obs_clust, theor_tracks, R_V, ext_coefs, st_dist_mass, N_fc, cmpl_rnd, - err_rnd, varIdxs, ranges, p_lst): + obs_clust, theor_tracks, R_V, ext_coefs, st_dist_mass, N_fc, m_ini, + cmpl_rnd, err_rnd, varIdxs, ranges, p_lst): ''' Evaluate each model in the objective function to obtain the fitness of each one. @@ -327,7 +328,8 @@ def evaluation( synth_clust = synthClust( fundam_params, varIdxs, model, ( theor_tracks, e_max, err_lst, completeness, max_mag_syn, - st_dist_mass, R_V, ext_coefs, N_fc, cmpl_rnd, err_rnd)) + st_dist_mass, R_V, ext_coefs, N_fc, m_ini, cmpl_rnd, + err_rnd)) # Call likelihood function for this model. # with timeblock(" Likelihood"): diff --git a/packages/best_fit/ptemcee_algor.py b/packages/best_fit/ptemcee_algor.py index 84b16d60..ed7425dd 100644 --- a/packages/best_fit/ptemcee_algor.py +++ b/packages/best_fit/ptemcee_algor.py @@ -12,7 +12,7 @@ def main( err_lst, completeness, e_max, max_mag_syn, obs_clust, ext_coefs, - st_dist_mass, N_fc, cmpl_rnd, err_rnd, lkl_method, fundam_params, + st_dist_mass, N_fc, m_ini, cmpl_rnd, err_rnd, lkl_method, fundam_params, theor_tracks, R_V, ntemps, nsteps_pt, nwalkers_pt, nburn_pt, pt_adapt, tmax_pt, priors_pt, hmax, **kwargs): """ @@ -22,7 +22,7 @@ def main( # Pack synthetic cluster arguments. synthcl_args = [ theor_tracks, e_max, err_lst, completeness, max_mag_syn, st_dist_mass, - R_V, ext_coefs, N_fc, cmpl_rnd, err_rnd] + R_V, ext_coefs, N_fc, m_ini, cmpl_rnd, err_rnd] if tmax_pt in ('n', 'none', 'None'): Tmax = None diff --git a/packages/best_fit/zaInterp.py b/packages/best_fit/zaInterp.py index 7867fcb4..9e0b3418 100644 --- a/packages/best_fit/zaInterp.py +++ b/packages/best_fit/zaInterp.py @@ -1,11 +1,10 @@ import numpy as np -import warnings -def main(theor_tracks, fundam_params, varIdxs, model): +def main(theor_tracks, fundam_params, varIdxs, model, m_ini): """ - Interpolate a new isochrone from the four closest points in the (z, a) + Average a new isochrone from the four closest points in the (z, a) grid. The mass value is not interpolated. @@ -21,109 +20,72 @@ def main(theor_tracks, fundam_params, varIdxs, model): bp: binary probabilities mb: binary masses m_ini,..., m_bol: six extra parameters. - This means that '-6' is the index of the initial masses. """ - # Select the proper values for the model (ie: that exist in the grid), and - # choose the (z, a) grid indexes to interpolate the isochrone. - model_proper, j = [], 0 - for i, par in enumerate(fundam_params): - # Check if this parameter is one of the 'free' parameters. - if i in varIdxs: - # If it is the parameter metallicity. - if i == 0: - # Select the closest value in the array of allowed values. - mh = min(len(par) - 1, np.searchsorted(par, model[i - j])) - ml = mh - 1 - model_proper.append(par[mh]) - # Define model's z value - z_model = model[i - j] - # If it is the parameter log(age). - elif i == 1: - # Select the closest value in the array of allowed values. - ah = min(len(par) - 1, np.searchsorted(par, model[i - j])) - al = ah - 1 - model_proper.append(par[ah]) - a_model = model[i - j] - else: - model_proper.append(model[i - j]) - else: - if i == 0: - ml = mh = 0 - z_model = fundam_params[0][0] - elif i == 1: - al = ah = 0 - a_model = fundam_params[1][0] - model_proper.append(par[0]) - j += 1 + # Define the 'proper' model with values for (z, a) taken from its grid, + # and filled values for those parameters that are fixed. + model_proper, z_model, a_model, ml, mh, al, ah = properModel( + fundam_params, model, varIdxs) - # TODO: IN PLACE FOR #439 - # # If (z, a) are both fixed, just return the single processed isochrone - # if ml == al == mh == ah == 0: - # return theor_tracks[ml][al], model_proper + # If (z, a) are both fixed, just return the single processed isochrone + if ml == al == mh == ah == 0: + return theor_tracks[ml][al], model_proper - # Values of the four points in the (z, age) grid that contain the model - # value (z_model, a_model) + # The four points in the (z, age) grid that define the box that contains + # the model value (z_model, a_model) z1, z2 = fundam_params[0][ml], fundam_params[0][mh] a1, a2 = fundam_params[1][al], fundam_params[1][ah] pts = np.array([(z1, a1), (z1, a2), (z2, a1), (z2, a2)]) + # Order: (z1, a1), (z1, a2), (z2, a1), (z2, a2) + isochs = np.array([ + theor_tracks[ml][al], theor_tracks[ml][ah], theor_tracks[mh][al], + theor_tracks[mh][ah]]) + # Define weights for the average, based on the inverse of the distance # to the grid (z, age) points. - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - - # Inverse distances between the (z, a) points in the 'model', and the - # four points in the (z, a) grid that contain the model point. - # Fast euclidean distance: https://stackoverflow.com/a/47775357/1391441 - a_min_b = np.array([(z_model, a_model)]) - pts - - # TODO: IN PLACE FOR #439 - # # If the model has a 0. distance in (z,a) to the closest isochrone, - # # then just return that isochrone. - # try: - # # Which one is faster? - # idx = np.nonzero(a_min_b == 0.)[0][0] - # idx = np.where(a_min_b == 0.)[0][0] - # return theor_tracks[pts[idx][0]][pts[idx][1]], model_proper - # except IndexError: - # pass - - # Inverse distance. - inv_d = 1. / np.sqrt(np.einsum('ij,ij->i', a_min_b, a_min_b)) - - # Order: (z1, a1), (z1, a2), (z2, a1), (z2, a2) - isochs = np.array([ - theor_tracks[ml][al], theor_tracks[ml][ah], theor_tracks[mh][al], - theor_tracks[mh][ah]]) - - # Sort according to smaller distance to the (z_model, a_model) point. - isoch_idx = np.argsort(inv_d)[::-1] - - # Masked weighted average. Source: - # https://stackoverflow.com/a/35758345/1391441 - x1, x2, x3, x4 = isochs[isoch_idx] - for x in (x2, x3, x4): - # Maximum mass difference allowed - msk = abs(x1[-6] - x[-6]) > .01 - # If the distance in this array is larger than the maximum allowed, - # mask with the values taken from 'x1'. - # x[:, msk] = x1[:, msk] - np.copyto(x, x1, where=msk) - - # # Weighted average with mass "alignment". - # isochrone = np.average( - # np.array([x1, x2, x3, x4]), weights=inv_d[isoch_idx], axis=0) - # Scale weights so they add up to 1, then add based on them - weights = inv_d[isoch_idx] / np.sum(inv_d[isoch_idx]) - # If the model has a 0. distance in (z,a) to the closest isochrone, - # then 'inv_d' will be 'inf', and 'weights' will be 'nan' (all other - # values will be zero) Replace that weight with 1. - weights[np.isnan(weights)] = 1. - isochrone = x1 * weights[0] + x2 * weights[1] + x3 * weights[2] +\ - x4 * weights[3] + # Inverse distances between the (z, a) points in the 'model', and the + # four points in the (z, a) grid that contain the model point. + # Fast euclidean distance: https://stackoverflow.com/a/47775357/1391441 + a_min_b = np.array([(z_model, a_model)]) - pts + dist = np.sqrt(np.einsum('ij,ij->i', a_min_b, a_min_b)) + + # If the model has a 0. distance in (z, a) to the closest isochrone, + # then just return that isochrone. + try: + idx = np.where(dist == 0.)[0][0] + return isochs[idx], model_proper + except IndexError: + pass + + # Inverse of the distance. + inv_d = 1. / dist + + # Sort according to smaller distances to the (z_model, a_model) point. + isoch_idx = np.argsort(inv_d)[::-1] + + # This block "aligns" the masses such that if the 'mass distance' between a + # star in the closest grid isochrone to the model (x1) and a star from any + # of the remaining three isochrones (x2, x3, x4) is larger than a given + # max value (.01, then we replace the star in the second isochrone with its + # 'aligned' star taken from the closest isochrone (x1). + x1, x2, x3, x4 = isochs[isoch_idx] + for x in (x2, x3, x4): + # Mask according to the maximum mass difference allowed + msk = abs(x1[m_ini] - x[m_ini]) > .01 + # Replace stars in 'x' with large mass differences, with stars in 'x1' + np.copyto(x, x1, where=msk) + + # Weighted average with mass "alignment". This way is faster than using + # 'np.average()'. + # Scale weights so they add up to 1. + weights = inv_d[isoch_idx] / np.sum(inv_d) + isochrone = x1 * weights[0] + x2 * weights[1] + x3 * weights[2] +\ + x4 * weights[3] + + ######################## # This way is *marginally* faster # wgts = D * inv_d[isoch_idx] # x1, x2, x3, x4 = isochs[isoch_idx] @@ -134,25 +96,15 @@ def main(theor_tracks, fundam_params, varIdxs, model): # isochrone = np.sum(np.array([ # x1 * wgts[0], x2 * wgts[1], x3 * wgts[2], x4 * wgts[3]]), 0) - # Weighted version, no mass alignment. - # isochrone = np.average(np.array([ - # theor_tracks[ml][al], theor_tracks[ml][ah], theor_tracks[mh][al], - # theor_tracks[mh][ah]]), weights=D * inv_d, axis=0) - - # Simpler mean version, no weights. - # isochrone = np.mean([ - # theor_tracks[ml][al], theor_tracks[ml][ah], theor_tracks[mh][al], - # theor_tracks[mh][ah]], axis=0) - # nn = np.random.randint(0, 100) - # if nn == 50: - # print(model) - # print(model_proper) + # if nn == 50: + # print(model, model_proper) + # print(model_proper[0] - model[0]) # import matplotlib.pyplot as plt # plt.subplot(131) # plt.scatter(*pts[isoch_idx][0], c='r') # plt.scatter(*pts[isoch_idx][1:].T, c='g') - # plt.scatter(model[0], model[1], marker='x') + # plt.scatter(model[0], model[1], marker='x', c='g') # # First color # plt.subplot(132) # plt.plot(theor_tracks[ml][al][1], theor_tracks[ml][al][0], c='b') @@ -172,3 +124,53 @@ def main(theor_tracks, fundam_params, varIdxs, model): # plt.show() return isochrone, model_proper + + +def properModel(fundam_params, model, varIdxs): + """ + + Returns + ------- + model_proper : list + Stores the closest (z, a) values in the grid for the parameters in + 'model', and add the fixed parameters that are missing from 'model'. + z_model, a_model : floats + The (z, a) values for this model's isochrone. + ml, mh, al, ah : ints + Indexes of the (z, a) values in the grid that define the box that enclose + the (z_model, a_model) values. + + """ + + model_proper, j = [], 0 + for i, par in enumerate(fundam_params): + # Check if this parameter is one of the 'free' parameters. + if i in varIdxs: + # If it is the parameter metallicity. + if i == 0: + # Select the closest value in the array of allowed values. + mh = min(len(par) - 1, np.searchsorted(par, model[i - j])) + ml = mh - 1 + model_proper.append(par[mh]) + # Define model's z value + z_model = model[i - j] + # If it is the parameter log(age). + elif i == 1: + # Select the closest value in the array of allowed values. + ah = min(len(par) - 1, np.searchsorted(par, model[i - j])) + al = ah - 1 + model_proper.append(par[ah]) + a_model = model[i - j] + else: + model_proper.append(model[i - j]) + else: + if i == 0: + ml = mh = 0 + z_model = fundam_params[0][0] + elif i == 1: + al = ah = 0 + a_model = fundam_params[1][0] + model_proper.append(par[0]) + j += 1 + + return model_proper, z_model, a_model, ml, mh, al, ah diff --git a/packages/synth_clust/synth_cluster.py b/packages/synth_clust/synth_cluster.py index ebda7d8f..acc8f463 100644 --- a/packages/synth_clust/synth_cluster.py +++ b/packages/synth_clust/synth_cluster.py @@ -8,8 +8,9 @@ from . import add_errors -def main(err_max, err_lst, completeness, max_mag_syn, st_dist_mass, - isochrone, R_V, ext_coefs, N_fc, cmpl_rnd, err_rnd, synth_cl_params): +def main( + err_max, err_lst, completeness, max_mag_syn, st_dist_mass, isochrone, + R_V, ext_coefs, N_fc, m_ini, cmpl_rnd, err_rnd, synth_cl_params): """ Takes an isochrone and returns a synthetic cluster created according to a certain mass distribution. @@ -72,9 +73,6 @@ def main(err_max, err_lst, completeness, max_mag_syn, st_dist_mass, # t3 = time.clock() - s # s = time.clock() - # Index of m_ini (theoretical initial mass), stored in the theoretical - # isochrones. - m_ini = 2 * N_fc[0] + 2 * N_fc[1] + 2 # Interpolate masses in mass_dist into the isochrone rejecting those # masses that fall outside of the isochrone's mass range. isoch_mass = mass_interp.main(isoch_cut, mass_dist, m_ini) From f2abc4a06fc1c3c6a390df09cffb7106a93406d9 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Wed, 2 Oct 2019 15:33:42 -0300 Subject: [PATCH 18/27] minor minor --- packages/best_fit/genetic_algorithm.py | 1 - 1 file changed, 1 deletion(-) diff --git a/packages/best_fit/genetic_algorithm.py b/packages/best_fit/genetic_algorithm.py index 4f44fa3c..4c47172a 100644 --- a/packages/best_fit/genetic_algorithm.py +++ b/packages/best_fit/genetic_algorithm.py @@ -4,7 +4,6 @@ import random import numpy as np from .bf_common import varPars, rangeCheck, synthClust, random_population -# from ..synth_clust import synth_cluster from . import likelihood ############################################################# From 7e34c832ce43c0502e09c50129276138a4e4ebf8 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Wed, 2 Oct 2019 16:59:17 -0300 Subject: [PATCH 19/27] intermediate commit to save the changes mainly in isoch_params() --- packages/best_fit/bf_common.py | 4 +- .../best_fit/{zaInterp.py => zaWAverage.py} | 108 +++++------ packages/inp/isoch_params.py | 168 ++++++++++++------ 3 files changed, 170 insertions(+), 110 deletions(-) rename packages/best_fit/{zaInterp.py => zaWAverage.py} (64%) diff --git a/packages/best_fit/bf_common.py b/packages/best_fit/bf_common.py index 2d050ae0..425e101a 100644 --- a/packages/best_fit/bf_common.py +++ b/packages/best_fit/bf_common.py @@ -5,7 +5,7 @@ from scipy import stats import warnings from ..synth_clust import synth_cluster -from . import likelihood, zaInterp +from . import likelihood, zaWAverage from .. import update_progress @@ -127,7 +127,7 @@ def synthClust(fundam_params, varIdxs, model, synthcl_args): R_V, ext_coefs, N_fc, m_ini, cmpl_rnd, err_rnd = synthcl_args # Interpolate (z, a) data - isochrone, model_proper = zaInterp.main( + isochrone, model_proper = zaWAverage.main( theor_tracks, fundam_params, varIdxs, model, m_ini) # Generate synthetic cluster. diff --git a/packages/best_fit/zaInterp.py b/packages/best_fit/zaWAverage.py similarity index 64% rename from packages/best_fit/zaInterp.py rename to packages/best_fit/zaWAverage.py index 9e0b3418..627edfcb 100644 --- a/packages/best_fit/zaInterp.py +++ b/packages/best_fit/zaWAverage.py @@ -4,10 +4,8 @@ def main(theor_tracks, fundam_params, varIdxs, model, m_ini): """ - Average a new isochrone from the four closest points in the (z, a) - grid. - - The mass value is not interpolated. + Average a new (weighted) isochrone from the four closest points in the + (z, a) grid. theor_tracks = [m1, m2, .., mN] mX = [age1, age2, ..., ageM] @@ -21,6 +19,9 @@ def main(theor_tracks, fundam_params, varIdxs, model, m_ini): mb: binary masses m_ini,..., m_bol: six extra parameters. + WARNING: currently only reading the 'm_ini' extra parameter in + 'isochs_format()' + """ # Define the 'proper' model with values for (z, a) taken from its grid, @@ -39,9 +40,9 @@ def main(theor_tracks, fundam_params, varIdxs, model, m_ini): pts = np.array([(z1, a1), (z1, a2), (z2, a1), (z2, a2)]) # Order: (z1, a1), (z1, a2), (z2, a1), (z2, a2) - isochs = np.array([ + isochs = [ theor_tracks[ml][al], theor_tracks[ml][ah], theor_tracks[mh][al], - theor_tracks[mh][ah]]) + theor_tracks[mh][ah]] # Define weights for the average, based on the inverse of the distance # to the grid (z, age) points. @@ -64,64 +65,69 @@ def main(theor_tracks, fundam_params, varIdxs, model, m_ini): inv_d = 1. / dist # Sort according to smaller distances to the (z_model, a_model) point. - isoch_idx = np.argsort(inv_d)[::-1] + # isoch_idx = np.argsort(inv_d)[::-1] + + # Indexes that identify which of the four arrays in 'x-all' contain the + # the largest minimum and the smallest maximum in their '-1' sub-array + vmin = np.argmax([isochs[_][m_ini][0] for _ in range(4)]) + vmax = np.argmin([isochs[_][m_ini][-1] for _ in range(4)]) + # Identify lower and upper range limits for the '-1' sub-array for all + # arrays in 'isochs' + xmin, xmax = isochs[vmin][m_ini][0], isochs[vmax][m_ini][-1] + + # Align from bottom and top + for i in [0, 1, 2, 3]: + imin = np.searchsorted(isochs[i][m_ini], xmin) + imax = np.searchsorted(isochs[i][m_ini], xmax) + isochs[i] = isochs[i][:, imin:imax] # This block "aligns" the masses such that if the 'mass distance' between a # star in the closest grid isochrone to the model (x1) and a star from any # of the remaining three isochrones (x2, x3, x4) is larger than a given # max value (.01, then we replace the star in the second isochrone with its # 'aligned' star taken from the closest isochrone (x1). - x1, x2, x3, x4 = isochs[isoch_idx] - for x in (x2, x3, x4): - # Mask according to the maximum mass difference allowed - msk = abs(x1[m_ini] - x[m_ini]) > .01 - # Replace stars in 'x' with large mass differences, with stars in 'x1' - np.copyto(x, x1, where=msk) + # x1, x2, x3, x4 = isochs[isoch_idx] + # for x in (x2, x3, x4): + # # Mask according to the maximum mass difference allowed + # msk = abs(x1[m_ini] - x[m_ini]) > .01 + # # Replace stars in 'x' with large mass differences, with stars in 'x1' + # np.copyto(x, x1, where=msk) # Weighted average with mass "alignment". This way is faster than using # 'np.average()'. # Scale weights so they add up to 1. - weights = inv_d[isoch_idx] / np.sum(inv_d) - isochrone = x1 * weights[0] + x2 * weights[1] + x3 * weights[2] +\ - x4 * weights[3] - - ######################## - # This way is *marginally* faster - # wgts = D * inv_d[isoch_idx] - # x1, x2, x3, x4 = isochs[isoch_idx] - # for x in (x2, x3, x4): - # # Maximum mass difference allowed - # msk = abs(x1[-6] - x[-6]) > .01 - # x[:, msk] = x1[:, msk] - # isochrone = np.sum(np.array([ - # x1 * wgts[0], x2 * wgts[1], x3 * wgts[2], x4 * wgts[3]]), 0) + weights = inv_d / np.sum(inv_d) + isochrone = isochs[0] * weights[0] + isochs[1] * weights[1] +\ + isochs[2] * weights[2] + isochs[3] * weights[3] # nn = np.random.randint(0, 100) # if nn == 50: - # print(model, model_proper) - # print(model_proper[0] - model[0]) - # import matplotlib.pyplot as plt - # plt.subplot(131) - # plt.scatter(*pts[isoch_idx][0], c='r') - # plt.scatter(*pts[isoch_idx][1:].T, c='g') - # plt.scatter(model[0], model[1], marker='x', c='g') - # # First color - # plt.subplot(132) - # plt.plot(theor_tracks[ml][al][1], theor_tracks[ml][al][0], c='b') - # plt.plot(theor_tracks[ml][ah][1], theor_tracks[ml][ah][0], c='r') - # plt.plot(theor_tracks[mh][al][1], theor_tracks[mh][al][0], c='cyan') - # plt.plot(theor_tracks[mh][ah][1], theor_tracks[mh][ah][0], c='orange') - # plt.plot(isochrone[1], isochrone[0], c='g', ls='--') - # plt.gca().invert_yaxis() - # # Second color - # plt.subplot(133) - # plt.plot(theor_tracks[ml][al][2], theor_tracks[ml][al][0], c='b') - # plt.plot(theor_tracks[ml][ah][2], theor_tracks[ml][ah][0], c='r') - # plt.plot(theor_tracks[mh][al][2], theor_tracks[mh][al][0], c='cyan') - # plt.plot(theor_tracks[mh][ah][2], theor_tracks[mh][ah][0], c='orange') - # plt.plot(isochrone[2], isochrone[0], c='g', ls='--') - # plt.gca().invert_yaxis() - # plt.show() + print(model) + print(model_proper) + import matplotlib.pyplot as plt + plt.subplot(131) + plt.scatter(*pts.T, c='r') + # plt.scatter(*pts[isoch_idx][1:].T, c='g') + plt.scatter(model[0], model[1], marker='x', c='g') + # First color + plt.subplot(132) + plt.plot(theor_tracks[ml][al][1], theor_tracks[ml][al][0], c='b') + plt.plot(theor_tracks[ml][ah][1], theor_tracks[ml][ah][0], c='r') + plt.plot(theor_tracks[mh][al][1], theor_tracks[mh][al][0], c='cyan') + plt.plot(theor_tracks[mh][ah][1], theor_tracks[mh][ah][0], c='orange') + plt.plot(isochrone[1], isochrone[0], c='g', ls='--') + plt.gca().invert_yaxis() + # Second color + plt.subplot(133) + plt.plot(theor_tracks[ml][al][-1], theor_tracks[ml][al][0], c='b') + plt.plot(theor_tracks[ml][ah][-1], theor_tracks[ml][ah][0], c='r') + plt.plot(theor_tracks[mh][al][-1], theor_tracks[mh][al][0], c='cyan') + plt.plot(theor_tracks[mh][ah][-1], theor_tracks[mh][ah][0], c='orange') + plt.plot(isochrone[-1], isochrone[0], c='g', ls='--') + plt.gca().invert_yaxis() + # plt.show() + import pdb; pdb.set_trace() # breakpoint d9c1b180 // + return isochrone, model_proper diff --git a/packages/inp/isoch_params.py b/packages/inp/isoch_params.py index a3ba6d52..cf82dd10 100644 --- a/packages/inp/isoch_params.py +++ b/packages/inp/isoch_params.py @@ -1,6 +1,7 @@ import sys import numpy as np +from scipy.interpolate import interp1d from . import read_isochs from ..synth_clust import binarity from .. import update_progress @@ -39,27 +40,30 @@ def main(met_f_filter, age_values, cmd_evol_tracks, evol_track, bin_mr, # masses to be more accurately interpolated into the theoretical # isochrones. - # Find the maximum number of points in all the read ages for all the - # metallicities. - N_pts_max = 0 - for z in mags_theor: - for a in z: - N_pts_max = max(N_pts_max, len(a[0])) + # # Find the maximum number of points in all the read ages for all the + # # metallicities. + # N_pts_max, m_range = 0, 0. + # for z in extra_pars: + # for j, m_ini in enumerate(z): + # m_min, m_max = m_ini[0][0], m_ini[0][-1] + # N_pts_max = max(N_pts_max, len(m_ini[0])) + # m_range = max(m_range, m_max - m_min) + # # print("{:.3f} {} {:.2f} {:.3f} {:.3f}".format( + # # age_values[j], len(m_ini[0]), m_max-m_min, m_min, m_max)) + # print(N_pts_max, m_range) - if N_mass_interp > N_pts_max: - print("Interpolating extra points ({}) into the isochrones".format( - N_mass_interp)) - interp_data = [] - for i, data in enumerate([ - mags_theor, cols_theor, mags_cols_theor, extra_pars]): - interp_data.append(interp_isoch_data(data, N_mass_interp)) - update_progress.updt(4, i + 1) - a, b, c, d = interp_data - else: - sys.exit( - ("ERROR: N_interp={} must be larger than the maximum\n" + - "number of points in any isochrone read: {}").format( - N_mass_interp, N_pts_max)) + print("Interpolating masses into the isochrones") + interp_data = interp_isoch_data( + mags_theor, cols_theor, mags_cols_theor, extra_pars) + + # interp_data = [] + # for i, data in enumerate([ + # mags_theor, cols_theor, mags_cols_theor, extra_pars]): + # interp_data.append(interp_isoch_data(data, N_mass_interp)) + # update_progress.updt(4, i + 1) + + # a, b, c, d = mags_theor, cols_theor, mags_cols_theor, extra_pars + a, b, c, d = interp_data # # Size of arrays in memory # sz = 0. @@ -71,7 +75,7 @@ def main(met_f_filter, age_values, cmd_evol_tracks, evol_track, bin_mr, # discarded after the colors (and magnitudes) with binarity assignment # are obtained. mags_binar, cols_binar, probs_binar, mass_binar = binarity.binarGen( - fundam_params[5], N_mass_interp, a, b, c, d, bin_mr) + fundam_params[5], a, b, c, d, bin_mr) # Create list structured as: # theor_tracks = [m1, m2, .., mN] @@ -85,29 +89,46 @@ def main(met_f_filter, age_values, cmd_evol_tracks, evol_track, bin_mr, # bp: binary probabilities # mb: binary masses # Mini,...: extra parameters. + # theor_tracks = [[[] for _ in a[0]] for _ in a] + theor_tracks = [] + for i, mx in enumerate(a): + m = [] + for j, ax in enumerate(mx): + m.append([ax, b[i][j], mags_binar[i][j], cols_binar[i][j], + probs_binar[i][j], mass_binar[i][j], d[i][j]]) + theor_tracks.append(m) + + dd = [] + for m in theor_tracks: + m1 = [] + for a in m: + m1.append(np.array(a)[:,0,:]) + dd.append(m1) + + theor_tracks = dd # Combine all data into a single array of shape: # (N_z, N_age, N_data, N_IMF_interp), where 'N_data' is the number of # sub-arrays in each array. - comb_data = np.concatenate( - (a, b, mags_binar, cols_binar, probs_binar, mass_binar, d), axis=2) - - # Sort all isochrones according to the main magnitude (min to max). - # This is necessary so that the cut_max_mag() function does not need - # to do this every time a new synthetic cluster is generated. - theor_tracks = [[[] for _ in a[0]] for _ in a] - for i, mx in enumerate(comb_data): - for j, ax in enumerate(mx): - # TODO ordering the data according to magnitude is necessary - # if the cut_max_mag() function expects the data to be sorted - # this way. This however, prevents us from being able to - # interpolate new (z, a) values when the MCMC samples them - # because the mass order is not preserved anymore. So, in order - # to be able to generate that interpolation of values (which - # greatly improves the ptemcee performance), we're back to the - # 'old' cut_max_mag() function and no magnitude ordering. - # theor_tracks[i][j] = ax[:, ax[0].argsort(kind='mergesort')] - theor_tracks[i][j] = ax + # comb_data = np.concatenate( + # (a, b, mags_binar, cols_binar, probs_binar, mass_binar, d), axis=2) + + # # Sort all isochrones according to the main magnitude (min to max). + # # This is necessary so that the cut_max_mag() function does not need + # # to do this every time a new synthetic cluster is generated. + # theor_tracks = [[[] for _ in a[0]] for _ in a] + # for i, mx in enumerate(comb_data): + # for j, ax in enumerate(mx): + # # TODO ordering the data according to magnitude is necessary + # # if the cut_max_mag() function expects the data to be sorted + # # this way. This however, prevents us from being able to + # # interpolate new (z, a) values when the MCMC samples them + # # because the mass order is not preserved anymore. So, in order + # # to be able to generate that interpolation of values (which + # # greatly improves the ptemcee performance), we're back to the + # # 'old' cut_max_mag() function and no magnitude ordering. + # # theor_tracks[i][j] = ax[:, ax[0].argsort(kind='mergesort')] + # theor_tracks[i][j] = ax # The above sorting destroys the original order of the isochrones. This # results in messy plots for the "best fit isochrone" at the end. @@ -203,23 +224,56 @@ def arrange_filters(isoch_list, all_syst_filters, filters, colors): return mags_theor, cols_theor, mags_cols_theor -def interp_isoch_data(data, N): - ''' +def interp_isoch_data(mags_theor, cols_theor, mags_cols_theor, extra_pars): + """ Interpolate extra values for all the parameters in the theoretic isochrones. - ''' - interp_data = [] - # For each metallicity value list. - for met in data: - m = [] - # For each age value list. - for age in met: - a = [] - # For each filter/color/extra parameter in list. - for fce in age: - t, xp = np.linspace(0., 1., N), np.linspace(0, 1, len(fce)) - a.append(np.interp(t, xp, fce)) - m.append(a) - interp_data.append(m) + """ + # mags_theor, cols_theor, mags_cols, extra_pars --> list + # list = [m0, m1, .., mN] <-- N metallicities + # mx = [a0, a1, .., aM] <-- M log(ages) + # ax = [[s1, s2, .., sQ]] <-- Q stars + # list[i][j][0] --> i metallicity, j age + + mass_step = 0.01 + mass_vals1 = np.arange(.01, 5., mass_step) + mass_step = 0.001 + mass_vals2 = np.arange(5., 18., mass_step) + mass_step = 0.0005 + mass_vals3 = np.arange(18., 150., mass_step) + mass_vals = np.array( + mass_vals1.tolist() + mass_vals2.tolist() + mass_vals3.tolist()) + + all_interp_data = [] + for data in (mags_theor, cols_theor, mags_cols_theor, extra_pars): + interp_data = [] + # For each metallicity value. + for i, met in enumerate(data): + m = [] + # For each age value. + for j, age in enumerate(met): + a = [] + # For each filter/color/extra parameter in list. + for fce in age: + # Initial masses for this isochrone. + mass_ini = np.array(extra_pars[i][j][0]) + # mass values must be in the range of the initial masses + mmin = np.searchsorted(mass_vals, mass_ini.min()) + mmax = np.searchsorted(mass_vals, mass_ini.max()) + mass_vals_cut = mass_vals[mmin:mmax] + + y = np.array(data[i][j][0]) + mfunc = interp1d(mass_ini, y) + ynew = mfunc(mass_vals_cut) + a.append(ynew) + + m.append(a) + interp_data.append(m) + all_interp_data.append(interp_data) + + # np.shape(all_interp_data): (4, Nz, Na) + # 4 : ((mags, cols, mags_cols, extra_pars)) + # Nz : metallicity files (values) + # Na : log(age) values per metallicity file - return interp_data + return np.array(all_interp_data) From 9129c9773bd83820c70a305ac1d4e8ecb7c5fa41 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Wed, 2 Oct 2019 17:10:51 -0300 Subject: [PATCH 20/27] just re-arranging, and change np.zeros with np.empty --- packages/synth_clust/binarity.py | 85 ++++++++++++++++---------------- 1 file changed, 43 insertions(+), 42 deletions(-) diff --git a/packages/synth_clust/binarity.py b/packages/synth_clust/binarity.py index 019d9229..5493e3e6 100644 --- a/packages/synth_clust/binarity.py +++ b/packages/synth_clust/binarity.py @@ -4,18 +4,38 @@ from .. import update_progress -def mag_combine(m1, m2): - """ - Combine two magnitudes. This is a faster re-ordering of the standard - formula: +def main(isoch_mass, bin_frac, m_ini, N_fc): + ''' + Update the randomly selected fraction of binary stars. + ''' - -2.5 * np.log10(10 ** (-0.4 * m1) + 10 ** (-0.4 * m2)) + # If the binary fraction is zero, skip the whole process. + if bin_frac > 0.: - """ - c = 10 ** -.4 - mbin = -2.5 * (-.4 * m1 + np.log10(1. + c ** (m2 - m1))) + # Select a fraction of stars to be binaries, according to the random + # probabilities assigned before. + bin_indxs = isoch_mass[m_ini - 2] <= bin_frac - return mbin + # Index of the first binary magnitude, stored in the theoretical + # isochrones list. + mag_ini = N_fc[0] + N_fc[1] + + # Update array with new values of magnitudes, colors, and masses. + # New magnitudes. + for i in range(N_fc[0]): + isoch_mass[i][bin_indxs] = isoch_mass[mag_ini + i][bin_indxs] + # New colors. + for i in range(N_fc[1]): + isoch_mass[N_fc[0] + i][bin_indxs] =\ + isoch_mass[mag_ini + N_fc[0] + i][bin_indxs] + + # Update masses. + isoch_mass[m_ini][bin_indxs] = isoch_mass[m_ini - 1][bin_indxs] + + # Update binary systems to a '2.'. + isoch_mass[m_ini - 2][bin_indxs] = 2. + + return isoch_mass def binarGen( @@ -57,7 +77,8 @@ def binarGen( # For each age defined. for ax, isoch in enumerate(_): - # Extract initial masses for this isochrone. + # Extract initial masses for this isochrone. Assumes that + # the initial masses are in the '0' index. mass_ini = extra_pars[mx][ax][0] # Calculate random secondary masses of these binary stars @@ -121,43 +142,23 @@ def binarGen( update_progress.updt(N_mags_theor, mx + 1) else: mags_binar, cols_binar, mass_binar =\ - np.zeros(np.shape(mags_theor)),\ - np.zeros(np.shape(cols_theor)),\ + np.empty(np.shape(mags_theor)),\ + np.empty(np.shape(cols_theor)),\ np.array(extra_pars)[:, :, :1, :] - probs_binar = np.zeros(np.shape(mass_binar)) + probs_binar = np.empty(np.shape(mass_binar)) return mags_binar, cols_binar, probs_binar, mass_binar -def main(isoch_mass, bin_frac, m_ini, N_fc): - ''' - Update the randomly selected fraction of binary stars. - ''' - - # If the binary fraction is zero, skip the whole process. - if bin_frac > 0.: - - # Select a fraction of stars to be binaries, according to the random - # probabilities assigned before. - bin_indxs = isoch_mass[m_ini - 2] <= bin_frac - - # Index of the first binary magnitude, stored in the theoretical - # isochrones list. - mag_ini = N_fc[0] + N_fc[1] - - # Update array with new values of magnitudes, colors, and masses. - # New magnitudes. - for i in range(N_fc[0]): - isoch_mass[i][bin_indxs] = isoch_mass[mag_ini + i][bin_indxs] - # New colors. - for i in range(N_fc[1]): - isoch_mass[N_fc[0] + i][bin_indxs] =\ - isoch_mass[mag_ini + N_fc[0] + i][bin_indxs] +def mag_combine(m1, m2): + """ + Combine two magnitudes. This is a faster re-ordering of the standard + formula: - # Update masses. - isoch_mass[m_ini][bin_indxs] = isoch_mass[m_ini - 1][bin_indxs] + -2.5 * np.log10(10 ** (-0.4 * m1) + 10 ** (-0.4 * m2)) - # Update binary systems to a '2.'. - isoch_mass[m_ini - 2][bin_indxs] = 2. + """ + c = 10 ** -.4 + mbin = -2.5 * (-.4 * m1 + np.log10(1. + c ** (m2 - m1))) - return isoch_mass + return mbin From fb4dfab9707a62e868a67ca3b39fa6e432bc2cc8 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Wed, 2 Oct 2019 17:56:07 -0300 Subject: [PATCH 21/27] revert to previous working version --- packages/inp/isoch_params.py | 168 ++++++++++++----------------------- 1 file changed, 57 insertions(+), 111 deletions(-) diff --git a/packages/inp/isoch_params.py b/packages/inp/isoch_params.py index cf82dd10..a3ba6d52 100644 --- a/packages/inp/isoch_params.py +++ b/packages/inp/isoch_params.py @@ -1,7 +1,6 @@ import sys import numpy as np -from scipy.interpolate import interp1d from . import read_isochs from ..synth_clust import binarity from .. import update_progress @@ -40,30 +39,27 @@ def main(met_f_filter, age_values, cmd_evol_tracks, evol_track, bin_mr, # masses to be more accurately interpolated into the theoretical # isochrones. - # # Find the maximum number of points in all the read ages for all the - # # metallicities. - # N_pts_max, m_range = 0, 0. - # for z in extra_pars: - # for j, m_ini in enumerate(z): - # m_min, m_max = m_ini[0][0], m_ini[0][-1] - # N_pts_max = max(N_pts_max, len(m_ini[0])) - # m_range = max(m_range, m_max - m_min) - # # print("{:.3f} {} {:.2f} {:.3f} {:.3f}".format( - # # age_values[j], len(m_ini[0]), m_max-m_min, m_min, m_max)) - # print(N_pts_max, m_range) + # Find the maximum number of points in all the read ages for all the + # metallicities. + N_pts_max = 0 + for z in mags_theor: + for a in z: + N_pts_max = max(N_pts_max, len(a[0])) - print("Interpolating masses into the isochrones") - interp_data = interp_isoch_data( - mags_theor, cols_theor, mags_cols_theor, extra_pars) - - # interp_data = [] - # for i, data in enumerate([ - # mags_theor, cols_theor, mags_cols_theor, extra_pars]): - # interp_data.append(interp_isoch_data(data, N_mass_interp)) - # update_progress.updt(4, i + 1) - - # a, b, c, d = mags_theor, cols_theor, mags_cols_theor, extra_pars - a, b, c, d = interp_data + if N_mass_interp > N_pts_max: + print("Interpolating extra points ({}) into the isochrones".format( + N_mass_interp)) + interp_data = [] + for i, data in enumerate([ + mags_theor, cols_theor, mags_cols_theor, extra_pars]): + interp_data.append(interp_isoch_data(data, N_mass_interp)) + update_progress.updt(4, i + 1) + a, b, c, d = interp_data + else: + sys.exit( + ("ERROR: N_interp={} must be larger than the maximum\n" + + "number of points in any isochrone read: {}").format( + N_mass_interp, N_pts_max)) # # Size of arrays in memory # sz = 0. @@ -75,7 +71,7 @@ def main(met_f_filter, age_values, cmd_evol_tracks, evol_track, bin_mr, # discarded after the colors (and magnitudes) with binarity assignment # are obtained. mags_binar, cols_binar, probs_binar, mass_binar = binarity.binarGen( - fundam_params[5], a, b, c, d, bin_mr) + fundam_params[5], N_mass_interp, a, b, c, d, bin_mr) # Create list structured as: # theor_tracks = [m1, m2, .., mN] @@ -89,46 +85,29 @@ def main(met_f_filter, age_values, cmd_evol_tracks, evol_track, bin_mr, # bp: binary probabilities # mb: binary masses # Mini,...: extra parameters. - # theor_tracks = [[[] for _ in a[0]] for _ in a] - theor_tracks = [] - for i, mx in enumerate(a): - m = [] - for j, ax in enumerate(mx): - m.append([ax, b[i][j], mags_binar[i][j], cols_binar[i][j], - probs_binar[i][j], mass_binar[i][j], d[i][j]]) - theor_tracks.append(m) - - dd = [] - for m in theor_tracks: - m1 = [] - for a in m: - m1.append(np.array(a)[:,0,:]) - dd.append(m1) - - theor_tracks = dd # Combine all data into a single array of shape: # (N_z, N_age, N_data, N_IMF_interp), where 'N_data' is the number of # sub-arrays in each array. - # comb_data = np.concatenate( - # (a, b, mags_binar, cols_binar, probs_binar, mass_binar, d), axis=2) - - # # Sort all isochrones according to the main magnitude (min to max). - # # This is necessary so that the cut_max_mag() function does not need - # # to do this every time a new synthetic cluster is generated. - # theor_tracks = [[[] for _ in a[0]] for _ in a] - # for i, mx in enumerate(comb_data): - # for j, ax in enumerate(mx): - # # TODO ordering the data according to magnitude is necessary - # # if the cut_max_mag() function expects the data to be sorted - # # this way. This however, prevents us from being able to - # # interpolate new (z, a) values when the MCMC samples them - # # because the mass order is not preserved anymore. So, in order - # # to be able to generate that interpolation of values (which - # # greatly improves the ptemcee performance), we're back to the - # # 'old' cut_max_mag() function and no magnitude ordering. - # # theor_tracks[i][j] = ax[:, ax[0].argsort(kind='mergesort')] - # theor_tracks[i][j] = ax + comb_data = np.concatenate( + (a, b, mags_binar, cols_binar, probs_binar, mass_binar, d), axis=2) + + # Sort all isochrones according to the main magnitude (min to max). + # This is necessary so that the cut_max_mag() function does not need + # to do this every time a new synthetic cluster is generated. + theor_tracks = [[[] for _ in a[0]] for _ in a] + for i, mx in enumerate(comb_data): + for j, ax in enumerate(mx): + # TODO ordering the data according to magnitude is necessary + # if the cut_max_mag() function expects the data to be sorted + # this way. This however, prevents us from being able to + # interpolate new (z, a) values when the MCMC samples them + # because the mass order is not preserved anymore. So, in order + # to be able to generate that interpolation of values (which + # greatly improves the ptemcee performance), we're back to the + # 'old' cut_max_mag() function and no magnitude ordering. + # theor_tracks[i][j] = ax[:, ax[0].argsort(kind='mergesort')] + theor_tracks[i][j] = ax # The above sorting destroys the original order of the isochrones. This # results in messy plots for the "best fit isochrone" at the end. @@ -224,56 +203,23 @@ def arrange_filters(isoch_list, all_syst_filters, filters, colors): return mags_theor, cols_theor, mags_cols_theor -def interp_isoch_data(mags_theor, cols_theor, mags_cols_theor, extra_pars): - """ +def interp_isoch_data(data, N): + ''' Interpolate extra values for all the parameters in the theoretic isochrones. - """ - # mags_theor, cols_theor, mags_cols, extra_pars --> list - # list = [m0, m1, .., mN] <-- N metallicities - # mx = [a0, a1, .., aM] <-- M log(ages) - # ax = [[s1, s2, .., sQ]] <-- Q stars - # list[i][j][0] --> i metallicity, j age - - mass_step = 0.01 - mass_vals1 = np.arange(.01, 5., mass_step) - mass_step = 0.001 - mass_vals2 = np.arange(5., 18., mass_step) - mass_step = 0.0005 - mass_vals3 = np.arange(18., 150., mass_step) - mass_vals = np.array( - mass_vals1.tolist() + mass_vals2.tolist() + mass_vals3.tolist()) - - all_interp_data = [] - for data in (mags_theor, cols_theor, mags_cols_theor, extra_pars): - interp_data = [] - # For each metallicity value. - for i, met in enumerate(data): - m = [] - # For each age value. - for j, age in enumerate(met): - a = [] - # For each filter/color/extra parameter in list. - for fce in age: - # Initial masses for this isochrone. - mass_ini = np.array(extra_pars[i][j][0]) - # mass values must be in the range of the initial masses - mmin = np.searchsorted(mass_vals, mass_ini.min()) - mmax = np.searchsorted(mass_vals, mass_ini.max()) - mass_vals_cut = mass_vals[mmin:mmax] - - y = np.array(data[i][j][0]) - mfunc = interp1d(mass_ini, y) - ynew = mfunc(mass_vals_cut) - a.append(ynew) - - m.append(a) - interp_data.append(m) - all_interp_data.append(interp_data) - - # np.shape(all_interp_data): (4, Nz, Na) - # 4 : ((mags, cols, mags_cols, extra_pars)) - # Nz : metallicity files (values) - # Na : log(age) values per metallicity file + ''' + interp_data = [] + # For each metallicity value list. + for met in data: + m = [] + # For each age value list. + for age in met: + a = [] + # For each filter/color/extra parameter in list. + for fce in age: + t, xp = np.linspace(0., 1., N), np.linspace(0, 1, len(fce)) + a.append(np.interp(t, xp, fce)) + m.append(a) + interp_data.append(m) - return np.array(all_interp_data) + return interp_data From 91c80ccf11358fb26c6353dea0d0a925b8c69b58 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Wed, 2 Oct 2019 17:56:49 -0300 Subject: [PATCH 22/27] use no mass alignment, and re-format model_proper --- packages/best_fit/zaWAverage.py | 102 +++++++++++--------------------- 1 file changed, 35 insertions(+), 67 deletions(-) diff --git a/packages/best_fit/zaWAverage.py b/packages/best_fit/zaWAverage.py index 627edfcb..cacac9b8 100644 --- a/packages/best_fit/zaWAverage.py +++ b/packages/best_fit/zaWAverage.py @@ -2,7 +2,7 @@ import numpy as np -def main(theor_tracks, fundam_params, varIdxs, model, m_ini): +def main(theor_tracks, fundam_params, varIdxs, model): """ Average a new (weighted) isochrone from the four closest points in the (z, a) grid. @@ -40,15 +40,12 @@ def main(theor_tracks, fundam_params, varIdxs, model, m_ini): pts = np.array([(z1, a1), (z1, a2), (z2, a1), (z2, a2)]) # Order: (z1, a1), (z1, a2), (z2, a1), (z2, a2) - isochs = [ + isochs = np.array([ theor_tracks[ml][al], theor_tracks[ml][ah], theor_tracks[mh][al], - theor_tracks[mh][ah]] + theor_tracks[mh][ah]]) - # Define weights for the average, based on the inverse of the distance - # to the grid (z, age) points. - - # Inverse distances between the (z, a) points in the 'model', and the - # four points in the (z, a) grid that contain the model point. + # Distances between the (z, a) points in the 'model', and the four + # points in the (z, a) grid that contain the model point. # Fast euclidean distance: https://stackoverflow.com/a/47775357/1391441 a_min_b = np.array([(z_model, a_model)]) - pts dist = np.sqrt(np.einsum('ij,ij->i', a_min_b, a_min_b)) @@ -64,70 +61,40 @@ def main(theor_tracks, fundam_params, varIdxs, model, m_ini): # Inverse of the distance. inv_d = 1. / dist - # Sort according to smaller distances to the (z_model, a_model) point. - # isoch_idx = np.argsort(inv_d)[::-1] - - # Indexes that identify which of the four arrays in 'x-all' contain the - # the largest minimum and the smallest maximum in their '-1' sub-array - vmin = np.argmax([isochs[_][m_ini][0] for _ in range(4)]) - vmax = np.argmin([isochs[_][m_ini][-1] for _ in range(4)]) - # Identify lower and upper range limits for the '-1' sub-array for all - # arrays in 'isochs' - xmin, xmax = isochs[vmin][m_ini][0], isochs[vmax][m_ini][-1] - - # Align from bottom and top - for i in [0, 1, 2, 3]: - imin = np.searchsorted(isochs[i][m_ini], xmin) - imax = np.searchsorted(isochs[i][m_ini], xmax) - isochs[i] = isochs[i][:, imin:imax] - - # This block "aligns" the masses such that if the 'mass distance' between a - # star in the closest grid isochrone to the model (x1) and a star from any - # of the remaining three isochrones (x2, x3, x4) is larger than a given - # max value (.01, then we replace the star in the second isochrone with its - # 'aligned' star taken from the closest isochrone (x1). - # x1, x2, x3, x4 = isochs[isoch_idx] - # for x in (x2, x3, x4): - # # Mask according to the maximum mass difference allowed - # msk = abs(x1[m_ini] - x[m_ini]) > .01 - # # Replace stars in 'x' with large mass differences, with stars in 'x1' - # np.copyto(x, x1, where=msk) - # Weighted average with mass "alignment". This way is faster than using # 'np.average()'. # Scale weights so they add up to 1. weights = inv_d / np.sum(inv_d) + isochrone = isochs[0] * weights[0] + isochs[1] * weights[1] +\ isochs[2] * weights[2] + isochs[3] * weights[3] # nn = np.random.randint(0, 100) - # if nn == 50: - print(model) - print(model_proper) - import matplotlib.pyplot as plt - plt.subplot(131) - plt.scatter(*pts.T, c='r') - # plt.scatter(*pts[isoch_idx][1:].T, c='g') - plt.scatter(model[0], model[1], marker='x', c='g') - # First color - plt.subplot(132) - plt.plot(theor_tracks[ml][al][1], theor_tracks[ml][al][0], c='b') - plt.plot(theor_tracks[ml][ah][1], theor_tracks[ml][ah][0], c='r') - plt.plot(theor_tracks[mh][al][1], theor_tracks[mh][al][0], c='cyan') - plt.plot(theor_tracks[mh][ah][1], theor_tracks[mh][ah][0], c='orange') - plt.plot(isochrone[1], isochrone[0], c='g', ls='--') - plt.gca().invert_yaxis() - # Second color - plt.subplot(133) - plt.plot(theor_tracks[ml][al][-1], theor_tracks[ml][al][0], c='b') - plt.plot(theor_tracks[ml][ah][-1], theor_tracks[ml][ah][0], c='r') - plt.plot(theor_tracks[mh][al][-1], theor_tracks[mh][al][0], c='cyan') - plt.plot(theor_tracks[mh][ah][-1], theor_tracks[mh][ah][0], c='orange') - plt.plot(isochrone[-1], isochrone[0], c='g', ls='--') - plt.gca().invert_yaxis() - # plt.show() - import pdb; pdb.set_trace() # breakpoint d9c1b180 // - + # if nn == 50: + # print(model) + # print(model_proper) + # import matplotlib.pyplot as plt + # plt.subplot(131) + # plt.scatter(*pts.T, c='r') + # # plt.scatter(*pts[isoch_idx][1:].T, c='g') + # plt.scatter(model[0], model[1], marker='x', c='g') + # # First color + # plt.subplot(132) + # plt.plot(theor_tracks[ml][al][1], theor_tracks[ml][al][0], c='b') + # plt.plot(theor_tracks[ml][ah][1], theor_tracks[ml][ah][0], c='r') + # plt.plot(theor_tracks[mh][al][1], theor_tracks[mh][al][0], c='cyan') + # plt.plot(theor_tracks[mh][ah][1], theor_tracks[mh][ah][0], c='orange') + # plt.plot(isochrone[1], isochrone[0], c='g', ls='--') + # plt.gca().invert_yaxis() + # # Second color + # plt.subplot(133) + # plt.plot(theor_tracks[ml][al][-1], theor_tracks[ml][al][0], c='b') + # plt.plot(theor_tracks[ml][ah][-1], theor_tracks[ml][ah][0], c='r') + # plt.plot(theor_tracks[mh][al][-1], theor_tracks[mh][al][0], c='cyan') + # plt.plot(theor_tracks[mh][ah][-1], theor_tracks[mh][ah][0], c='orange') + # plt.plot(isochrone[-1], isochrone[0], c='g', ls='--') + # plt.gca().invert_yaxis() + # plt.show() return isochrone, model_proper @@ -157,7 +124,7 @@ def properModel(fundam_params, model, varIdxs): # Select the closest value in the array of allowed values. mh = min(len(par) - 1, np.searchsorted(par, model[i - j])) ml = mh - 1 - model_proper.append(par[mh]) + # model_proper.append(par[mh]) # Define model's z value z_model = model[i - j] # If it is the parameter log(age). @@ -165,7 +132,7 @@ def properModel(fundam_params, model, varIdxs): # Select the closest value in the array of allowed values. ah = min(len(par) - 1, np.searchsorted(par, model[i - j])) al = ah - 1 - model_proper.append(par[ah]) + # model_proper.append(par[ah]) a_model = model[i - j] else: model_proper.append(model[i - j]) @@ -176,7 +143,8 @@ def properModel(fundam_params, model, varIdxs): elif i == 1: al = ah = 0 a_model = fundam_params[1][0] - model_proper.append(par[0]) + else: + model_proper.append(par[0]) j += 1 return model_proper, z_model, a_model, ml, mh, al, ah From 39463ae6dd3d1beb4462648579d8ef765b3bb882 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 3 Oct 2019 10:07:03 -0300 Subject: [PATCH 23/27] close #441 and several other minor issues mentioned there --- packages/best_fit/bf_common.py | 51 +++++++++++------------- packages/best_fit/bootstrap.py | 13 ++++-- packages/best_fit/genetic_algorithm.py | 41 ++++++++++--------- packages/best_fit/ptemcee_algor.py | 2 +- packages/best_fit/zaWAverage.py | 30 +++++++++----- packages/check/read_met_files.py | 5 +-- packages/inp/isoch_params.py | 22 ++++++----- packages/out/synth_cl_file.py | 18 ++++++--- packages/synth_clust/synth_cl_plot.py | 55 ++++++++++---------------- packages/synth_clust/synth_cluster.py | 6 +-- 10 files changed, 128 insertions(+), 115 deletions(-) diff --git a/packages/best_fit/bf_common.py b/packages/best_fit/bf_common.py index 425e101a..204ac04b 100644 --- a/packages/best_fit/bf_common.py +++ b/packages/best_fit/bf_common.py @@ -46,21 +46,22 @@ def fillParams(fundam_params, varIdxs, model): return model_filled -def closeSol(fundam_params, model, pushidxs): - """ - Find the closest value in the parameters list for the discrete parameters - metallicity, age, and mass. - """ - model_proper = [] - for i, par in enumerate(fundam_params): - # If it is the parameter metallicity, age or mass. - if i in pushidxs: - # Select the closest value in the array of allowed values. - model_proper.append(min(par, key=lambda x: abs(x - model[i]))) - else: - model_proper.append(model[i]) +# DEPRECATED 02-10-2019 +# def closeSol(fundam_params, model, pushidxs): +# """ +# Find the closest value in the parameters list for the discrete parameters +# metallicity, age, and mass. +# """ +# model_proper = [] +# for i, par in enumerate(fundam_params): +# # If it is the parameter metallicity, age or mass. +# if i in pushidxs: +# # Select the closest value in the array of allowed values. +# model_proper.append(min(par, key=lambda x: abs(x - model[i]))) +# else: +# model_proper.append(model[i]) - return model_proper +# return model_proper def initPop( @@ -109,32 +110,28 @@ def random_population(fundam_params, varIdxs, n_ran): """ Generate a random set of parameter values to use as a random population. - Pick n_ran initial random solutions from each list storing all the - possible parameters values. + Pick n_ran initial random solutions from the valid ranges. """ p_lst = [] for i in varIdxs: - p_lst.append([random.choice(fundam_params[i]) for _ in range(n_ran)]) + p_lst.append( + np.random.uniform( + fundam_params[i][0], fundam_params[i][-1], n_ran)) return np.array(p_lst).T -def synthClust(fundam_params, varIdxs, model, synthcl_args): +def synthClust(fundam_params, varIdxs, synthcl_args, model): """ Generate a synthetic cluster. """ - theor_tracks, e_max, err_lst, completeness, max_mag_syn, st_dist_mass,\ - R_V, ext_coefs, N_fc, m_ini, cmpl_rnd, err_rnd = synthcl_args - - # Interpolate (z, a) data + # Average a new isochrone + # theor_tracks = synthcl_args[0] isochrone, model_proper = zaWAverage.main( - theor_tracks, fundam_params, varIdxs, model, m_ini) + synthcl_args[0], fundam_params, varIdxs, model) # Generate synthetic cluster. - return synth_cluster.main( - e_max, err_lst, completeness, max_mag_syn, st_dist_mass, isochrone, - R_V, ext_coefs, N_fc, m_ini, cmpl_rnd, err_rnd, model_proper) - + return synth_cluster.main(isochrone, model_proper, *synthcl_args[1:]) # DEPRECATED 24-09-2019 # def discreteParams(fundam_params, varIdxs, chains_nruns, pushidxs): diff --git a/packages/best_fit/bootstrap.py b/packages/best_fit/bootstrap.py index 65c3f7ca..ae55be91 100644 --- a/packages/best_fit/bootstrap.py +++ b/packages/best_fit/bootstrap.py @@ -27,6 +27,9 @@ def main( else: available_secs = max_secs * (1. - pd['hperc_btstrp']) + # Identify free parameters + varIdxs = varPars(pd['fundam_params'])[0] + # # TODO not yet fully implemented # pd['best_fit_algor'] = 'boot+DE' @@ -35,8 +38,12 @@ def main( # the observed data. Use a large number for 'N_gen' N_pop_init, N_gen = pd['N_pop'], int(5e5) # Use random initial population - init_pop = random_population( - pd['fundam_params'], (0, 1, 2, 3, 4, 5), N_pop_init).tolist() + init_pop_free = random_population( + pd['fundam_params'], varIdxs, N_pop_init) + init_pop = [] + for model in init_pop_free: + init_pop.append(fillParams(pd['fundam_params'], varIdxs, model)) + flag_print_perc = True argsOF = [ pd, clp, max_mag_syn, obs_clust, ext_coefs, st_dist_mass, N_fc, @@ -53,6 +60,7 @@ def main( # First run of the numerical optimizer function with the observed data. isoch_fit_params = optimizerFunc(pd['best_fit_algor'], argsOF) + isoch_fit_params['varIdxs'] = varIdxs # Check time available for bootstrap after the previous run. available_secs = max_secs - isoch_fit_params['OF_elapsed'] @@ -96,7 +104,6 @@ def main( # As array with (params, btstrp_runs) shape isoch_fit_params['params_boot'] = np.array(params_boot).T - isoch_fit_params['varIdxs'] = varPars(pd['fundam_params'])[0] # Remove fixed parameters from array. isoch_fit_params['params_boot'] = isoch_fit_params['params_boot'][ isoch_fit_params['varIdxs']] diff --git a/packages/best_fit/genetic_algorithm.py b/packages/best_fit/genetic_algorithm.py index 4c47172a..0066e09f 100644 --- a/packages/best_fit/genetic_algorithm.py +++ b/packages/best_fit/genetic_algorithm.py @@ -3,7 +3,8 @@ import textwrap import random import numpy as np -from .bf_common import varPars, rangeCheck, synthClust, random_population +from .bf_common import varPars, rangeCheck, synthClust, random_population,\ + fillParams from . import likelihood ############################################################# @@ -57,11 +58,15 @@ def main( varIdxs, ndim, ranges = varPars(fundam_params) + # Pack synthetic cluster arguments. + synthcl_args = [ + theor_tracks, e_max, err_lst, completeness, max_mag_syn, st_dist_mass, + R_V, ext_coefs, N_fc, m_ini, cmpl_rnd, err_rnd] + # Evaluate initial random solutions in the objective function. generation, lkl = evaluation( - lkl_method, e_max, err_lst, completeness, max_mag_syn, fundam_params, - obs_clust, theor_tracks, R_V, ext_coefs, st_dist_mass, N_fc, m_ini, - cmpl_rnd, err_rnd, varIdxs, ranges, ran_pop) + lkl_method, fundam_params, obs_clust, varIdxs, ranges, synthcl_args, + ran_pop) # Store best solution(s) for passing along in the 'Elitism' block. best_sol = generation[:N_el] @@ -117,10 +122,8 @@ def main( # according to the best solutions found. # with timeblock("Evaluation"): generation, lkl = evaluation( - lkl_method, e_max, err_lst, completeness, max_mag_syn, - fundam_params, obs_clust, theor_tracks, R_V, ext_coefs, - st_dist_mass, N_fc, m_ini, cmpl_rnd, err_rnd, varIdxs, ranges, - p_lst_e) + lkl_method, fundam_params, obs_clust, varIdxs, ranges, + synthcl_args, p_lst_e) if flag_print_perc: @@ -188,7 +191,7 @@ def main( ext_imm_count = 0 # Apply Extinction/Immigration operator. - generation = ext_imm(best_sol, fundam_params, N_popl) + generation = ext_imm(varIdxs, best_sol, fundam_params, N_popl) # Reset best solution counter. best_sol_count = 0 else: @@ -306,9 +309,8 @@ def decode(fundam_params, n_bin, p_delta, p_mins, mut_chrom): def evaluation( - lkl_method, e_max, err_lst, completeness, max_mag_syn, fundam_params, - obs_clust, theor_tracks, R_V, ext_coefs, st_dist_mass, N_fc, m_ini, - cmpl_rnd, err_rnd, varIdxs, ranges, p_lst): + lkl_method, + fundam_params, obs_clust, varIdxs, ranges, synthcl_args, p_lst): ''' Evaluate each model in the objective function to obtain the fitness of each one. @@ -325,10 +327,7 @@ def evaluation( # Generate synthetic cluster. synth_clust = synthClust( - fundam_params, varIdxs, model, ( - theor_tracks, e_max, err_lst, completeness, max_mag_syn, - st_dist_mass, R_V, ext_coefs, N_fc, m_ini, cmpl_rnd, - err_rnd)) + fundam_params, varIdxs, synthcl_args, model) # Call likelihood function for this model. # with timeblock(" Likelihood"): @@ -454,16 +453,20 @@ def elitism(best_sol, p_lst): return p_lst_r -def ext_imm(best_sol, fundam_params, N_popl): +def ext_imm(varIdxs, best_sol, fundam_params, N_popl): ''' Append a new random population to the best solution so far. ''' # Generate (N_popl-N_el) random solutions. n_ran = N_popl - len(best_sol) - p_lst_r = random_population(fundam_params, (0, 1, 2, 3, 4, 5), n_ran) + + p_lst_r_free = random_population(fundam_params, varIdxs, n_ran) + p_lst_r = [] + for model in p_lst_r_free: + p_lst_r.append(fillParams(fundam_params, varIdxs, model)) # Append immigrant random population to the best solution. - generation_ei = best_sol + p_lst_r.tolist() + generation_ei = best_sol + p_lst_r return generation_ei diff --git a/packages/best_fit/ptemcee_algor.py b/packages/best_fit/ptemcee_algor.py index ed7425dd..9410f7af 100644 --- a/packages/best_fit/ptemcee_algor.py +++ b/packages/best_fit/ptemcee_algor.py @@ -247,7 +247,7 @@ def loglkl( logpost = -1e9 if rangeFlag: # Generate synthetic cluster. - synth_clust = synthClust(fundam_params, varIdxs, model, synthcl_args) + synth_clust = synthClust(fundam_params, varIdxs, synthcl_args, model) # Call likelihood function for this model. lkl = likelihood.main(lkl_method, synth_clust, obs_clust) log_p = 0. diff --git a/packages/best_fit/zaWAverage.py b/packages/best_fit/zaWAverage.py index cacac9b8..17d89b5c 100644 --- a/packages/best_fit/zaWAverage.py +++ b/packages/best_fit/zaWAverage.py @@ -24,8 +24,6 @@ def main(theor_tracks, fundam_params, varIdxs, model): """ - # Define the 'proper' model with values for (z, a) taken from its grid, - # and filled values for those parameters that are fixed. model_proper, z_model, a_model, ml, mh, al, ah = properModel( fundam_params, model, varIdxs) @@ -61,14 +59,21 @@ def main(theor_tracks, fundam_params, varIdxs, model): # Inverse of the distance. inv_d = 1. / dist - # Weighted average with mass "alignment". This way is faster than using - # 'np.average()'. + # Weighted average. This way is faster than using 'np.average()'. # Scale weights so they add up to 1. weights = inv_d / np.sum(inv_d) - isochrone = isochs[0] * weights[0] + isochs[1] * weights[1] +\ isochs[2] * weights[2] + isochs[3] * weights[3] + # The above averaging assumes that the four isochrones are mass-aligned. + # This is strictly *not true*, but since the four isochrones have + # relatively similar values, it should be close enough to true. + # As more points are interpolated into the theoretical isochrones, this + # will work better. + # + # In my tests, attempting an alignment worsens the final averaged + # isochrone. + # nn = np.random.randint(0, 100) # if nn == 50: # print(model) @@ -101,12 +106,19 @@ def main(theor_tracks, fundam_params, varIdxs, model): def properModel(fundam_params, model, varIdxs): """ + Define the 'proper' model with values for (z, a) taken from its grid, + and filled values for those parameters that are fixed. + + Parameters + ---------- + model : array + Array of *free* fundamental parameters only (ie: not in varIdxs). Returns ------- model_proper : list - Stores the closest (z, a) values in the grid for the parameters in - 'model', and add the fixed parameters that are missing from 'model'. + Stores (E_BV, dm, Mass, b_fr) including the fixed parameters that are + missing from 'model'. z_model, a_model : floats The (z, a) values for this model's isochrone. ml, mh, al, ah : ints @@ -124,15 +136,13 @@ def properModel(fundam_params, model, varIdxs): # Select the closest value in the array of allowed values. mh = min(len(par) - 1, np.searchsorted(par, model[i - j])) ml = mh - 1 - # model_proper.append(par[mh]) - # Define model's z value + # Define the model's z value z_model = model[i - j] # If it is the parameter log(age). elif i == 1: # Select the closest value in the array of allowed values. ah = min(len(par) - 1, np.searchsorted(par, model[i - j])) al = ah - 1 - # model_proper.append(par[ah]) a_model = model[i - j] else: model_proper.append(model[i - j]) diff --git a/packages/check/read_met_files.py b/packages/check/read_met_files.py index cee124c3..0c2364d2 100644 --- a/packages/check/read_met_files.py +++ b/packages/check/read_met_files.py @@ -12,7 +12,7 @@ def check_get(pd): store the data and pass it. """ # Only read files if best fit method is set to run. Else pass empty list. - pd['fundam_params'], pd['theor_tracks'], pd['plot_isoch_data'] = [], [], [] + pd['fundam_params'], pd['theor_tracks'] = [], [] if pd['bf_flag']: # Read metallicity files' names, store proper ranges for all @@ -39,8 +39,7 @@ def check_get(pd): pd['fundam_params'] = [met_values[0], age_values] + params_values[2:] # Store all isochrones in all the metallicity files. - pd['theor_tracks'], pd['plot_isoch_data'] = isoch_params.main( - met_f_filter, age_values, **pd) + pd['theor_tracks'] = isoch_params.main(met_f_filter, age_values, **pd) return pd diff --git a/packages/inp/isoch_params.py b/packages/inp/isoch_params.py index a3ba6d52..79545234 100644 --- a/packages/inp/isoch_params.py +++ b/packages/inp/isoch_params.py @@ -109,14 +109,18 @@ def main(met_f_filter, age_values, cmd_evol_tracks, evol_track, bin_mr, # theor_tracks[i][j] = ax[:, ax[0].argsort(kind='mergesort')] theor_tracks[i][j] = ax - # The above sorting destroys the original order of the isochrones. This - # results in messy plots for the "best fit isochrone" at the end. - # (see: https://stackoverflow.com/q/35606712/1391441, - # https://stackoverflow.com/q/37742358/1391441) - # To avoid this, we also store the not-interpolated, not sorted original - # values for the magnitudes and colors; just for the purpose of plotting - # the final isochrones. - plot_isoch_data = np.concatenate((mags_theor, cols_theor), axis=2) + # DEPRECATED 02-10-2019: not needed anymore since there is no sorting, + # as stated above, and the (z, a) final values can now be averaged from + # grid values. + # + # # The above sorting destroys the original order of the isochrones. This + # # results in messy plots for the "best fit isochrone" at the end. + # # (see: https://stackoverflow.com/q/35606712/1391441, + # # https://stackoverflow.com/q/37742358/1391441) + # # To avoid this, we also store the not-interpolated, not sorted original + # # values for the magnitudes and colors; just for the purpose of plotting + # # the final isochrones. + # plot_isoch_data = np.concatenate((mags_theor, cols_theor), axis=2) print("\nGrid values: {} [z], {} [log(age)]".format( len(fundam_params[0]), len(fundam_params[1]))) @@ -127,7 +131,7 @@ def main(met_f_filter, age_values, cmd_evol_tracks, evol_track, bin_mr, # theor_tracks_mp = np.memmap(filename, dtype='float64', mode='w+', shape=sp) # theor_tracks_mp[:] = theor_tracks[:] - return theor_tracks, plot_isoch_data + return theor_tracks def arrange_filters(isoch_list, all_syst_filters, filters, colors): diff --git a/packages/out/synth_cl_file.py b/packages/out/synth_cl_file.py index 3d6ea14c..e48a533f 100644 --- a/packages/out/synth_cl_file.py +++ b/packages/out/synth_cl_file.py @@ -3,18 +3,26 @@ def main(clp, npd, bf_flag, best_fit_algor, fundam_params, filters, colors, - theor_tracks, plot_isoch_data, R_V, **kwargs): + theor_tracks, R_V, **kwargs): ''' Create output data file with stars in the best fit synthetic cluster found by the 'Best Fit' function. ''' clp['synth_clst'], clp['shift_isoch'] = [], [] if bf_flag: + + # Index of m_ini (theoretical initial mass), stored in the theoretical + # isochrones. + m_ini = 2 * clp['N_fc'][0] + 2 * clp['N_fc'][1] + 2 + # Pack synthetic cluster arguments. + synthcl_args = [ + theor_tracks, clp['em_float'], clp['err_lst'], clp['completeness'], + clp['max_mag_syn'], clp['st_dist_mass'], R_V, clp['ext_coefs'], + clp['N_fc'], m_ini, clp['cmpl_rnd'], clp['err_rnd']] + shift_isoch, synth_clst = synth_cl_plot.main( - best_fit_algor, fundam_params, theor_tracks, plot_isoch_data, R_V, - clp['em_float'], clp['isoch_fit_params'], clp['err_lst'], - clp['completeness'], clp['max_mag_syn'], clp['st_dist_mass'], - clp['ext_coefs'], clp['N_fc'], clp['cmpl_rnd'], clp['err_rnd']) + best_fit_algor, fundam_params, clp['isoch_fit_params'], + synthcl_args) # If cluster is not empty. if synth_clst: diff --git a/packages/synth_clust/synth_cl_plot.py b/packages/synth_clust/synth_cl_plot.py index facce82e..5196044e 100644 --- a/packages/synth_clust/synth_cl_plot.py +++ b/packages/synth_clust/synth_cl_plot.py @@ -1,51 +1,36 @@ +import numpy as np from . import move_isochrone from . import synth_cluster -from ..best_fit.bf_common import closeSol +# from ..best_fit.bf_common import closeSol +from ..best_fit import zaWAverage -def main( - best_fit_algor, fundam_params, theor_tracks, plot_isoch_data, R_V, - e_max, isoch_fit_params, err_lst, completeness, max_mag_syn, st_dist_mass, - ext_coefs, N_fc, cmpl_rnd, err_rnd): - ''' - # Generate shifted isochrone and synthetic cluster for plotting. - ''' +def main(best_fit_algor, fundam_params, isoch_fit_params, synthcl_args): + """ + Generate shifted isochrone and synthetic cluster for plotting. + """ if best_fit_algor == 'boot+GA': # Use ML fit values for all parameters. - synth_cl_params = isoch_fit_params['map_sol'] + model = isoch_fit_params['map_sol'] elif best_fit_algor == 'ptemcee': # Use mean fit values for all parameters. - synth_cl_params = isoch_fit_params['mean_sol'] - # Grid values for (z, a) - synth_cl_params_grid = closeSol(fundam_params, synth_cl_params, [0, 1]) - - # Find indexes for metallicity and age. If indexes are not found due - # to some difference in the significant figures, use the indexes - # [0, 0] to prevent the code from halting. - try: - m_i = fundam_params[0].index(synth_cl_params_grid[0]) - except Exception: - m_i = 0 - print(" WARNING: metallicity for the best match synthetic\n" - " cluster not found.") - try: - a_i = fundam_params[1].index(synth_cl_params_grid[1]) - except Exception: - a_i = 0 - print(" WARNING: age for the best match synthetic\n" - " cluster not found.") + model = isoch_fit_params['mean_sol'] + + theor_tracks = synthcl_args[0] + R_V, ext_coefs, N_fc = synthcl_args[6:9] + E_BV, dm = model[2], model[3] + + # Generate a model with the fitted parameters only. + model_var = np.array(model)[isoch_fit_params['varIdxs']] # Generate shifted best fit isochrone. - isochrone = plot_isoch_data[m_i][a_i] + isochrone, _ = zaWAverage.main( + theor_tracks, fundam_params, isoch_fit_params['varIdxs'], model_var) shift_isoch = move_isochrone.main( - isochrone, synth_cl_params[2], synth_cl_params[3], R_V, ext_coefs, - N_fc) + isochrone, E_BV, dm, R_V, ext_coefs, N_fc) # Generate best fit synthetic cluster. - isochrone = theor_tracks[m_i][a_i] - synth_clst = synth_cluster.main( - e_max, err_lst, completeness, max_mag_syn, st_dist_mass, isochrone, - R_V, ext_coefs, N_fc, cmpl_rnd, err_rnd, synth_cl_params_grid) + synth_clst = synth_cluster.main(isochrone, model[2:], *synthcl_args[1:]) return shift_isoch, synth_clst diff --git a/packages/synth_clust/synth_cluster.py b/packages/synth_clust/synth_cluster.py index acc8f463..779698ff 100644 --- a/packages/synth_clust/synth_cluster.py +++ b/packages/synth_clust/synth_cluster.py @@ -9,8 +9,8 @@ def main( - err_max, err_lst, completeness, max_mag_syn, st_dist_mass, isochrone, - R_V, ext_coefs, N_fc, m_ini, cmpl_rnd, err_rnd, synth_cl_params): + isochrone, model_proper, err_max, err_lst, completeness, max_mag_syn, + st_dist_mass, R_V, ext_coefs, N_fc, m_ini, cmpl_rnd, err_rnd): """ Takes an isochrone and returns a synthetic cluster created according to a certain mass distribution. @@ -36,7 +36,7 @@ def main( # Unpack synthetic cluster parameters. The first two elements are the # metallicity and the age, which are already incorporated in the selected # isochrone. - e, d, M_total, bin_frac = synth_cl_params[2:] + e, d, M_total, bin_frac = model_proper # import time # t1, t2, t3, t4, t5, t6, t7 = 0., 0., 0., 0., 0., 0., 0. From c44e981dfb4f6dd20d9746e1a748982345fbc6ed Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 3 Oct 2019 10:50:28 -0300 Subject: [PATCH 24/27] minor, recommend upgrading --- packages/check/update.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/packages/check/update.py b/packages/check/update.py index 705316f1..a0ae4455 100644 --- a/packages/check/update.py +++ b/packages/check/update.py @@ -32,8 +32,9 @@ def check(): if s[1][1:] != __version__.strip('-dev'): print("*******************************************") print(" IMPORTANT\n") - print("There is a new version of ASteCA available!") - print(" Get the latest version '{}' from:\n".format(s[1][1:])) + print("There is a new version of ASteCA available") + print(" Upgrading is STRONGLY recommended!\n") + print(" Get the latest version '{}' from:".format(s[1][1:])) print(" http://asteca.github.io/") print("*******************************************\n") From ffe666a5da1713f0927ab5717838239a454fcbcb Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 3 Oct 2019 10:51:04 -0300 Subject: [PATCH 25/27] close #438 --- packages/defvals/params_input.dat | 7 ++----- packages/inp/input_params.py | 5 +---- packages/inp/isoch_params.py | 35 +++++++++++++++---------------- 3 files changed, 20 insertions(+), 27 deletions(-) diff --git a/packages/defvals/params_input.dat b/packages/defvals/params_input.dat index 8b73d107..3f87ac1f 100644 --- a/packages/defvals/params_input.dat +++ b/packages/defvals/params_input.dat @@ -478,11 +478,8 @@ ET PAR12 No # * z, log(age): [floats] # Input [z, log(age)] steps to read the stored isochrones. # -# * N_interp: [int] -# Number of extra values to interpolate into the isochrones. -# -# z log(age) N_interp -PS 0.005 0.5 1000 +# z log(age) +PS 0.005 0.5 # IMF used to populate the synthetic clusters. # diff --git a/packages/inp/input_params.py b/packages/inp/input_params.py index b0a3a053..c452c8ad 100644 --- a/packages/inp/input_params.py +++ b/packages/inp/input_params.py @@ -217,7 +217,6 @@ def main(mypath, pars_f_path): colibri = str(reader[2]) elif reader[0] == 'PS': za_steps = list(map(float, reader[1:3])) - N_mass_interp = int(reader[3]) elif reader[0] == 'MF': IMF_name = str(reader[1]) elif reader[0] == 'BR': @@ -321,9 +320,7 @@ def main(mypath, pars_f_path): 'lkl_weight': lkl_weight, # Synthetic cluster parameters 'evol_track': evol_track, 'za_steps': za_steps, - 'max_mag': max_mag, 'IMF_name': IMF_name, - 'N_mass_interp': N_mass_interp, 'R_V': R_V, - 'bin_mr': bin_mr, + 'max_mag': max_mag, 'IMF_name': IMF_name, 'R_V': R_V, 'bin_mr': bin_mr, # parameters ranges 'par_ranges': par_ranges, # ptemcee algorithm parameters. diff --git a/packages/inp/isoch_params.py b/packages/inp/isoch_params.py index 79545234..19fb8715 100644 --- a/packages/inp/isoch_params.py +++ b/packages/inp/isoch_params.py @@ -6,9 +6,9 @@ from .. import update_progress -def main(met_f_filter, age_values, cmd_evol_tracks, evol_track, bin_mr, - all_syst_filters, cmd_systs, filters, colors, fundam_params, - N_mass_interp, **kwargs): +def main( + met_f_filter, age_values, cmd_evol_tracks, evol_track, bin_mr, + all_syst_filters, cmd_systs, filters, colors, fundam_params, **kwargs): ''' Read isochrones and parameters if best fit function is set to run. ''' @@ -45,21 +45,20 @@ def main(met_f_filter, age_values, cmd_evol_tracks, evol_track, bin_mr, for z in mags_theor: for a in z: N_pts_max = max(N_pts_max, len(a[0])) - - if N_mass_interp > N_pts_max: - print("Interpolating extra points ({}) into the isochrones".format( - N_mass_interp)) - interp_data = [] - for i, data in enumerate([ - mags_theor, cols_theor, mags_cols_theor, extra_pars]): - interp_data.append(interp_isoch_data(data, N_mass_interp)) - update_progress.updt(4, i + 1) - a, b, c, d = interp_data - else: - sys.exit( - ("ERROR: N_interp={} must be larger than the maximum\n" + - "number of points in any isochrone read: {}").format( - N_mass_interp, N_pts_max)) + # Fixed value. Will only change if some isochrone contains more than + # this number of stars (not likely) + N_mass_interp = 1500 + if N_mass_interp < N_pts_max: + N_mass_interp = N_pts_max + 100 + + print("Interpolating extra points ({}) into the isochrones".format( + N_mass_interp)) + interp_data = [] + for i, data in enumerate([ + mags_theor, cols_theor, mags_cols_theor, extra_pars]): + interp_data.append(interp_isoch_data(data, N_mass_interp)) + update_progress.updt(4, i + 1) + a, b, c, d = interp_data # # Size of arrays in memory # sz = 0. From 0704150305d87961011784b2df3044c0509be586 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 3 Oct 2019 11:24:53 -0300 Subject: [PATCH 26/27] update changelog --- CHANGELOG.md | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 86ee00a3..039bb305 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,16 @@ # Change Log +## [[v0.2.7]][239] - 2019-10-03 + +### Changed + +* Use inverse transform sampling to sample the IMF ([434][238]) +* Interpolation of (z,a) values uses wrong m_ini index ([440][237]) +* Interpolation of isochrone fails when (z,a) are both fixed ([439][236]) +* Mass 'alignment' in zaInterp() gives poor result ([441][235]) +* Select the N_mass_interp number automatically ([438][234]) + + ## [[v0.2.6]][233] - 2019-09-19 ### Changed @@ -708,3 +719,9 @@ ________________________________________________________________________________ [231]: https://github.com/asteca/ASteCA/issues/428 [232]: https://github.com/asteca/ASteCA/issues/426 [233]: https://github.com/asteca/asteca/releases/tag/v0.2.6 +[234]: https://github.com/asteca/ASteCA/issues/438 +[235]: https://github.com/asteca/ASteCA/issues/441 +[236]: https://github.com/asteca/ASteCA/issues/440 +[237]: https://github.com/asteca/ASteCA/issues/439 +[238]: https://github.com/asteca/ASteCA/issues/434 +[239]: https://github.com/asteca/asteca/releases/tag/v0.2.7 \ No newline at end of file From b7369e80a99f21a8d46215598d66fcd0e363c961 Mon Sep 17 00:00:00 2001 From: Gabriel Date: Thu, 3 Oct 2019 11:25:21 -0300 Subject: [PATCH 27/27] Bumped version number to v0.2.7 --- packages/_version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/_version.py b/packages/_version.py index abf19198..5f3a01b8 100644 --- a/packages/_version.py +++ b/packages/_version.py @@ -1 +1 @@ -__version__ = "v0.2.6-dev" \ No newline at end of file +__version__ = "v0.2.7" \ No newline at end of file