Merge branch 'release-026'

asteca · Sep 19, 2019 · 98ed626 · 98ed626
2 parents 9c64caa + d02ba18
commit 98ed626
Show file tree

Hide file tree

Showing 98 changed files with 1,275 additions and 2,048 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,13 +1,25 @@
 # Change Log
 
+## [[v0.2.6]][233] - 2019-09-19
 
-## [[v0.2.5]][225] - 2019-10-07
+### Changed
+
+* Fix normalization in Bayesian DA ([426][232])
+* Fix function to detect X11 that fails in Mac OS (Windows too?) ([428][231])
+* Merge `semi_input.dat` file into `params_input.dat` and copy input file as output ([427][230])
+* Remove modes ([429][229])
+* Use one photometric systems file instead of two identical ones ([421][228])
+* Fix Ext/Imm operator causing spurious points in the GA ([424][227])
+
+
+## [[v0.2.5]][226] - 2019-08-07
 
 ### Changed
 
-* Added the `ptemcee` method, and deprecated (for now) the BF ([367][224])
-* Accept a CMD/CCD from mixed photometric systems ([228][222], [229][223])
-* Add support for the new set of isochrones PARSEC+COLIBRI ([322][221])
+* Added the `ptemcee` method, and deprecated (for now) the BF ([367][225])
+* Accept a CMD/CCD from mixed photometric systems ([228][223], [229][224])
+* Add support for the new set of isochrones PARSEC+COLIBRI ([322][222])
+* Output all information obtained from the bootstrap ([279][221])
 * Mask stars with photometry outside of reasonable range ([414][220])
 * Add proper motions, parallax, and radial velocity support to Bayesian DA ([220][219])
 * Use stars with no complete data in the Bayesian equation ([377][218]).
@@ -683,8 +695,16 @@ ________________________________________________________________________________
 [218]: https://github.com/asteca/ASteCA/issues/377
 [219]: https://github.com/asteca/ASteCA/issues/220
 [220]: https://github.com/asteca/ASteCA/issues/414
