Merge pull request #17 from js51/develop

More distances, pytest now working with sage

.gitignore

@@ -135,3 +135,4 @@ examples/_output/.DS_Store
 examples/.DS_Store
 cgt/scratchpad.ipynb
 _output
+.vscode/settings.json

VERSION

@@ -1 +1 @@
-0.1.0
+0.1.1

cgt/__init__.py

@@ -1,12 +1,13 @@
 """
 """
-
 name = "cgt"
-from . import visualisation
-from . import models
-from . import position_paradigm
-from . import enums
-from . import rearrangements
+from sage.all_cmdline import gap
+gigs = 30
+old_command = gap._Expect__command
+s = old_command.index('-o ')
+e = old_command.index(' -', s)
+gap._Expect__command = gap._Expect__command.replace(gap._Expect__command[s:e], f'-o {gigs}G')
+from .examples import *
 from .enums import *
 from .position_paradigm import PositionParadigmFramework
-from .models import Model
+from .models import Model

cgt/distances.py

@@ -1,14 +1,107 @@
 """
 """
+from cgt.enums import ALGEBRA
+import numpy as np
+import networkx as nx
+from sage.all import ComplexDoubleField, UniversalCyclotomicField, matrix, Matrix, real, exp, round
+from scipy.optimize import minimize
+from warnings import warn
+import matplotlib.pyplot as plt
 
-def mle(framework, model, genomes):
-    """Return dictionary of MLEs (if they exist) for time elapsed rearranging ref->genome for each genome"""
+def mles(framework, model, genome_reps, plots=False):
+    """
+    Return dictionary of MLEs (if they exist) for time elapsed rearranging ref->genome for each genome
+    """
+    warn("This... might take a while!")
+    CDF, UCF = ComplexDoubleField(), UniversalCyclotomicField()
+    G = framework.genome_group()
+    Z = framework.symmetry_group()
+    z = framework.symmetry_element()
+    s = model.s_element(in_algebra=ALGEBRA.genome)
+    irreps_of_z, irreps_of_s = framework.irreps(z), framework.irreps(s)
+    irreps_of_zs = [ Matrix(CDF, matrix(UCF, irrep_z*irrep_s)) for irrep_z, irrep_s in zip(irreps_of_z, irreps_of_s) ]
+    for irrep in irreps_of_zs:
+	    irrep.set_immutable()
+    dims = [irrep_of_zs.nrows() for irrep_of_zs in irreps_of_zs]
+    eig_lists   = [_eigenvalues(irrep_zs, round_to=9, make_real=True, inc_repeated=False) for irrep_zs in irreps_of_zs]
+    projections = [_projection_operators(*vals) for vals in zip(irreps_of_zs, eig_lists)]
+    traces = {
+        r: {
+            genome_rep: { 
+                eigenvalue: {} for eigenvalue in eig_lists[r]
+            } for genome_rep in genome_reps
+        } for r in range(len(irreps_of_zs))
+    }
+    inv_sigmas = { instance : framework.cycles(instance.inverse()) for instance in genome_reps } # pre-compute for efficiency
+    irreps = framework.irreps()
+    for r, irrep in enumerate(irreps): # Iterate over irreducible representations
+        for instance in genome_reps:
+            sigd = Matrix(CDF, matrix(UCF, irrep(inv_sigmas[instance]))) * irreps_of_z[r]
+            for e, eigenvalue in enumerate(eig_lists[r]):
+                traces[r][instance][eigenvalue] = round(real((sigd*projections[r][e]).trace()), 6)
+    def likelihood(t, instance):
+        ans = 0
+        for r in range(len(irreps_of_zs)):
+            ans += Z.order()*(exp(-t)/G.order())*dims[r]*sum(exp(eigenvalue*CDF(t)) * traces[r][instance][eigenvalue] for eigenvalue in eig_lists[r])
+        return real(ans)
+    mle_dict = {}
+    bound = (0, 35)
+    for instance in genome_reps:
+        def L(time):
+            return -likelihood(time, instance)
+        mle = minimize(L, 0.1, method="TNC", bounds=[bound])
+        mle_dict[instance] = mle.x[0]
+        if plots:
+            plt.figure()
+            times = np.arange(bound[0], bound[1], 0.5)
+            likelihood_functional_values = [likelihood(time, instance) for time in times]
+            plt.plot(times, likelihood_functional_values)
+            plt.savefig(f'./cgt_{framework.one_row(instance)}')
+
+    return {framework.one_row(key) : value for key, value in mle_dict.items()}
+
+def _projection_operators(mat, eigs):
+        """Return projection operators for given matrix and its eigenvalues"""
+        dim = mat.nrows()
+        projections = [matrix.identity(dim) for _ in eigs]
+        for e1, eig1 in enumerate(eigs):
+            for eig2 in eigs:
+                if eig1 != eig2:
+                    projections[e1] *= (mat-(eig2*np.eye(dim)))*(1/(eig1-eig2))
+        return projections
+
+def _eigenvalues(mat, round_to=9, make_real=True, inc_repeated=False):
+    col = list if inc_repeated else set
+    collection = col(round(real(eig) if make_real else eig, round_to) for eig in mat.eigenvalues())
+    return sorted(collection)
+
+def likelihood_function():
+    """Return the likelihood function for a given genome"""
     pass
 
