Merge branch 'release-027'

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

packages/_version.py

@@ -1 +1 @@
-__version__ = "v0.2.6"
+__version__ = "v0.2.7"

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):

packages/best_fit/best_fit_synth_cl.py

@@ -34,12 +34,15 @@ 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'])]
 
+        # 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)
@@ -55,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)
 
@@ -67,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)
 

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
 
 
@@ -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,58 +110,55 @@ 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, cmpl_rnd, err_rnd = synthcl_args
-
-    # Interpolate (z, a) data
-    isochrone, model_proper = zaInterp.main(
-        theor_tracks, fundam_params, varIdxs, model)
+    # Average a new isochrone
+    # theor_tracks = synthcl_args[0]
+    isochrone, model_proper = zaWAverage.main(
+        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, 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
+    return synth_cluster.main(isochrone, model_proper, *synthcl_args[1:])
+
+#  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):

packages/best_fit/bootstrap.py

@@ -6,12 +6,12 @@
 from . import obs_clust_prepare
 from .. import update_progress
 from .bf_common import random_population, varPars, modeKDE, fillParams,\
-    r2Dist, closeSol
+    r2Dist
 
 
 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
@@ -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,23 +38,29 @@ 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,
-            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)
+    isoch_fit_params['varIdxs'] = varIdxs
 
     # Check time available for bootstrap after the previous run.
     available_secs = max_secs - isoch_fit_params['OF_elapsed']
@@ -77,14 +86,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))
@@ -95,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']]
@@ -120,8 +128,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(
@@ -136,40 +142,40 @@ 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
 
         if flag_print_perc:
             # 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
@@ -213,11 +219,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