From cf1976c1ff46219d061f5743afd70067b7d4b85a Mon Sep 17 00:00:00 2001 From: Johannes Buchner Date: Sun, 20 Jun 2021 21:49:53 +0200 Subject: [PATCH] Improve doc strings; remove UniformOrderAccumulator arg --- tests/test_ordertest.py | 6 ++-- ultranest/hotstart.py | 64 ++++++++++++++++++++++++++++++------- ultranest/integrator.py | 66 ++++++++++++++++++++------------------- ultranest/netiter.py | 60 +++++++++++++++++++++++++++-------- ultranest/ordertest.py | 18 ++++++----- ultranest/samplingpath.py | 5 +-- ultranest/store.py | 3 +- ultranest/utils.py | 23 ++++++-------- 8 files changed, 163 insertions(+), 82 deletions(-) diff --git a/tests/test_ordertest.py b/tests/test_ordertest.py index 8856c6e7..9301d345 100644 --- a/tests/test_ordertest.py +++ b/tests/test_ordertest.py @@ -3,7 +3,7 @@ from ultranest.ordertest import UniformOrderAccumulator, infinite_U_zscore def test_invalid_order(): - sample_acc = UniformOrderAccumulator(3) + sample_acc = UniformOrderAccumulator() sample_acc.add(2, 3) try: sample_acc.add(4, 3) @@ -12,7 +12,7 @@ def test_invalid_order(): pass def test_diff_expand(): - sample_acc = UniformOrderAccumulator(3) + sample_acc = UniformOrderAccumulator() sample_acc.add(1, 3) sample_acc.add(4, 5) sample_acc.add(5, 6) @@ -24,7 +24,7 @@ def test_order_correctness(): nruns = [] for frac in 1, 0.9: print("frac:", frac) - sample_acc = UniformOrderAccumulator(Nlive) + sample_acc = UniformOrderAccumulator() runlength = [] samples = [] for i in range(N): diff --git a/ultranest/hotstart.py b/ultranest/hotstart.py index 6b76875b..fcdd582d 100644 --- a/ultranest/hotstart.py +++ b/ultranest/hotstart.py @@ -4,11 +4,12 @@ import scipy.stats from .utils import vectorize, resample_equal + def get_auxiliary_problem(loglike, transform, ctr, invcov, enlargement_factor, df=1): """Return a new loglike and transform based on an auxiliary distribution. Given a likelihood and prior transform, and information about - the (expected) posterior peak, generates a auxiliary + the (expected) posterior peak, generates a auxiliary likelihood and prior transform that is identical but requires fewer nested sampling iterations. @@ -56,14 +57,13 @@ def get_auxiliary_problem(loglike, transform, ctr, invcov, enlargement_factor, d The first d return coordinates are identical to what ``transform`` would return. The final coordinate is the correction weight. """ - ndim, = ctr.shape assert invcov.shape == (ndim, ndim) assert df >= 1, ('Degrees of freedom must be above 1', df) l, v = np.linalg.eigh(invcov) rotation_matrix = np.dot(v, enlargement_factor * np.diag(1. / np.sqrt(l))) - + rv_auxiliary1d = scipy.stats.t(df) def aux_rotator(coords): @@ -78,19 +78,20 @@ def aux_loglikelihood(u): if not (x > 0).all() or not (x < 1).all(): return -1e300 # undo the effect of the auxiliary distribution - l = rv_auxiliary1d.logpdf(coords).sum() - return loglike(transform(x)) - l - + loglike = rv_auxiliary1d.logpdf(coords).sum() + return loglike(transform(x)) - loglike + def aux_aftertransform(u): return transform(aux_rotator(rv_auxiliary1d.ppf(u))) return aux_loglikelihood, aux_aftertransform + def get_extended_auxiliary_problem(loglike, transform, ctr, invcov, enlargement_factor, df=1): """Return a new loglike and transform based on an auxiliary distribution. Given a likelihood and prior transform, and information about - the (expected) posterior peak, generates a auxiliary + the (expected) posterior peak, generates a auxiliary likelihood and prior transform that is identical but requires fewer nested sampling iterations. @@ -135,13 +136,13 @@ def get_extended_auxiliary_problem(loglike, transform, ctr, invcov, enlargement_ ndim, = ctr.shape assert invcov.shape == (ndim, ndim) assert df >= 1, ('Degrees of freedom must be above 1', df) - + l, v = np.linalg.eigh(invcov) rotation_matrix = np.dot(v, enlargement_factor * np.diag(1. / np.sqrt(l))) rv_auxiliary1d = scipy.stats.t(df) weight_ref = rv_auxiliary1d.logpdf(0) * ndim - + def aux_transform(u): # get uniform gauss/t distributed values: coords = rv_auxiliary1d.ppf(u) @@ -166,12 +167,53 @@ def aux_loglikelihood(x): return aux_loglikelihood, aux_transform + def reuse_samples( - param_names, loglike, points, logl, logw=None, + param_names, loglike, points, logl, logw=None, logz=0.0, logzerr=0.0, upoints=None, batchsize=128, vectorized=False, log_weight_threshold=-10, **kwargs ): + """ + Reweight existing nested sampling run onto a new loglikelihood. + + Parameters + ------------ + param_names: list of strings + Names of the parameters + loglike: function + New likelihood function + points: np.array of shape (npoints, ndim) + Equally weighted (unless logw is passed) posterior points + logl: np.array(npoints) + Previously likelihood values of points + logw: np.array(npoints) + Log-weights of existing points. + logz: float + Previous evidence / marginal likelihood value. + logzerr: float + Previous evidence / marginal likelihood uncertainty. + upoints: np.array of shape (npoints, ndim) + Posterior points before transformation. + vectorized: bool + Whether loglike function is vectorized + batchsize: int + Number of points simultaneously passed to vectorized loglike function + log_weight_threshold: float + Lowest log-weight to consider + + Returns: + --------- + results: dict + All information of the run. Important keys: + Number of nested sampling iterations (niter), + Evidence estimate (logz), + Effective Sample Size (ess), + weighted samples (weighted_samples), + equally weighted samples (samples), + best-fit point information (maximum_likelihood), + posterior summaries (posterior). + """ if not vectorized: loglike = vectorize(loglike) @@ -188,7 +230,7 @@ def reuse_samples( indices = np.argsort(logl + logw)[::-1] ncall = 0 for i in range(int(np.ceil(Npoints / batchsize))): - batch = indices[i * batchsize : (i + 1) * batchsize] + batch = indices[i * batchsize:(i + 1) * batchsize] logl_new[batch] = loglike(points[batch,:]) logw_new[batch] = logw[batch] + logl_new[batch] ncall += len(batch) diff --git a/ultranest/integrator.py b/ultranest/integrator.py index 0f3f5b45..59c55fb8 100644 --- a/ultranest/integrator.py +++ b/ultranest/integrator.py @@ -16,7 +16,8 @@ from numpy import log, exp, logaddexp import numpy as np -from .utils import create_logger, make_run_dir, resample_equal, vol_prefactor, vectorize, listify as _listify, is_affine_transform, normalised_kendall_tau_distance +from .utils import create_logger, make_run_dir, resample_equal, vol_prefactor, vectorize, listify as _listify +from .utils import is_affine_transform, normalised_kendall_tau_distance from ultranest.mlfriends import MLFriends, AffineLayer, ScalingLayer, find_nearby, WrappingEllipsoid, RobustEllipsoidRegion from .store import HDF5PointStore, TextPointStore, NullPointStore from .viz import get_default_viz_callback, nicelogger @@ -127,9 +128,10 @@ def _explore_iterator_batch(explorer, pop, x_dim, num_params, pointpile, batchsi yield batch - -def resume_from_similar_file(log_dir, x_dim, loglikelihood, transform, - max_tau=0, verbose=False, ndraw=400): +def resume_from_similar_file( + log_dir, x_dim, loglikelihood, transform, + max_tau=0, verbose=False, ndraw=400 +): """ Change a stored UltraNest run to a modified loglikelihood/transform. @@ -150,7 +152,7 @@ def resume_from_similar_file(log_dir, x_dim, loglikelihood, transform, max_tau: float Allowed dissimilarity in the live point ordering, quantified as normalised Kendall tau distance. - + max_tau=0 is the very conservative choice of stopping the warm start when the live point order differs. Near 1 are completely different live point orderings. @@ -191,7 +193,7 @@ def resume_from_similar_file(log_dir, x_dim, loglikelihood, transform, pointpile2 = PointPile(x_dim, num_params) def pop(Lmin): - """ find matching sample from points file """ + """Find matching sample from points file.""" # look forward to see if there is an exact match # if we do not use the exact matches # this causes a shift in the loglikelihoods @@ -273,8 +275,7 @@ def pop(Lmin): print("stopping, number of live points differ (%d vs %d)" % (len(active_values), len(active_values2))) good_state = False break - - + if len(active_values) != len(indices1): indices1, indices2 = np.meshgrid(np.arange(len(active_values)), np.arange(len(active_values2))) tau = normalised_kendall_tau_distance(active_values, active_values2, indices1, indices2) @@ -289,8 +290,8 @@ def pop(Lmin): if verbose == 2: print(niter, tau) if good_state: - #print(" (%.1e) L=%.1f" % (last_good_like, Lmin2)) - #assert last_good_like < Lmin2, (last_good_like, Lmin2) + # print(" (%.1e) L=%.1f" % (last_good_like, Lmin2)) + # assert last_good_like < Lmin2, (last_good_like, Lmin2) last_good_like = Lmin2 last_good_state = niter else: @@ -306,7 +307,7 @@ def pop(Lmin): logl_new = logls_new[j] j += 1 - #print(j, Lmin2, '->', logl_new, 'instead of', Lmin, '->', [c.value for c in node2.children]) + # print(j, Lmin2, '->', logl_new, 'instead of', Lmin, '->', [c.value for c in node2.children]) child2 = pointpile2.make_node(logl_new, u, v) node2.children.append(child2) if logl_new > Lmin2: @@ -314,8 +315,8 @@ def pop(Lmin): else: if verbose == 2: print("cannot use new point because it would decrease likelihood (%.1f->%.1f)" % (Lmin2, logl_new)) - #good_state = False - #break + # good_state = False + # break main_iterator2.passing_node(node2, active_nodes2) @@ -352,7 +353,8 @@ def pop(Lmin): def _update_region_bootstrap(region, nbootstraps, minvol=0., comm=None, mpi_size=1): """ - update *region* with *nbootstraps* rounds of excluding points randomly. + Update *region* with *nbootstraps* rounds of excluding points randomly. + Stiffen ellipsoid size using the minimum volume *minvol*. If the mpi communicator *comm* is not None, use MPI to distribute @@ -379,7 +381,7 @@ def _update_region_bootstrap(region, nbootstraps, minvol=0., comm=None, mpi_size recv_enlarge = comm.gather(f, root=0) recv_enlarge = comm.bcast(recv_enlarge, root=0) f = np.max(recv_enlarge[:nbootstraps]) - + if not np.isfinite(r) and not np.isfinite(r): # reraise error if needed if e is None: @@ -472,7 +474,7 @@ def __init__(self, assert np.isfinite(logl).all(), ("Error in loglikelihood function: returned non-finite number: %s for input u=%s p=%s" % (logl, u, p)) def safe_loglike(x): - """ wrapped likelihood function """ + """Call likelihood function safely wrapped to avoid non-finite values.""" x = np.asarray(x) logl = loglike(x) assert np.isfinite(logl).all(), ( @@ -1004,7 +1006,7 @@ def __init__(self, warmstart_max_tau: float Maximum disorder to accept when resume='resume-similar'; - Live points are reused as long as the live point order + Live points are reused as long as the live point order is below this normalised Kendall tau distance. Values from 0 (highly conservative) to 1 (extremely negligent). """ @@ -1129,7 +1131,8 @@ def _setup_distributed_seeds(self): np.random.seed(seed) def _check_likelihood_function(self, transform, loglike, num_test_samples): - """Tests the `transform` and `loglike`lihood functions. + """Test the `transform` and `loglike`lihood functions. + `num_test_samples` samples are used to check whether they work and give the correct output. returns whether the most recently stored point (if any) @@ -1195,7 +1198,7 @@ def _set_likelihood_function(self, transform, loglike, num_test_samples, make_sa """ def safe_loglike(x): - """ safe wrapper of likelihood function """ + """Safe wrapper of likelihood function.""" x = np.asarray(x) if len(x.shape) == 1: assert x.shape[0] == self.x_dim @@ -1215,7 +1218,7 @@ def safe_loglike(x): self.transform = lambda x: x elif make_safe: def safe_transform(x): - """ safe wrapper of transform function """ + """Safe wrapper of transform function.""" x = np.asarray(x) if len(x.shape) == 1: assert x.shape[0] == self.x_dim @@ -1672,7 +1675,7 @@ def _create_point(self, Lmin, ndraw, active_u, active_values): else: u, v, logl, nc, quality = self._refill_samples(Lmin, ndraw, nit) nit += 1 - + if logl is None: u = np.empty((0, self.x_dim)) v = np.empty((0, self.num_params)) @@ -1724,8 +1727,8 @@ def _update_region( bootstrap_rootids=None, active_rootids=None, nbootstraps=30, minvol=0., active_p=None ): - """Build a new MLFriends region from `active_u`, - and wrapping ellipsoid. + """Build a new MLFriends region from `active_u`, and wrapping ellipsoid. + Both are safely built using bootstrapping, so that the region can be used for sampling and rejecting points. If MPI is enabled, this computation is parallelised. @@ -2058,7 +2061,7 @@ def run( min_num_live_points=400, cluster_num_live_points=40, insertion_test_window=10, - insertion_test_zscore_threshold=2, + insertion_test_zscore_threshold=4, region_class=MLFriends, ): """Run until target convergence criteria are fulfilled. @@ -2125,14 +2128,13 @@ def run( z-score used as a threshold for the insertion order test. Set to infinity to disable. - insertion_test_window: float + insertion_test_window: int Number of iterations after which the insertion order test is reset. region_class: MLFriends or RobustEllipsoidRegion Whether to use MLFriends+ellipsoidal+tellipsoidal region (better for multi-modal problems) or just ellipsoidal sampling (faster for high-dimensional, gaussian-like problems). """ - for result in self.run_iter( update_interval_volume_fraction=update_interval_volume_fraction, update_interval_ncall=update_interval_ncall, @@ -2174,7 +2176,7 @@ def run_iter( cluster_num_live_points=40, show_status=True, viz_callback='auto', - insertion_test_window=10, + insertion_test_window=10000, insertion_test_zscore_threshold=2, region_class=MLFriends ): @@ -2262,7 +2264,7 @@ def run_iter( nbootstraps=max(1, self.num_bootstraps // self.mpi_size), random=False, check_insertion_order=False) main_iterator.Lmax = max(Lmax, max(n.value for n in roots)) - insertion_test = UniformOrderAccumulator(nroots) + insertion_test = UniformOrderAccumulator() insertion_test_runs = [] insertion_test_quality = np.inf insertion_test_direction = 0 @@ -2394,7 +2396,7 @@ def run_iter( insertion_test_quality = insertion_test.N insertion_test_direction = np.sign(insertion_test.zscore) insertion_test.reset() - elif insertion_test.N > nlive * insertion_test_window: + elif insertion_test.N > insertion_test_window: insertion_test_quality = np.inf insertion_test_direction = 0 insertion_test.reset() @@ -2787,7 +2789,7 @@ def read_file(log_dir, x_dim, num_bootstraps=20, random=True, verbose=False, che pointpile = PointPile(x_dim, num_params) def pop(Lmin): - """ find matching sample from points file """ + """Find matching sample from points file.""" # look forward to see if there is an exact match # if we do not use the exact matches # this causes a shift in the loglikelihoods @@ -2812,11 +2814,11 @@ def pop(Lmin): root = TreeNode(id=-1, value=-np.inf, children=roots) def onNode(node, main_iterator): - """ insert (single) child of node if available """ + """Insert (single) child of node if available.""" while True: _, row = pop(node.value) if row is None: - break + break if row is not None: logl = row[1] u = row[3:3 + x_dim] diff --git a/ultranest/netiter.py b/ultranest/netiter.py index d3785ef1..035c9cae 100644 --- a/ultranest/netiter.py +++ b/ultranest/netiter.py @@ -23,7 +23,6 @@ import numpy as np from numpy import log, log1p, exp, logaddexp import math -import operator import sys from .utils import resample_equal from .ordertest import UniformOrderAccumulator @@ -121,8 +120,7 @@ def drop_next_node(self): assert len(self.active_nodes) == len(self.active_node_values) def expand_children_of(self, rootid, node): - """Replace the current node with its children in the iterators - list of active nodes. + """Replace the current node with its children in the iterators list of active nodes. Parameters ---------- @@ -164,7 +162,7 @@ def expand_children_of(self, rootid, node): def _stringify_lanes(lanes, char='║'): - """ unicode-draw lanes, fill with vertical stripes or spaces """ + """unicode-draw lanes, fill with vertical stripes or spaces.""" return ''.join([' ' if n is None else char for n in lanes]) @@ -286,7 +284,9 @@ def count_tree(roots): def count_tree_between(roots, lo, hi): - """Return the total number of nodes and maximum number of parallel edges, + """Compute basic statistics about a tree. + + Return the total number of nodes and maximum number of parallel edges, but only considering a interval of the tree. Parameters @@ -621,7 +621,7 @@ def __init__(self, nroots, nbootstraps=10, random=False, check_insertion_order=F self.check_insertion_order = check_insertion_order self.insertion_order_threshold = 4 - self.insertion_order_accumulator = UniformOrderAccumulator(self.rootids.shape[1]) + self.insertion_order_accumulator = UniformOrderAccumulator() self.reset(len(self.rootids)) @@ -658,6 +658,22 @@ def logZerr_bs(self): @property def insertion_order_runlength(self): + """Get shortest insertion order test run. + + Returns + -------- + shortest_run_length: int + Shortest insertion order test run length. + + The MWW (U-test) statistic is considered at each iteration. + When it exceeds a threshold (4 sigma by default, `insertion_order_threshold`), + the statistic is reset. The run length is recorded. + This property returns the shortest run length of all recorded + so far, or infinity otherwise. + + At 4 sigma, run lengths no shorter than 10^5.5 are expected + in unbiased runs. + """ runs = self.insertion_order_runs if len(runs) == 0: return np.inf @@ -666,6 +682,27 @@ def insertion_order_runlength(self): @property def insertion_order_converged(self): + """Check convergence. + + Returns + -------- + converged: bool + Whether the run is unbiased according to a U-test. + + The MWW (U-test) statistic is considered at each iteration. + When it exceeds a threshold (4 sigma by default, `insertion_order_threshold`), + the statistic is reset. The run length is recorded. + This property returns the shortest run length of all recorded + so far, or infinity otherwise. + + At 4 sigma, run lengths no shorter than 10^5.5 are expected + in unbiased runs. If the number of runs exceeds the number + of iterations divided by 10^5.5, the run is likely unbiased + and thus not converged. + + If not converged, the step sampler may need to use more steps, + or the problem needs to be reparametrized. + """ # we expect run lengths not shorter than 300000 for 4sigma # if we get many more than expected from the number of iterations # the run has not converged @@ -765,7 +802,6 @@ def passing_node(self, rootid, node, rootids, parallel_values): self.logVolremaining = self.all_logVolremaining[0] if self.check_insertion_order and len(np.unique(parallel_values)) == len(parallel_values): - order_max = nlive.max() + 1 acc = self.insertion_order_accumulator parallel_values_here = parallel_values[self.rootids[0, rootids]] for child in node.children: @@ -826,14 +862,14 @@ def combine_results(saved_logl, saved_nodeids, pointpile, main_iterator, mpi_com iterator used mpi_comm: MPI communicator object, or None if MPI is not used. - + Returns -------- results: dict All information of the run. Important keys: - Number of nested sampling iterations (niter), - Evidence estimate (logz), - Effective Sample Size (ess), + Number of nested sampling iterations (niter), + Evidence estimate (logz), + Effective Sample Size (ess), H (information gain), weighted samples (weighted_samples), equally weighted samples (samples), @@ -842,7 +878,6 @@ def combine_results(saved_logl, saved_nodeids, pointpile, main_iterator, mpi_com The rank order test score (insertion_order_MWW_test) is included if the iterator has it. """ - assert np.shape(main_iterator.logweights) == (len(saved_logl), len(main_iterator.all_logZ)), ( np.shape(main_iterator.logweights), np.shape(saved_logl), @@ -1005,7 +1040,6 @@ def logz_sequence(root, pointpile, nbootstraps=12, random=True, onNode=None, ver # first time they are all the same logzerr.append(main_iterator.logZerr_bs) - insert_order_value = 0.0 if len(np.unique(active_values)) == len(active_values) and len(node.children) > 0: child_insertion_order = (active_values > node.children[0].value).sum() insert_order.append(2 * (child_insertion_order + 1.) / len(active_values)) diff --git a/ultranest/ordertest.py b/ultranest/ordertest.py index 9a189ed6..06df14de 100644 --- a/ultranest/ordertest.py +++ b/ultranest/ordertest.py @@ -1,6 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- """ +Mann-Whitney-Wilcoxon U test for a uniform distribution of integers. This implements the same idea as https://arxiv.org/abs/2006.03371 except their KS test is problematic because the variable (insertion order) @@ -25,8 +26,7 @@ def infinite_U_zscore(sample, B): """ - Compute Mann-Whitney-Wilcoxon U test for a - *sample* of integers to be uniformly distributed between 0 and *B*. + Compute Mann-Whitney-Wilcoxon U test for a *sample* of integers to be uniformly distributed between 0 and *B*. Parameters ---------- @@ -44,15 +44,18 @@ def infinite_U_zscore(sample, B): class UniformOrderAccumulator(): + """Mann-Whitney-Wilcoxon U test accumulator. + + Stores rank orders (1 to N), for comparison with a uniform order. """ - Store orders (1 to N), for comparison with a uniform order. - """ - def __init__(self, N): + + def __init__(self): + """Initiate empty accumulator.""" self.N = 0 self.U = 0.0 def reset(self): - """Set all counts to zero. """ + """Set all counts to zero.""" self.N = 0 self.U = 0.0 @@ -74,7 +77,7 @@ def add(self, order, N): @property def zscore(self): - """ Mann-Whitney-Wilcoxon U test z-score, against a uniform distribution. """ + """z-score of the null hypothesis (uniform distribution) probability.""" N = self.N if N == 0: return 0.0 @@ -83,4 +86,5 @@ def zscore(self): return (self.U - m_U) / sigma_U_corr def __len__(self): + """Return number of samples accumulated so far.""" return self.N diff --git a/ultranest/samplingpath.py b/ultranest/samplingpath.py index c8cf469b..0a800a44 100644 --- a/ultranest/samplingpath.py +++ b/ultranest/samplingpath.py @@ -10,7 +10,7 @@ def nearest_box_intersection_line(ray_origin, ray_direction, fwd=True): - """Compute intersection of a line (ray) and a unit box (0:1 in all axes). + r"""Compute intersection of a line (ray) and a unit box (0:1 in all axes). Based on http://www.iquilezles.org/www/articles/intersectors/intersectors.htm @@ -463,7 +463,8 @@ def add(self, i, x, v, L): def interpolate(self, i): """Interpolate point with index `i` on path.""" - return interpolate(i, self.samplingpath.points, + return interpolate( + i, self.samplingpath.points, fwd_possible=self.samplingpath.fwd_possible, rwd_possible=self.samplingpath.rwd_possible, contourpath=self) diff --git a/ultranest/store.py b/ultranest/store.py index e1640835..50311515 100644 --- a/ultranest/store.py +++ b/ultranest/store.py @@ -157,6 +157,7 @@ class HDF5PointStore(FilePointStore): The format is a HDF5 file, which grows as needed. """ + FILES_OPENED = [] def __init__(self, filepath, ncols, **h5_file_args): @@ -169,7 +170,7 @@ def __init__(self, filepath, ncols, **h5_file_args): self.ncols = int(ncols) self.stack_empty = True h5_file_args['mode'] = h5_file_args.get('mode', 'a') - + # An annoying part of jupyter notebooks is that they keep all the variables # This means a old pointstore can survive, as we don't usually close them # Opening a new one with the same path will then fail with diff --git a/ultranest/utils.py b/ultranest/utils.py index 56888fd6..0cc5dbc4 100644 --- a/ultranest/utils.py +++ b/ultranest/utils.py @@ -57,7 +57,7 @@ def create_logger(module_name, log_dir=None, level=logging.INFO): def _makedirs(name): - """python2-compatible makedir """ + """python2-compatible makedir.""" # for Python2 compatibility: try: os.makedirs(name) @@ -118,7 +118,7 @@ def make_run_dir(log_dir, run_num=None, append_run_num=True): def vectorize(function): """Vectorize likelihood or prior_transform function.""" def vectorized(args): - """ vectorized version of function""" + """Vectorized version of function.""" return np.asarray([function(arg) for arg in args]) vectorized.__name__ = function.__name__ @@ -198,7 +198,7 @@ def resample_equal(samples, weights, rstate=None): def listify(*args): """ - concatenate args, which are (made to be) lists + Concatenate args, which are (made to be) lists. Parameters ---------- @@ -291,8 +291,7 @@ def vol_prefactor(n): def is_affine_transform(a, b): """ - checks if one group of points *a* is an affine transform - of another group of points *b*. + Check if one points *a* and *b* are related by an affine transform. The implementation currently returns False for rotations. @@ -312,7 +311,7 @@ def is_affine_transform(a, b): nb, db = b.shape assert n == nb assert db >= da - + n = (n // 2) * 2 a1 = a[0:n:2] a2 = a[1:n:2] @@ -329,9 +328,8 @@ def is_affine_transform(a, b): def normalised_kendall_tau_distance(values1, values2, i=None, j=None): """ - Normalised Kendall tau distance between two arrays, - *values1* and *values2*, both of length N. - + Normalised Kendall tau distance between two equally sized arrays. + see https://en.wikipedia.org/wiki/Kendall_tau_distance You can optionally pass precomputed indices:: @@ -364,10 +362,9 @@ def normalised_kendall_tau_distance(values1, values2, i=None, j=None): return ndisordered / (N * (N - 1)) - def _merge_transform_loglike_gradient_function(transform, loglike, gradient): def transform_loglike_gradient(u): - """combine transform, likelihood and gradient function""" + """Combine transform, likelihood and gradient function.""" p = transform(u.reshape((1, -1))) return p[0], loglike(p)[0], gradient(u) return transform_loglike_gradient @@ -375,9 +372,9 @@ def transform_loglike_gradient(u): def verify_gradient(ndim, transform, loglike, gradient, verbose=False, combination=False): """ - check with numerical differentiation if the gradient function - is plausibly correct. Raises AssertError otherwise. + Check with numerical differentiation if gradient function is plausibly correct. + Raises AssertError if not fulfilled. All functions are vectorized. Parameters