-def min_distance(framework, model, genomes, weighted=False):
+def min_distance(framework, model, weighted=False):
     """Return dictionary of minimum distances ref->genome, using inverse of model probabilities as weights (or not)"""
-    pass
+    matrix = model.reg_rep_of_zs().toarray()
+    if weighted:
+        matrix = np.reciprocal(matrix, where=matrix!=0)
+    graph = nx.Graph(matrix)
+    genome_reps = framework.genomes().keys()
+    return { 
+        rep : nx.shortest_path_length(graph, source=0, target=r, weight='weighted' if weighted else None) 
+        for r, rep in enumerate(genome_reps) 
+    }
 
-def MFPT(framework, model, genomes, mean_inter_arrival_time=1):
+def MFPT(framework, model, scale_by=1):
     """Return the mean time elapsed rearranging ref->genome where the target is an absorbing state"""
-    pass
+    genomes = framework.genomes()
+    P = model.reg_rep_of_zs().toarray()
+    P_0 = P.copy()
+    for i in range(len(P)):
+        P_0[0,i] = 0 
+    Z = np.linalg.inv(np.eye(len(P)) - P_0)
+    def MFTP_dists():
+        return (Z@P_0@Z)[:,0]
+    reps = list(genomes.keys())
+    MFTP_distances = list(MFTP_dists())
+    MFTP_distances = { rep : scale_by * MFTP_distances[r] for r, rep in enumerate(reps) }
+    return MFTP_distances

cgt/enums.py

@@ -22,17 +22,27 @@ class FORMAT(Enum):
     only_reps      = auto()
 
 class SET(Enum):
-	signed_cycles 	= auto()
-	unsigned_cycles = auto()
-	one_row			= auto()
-	wreath			= auto()
-	wreath_s2		= auto()
-	
+    signed_cycles 	= auto()
+    unsigned_cycles = auto()
+    one_row			= auto()
+    wreath			= auto()
+    wreath_s2		= auto()
+    
 class SYMMETRY(Enum):
     circular = dihedral = D_n              = auto()
     linear = S_2 = C_2 = flip = reflection = auto()
 
 class CLASSES(Enum):
     conjugacy_classes = auto()
     cosets            = auto()
-    double_cosets     = auto()
+    double_cosets     = auto()
+
+class DISPLAY(Enum):
+    arrows  = auto()
+    one_row = auto()
+    cycles  = auto()
+
+class ALGEBRA(Enum):
+    group        = auto()
+    genome       = auto()
+    genome_class = auto()

cgt/examples.py

@@ -0,0 +1,14 @@
+"""
+"""
+from sage.modular.modform.element import ModularFormElement
+from .enums import MODEL, SYMMETRY
+from .models import Model
+from .position_paradigm import PositionParadigmFramework
+
+def example(n):
+    framework = PositionParadigmFramework(n, symmetry=SYMMETRY.circular)
+    model = Model.named_model_with_relative_probs(framework, {
+        MODEL.one_region_inversions: 0.5, 
+        MODEL.two_region_inversions: 0.5
+    })
+    return framework, model

cgt/models.py

@@ -1,60 +1,72 @@
 """
 """
 
-from enum import Enum, auto
 from sage.all_cmdline import *
 from .enums import *
 from . import rearrangements
-from . import structures
 
 class Model:
-	"""Defines a model. A model consists of some collection of permutations and a map from these permutations to probabilities [0,1]"""
-	def __init__(self, framework, generating_dictionary):
-		"""Define a model from a dictionary of single permutations, with their probabilities as the values."""
-		self.framework = framework
-		self.generating_dictionary = generating_dictionary
-		# TODO: Implement checks for certain model properties, for example time reversibility, symmetry and missing rearrangements
+    """Defines a model. A model consists of some collection of permutations and a map from these permutations to probabilities [0,1]"""
+    def __init__(self, framework, generating_dictionary):
+        """Define a model from a dictionary of single permutations, with their probabilities as the values."""
+        self.framework = framework
+        self.generating_dictionary = generating_dictionary
+        self.s = {} # model elements
+        # TODO: Implement checks for certain model properties, for example time reversibility, symmetry and missing rearrangements
 
