Issue/83/referee

README.md

@@ -30,6 +30,6 @@ As a disclaimer, it will need a major upgrade for flexibility and computational
 ## License, Contributing etc
 
 The code in this repo is available for re-use under the MIT license, which means that you can do whatever you like with it, just don't blame me.
-If you end up using any of the code or ideas you find here in your academic research, please cite me as `Malz et al, in preparation\footnote{\texttt{https://github.com/aimalz/chippr}}`.
+If you end up using any of the code or ideas you find here in your academic research, please cite me as `Malz & Hogg, 2020 \footnote{[arXiv:2007.12178](https://arxiv.org/abs/2007.12178)\texttt{https://github.com/aimalz/chippr}}`.
 If you are interested in this project, please do drop me a line via the hyperlinked contact name above, or by [writing me an issue](https://github.com/aimalz/chippr/issues/new).
 To get started contributing to the `chippr` project, just fork the repo -- pull requests are always welcome!

chippr/__init__.py

@@ -1,17 +1,17 @@
-from defaults import *
+from .defaults import *
 
-from utils import *
-from stat_utils import *
-from plot_utils import *
+from .utils import *
+from .stat_utils import *
+from .plot_utils import *
 
-from sim_utils import *
-from discrete import *
-from gauss import *
-from gamma import *
-from gmix import *
-from mvn import *
-from multi_dist import *
-from catalog import *
+from .sim_utils import *
+from .discrete import *
+from .gauss import *
+from .gamma import *
+from .gmix import *
+from .mvn import *
+from .multi_dist import *
+from .catalog import *
 
-from log_z_dens import *
-from log_z_dens_plots import *
+from .log_z_dens import *
+from .log_z_dens_plots import *

chippr/catalog.py

@@ -5,7 +5,7 @@
 import pickle as pkl
 
 import matplotlib as mpl
-mpl.use('PS')
+# mpl.use('PS')
 import matplotlib.pyplot as plt
 
 import chippr
@@ -45,7 +45,7 @@ def __init__(self, params={}, vb=True, loc='.', prepend=''):
         self.params = d.check_sim_params(self.params)
 
         if vb:
-            print self.params
+            print(self.params)
 
         np.random.seed(d.seed)
         self.cat = {}
@@ -427,7 +427,7 @@ def evaluate_lfs(self, pspace,  vb=True):
         lfs = []
         for n in self.N_range:
             points = zip(self.z_fine, [self.samps[n][1]] * self.n_tot)
-            cur=pspace.pdf(np.array(points))
+            cur = pspace.pdf(np.array([p for p in points]))
             lfs.append(cur)
         lfs = np.array(lfs)
         lfs /= np.sum(lfs, axis=-1)[:, np.newaxis] * self.dz_fine

chippr/catalog_plots.py

@@ -2,7 +2,7 @@
 import os
 
 import matplotlib as mpl
-mpl.use('PS')
+# mpl.use('PS')
 import matplotlib.pyplot as plt
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 
@@ -31,8 +31,8 @@ def plot_true_histogram(true_samps, n_bins=(10, 50), plot_loc='', prepend='', pl
     pu.set_up_plot()
     f = plt.figure(figsize=(5, 5))
     sps = f.add_subplot(1, 1, 1)
-    sps.hist(true_samps, bins=n_bins[1], normed=1, color='k', alpha=0.5, log=True)
-    sps.hist(true_samps, bins=n_bins[0], normed=1, color='y', alpha=0.5, log=True)
+    sps.hist(true_samps, bins=n_bins[1], color='k', alpha=0.5, log=True)#, normed=1
+    sps.hist(true_samps, bins=n_bins[0], color='y', alpha=0.5, log=True)#, normed=1
     sps.set_xlabel(r'$z_{true}$')
     sps.set_ylabel(r'$n(z_{true})$')
     f.savefig(os.path.join(plot_loc, prepend+plot_name), bbox_inches='tight', pad_inches = 0, dpi=d.dpi)
@@ -105,6 +105,9 @@ def plot_mega_scatter(zs, pfs, z_grid, grid_ends, truth=None, plot_loc='', prepe
     int_pr: numpy.ndarray, float, optional
         plot the interim prior with the histograms?
     """
+    print('debug demo')
+    print((truth, np.shape(zs), np.shape(pfs), np.shape(z_grid), np.shape(grid_ends)))
+
     n = len(zs)
     zs = zs.T
     true_zs = zs[0]
@@ -154,10 +157,10 @@ def plot_mega_scatter(zs, pfs, z_grid, grid_ends, truth=None, plot_loc='', prepe
     histy.hist(obs_zs, bins=grid_ends, orientation='horizontal', alpha=0.25, color='k', stacked=False, density=True)
     if truth is not None:
         histx.step(truth[0], truth[1], c='k', where='mid')
-        pu.plot_step(histy, np.pad(truth[1], (1, 1), 'constant', constant_values=(0, 0)), grid_ends, c='k', w=1.5)
+        pu.plot_step(histy, np.pad(truth[1], (1, 1), 'constant', constant_values=(0, 0)), u.ends_from_mids(truth[0]), c='k', w=1.5)
     if int_pr is not None:
         histx.step(int_pr[0], int_pr[1], c='k', alpha=0.5, where='mid')
-        pu.plot_step(histy, np.pad(int_pr[1], (1, 1), 'constant', constant_values=(0, 0)), grid_ends, c='k', a=0.5, w=1.5)
+        pu.plot_step(histy, np.pad(int_pr[1], (1, 1), 'constant', constant_values=(0, 0)), u.ends_from_mids(int_pr[0]), c='k', a=0.5, w=1.5)
     histx.xaxis.set_tick_params(labelbottom=False)
     histy.yaxis.set_tick_params(labelleft=False)
     histx.set_yticks([])

chippr/discrete.py

@@ -30,7 +30,7 @@ def __init__(self, bin_ends, weights):
         self.bin_ends = bin_ends
         self.dbins = self.bin_ends[1:] - self.bin_ends[:-1]
         self.n_bins = len(self.bin_ends)-1
-        self.bin_range = range(self.n_bins)
+        self.bin_range = list(range(self.n_bins))
 
         self.weights = weights
         # print('dbins: '+str(self.dbins))

chippr/log_z_dens.py

@@ -2,11 +2,11 @@
 import scipy as sp
 import os
 import scipy.optimize as op
-import cPickle as cpkl
+import pickle as cpkl
 import emcee
 
 import matplotlib as mpl
-mpl.use('PS')
+# mpl.use('PS')
 import matplotlib.pyplot as plt
 
 import chippr
@@ -61,7 +61,7 @@ def __init__(self, catalog, hyperprior, truth=None, loc='.', prepend='', vb=True
         self.info['log_interim_posteriors'] = self.log_pdfs
 
         if vb:
-            print(str(self.n_bins) + ' bins, ' + str(len(self.log_pdfs)) + ' interim posterior PDFs')
+            print((str(self.n_bins) + ' bins, ' + str(len(self.log_pdfs)) + ' interim posterior PDFs'))
 
         self.hyper_prior = hyperprior
 
@@ -214,12 +214,12 @@ def _objective(log_nz):
                 return -2. * self.evaluate_log_hyper_posterior(log_nz)
 
         if vb:
-            print(self.dir + ' starting at ', start, _objective(start))
+            print((self.dir + ' starting at ', start, _objective(start)))
 
         res = op.minimize(_objective, start, method="Nelder-Mead", options={"maxfev": 1e5, "maxiter":1e5})
 
         if vb:
-            print(self.dir + ': ' + str(res))
+            print((self.dir + ': ' + str(res)))
         return res.x
 
     def calculate_mmle(self, start, vb=True, no_data=0, no_prior=0):
@@ -428,7 +428,7 @@ def distribution(log_nz):
             full_chain = np.array([[vals[w]] for w in range(self.n_walkers)])
             while self.burning_in:
                 if vb:
-                    print('beginning sampling '+str(self.burn_ins))
+                    print(('beginning sampling '+str(self.burn_ins)))
                 burn_in_mcmc_outputs = self.sample(vals, 10**n_burned)
                 chain = burn_in_mcmc_outputs['chains']
                 burn_in_mcmc_outputs['chains'] -= u.safe_log(np.sum(np.exp(chain) * self.bin_difs[np.newaxis, np.newaxis, :], axis=2))[:, :, np.newaxis]
@@ -489,10 +489,10 @@ def compare(self, vb=True):
                 self.info['stats']['rms']['true_nz' + '__' + key[4:]] = s.calculate_rms(np.exp(self.info['log_tru_nz']), np.exp(self.info['estimators'][key]))
                 self.info['stats']['log_rms']['log_true_nz'+ '__' + key] = s.calculate_rms(self.info['log_tru_nz'], self.info['estimators'][key])
 
-        for i in range(len(self.info['estimators'].keys())):
-            key_1 = self.info['estimators'].keys()[i]
-            for j in range(len(self.info['estimators'].keys()[:i])):
-                key_2 = self.info['estimators'].keys()[j]
+        for i in range(len(list(self.info['estimators'].keys()))):
+            key_1 = list(self.info['estimators'].keys())[i]
+            for j in range(len(list(self.info['estimators'].keys())[:i])):
+                key_2 = list(self.info['estimators'].keys())[j]
                 # print(((i,j), (key_1, key_2)))
                 self.info['stats']['log_rms'][key_1 + '__' + key_2] = s.calculate_rms(self.info['estimators'][key_1], self.info['estimators'][key_2])
                 self.info['stats']['rms'][key_1[4:] + '__' + key_2[4:]] = s.calculate_rms(np.exp(self.info['estimators'][key_1]), np.exp(self.info['estimators'][key_2]))
@@ -537,11 +537,11 @@ def read(self, read_loc, style='pickle', vb=True):
         with open(os.path.join(self.res_dir, read_loc), 'rb') as file_location:
             self.info = cpkl.load(file_location)
         if vb:
-            print('The following quantities were read from '+read_loc+' in the '+style+' format:')
+            print(('The following quantities were read from '+read_loc+' in the '+style+' format:'))
             for key in self.info:
                 print(key)
             if 'estimators' in self.info:
-                print(self.info['estimators'].keys())
+                print((list(self.info['estimators'].keys())))
         return self.info
 
     def write(self, write_loc, style='pickle', vb=True):
@@ -560,7 +560,7 @@ def write(self, write_loc, style='pickle', vb=True):
         with open(os.path.join(self.res_dir, write_loc), 'wb') as file_location:
             cpkl.dump(self.info, file_location)
         if vb:
-            print('The following quantities were written to '+write_loc+' in the '+style+' format:')
+            print(('The following quantities were written to '+write_loc+' in the '+style+' format:'))
             for key in self.info:
                 print(key)
         return

chippr/log_z_dens_plots.py

@@ -3,7 +3,7 @@
 import scipy as sp
 
 import matplotlib as mpl
-mpl.use('PS')
+# mpl.use('PS')
 import matplotlib.pyplot as plt
 from matplotlib import gridspec
 

chippr/mvn.py

@@ -29,7 +29,7 @@ def __init__(self, mean, var):
         self.sigma = self.norm_var()
         self.invvar = self.invert_var()
 
-        assert np.linalg.eig(self.var) > 0.
+        assert(np.all(np.linalg.eig(self.var)[0] > 0.))
 
         self.dist = MGD(self.mean, self.var)
 

chippr/plot_utils.py

@@ -1,7 +1,7 @@
 import numpy as np
 
 import matplotlib as mpl
-mpl.use('PS')
+# mpl.use('PS')
 import matplotlib.pyplot as plt
 import matplotlib.cm as cm
 

chippr/stat_utils.py

@@ -151,7 +151,7 @@ def multi_parameter_gr_stat(sample):
     """
     dims = np.shape(sample)
     (n_walkers, n_iterations, n_params) = dims
-    n_burn_ins = n_iterations / 2
+    n_burn_ins = int(n_iterations / 2)
     chain_ensemble = np.swapaxes(sample, 0, 1)
     chain_ensemble = chain_ensemble[n_burn_ins:, :]
     Rs = np.zeros((n_params))
@@ -177,7 +177,7 @@ def gr_test(sample, threshold=d.gr_threshold):
         True if burning in, False if post-burn in
     """
     gr = multi_parameter_gr_stat(sample)
-    print('Gelman-Rubin test statistic = '+str(gr))
+    print(('Gelman-Rubin test statistic = '+str(gr)))
     test_result = np.max(gr) > threshold
     return test_result
 
@@ -198,7 +198,7 @@ def cft(xtimes, lag):#xtimes has ntimes elements
         autocorrelation time for one time lag for one parameter of one walker
     """
     lent = len(xtimes) - lag
-    allt = xrange(lent)
+    allt = range(lent)
     ans = np.array([xtimes[t+lag] * xtimes[t] for t in allt])
     return ans
 
@@ -217,8 +217,8 @@ def cf(xtimes):#xtimes has ntimes elements
         autocorrelation time over all time lags for one parameter of one walker
     """
     cf0 = np.dot(xtimes, xtimes)
-    allt = xrange(len(xtimes) / 2)
-    cf = np.array([sum(cft(xtimes,lag)[len(xtimes) / 2:]) for lag in allt]) / cf0
+    allt = range(int(len(xtimes) / 2))
+    cf = np.array([sum(cft(xtimes, lag)[int(len(xtimes) / 2):]) for lag in allt]) / cf0
     return cf
 
 def cfs(x, mode):#xbinstimes has nbins by ntimes elements

chippr/utils.py

@@ -56,3 +56,39 @@ def ingest(in_info):
     else:
         in_dict = in_info
     return in_dict
+
+def mids_from_ends(inarr):
+    """
+    Function to make midpoints from grid of endpoints
+
+    Parameters
+    ----------
+    inarr: array, shape=(N)
+        array of grid endpoints
+
+    Returns
+    -------
+    outarr: array, shape=(N-1)
+        array of corresponding midpoints
+    """
+    outarr = (inarr[1:] + inarr[:-1]) / 2.
+    return outarr
+
+def ends_from_mids(inarr):
+    """
+    Function to make endpoints from grid of midpoints, assuming equal spacing
+
+    Parameters
+    ----------
+    inarr: array, shape=(N)
+        array of grid midpoints
+
+    Returns
+    -------
+    outarr: array, shape=(N+1)
+        array of corresponding endpoints
+    """
+    dif = np.mean(inarr[1:] - inarr[:-1])
+    int_ends = mids_from_ends(inarr)
+    outarr = np.concatenate((np.array([int_ends[0] - dif]), int_ends, np.array([int_ends[-1] + dif])))
+    return outarr

docs/notebooks/demo2.ipynb

@@ -8,7 +8,9 @@
     "\n",
     "This notebook demonstrates the use of the Cosmological Hierarchical Inference with Probabilistic Photometric Redshifts (CHIPPR) package to estimate population distributions based on a catalog of probability distributions.\n",
     "\n",
-    "The package supports two primary objectives: simulation of catalogs and inference of posterior distributions over parameters defining population distributions."
+    "The package supports two primary objectives: simulation of catalogs and inference of posterior distributions over parameters defining population distributions.\n",
+    "\n",
+    "*Note that this notebook dates back to the [Python 2 tagged release](https://github.com/aimalz/chippr/releases/tag/v0.1-beta) and may not work under the latest version.*"
    ]
   },
   {
@@ -451,14 +453,14 @@
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2",
-   "version": "2.7.15rc1"
+   "pygments_lexer": "ipython3",
+   "version": "3.8.5"
   }
  },
  "nbformat": 4,