Skip to content

Commit

Permalink
fast min dist + misc changes
Browse files Browse the repository at this point in the history
  • Loading branch information
js51 committed Jul 16, 2024
1 parent 0c930ad commit d2a1bfe
Show file tree
Hide file tree
Showing 5 changed files with 29 additions and 2 deletions.
5 changes: 4 additions & 1 deletion cgt/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,4 +18,7 @@
from .position_paradigm import GenomeFramework, PositionParadigmFramework, Framework
from .models import Model
from . import distances
from . import plotting
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")
21 changes: 21 additions & 0 deletions cgt/distances.py
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
1 change: 1 addition & 0 deletions cgt/enums.py
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Expand Down
4 changes: 3 additions & 1 deletion cgt/simulations.py
Original file line number Diff line number Diff line change
Expand Up @@ -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":
Expand All @@ -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:
Expand Down
Binary file modified plot.pdf
Binary file not shown.

0 comments on commit d2a1bfe

Please sign in to comment.