-	@classmethod
-	def named_model_with_relative_probs(cls, framework, named_model_dictionary):
-		if not framework.oriented:
-			raise NotImplementedError(f"not yet implemented for {str(framework)}")
-		# look through all specified model names and handle them appropriately
-		model = {}
-		if abs(sum(named_model_dictionary.values()) - 1) > 0.00001: 
-			raise ValueError("supplied probabilities do not sum to 1")
-		for model_name, relative_prob in named_model_dictionary.items():
-			if model_name is MODEL.one_region_inversions:
-				gens = rearrangements.all_inversions_representatives(framework, num_regions=1)
-				for generator in gens:
-					model[generator] = relative_prob/len(gens)
-			if model_name is MODEL.two_region_inversions:
-				gens = rearrangements.all_inversions_representatives(framework, num_regions=2)
-				for generator in gens:
-					model[generator] = relative_prob/len(gens)
-			if model_name is MODEL.all_inversions:
-				gens = rearrangements.all_inversions_representatives(framework)
-				for generator in gens:
-					model[generator] = relative_prob/len(gens)
-		return cls(framework, model)
+    @classmethod
+    def named_model_with_relative_probs(cls, framework, named_model_dictionary):
+        if not framework.oriented:
+            raise NotImplementedError(f"not yet implemented for {str(framework)}")
+        # look through all specified model names and handle them appropriately
+        model = {}
+        if abs(sum(named_model_dictionary.values()) - 1) > 0.00001: 
+            raise ValueError("supplied probabilities do not sum to 1")
+        for model_name, relative_prob in named_model_dictionary.items():
+            if model_name is MODEL.one_region_inversions:
+                gens = rearrangements.all_inversions_representatives(framework, num_regions=1)
+                for generator in gens:
+                    model[generator] = relative_prob/len(gens)
+            if model_name is MODEL.two_region_inversions:
+                gens = rearrangements.all_inversions_representatives(framework, num_regions=2)
+                for generator in gens:
+                    model[generator] = relative_prob/len(gens)
+            if model_name is MODEL.all_inversions:
+                gens = rearrangements.all_inversions_representatives(framework)
+                for generator in gens:
+                    model[generator] = relative_prob/len(gens)
+        return cls(framework, model)
 
-	def generators(self):
-		return self.generating_dictionary
+    def reg_rep_of_zs(self):
+        """Return the regular representation of zs as comptued by PositionParadigmFramework.reg_rep_zs, but store the sparse result"""
+        try:
+            return self._reg_rep_of_zs
+        except AttributeError:
+            self._reg_rep_of_zs = self.framework.reg_rep_of_zs(self, sparse=True)
+            return self._reg_rep_of_zs
 
-	def element(self):
-		try:
-			return self.s
-		except AttributeError:
-			A = self.framework.group_algebra()
-			s = A(0) # Zero element in algebra
-			gens_dict = self.generators()
-			gens = {x for x in self.generators().keys()}
-			while len(gens) > 0:
-				gen = gens.pop()
-				prob = gens_dict[gen]
-				conj_class = rearrangements.conjugacy_class(self.framework, gen)
-				for perm in conj_class:
-					s += (prob/len(conj_class)) * A(perm)
-				gens -= conj_class
-			self.s = s
-			return s
+    def generators(self):
+        return self.generating_dictionary
+
+    def s_element(self, in_algebra=ALGEBRA.group):
+        if in_algebra not in {ALGEBRA.group, ALGEBRA.genome}:
+            raise NotImplementedError(f"Model element for {str(in_algebra)} algebra not yet implemented")
+        try:
+            return self.s[in_algebra]
+        except (AttributeError, KeyError):
+            A = self.framework.group_algebra()
+            s = A(0) # Zero element in algebra
+            gens_dict = self.generators()
+            gens = {x for x in self.generators().keys()}
+            while len(gens) > 0:
+                gen = gens.pop()
+                prob = gens_dict[gen]
+                if in_algebra is ALGEBRA.group:
+                    conj_class = rearrangements.conjugacy_class(self.framework, gen)
+                elif in_algebra is ALGEBRA.genome:
+                    conj_class = {gen}
+                for perm in conj_class:
+                    s += QQ(prob/len(conj_class)) * A(perm)
+                gens -= conj_class
+            self.s[in_algebra] = s
+            return s