How to use dependency_sequence method in molecule

Best Python code snippet using molecule_python

test_diff.py

Source:test_diff.py

1import itertools2import numpy as np3import pytest4import sage.all5from scipy.special import factorial6import gf.gf as gflib7import gf.mutations as mutations8#import gf.togimble as tg9import gf.matrix_representation as gfmat10import gf.diff as gfdiff11@pytest.mark.taylor12class Test_taylor_series_coefficients:13	def test_taylor_polynomial(self):14		denom = np.arange(1,5, dtype=int)15		var_array = np.arange(1,5, dtype=float)/1016		diff_array = (2,2)17		num_mutypes = 218		result = gfdiff.taylor_coeff_inverse_polynomial(denom, var_array, diff_array, num_mutypes)19		expected = 3.5555620		assert np.isclose(expected, result)21	def test_taylor_exponential(self):22		denom = np.arange(1,5, dtype=int)23		var_array = np.arange(1,5, dtype=float)/1024		diff_array = (2,2)25		num_mutypes = 226		time = 227		exponential_part = np.exp(-2*np.dot(var_array, denom))28		result = gfdiff.taylor_coeff_exponential(-time, denom, exponential_part, diff_array, num_mutypes)29		expected = 1.4277630		assert np.isclose(expected, result)31	def test_prod_of_polynomials(self):32		eq_matrix = np.array([33    		[[1,0,0,0],[1,1,1,1]],34    		[[0,1,0,0],[1,1,1,1]],35    		[[1,0,0,0],[1,0,1,0]]],36    		dtype=np.int8)37		shape = (2,2)38		var_array = np.array([0.1,0.2,0.3,0.3], dtype=np.float64)39		symbolic_var_array = np.array(40			[sage.all.Rational(0.1), sage.all.Rational(0.2), sage.all.SR.var('m1'), sage.all.SR.var('m2')], 41			dtype=object42			)43		result = gfdiff.compile_non_inverted_eq(eq_matrix, shape, shape)(var_array)44		subsetdict = gfdiff.product_subsetdict(shape)45		combined_result = gfdiff.product_f(subsetdict, result)46		#from symbolic eq47		subs_dict = {b:var_array[-1] for b in symbolic_var_array[-2:]}48		symbolic_eq = np.prod(gfmat.equations_from_matrix(eq_matrix, symbolic_var_array))49		expected = return_symbolic_result(symbolic_eq, subs_dict, symbolic_var_array[-2:], shape)50		assert np.allclose(expected, combined_result)51	@pytest.mark.parametrize("eq_matrix, delta_in_nom",52		[(np.array([53    		[[1,0,0,0],[1,1,1,1]],54    		[[0,1,0,0],[2,1,0,2]],55    		[[1,0,0,0],[1,1,1,0]]],56    		dtype=np.int8), True),57		(np.array([58    		[[1,0,0,0],[1,1,1,1]],59    		[[1,0,0,0],[2,1,0,2]],60    		[[1,0,0,0],[1,1,1,0]]],61    		dtype=np.int8), False)])62	def test_diff_inverse_laplace(self, eq_matrix, delta_in_nom):63		dummy_variable_idx = 164		dummy_variable = sage.all.SR.var('E')65		eq_matrix_no_delta = np.delete(eq_matrix, dummy_variable_idx, axis=2)66		shape = (2,2)67		var_array = np.array([0.1,0.3,0.3], dtype=np.float64)68		symbolic_var_array = np.array(69			[sage.all.Rational(0.1), dummy_variable, sage.all.SR.var('m1'), sage.all.SR.var('m2')], 70			dtype=object71			)72		time = 1.073		subsetdict = gfdiff.product_subsetdict(shape)74		result = gfdiff.compile_inverted_eq(eq_matrix_no_delta, shape, subsetdict, delta_in_nom, shape)(var_array, time)75		print(result)76		#from symbolic eq77		subs_dict = {b:var_array[-1] for b in symbolic_var_array[-2:]}78		subs_dict[sage.all.SR.var('T')] = time79		symbolic_eq = np.prod(gfmat.equations_from_matrix(eq_matrix, symbolic_var_array))80		inverted_symbolic_eq = gflib.return_inverse_laplace(symbolic_eq, dummy_variable)81		expected = return_symbolic_result(inverted_symbolic_eq, subs_dict, symbolic_var_array[-2:], shape)82		assert np.allclose(expected, result)83	@pytest.mark.parametrize("eq_matrix, delta_in_denom, delta_in_nom",84		[(np.array([85    	[[1,0,0,0],[1,0,1,1]],86    	[[0,1,0,0],[2,1,0,2]],87    	[[1,0,0,0],[1,0,1,0]]],88    	dtype=np.int8), np.array([False, True, False], dtype=bool), True),89    	(np.array([90    		[[1,0,0,0],[1,0,1,1]],91    		[[1,0,0,0],[2,1,0,2]],92    		[[1,0,0,0],[1,0,1,0]]],93    		dtype=np.int8), np.array([False, True, False], dtype=bool), False)94		])95	def test_diff_inverse_laplace2(self, eq_matrix, delta_in_denom, delta_in_nom):96		#eq in which some factors have dummy variable, others don't97		dummy_variable_idx = 198		dummy_variable = sage.all.SR.var('E')99		eq_matrix_no_delta = np.delete(eq_matrix, dummy_variable_idx, axis=2)100		shape = (2,2)101		var_array = np.array([0.1,0.3,0.3], dtype=np.float64)102		symbolic_var_array = np.array(103			[sage.all.Rational(0.1), dummy_variable, sage.all.SR.var('m1'), sage.all.SR.var('m2')], 104			dtype=object105			)106		time = 1.0107		subsetdict = gfdiff.product_subsetdict(shape)108		result_inverted_part = gfdiff.compile_inverted_eq(eq_matrix_no_delta[delta_in_denom], shape, subsetdict, delta_in_nom, shape)(var_array, time)109		result_non_inverted_part = gfdiff.compile_non_inverted_eq(eq_matrix_no_delta[~delta_in_denom], shape, shape)(var_array)110		result = gfdiff.product_f(subsetdict, np.vstack((result_inverted_part[None, :], result_non_inverted_part)))111		print('result')112		print(result)113		#from symbolic eq114		subs_dict = {b:var_array[-1] for b in symbolic_var_array[-2:]}115		subs_dict[sage.all.SR.var('T')] = time116		symbolic_eq = np.prod(gfmat.equations_from_matrix(eq_matrix, symbolic_var_array))117		inverted_symbolic_eq = gflib.return_inverse_laplace(symbolic_eq, dummy_variable)118		expected = return_symbolic_result(inverted_symbolic_eq, subs_dict, symbolic_var_array[-2:], shape)119		print('expected')120		print(expected)121		assert np.allclose(expected, result)122def return_symbolic_result(eq, subs_dict, branchtypes, shape, root=(0,0)):123	symbolic_diff = np.zeros(shape, dtype=np.float64)124	for mutype in itertools.product(*(range(i) for i in shape)):125		if mutype == root:126			symbolic_diff[mutype] = eq.subs(subs_dict)127		else:128			how_to_diff = mutations.single_partial(branchtypes, mutype)129			symbolic_diff[mutype] = 1/np.prod(factorial(mutype)) * sage.all.diff(eq, *how_to_diff).subs(subs_dict)130			131	return symbolic_diff132@pytest.mark.collapse133class Test_collapse_graph:134	def test_collapse_graph(self):135		graph_array = ((1,2,4,5),(2,),(3,),(6,),(3,),(3,),tuple())136		eq_matrix = np.array([137			[[0,0],[0,1]],138			[[0,1],[0,1]],139			[[0,0],[0,0]]],140			dtype=np.uint8141			)142		adjacency_matrix = np.full((len(graph_array), len(graph_array)), fill_value=255, dtype=np.int8)143		adjacency_matrix[np.array([0,0,1]), np.array([1,2,2])] = 0144		adjacency_matrix[np.array([2]), np.array([3])] = 1145		adjacency_matrix[np.array([0,0,3,4,5]), np.array([4,5,6,3,3])] = 2146		sample_list = [('a', 'a'),('b', 'b')]147		coalescence_rates = (0,1)148		branchtype_dict = {}149		exodus_rate = 1150		exodus_direction = [(1,0)]151		gfObj = gflib.GFMatrixObject(sample_list, coalescence_rates, branchtype_dict, exodus_rate=exodus_rate, exodus_direction=exodus_direction)152		collapsed_graph_array, adjacency_matrix_b, eq_array, to_invert_array = gfObj.collapse_graph(graph_array, adjacency_matrix, eq_matrix)153		expected_graph_array= ((1, 2, 5, 6), (3,), (3,), (4,), (), (4,), (4,))154		expected_to_invert_array = np.zeros(9, dtype=bool)155		expected_to_invert_array[-2:] = 1156		expected_eq_array = ((2,), (2,), (2,), (2,), (2,), (2,), (2,), (0, 1), (0, 0, 1))157		def compare(ar1, ar2):158			for a, b in zip(ar1, ar2):159				assert np.array_equal(np.array(a), np.array(b))160		compare(expected_graph_array, collapsed_graph_array)161		compare(expected_eq_array, eq_array)162		assert np.array_equal(expected_to_invert_array, to_invert_array)163	@pytest.mark.parametrize(164		'sample_list, k_max, branchtype_dict, exp_graph_array',165		[166		([(), ('a', 'a')], {'m_1':2}, {'a':0}, [(1,3), (2,), (),()]),167		([(), ('a', 'a', 'a')], {'m_1':2, 'm_2':2}, {'a':0, 'aa':1}, [(1, 4, 5), (2,), (3,), (), (3,), ()]),168		])169	def test_graph_with_multiple_endpoints(self, sample_list, k_max, branchtype_dict, exp_graph_array):170		gfobj = self.get_gf_no_mutations(sample_list, k_max, branchtype_dict)171		delta_idx = gfobj.exodus_rate172		graph_array, adjacency_matrix, eq_matrix = gfobj.make_graph()173		collapsed_graph_array, *_ = gfobj.collapse_graph(graph_array, adjacency_matrix, eq_matrix)		174		print(exp_graph_array)175		print(collapsed_graph_array)176		for o,e in zip(collapsed_graph_array, exp_graph_array):177			assert o==e178	def get_gf_no_mutations(self, sample_list, k_max, branchtype_dict):179		coalescence_rate_idxs = (0, 1)180		exodus_rate_idx = 2 181		exodus_direction = [(1,0),]182		mutype_labels, max_k = zip(*sorted(k_max.items()))183		gfobj = gflib.GFMatrixObject(184			sample_list, 185			coalescence_rate_idxs, 186			branchtype_dict,187			exodus_direction=exodus_direction,188			exodus_rate=exodus_rate_idx189			)190		return gfobj191@pytest.mark.taylor2192class Test_taylor2:193	@pytest.mark.parametrize('size', [2, 3])194	def test_combining_probabilities(self, size):195		gfobj = self.get_gf_no_mutations(size)196		max_k = np.full(size-1,fill_value=2, dtype=int)197		shape = tuple(max_k+1)198		variable_array = np.array([1., 2.], dtype=np.float64)199		theta = .5200		theta_array = np.full(size-1, fill_value=theta)201		variable_array = np.hstack((variable_array, theta_array))202		time = 1.5203		ordered_mutype_list = [sage.all.SR.var(f'm_{idx}') for idx in range(1,size)]204		num_mutypes = len(ordered_mutype_list)205		alt_variable_array = np.hstack((variable_array[:2], np.array(ordered_mutype_list)))206		result = self.evaluate_graph2(gfobj, shape, theta, variable_array, time)207		result_with_marginals = self.evaluate_graph_marginals(gfobj, max_k, theta, variable_array, time)208		print(result_with_marginals)209		exp_result = self.evaluate_symbolic_equation(gfobj, ordered_mutype_list, max_k, theta, alt_variable_array, time)210		subidx = tuple([slice(0,s) for s in shape])211		print(exp_result[subidx])212		assert np.allclose(exp_result[subidx], result)213		assert np.allclose(exp_result, result_with_marginals)214	def test_IM_models(self, get_IM_gfobject):215		gfobj, variable_array = get_IM_gfobject216		num_variables = gfobj.num_variables if gfobj.exodus_rate is None else gfobj.num_variables-1217		max_k = np.array([2,2,2,2], dtype=int)218		ordered_mutype_list = [sage.all.SR.var(f'm_{idx}') for idx in range(1,len(max_k)+1)]219		shape = tuple(max_k+1)220		#variables depending on model: c0, c1, c2, M, E221		#variable_array = (np.arange(1,num_variables+1)/10).astype(np.float64)222		#variable_array = np.array([1.0, 0.5, 0.9, 0.001], dtype=np.float64)[:num_variables]223		theta = .51224		theta_array = np.full(len(max_k), fill_value=theta)225		var = np.hstack((variable_array, theta_array))226		#var_symbolic = np.hstack((variable_array, ordered_mutype_list))227		if gfobj.exodus_rate is not None:228			var_sage = np.zeros(gfobj.num_variables, dtype=object)229			var_sage[:-1] = [sage.all.Rational(v) for v in variable_array]230			var_sage[-1] = sage.all.SR('E')231			var_symbolic = np.hstack((var_sage, ordered_mutype_list))232		else:233			var_symbolic = np.hstack((variable_array, ordered_mutype_list))234		time = 1.5235		result = self.evaluate_graph2(gfobj, shape, theta, var, time)236		#print(result)237		result_with_marginals = self.evaluate_graph_marginals(gfobj, max_k, theta, var, time)238		#print(result_with_marginals)239		expected_result = self.evaluate_symbolic_equation(gfobj, ordered_mutype_list, max_k, theta, var_symbolic, time, sage_inverse=True)240		#print(expected_result)241		subidx = tuple([slice(0,s) for s in shape])242		assert np.allclose(expected_result[subidx], result)243		assert np.allclose(expected_result, result_with_marginals)244		245	def get_gf_no_mutations(self, size):246		sample_list = [(), ('a',)*size]247		coalescence_rate_idxs = (0, 1)248		exodus_rate_idx = 2 249		exodus_direction = [(1,0),]250		k_max = {f'm_{idx}':2 for idx in range(1,size)}251		mutype_labels, max_k = zip(*sorted(k_max.items()))252		branchtype_dict_mat = {'a'*idx:idx-1 for idx in range(1,size)} 253		gfobj = gflib.GFMatrixObject(254			sample_list, 255			coalescence_rate_idxs, 256			branchtype_dict_mat,257			exodus_direction=exodus_direction,258			exodus_rate=exodus_rate_idx259			)260		return gfobj261	def evaluate_graph(self, gfobj, shape, theta, var, time):262		delta_idx = gfobj.exodus_rate263		graph_array, adjacency_matrix, eq_matrix = gfobj.make_graph()264		collapsed_graph_array, adjacency_matrix, eq_array, to_invert = gfobj.collapse_graph(graph_array, adjacency_matrix, eq_matrix)		265		dependency_sequence = gfdiff.resolve_dependencies(collapsed_graph_array)266		#print(dependency_sequence)267		subsetdict = gfdiff.product_subsetdict(shape)268		f_non_inverted, f_inverted = gfdiff.prepare_graph_evaluation(eq_matrix, to_invert, eq_array, shape, delta_idx, subsetdict, shape)269		evaluator = gfdiff.evaluate_single_point(shape, f_non_inverted, *f_inverted)270		results = evaluator(var, time)271		final_result = gfdiff.iterate_graph(dependency_sequence, collapsed_graph_array, adjacency_matrix, results, subsetdict)272		final_result = final_result[0]273		multiplier_matrix = gfdiff.taylor_to_probability(gfdiff.taylor_to_probability_coeffs(shape, include_marginals=False), theta)274		assert final_result.shape==multiplier_matrix.shape275		return multiplier_matrix * final_result276	def evaluate_graph2(self, gfobj, shape, theta, var, time):277		delta_idx = gfobj.exodus_rate278		eq_graph_array, eq_array, to_invert, eq_matrix = gfobj.equations_graph()		279		dependency_sequence = gfdiff.resolve_dependencies(eq_graph_array)280		subsetdict = gfdiff.product_subsetdict(shape)281		f_non_inverted, f_inverted = gfdiff.prepare_graph_evaluation(eq_matrix, to_invert, eq_array, shape, delta_idx, subsetdict, shape)282		evaluator = gfdiff.evaluate_single_point(shape, f_non_inverted, *f_inverted)283		results = evaluator(var, time)284		final_result = gfdiff.iterate_eq_graph(dependency_sequence, eq_graph_array, results, subsetdict)285		#final_result = final_result[0]286		multiplier_matrix = gfdiff.taylor_to_probability(gfdiff.taylor_to_probability_coeffs(shape, include_marginals=False), theta)287		assert final_result.shape==multiplier_matrix.shape288		return multiplier_matrix * final_result289	def evaluate_graph_marginals(self, gfobj, k_max, theta, var, time):290		delta_idx = gfobj.exodus_rate291		eq_graph_array, eq_array, to_invert, eq_matrix = gfobj.equations_graph()		292		dependency_sequence = gfdiff.resolve_dependencies(eq_graph_array)293		final_result_shape = k_max+2	294		marg_iterator = gfdiff.marginals_nuissance_objects(k_max)295		marg_boolean, shapes, mutype_shapes, subsetdicts, slices = marg_iterator296		f_array = gfdiff.prepare_graph_evaluation_with_marginals(297			eq_matrix, 298			to_invert, 299			eq_array,300			marg_iterator, 301			delta_idx302			)303		num_eq_non_inverted = np.sum(to_invert==0) 304		num_eq_tuple = (num_eq_non_inverted, to_invert.size - num_eq_non_inverted)305		evaluator = gfdiff.evaluate_single_point_with_marginals(306			k_max, 307			f_array,308			num_eq_tuple,309			slices310			)311		results = evaluator(var, time)312		subsetdict_with_marginals = gfdiff.product_subsetdict_marg(tuple(final_result_shape))313		final_result = gfdiff.iterate_eq_graph(dependency_sequence, eq_graph_array, results, subsetdict_with_marginals)314		#final_result = gfdiff.iterate_eq_graph_with_marginals(dependency_sequence, eq_graph_array, results, subsetdicts, slices, shapes, final_result_shape)315		multiplier_matrix = gfdiff.taylor_to_probability(gfdiff.taylor_to_probability_coeffs(k_max+1, include_marginals=True), theta)316		assert final_result.shape==multiplier_matrix.shape317		return multiplier_matrix * final_result	318	def evaluate_symbolic_equation(self, gfobj, ordered_mutype_list, max_k, theta, var, time, sage_inverse=False):319		theta = sage.all.Rational(theta)320		rate_dict = {b:theta for b in ordered_mutype_list}321		paths, eq_matrix = gfobj.make_gf()322		if not sage_inverse:323			alt_eqs = equations_from_matrix_with_inverse(eq_matrix, paths, var, time, gfobj.exodus_rate)324		else:325			alt_eqs = equations_with_sage(eq_matrix, paths, var, sage.all.Rational(time), gfobj.exodus_rate)326		gf_alt = sum(alt_eqs)327		result = mutations.depth_first_mutypes(max_k, ordered_mutype_list, gf_alt, theta, rate_dict)328		return result.astype(np.float64)329	@pytest.fixture(330	scope='class', 331	params=[332		([(1,2,0)], sage.all.SR.var('E'), None, None, np.array([.1, .2, .3], dtype=np.float64)),333		(None, None, [(2,1)], sage.all.SR.var('M'), np.array([.1, .2, .3, .4], dtype=np.float64)),334		([(1,2,0)], sage.all.SR.var('E'), [(2,1)], sage.all.SR.var('M'), np.array([1.0, 0.5, 0.9, 0.001], dtype=np.float64))335		],336	ids=[337		'DIV',338		'MIG', 339		'IM'340		],341	)342	def get_IM_gfobject(self, request):343		sample_list = [(),('a','a'),('b','b')]344		ancestral_pop = 0345		coalescence_rates = (sage.all.SR.var('c0'), sage.all.SR.var('c1'), sage.all.SR.var('c2'))346		coalescence_rate_idxs = (0, 1, 2)347		k_max = {'m_1':2, 'm_2':2, 'm_3':2, 'm_4':2}348		mutype_labels, max_k = zip(*sorted(k_max.items()))349		branchtype_dict_mat = gfmat.make_branchtype_dict_idxs(sample_list, mapping='unrooted', labels=mutype_labels)350		exodus_direction, exodus_rate, migration_direction, migration_rate, variable_array = request.param351		352		variables_array = list(coalescence_rates)353		migration_rate_idx, exodus_rate_idx = None, None354		if migration_rate!=None:355			migration_rate_idx = len(variables_array)356			variables_array.append(migration_rate)357		if exodus_rate!=None:358			exodus_rate_idx = len(variables_array)359			variables_array.append(exodus_rate)360		variables_array += [sage.all.SR.var(m) for m in mutype_labels]361		variables_array = np.array(variables_array, dtype=object)362		363		gfobj = gflib.GFMatrixObject(364			sample_list,365			coalescence_rate_idxs, 366			branchtype_dict_mat,367			exodus_rate=exodus_rate_idx,368			exodus_direction=exodus_direction,369			migration_rate=migration_rate_idx,370			migration_direction=migration_direction371			)372		return gfobj, variable_array373def equations_from_matrix_with_inverse(multiplier_array, paths, var_array, time, delta_idx):374	split_paths = gfmat.split_paths_laplace(paths, multiplier_array, delta_idx)375	delta_in_nom_all = multiplier_array[:, 0, delta_idx]==1376	results = np.zeros(len(split_paths), dtype=object)377	subset_no_delta = np.arange(multiplier_array.shape[-1])!=delta_idx378	multiplier_array_no_delta = multiplier_array[:,:,subset_no_delta] 379	for idx, (no_delta, with_delta) in enumerate(split_paths):380		delta_in_nom_list = delta_in_nom_all[with_delta]381		inverse = gflib.inverse_laplace_single_event(multiplier_array_no_delta[with_delta], var_array, time, delta_in_nom_list)382		if isinstance(inverse, np.ndarray):383			inverse = np.sum(inverse)384		no_inverse = np.prod(gfmat.equations_from_matrix(multiplier_array_no_delta[no_delta], var_array))385		results[idx] = np.prod((inverse, no_inverse))386	return results387def equations_with_sage(multiplier_array, paths, var_array, time, delta_idx):388	delta = var_array[delta_idx] if delta_idx is not None else None389	eqs = np.zeros(len(paths), dtype=object)390	for i, path in enumerate(paths):391		ma = multiplier_array[np.array(path, dtype=int)]392		temp = np.prod(gfmat.equations_from_matrix(ma, var_array))393		eqs[i] = gflib.return_inverse_laplace(temp, delta).subs({sage.all.SR.var('T'):time})...

