New upstream version 0.5.3 (0945a217) · Commits · Debian Med / python-treetime

treetime/init.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -6,4 +6,4 @@ from .treetime import ttconf as treetime_conf
		from .gtr import GTR
		from .merger_models import Coalescent
		from .treeregression import TreeRegression
		version="0.5.2"
		version="0.5.3"

treetime/branch_len_interpolator.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -118,7 +118,7 @@ class BranchLenInterpolator (Distribution):

		@gamma.setter
		def gamma(self, value):
		self._gamma = value
		self._gamma = max(ttconf.TINY_NUMBER, value)

		@property
		def merger_cost(self):

treetime/clock_tree.py

+5 −4

Original line number	Diff line number	Diff line
		@@ -720,12 +720,13 @@ class ClockTree(TreeAnc):
		return ttconf.ERROR
		rate_std = np.sqrt(self.clock_model['cov'][0,0])

		current_rate = self.clock_model['slope']
		current_rate = np.abs(self.clock_model['slope'])
		upper_rate = self.clock_model['slope'] + rate_std
		lower_rate = max(0.1*current_rate, self.clock_model['slope'] - rate_std)
		for n in self.tree.find_clades():
		if n.up:
		n.branch_length_interpolator.gamma*=upper_rate/current_rate
		n._orig_gamma = n.branch_length_interpolator.gamma
		n.branch_length_interpolator.gamma = n._orig_gamma*upper_rate/current_rate

		self.logger("###ClockTree.calc_rate_susceptibility: run with upper bound of rate estimate", 1)
		self.make_time_tree(**params)
		@@ -733,7 +734,7 @@ class ClockTree(TreeAnc):
		for n in self.tree.find_clades():
		n.numdate_rate_variation = [(upper_rate, n.numdate)]
		if n.up:
		n.branch_length_interpolator.gamma*=lower_rate/upper_rate
		n.branch_length_interpolator.gamma = n._orig_gamma*lower_rate/current_rate

		self.logger("###ClockTree.calc_rate_susceptibility: run with lower bound of rate estimate", 1)
		self.make_time_tree(**params)
		@@ -741,7 +742,7 @@ class ClockTree(TreeAnc):
		for n in self.tree.find_clades():
		n.numdate_rate_variation.append((lower_rate, n.numdate))
		if n.up:
		n.branch_length_interpolator.gamma*=current_rate/lower_rate
		n.branch_length_interpolator.gamma = n._orig_gamma

		self.logger("###ClockTree.calc_rate_susceptibility: run with central rate estimate", 1)
		self.make_time_tree(**params)

treetime/config.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -27,7 +27,7 @@ MIN_INTEGRATION_PEAK = 0.001

		# clocktree parameters
		BRANCH_LEN_PENALTY = 0
		MAX_BRANCH_LENGTH = 1.5 # only relevant for time trees - upper boundary of interpolator objects
		MAX_BRANCH_LENGTH = 4.0 # only relevant for branch length optimization and time trees - upper boundary of interpolator objects
		NINTEGRAL = 300
		REL_TOL_PRUNE = 0.01
		REL_TOL_REFINE = 0.05

treetime/distribution.py

+27 −15

Original line number	Diff line number	Diff line
		@@ -17,7 +17,7 @@ class Distribution(object):
		"""

		@staticmethod
		def calc_fwhm(distribution, is_log=True):
		def calc_fwhm(distribution, is_neg_log=True):
		"""
		Assess the width of the probability distribution. This returns
		full-width-half-max
		@@ -25,29 +25,34 @@ class Distribution(object):

		if isinstance(distribution, interp1d):

		if is_log:
		if is_neg_log:
		ymin = distribution.y.min()
		prob = np.exp(-(distribution.y-ymin))
		log_prob = distribution.y-ymin
		else:
		prob = distribution.y
		log_prob = -np.log(distribution.y)
		log_prob -= log_prob.min()

		xvals = distribution.x

		elif isinstance(distribution, Distribution):
		# Distribution always stores log-prob
		# Distribution always stores neg log-prob with the peak value subtracted
		xvals = distribution._func.x
		prob = distribution.prob_relative(xvals)
		log_prob = distribution._func.y
		else:
		raise TypeError("Error in computing the FWHM for the distribution. "
		" The input should be either Distribution or interpolation object");

		x_idxs = binary_dilation(prob >= 0.4*(prob.max() - prob.min())+prob.min(), iterations=1)
		xs = xvals[x_idxs]
		if xs.shape[0] < 2:
		L = xvals.shape[0]
		# 0.69... is log(2), there is always one value for which this is true since
		# the minimum is subtracted
		tmp = np.where(log_prob < 0.693147)[0]
		x_l, x_u = tmp[0], tmp[-1]
		if L < 2:
		print ("Not enough points to compute FWHM: returning zero")
		return min(TINY_NUMBER, distribution.xmax - distribution.xmin)
		else:
		return xs.max() - xs.min()
		# need to guard against out-of-bounds errors
		return max(TINY_NUMBER, xvals[min(x_u+1,L-1)] - xvals[max(0,x_l-1)])


		@classmethod
		@@ -60,10 +65,12 @@ class Distribution(object):
		distribution.weight = weight
		return distribution


		@classmethod
		def shifted_x(cls, dist, delta_x):
		return Distribution(dist.x+delta_x, dist.y, kind=dist.kind)


		@staticmethod
		def multiply(dists):
		'''
		@@ -102,12 +109,13 @@ class Distribution(object):
		res = Distribution.delta_function(x_vals[0])
		else:
		res = Distribution(x_vals[ind], y_vals[ind], is_log=True,
		min_width=min_width, kind='linear')
		min_width=min_width, kind='linear', assume_sorted=True)

		return res


		def __init__(self, x, y, is_log=True, min_width = MIN_INTEGRATION_PEAK, kind='linear'):
		def __init__(self, x, y, is_log=True, min_width = MIN_INTEGRATION_PEAK,
		kind='linear', assume_sorted=False):

		"""
		Create Distribution instance
		@@ -117,8 +125,11 @@ class Distribution(object):
		if isinstance(x, Iterable) and isinstance (y, Iterable):

		self._delta = False # NOTE in classmethod this value is set explicitly to True.
		xvals, yvals = np.array(sorted(zip(x,y))).T
		# first, prepare x, y values
		if assume_sorted:
		xvals, yvals = x,y
		else:
		xvals, yvals = np.array(sorted(zip(x,y))).T
		if not is_log:
		yvals = -np.log(yvals)
		# just for safety
		@@ -135,7 +146,8 @@ class Distribution(object):
		yvals -= self._peak_val
		self._ymax = yvals.max()
		# store the interpolation object
		self._func= interp1d(xvals, yvals, kind=kind, fill_value=BIG_NUMBER, bounds_error=False)
		self._func= interp1d(xvals, yvals, kind=kind, fill_value=BIG_NUMBER,
		bounds_error=False, assume_sorted=True)
		self._fwhm = Distribution.calc_fwhm(self)

		elif np.isscalar(x):