diff --git a/cgt/__init__.py b/cgt/__init__.py index ddcdb42..c1c5a3e 100644 --- a/cgt/__init__.py +++ b/cgt/__init__.py @@ -18,4 +18,7 @@ from .position_paradigm import GenomeFramework, PositionParadigmFramework, Framework from .models import Model from . import distances -from . import plotting \ No newline at end of file +from . import plotting +import warnings +# This is required because we have complex numbers arising from numerical error that we want to ignore +warnings.filterwarnings("ignore", message="ComplexWarning") \ No newline at end of file diff --git a/cgt/distances.py b/cgt/distances.py index f5403e6..813bff0 100644 --- a/cgt/distances.py +++ b/cgt/distances.py @@ -516,6 +516,24 @@ def min_distance(framework, model, genome_instances=None, weighted=False): } +def fast_min_distance(framework, model, genome_instances=None): + num_genomes = int( + framework.genome_group().order() / framework.symmetry_group().order() + ) + if genome_instances is None: + genomes = list(framework.genomes().keys()) + genome_instances = genomes + dists = {} + for rep in genome_instances: + stepwise_dists = fast_step_probabilities(framework, model)[rep] + def a(k): + return stepwise_dists[k] + for k in range(0, num_genomes): + if a(k) > 10 ** (-10): + dists[rep] = k + break + return dists + def min_distance_using_irreps(framework, model, genome_instances=None): num_genomes = int( framework.genome_group().order() / framework.symmetry_group().order() @@ -646,6 +664,7 @@ def get_distance_function(distance_type): case DISTANCE.MFPT: return fast_MFPT case DISTANCE.discrete_MFPT: return discrete_MFPT case DISTANCE.min: return min_distance_using_irreps + case DISTANCE.min_exact: return fast_min_distance case DISTANCE.MLE: return mles case DISTANCE.min_weighted: return min_distance @@ -682,6 +701,8 @@ def distance_matrix( D[i, j] = distances[need_distances[p]] if replace_nan_with != np.nan: + if replace_nan_with == 'max': + replace_nan_with = D.max() + 0.1 D[np.isnan(D)] = replace_nan_with return D diff --git a/cgt/enums.py b/cgt/enums.py index 1fa91b6..f22de93 100644 --- a/cgt/enums.py +++ b/cgt/enums.py @@ -73,6 +73,7 @@ class DISTANCE(Enum): MFPT = mean_first_passage_time = auto() MLE = maximum_likelihood_distance = maximum_likelihood_estimate = auto() discrete_MFPT = DMFPT = discrete_mean_first_passage_time = auto() + min_exact = auto() class DATA(Enum): eig_data = 'eigen_data' diff --git a/cgt/simulations.py b/cgt/simulations.py index 5b01d15..22bc304 100644 --- a/cgt/simulations.py +++ b/cgt/simulations.py @@ -19,7 +19,7 @@ def newick_to_tree(newick_string): return tree -def evolve_on_tree(tree, framework, model, root="random", at_least_one_change=False, exactly_n_changes=None): +def evolve_on_tree(tree, framework, model, root="random", at_least_one_change=False, exactly_n_changes=None, exactly_bl_changes=False): # The tree should be a networkx DiGraph # Choose a genome at the root of the tree if root == "random": @@ -43,6 +43,8 @@ def evolve_on_tree(tree, framework, model, root="random", at_least_one_change=Fa # Decide how many rearrangements to apply if exactly_n_changes is not None: num_rearrangements = exactly_n_changes + elif exactly_bl_changes: + num_rearrangements = int(branch_length) else: num_rearrangements = np.random.poisson(branch_length) if at_least_one_change: diff --git a/plot.pdf b/plot.pdf index b007443..19af546 100644 Binary files a/plot.pdf and b/plot.pdf differ