-[221]: https://github.com/asteca/ASteCA/issues/322
-[222]: https://github.com/asteca/ASteCA/issues/228
-[223]: https://github.com/asteca/ASteCA/issues/229
-[224]: https://github.com/asteca/ASteCA/issues/367
-[225]: https://github.com/asteca/asteca/releases/tag/v0.2.5
+[221]: https://github.com/asteca/ASteCA/issues/279
+[222]: https://github.com/asteca/ASteCA/issues/322
+[223]: https://github.com/asteca/ASteCA/issues/228
+[224]: https://github.com/asteca/ASteCA/issues/229
+[225]: https://github.com/asteca/ASteCA/issues/367
+[226]: https://github.com/asteca/asteca/releases/tag/v0.2.5
+[227]: https://github.com/asteca/ASteCA/issues/424
+[228]: https://github.com/asteca/ASteCA/issues/421
+[229]: https://github.com/asteca/ASteCA/issues/429
+[230]: https://github.com/asteca/ASteCA/issues/427
+[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
diff --git a/packages/_version.py b/packages/_version.py
@@ -1 +1 @@
-__version__ = "v0.2.5"
+__version__ = "v0.2.6"
diff --git a/packages/best_fit/best_fit_synth_cl.py b/packages/best_fit/best_fit_synth_cl.py
@@ -45,12 +45,12 @@ def main(clp, pd):
         cmpl_rnd = np.random.uniform(0., 1., 1000000)
         err_rnd = np.random.normal(0., 1., 1000000)
 
-        print('Searching for optimal parameters.')
+        print("Searching for optimal parameters")
 
         # Calculate the best fitting parameters.
         if pd['best_fit_algor'] == 'boot+GA':
 
-            print('Using bootstrap + Genetic Algorithm ({}).'.format(
+            print("Using bootstrap + Genetic Algorithm ({})".format(
                 pd['lkl_method'] + '; ' + pd['lkl_binning'] if
                 pd['lkl_method'] == 'dolphin' else pd['lkl_method']))
             isoch_fit_params = bootstrap.main(
@@ -60,7 +60,7 @@ def main(clp, pd):
             isoch_fit_errors = params_errors(pd, isoch_fit_params)
 
         elif pd['best_fit_algor'] == 'ptemcee':
-            print('Using ptemcee algorithm ({}).'.format(
+            print("Using ptemcee algorithm ({})".format(
                 pd['lkl_method'] + '; ' + pd['lkl_binning'] if
                 pd['lkl_method'] == 'dolphin' else pd['lkl_method']))
 
@@ -99,14 +99,14 @@ def main(clp, pd):
         #     isoch_fit_errors, _ = params_errors(
         #         best_fit_algor, isoch_fit_params)
 
-        print("Best fit parameters obtained.")
+        print("Best fit parameters obtained")
 
         clp['max_mag_syn'], clp['ext_coefs'], clp['st_dist_mass'], \
             clp['N_fc'], clp['cmpl_rnd'], clp['err_rnd'], =\
             max_mag_syn, ext_coefs, st_dist_mass, N_fc, cmpl_rnd, err_rnd
 
     else:
-        print('Skip parameters fitting process.')
+        print("Skip parameters fitting process")
         # Pass dummy data used by data output and some plot functions.
         cl_max_mag = []
         isoch_fit_params = {

diff --git a/packages/best_fit/bf_common.py b/packages/best_fit/bf_common.py
@@ -73,11 +73,12 @@ def initPop(
 
     p0 = []
     if init_mode == 'random':
-        print("Random initial population")
+        # print("Random initial population")
         for _ in range(ntemps):
             p0.append(random_population(fundam_params, varIdxs, nwalkers_ptm))
 
     elif init_mode == 'diffevol':
+        # DEPRECATED 05/09/19 when #425 was implemented
         print("     DE init pop")
 
         with warnings.catch_warnings():

diff --git a/packages/best_fit/bootstrap.py b/packages/best_fit/bootstrap.py
@@ -60,7 +60,7 @@ def main(
     params_boot, N_steps_availb, btstrp_start_t = [], 0, t.time()
     if pd['hperc_btstrp'] > 0.:
 
-        print('Begin bootstrap process ({} | {}).'.format(
+        print("Begin bootstrap process ({} | {})".format(
             pd['N_pop_btstrp'], pd['N_step_btstrp']))
 
         btstrp_init, btstrp_t, btstrp_runs = t.time(), 0., 0
@@ -112,7 +112,7 @@ def main(
         isoch_fit_params['N_bootstrap'] = btstrp_runs
 
     else:
-        print('Skip bootstrap process.')
+        print("Skip bootstrap process")
         isoch_fit_params['N_bootstrap'] = 0
         isoch_fit_params['params_boot'] = np.array(params_boot).T
         isoch_fit_params['mean_sol'] = [np.nan] * 6

diff --git a/packages/best_fit/genetic_algorithm.py b/packages/best_fit/genetic_algorithm.py
@@ -122,6 +122,40 @@ def main(
             fundam_params, obs_clust, theor_tracks, R_V, ext_coefs,
             st_dist_mass, N_fc, cmpl_rnd, err_rnd, varIdxs, ranges, p_lst_e)
 
+        if flag_print_perc:
+
+            # For plotting purposes.
+            models_GA[N_popl * step:N_popl * (step + 1)] = generation
+            lkls_GA[N_popl * step:N_popl * (step + 1)] = lkl
+            lkl_best[step] = lkl[0]
+            # Discard large values associated with empty arrays from mean.
+            lkl_mean[step] = np.mean(
+                np.asarray(lkl)[np.asarray(lkl) < 9.9e08])
+
+            # Time used to check how fast the sampler is advancing.
+            elapsed_in = t.time() - start_in
+            percentage_complete = (100. * elapsed_in / available_secs)
+            if len(milestones) > 0 and percentage_complete >= milestones[0]:
+                m, s = divmod(max(1., available_secs - elapsed_in), 60)
+                h, m = divmod(m, 60)
+                # m += s / 60.
+                print("{:>3}% LP={:.1f} ({:.5f}, {:.3f}, {:.3f}, "
+                      "{:.2f}, {:.0f}, {:.2f})".format(
+                          milestones[0], lkl[0], *generation[0]) +
+                      " [{:.0f} m/s | {:.0f}h{:.0f}m]".format(
+                          (N_popl * step) / elapsed_in, h, m))
+                # Remove that milestone from the list.
+                milestones = milestones[1:]
+
+            # If the hardcoded maximum number of generations is reached.
+            if step == N_gener:
+                print(" WARNING: maximum allowed number of " +
+                      "generations ({}) reached.".format(N_gener))
+                break
+        else:
+            if step == N_gener:
+                break
+
         # *** Extinction/Immigration ***
         # If the best solution has remained unchanged for N_ei
         # generations, remove all chromosomes but the best ones (extinction)
@@ -144,7 +178,7 @@ def main(
                     ext_imm_count += 1
                     if ext_imm_count == N_es:
                         if flag_print_perc:
-                            print("  GA exit switch applied.")
+                            print("  GA exit switch applied")
                         elapsed_in = t.time() - start_in
                         break
                 else:
@@ -157,7 +191,6 @@ def main(
                 generation = ext_imm(best_sol, fundam_params, N_popl)
                 # Reset best solution counter.
                 best_sol_count = 0
-
         else:
             lkl_best_old = lkl[0]
             # For plotting purposes. Save index where a new best solution
@@ -168,40 +201,6 @@ def main(
             # Reset counter.
             best_sol_count = 0
 
-        if flag_print_perc:
-
-            # For plotting purposes.
-            models_GA[N_popl * step:N_popl * (step + 1)] = generation
-            lkls_GA[N_popl * step:N_popl * (step + 1)] = lkl
-            lkl_best[step] = lkl[0]
-            # Discard large values associated with empty arrays from mean.
-            lkl_mean[step] = np.mean(
-                np.asarray(lkl)[np.asarray(lkl) < 9.9e08])
-
-            # Time used to check how fast the sampler is advancing.
-            elapsed_in = t.time() - start_in
-            percentage_complete = (100. * elapsed_in / available_secs)
-            if len(milestones) > 0 and percentage_complete >= milestones[0]:
-                m, s = divmod(max(1., available_secs - elapsed_in), 60)
-                h, m = divmod(m, 60)
-                # m += s / 60.
-                print("{:>3}% LP={:.1f} ({:.5f}, {:.3f}, {:.3f}, "
-                      "{:.2f}, {:.0f}, {:.2f})".format(
-                          milestones[0], lkl[0], *generation[0]) +
-                      " [{:.0f} m/s | {:.0f}h{:.0f}m]".format(
-                          (N_popl * step) / elapsed_in, h, m))
-                # Remove that milestone from the list.
-                milestones = milestones[1:]
-
-            # If the hardcoded maximum number of generations is reached.
-            if step == N_gener:
-                print(" WARNING: maximum allowed number of " +
-                      "generations ({}) reached.".format(N_gener))
-                break
-
-        else:
-            if step == N_gener:
-                break
         step += 1
 
     # If this is a bootstrap run, return the best model found only.

diff --git a/packages/best_fit/likelihood.py b/packages/best_fit/likelihood.py
@@ -1,6 +1,6 @@
 
 import numpy as np
-from scipy.misc import logsumexp
+from scipy.special import logsumexp
 from scipy.stats import gaussian_kde, entropy
 
 # ############################################################

diff --git a/packages/best_fit/mcmc_convergence.py b/packages/best_fit/mcmc_convergence.py
@@ -301,13 +301,27 @@ def geweke(x, first=.1, last=.5, intervals=20):
 
 #     return mcmc_halves
 
-def convergenceVals(algor, ndim, varIdxs, N_conv, chains_nruns):
+def convergenceVals(algor, ndim, varIdxs, N_conv, chains_nruns, bi_steps):
     """
     Convergence statistics.
     """
     with warnings.catch_warnings():
         warnings.simplefilter("ignore")
 
+        # Mean Tau across chains, shape: (post-bi steps, ndims)
+        x = np.mean(chains_nruns.T, axis=1).T
+        tau_autocorr = []
+        for j in np.arange(50, x.shape[0], 50):
+            # tau.shape: ndim
+            tau = util.autocorr_integrated_time(x[:j])
+            # Autocorrelation time. Mean across dimensions.
+            tau_autocorr.append([bi_steps + j, np.mean(tau)])
+        # Add one last point with the entire chain.
+        if j < x.shape[0]:
+            tau = util.autocorr_integrated_time(x)
+            tau_autocorr.append([bi_steps + x.shape[0], np.mean(tau)])
+        tau_autocorr = np.array(tau_autocorr).T
+
         # Autocorrelation time for each parameter, mean across chains.
         if algor == 'emcee':
             acorr_t = autocorr.integrated_time(
@@ -359,10 +373,16 @@ def convergenceVals(algor, ndim, varIdxs, N_conv, chains_nruns):
         # Mean across chains
         geweke_z = np.nanmean(geweke_z, axis=1)
         acorr_function = np.nanmean(acorr_function, axis=1)
-        # Cut the autocorrelation function a bit after *all* the parameters
-        # have crossed the zero line.
-        lag_zero = max([np.where(_ < 0)[0][0] for _ in acorr_function])
-        acorr_function = acorr_function[:, :int(lag_zero + .2 * lag_zero)]
+
+        # # Cut the autocorrelation function just after *all* the parameters
+        # # have crossed the zero line.
+        # try:
+        #     lag_zero = max([np.where(_ < 0)[0][0] for _ in acorr_function])
+        # except IndexError:
+        #     # Could not obtain zero lag
+        #     lag_zero = acorr_function.shape[-1]
+        # acorr_function = acorr_function[:, :int(lag_zero + .2 * lag_zero)]
+
         # # Approx IAT
         # lag_iat = 1. + 2. * np.sum(acorr_function, axis=1)
         # print("  Approx (zero lag) IAT: ", lag_iat)
@@ -380,5 +400,5 @@ def convergenceVals(algor, ndim, varIdxs, N_conv, chains_nruns):
         #     mESS_epsilon[1].append(fminESS(ndim, alpha=alpha, ess=minESS))
         #     mESS_epsilon[2].append(fminESS(ndim, alpha=alpha, ess=mESS))
 
-    return acorr_t, med_at_c, all_taus, geweke_z, lag_zero, acorr_function,\
-        mcmc_ess
+    return tau_autocorr, acorr_t, med_at_c, all_taus, geweke_z,\
+        acorr_function, mcmc_ess
diff --git a/packages/best_fit/ptemcee/__init__.py b/packages/best_fit/ptemcee/__init__.py
@@ -4,8 +4,8 @@
 from __future__ import (division, print_function, absolute_import, unicode_literals)
 
 from .sampler import *
-from .interruptible_pool import InterruptiblePool
-from .mpi_pool import MPIPool
+# from .interruptible_pool import InterruptiblePool
+# from .mpi_pool import MPIPool
 from . import util
 
 __version__ = '1.0.0'