evaluate.py

Source:evaluate.py

1import numpy as np2import numba3import agemo.diff as gfdiff4import agemo.mutations as gfmuts5# getting out bSFS6class GfEvaluatorLegacy:7    def __init__(self, gfobj, k_max, mutype_array):8        delta_idx = gfobj.exodus_rate  # what value if None9        # only works with single delta_idx!10        (11            self.eq_graph_array,12            eq_array,13            to_invert,14            eq_matrix,15        ) = gfobj.equations_graph()16        self.dependency_sequence = gfdiff.resolve_dependencies(self.eq_graph_array)17        self.num_branchtypes = len(k_max)18        self.final_result_shape = k_max + 219        size = len(mutype_array)20        # marg_iterator = gfdiff.marginals_nuissance_objects(k_max)21        # slices = marg_iterator[-1]22        # prepare_graph_evaluation_with_marginals_alt(eq_matrix, to_invert_array,23        # eq_array, size, delta_idx, subsetdict, mutype_array, mutype_shape)24        self.subsetdict = gfdiff.product_subsetdict_marg(25            tuple(self.final_result_shape), mutype_array26        )27        f_tuple = gfdiff.prepare_graph_evaluation_with_marginals(28            eq_matrix,29            to_invert,30            eq_array,31            size,32            delta_idx,33            self.subsetdict,34            mutype_array,35            self.final_result_shape,36        )37        num_eq_non_inverted = np.sum(to_invert == 0)38        num_eq_tuple = (39            num_eq_non_inverted,40            to_invert.size - num_eq_non_inverted,41        )42        self.evaluator = gfdiff.evaluate_single_point_with_marginals(43            size, num_eq_tuple, f_tuple44        )45        self.multiplier_matrix = gfdiff.taylor_to_probability_coeffs(46            mutype_array, self.final_result_shape, include_marginals=True47        )48    def evaluate(self, theta, var, time):49        try:50            results = self.evaluator(var, time)51        except ZeroDivisionError:52            var = self.adjust_parameters(var)53            results = self.evaluator(var, time)54        final_result_flat = gfdiff.iterate_eq_graph(55            self.dependency_sequence,56            self.eq_graph_array,57            results,58            self.subsetdict,59        )60        theta_multiplier_matrix = gfdiff.taylor_to_probability(61            self.multiplier_matrix, theta62        )63        no_marginals = (theta_multiplier_matrix * final_result_flat).reshape(64            self.final_result_shape65        )66        # filling of array needed here!!!67        # final_result = np.zeros(self.final_result_shape, dtype=np.float64)68        return gfmuts.adjust_marginals_array(no_marginals, self.num_branchtypes)69    def adjust_parameters(self, var, factor=1e-5):70        epsilon = (71            np.random.randint(low=-100, high=100, size=len(var) - self.num_branchtypes)72            * factor73        )74        var[: -self.num_branchtypes] += epsilon75        return var76class BSFSEvaluator:77    """78    Class containing all methods needed to extract the block-wise site frequency spectrum from the generating function.79    :param gf: object defining the structured coalescent model for which the bSFS will be calculated.80    :type gf: class `agemo.GfMatrixObject`81    :param mutation_types: MutationTypeCounter object containing all information on the mutation types to be included in the bSFS,82    :type mutation_types: class `MutationTypeCounter`83    """84    def __init__(self, gfObj, MutationTypeCounter):85        if MutationTypeCounter.phased:86            raise NotImplementedError(87                "Calculating the bSFS for the fully phased case is still under development."88            )89        num_discrete_events = len(gfObj.discrete_events)90        if num_discrete_events == 0:91            delta_idx = None92        elif num_discrete_events == 1:93            delta_idx = gfObj.discrete_events[0]94        else:95            raise NotImplementedError(96                "BSFSEvaluator can only deal with 1 discrete event."97            )98        # only works with single delta_idx!99        (100            self.eq_graph_array,101            eq_array,102            to_invert,103            eq_matrix,104        ) = gfObj.equations_graph()105        self.dependency_sequence = gfdiff.resolve_dependencies(self.eq_graph_array)106        self.eq_graph_array = numba.typed.List(self.eq_graph_array)107        self.final_result_shape = MutationTypeCounter.mutype_shape108        size, self.num_branchtypes = MutationTypeCounter.all_mutypes.shape109        self.simple_reshape = np.prod(self.final_result_shape) == size110        self.subsetdict = self.make_product_subsetdict(MutationTypeCounter)111        self.all_mutypes_ravel = MutationTypeCounter.all_mutypes_ravel112        # final step:113        # eventually also: mapping mutypes with same probability (part of MutationTypeCounter obj)114        f_tuple = gfdiff.prepare_graph_evaluation_with_marginals(115            eq_matrix,116            to_invert,117            eq_array,118            size,119            delta_idx,120            self.subsetdict,121            MutationTypeCounter.all_mutypes,122            self.final_result_shape,123        )124        num_eq_non_inverted = np.sum(to_invert == 0)125        num_eq_tuple = (126            num_eq_non_inverted,127            to_invert.size - num_eq_non_inverted,128        )129        self.evaluator = gfdiff.evaluate_single_point_with_marginals(130            size, num_eq_tuple, f_tuple131        )132        self.multiplier_matrix = gfdiff.taylor_to_probability_coeffs(133            MutationTypeCounter.all_mutypes,134            self.final_result_shape,135            include_marginals=True,136        )137    def evaluate(self, theta, var, time=0):138        """139        Calculate the bSFS for a single point in paremeter space.140        :param theta: Scaled mutation rate (units 2Ne).141        :type theta: float142        :param var: array of floats containing the parameter values for ,143            First n entries are the coalescence rates for the populations as provided to the MutationTypeCounter.144            Order of the next elements determined by the respective indices of each of the events.145            Last m entries should equal theta. With m the number of branch types.146        :type var: class `np.ndarray(float)`147        :param time: Time to the discrete event. Only relevant in the presence of a discrete event.148        :type time: float149        :return bSFS: Shape is defined by k_{max} as defined by MutationTypeCounter for each branch type.150        :rtype bSFS: class `np.ndarray(float)`151        """152        try:153            results = self.evaluator(var, time)154        except ZeroDivisionError:155            var = self.adjust_parameters(var)156            results = self.evaluator(var, time)157        final_result_flat = gfdiff.iterate_eq_graph(158            self.dependency_sequence,159            self.eq_graph_array,160            results,161            self.subsetdict,162        )163        theta_multiplier_matrix = gfdiff.taylor_to_probability(164            self.multiplier_matrix, theta165        )166        no_marginals = theta_multiplier_matrix * final_result_flat167        if self.simple_reshape:168            final_result = no_marginals.reshape(self.final_result_shape)169        else:170            final_result = np.zeros(self.final_result_shape, dtype=np.float64)171            final_result.flat[self.all_mutypes_ravel] = no_marginals172        return gfmuts.adjust_marginals_array(final_result, self.num_branchtypes)173    def adjust_parameters(self, var, factor=1e-5):174        """175        Adjusts the var array in case we run into a zero-division error and multiplies each entry176        by a small factor epsilon.177        :param var: Array containing variables values.178        :type var: class `np.ndarray(float)`179        :param factor: multiplying factor used to determine epsilon.180        :type factor: float181        """182        epsilon = (183            np.random.randint(low=-100, high=100, size=len(var) - self.num_branchtypes)184            * factor185        )186        var[: -self.num_branchtypes] += epsilon187        return var188    def make_product_subsetdict(self, MutationTypeCounter):189        """190        Method returns dict containing the indices needed to calculate the product of two polynomials.191        """192        product_subsetdict_with_gaps = gfdiff.product_subsetdict_marg(193            self.final_result_shape, MutationTypeCounter.all_mutypes194        )195        if self.simple_reshape:196            return product_subsetdict_with_gaps197        else:198            reverse_mapping = {199                value: idx200                for idx, value in enumerate(MutationTypeCounter.all_mutypes_ravel)201            }202            product_subsetdict_no_gaps = numba.typed.Dict()203            for key, value in product_subsetdict_with_gaps.items():204                product_subsetdict_no_gaps[reverse_mapping[key]] = np.array(205                    [reverse_mapping[v] for v in value], dtype=np.uint64206                )...

dependencies.py

Source:dependencies.py

...48    def show_routes(self, current):49        for route in self._dependency_routes(current):50            print route51    @property52    def dependency_sequence(self):53        routes_dict = dict()54        end_nodes = list()55        dependency_sequence_list = list()56        for start in self.iterkeys():57            dependency_routes = self._dependency_routes(start)58            if len(dependency_routes) == 0:59                end_nodes.append(start)60            for route in dependency_routes:61                routes_dict.setdefault(len(route), []).append(route)62        for route_index in sorted(routes_dict, reverse=True):63            for route in routes_dict[route_index]:64                last_dependency = route[-1]65                if last_dependency not in dependency_sequence_list:66                    dependency_sequence_list.append(last_dependency)...

Automation Testing Tutorials

Learn to execute automation testing from scratch with LambdaTest Learning Hub. Right from setting up the prerequisites to run your first automation test, to following best practices and diving deeper into advanced test scenarios. LambdaTest Learning Hubs compile a list of step-by-step guides to help you be proficient with different test automation frameworks i.e. Selenium, Cypress, TestNG etc.