KaneAI - World’s First E2E Software Testing Agent.Book a Demo 
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.
LambdaTest Learning Hubs:
- JUnit Tutorial
- TestNG Tutorial
- Webdriver Tutorial
- WebDriverIO Tutorial
- Protractor Tutorial
- Selenium 4 Tutorial
- Jenkins Tutorial
- NUnit Tutorial
- Jest Tutorial
- Playwright Tutorial
- Cypress Tutorial
- PyTest Tutorial
YouTube
You could also refer to video tutorials over LambdaTest YouTube channel to get step by step demonstration from industry experts.
Run molecule automation tests on LambdaTest cloud grid
Perform automation testing on 3000+ real desktop and mobile devices online.
Try LambdaTest Now !!
Get 100 minutes of automation test minutes FREE!!
Was this article helpful?
