Adding Graphical Models in Experimental

experimental/GraphicalModels/README.md

@@ -0,0 +1,13 @@
+# An example template for the experimental section
+
+## Aims
+
+This is an example for a file structure to set up a new package 
+in the experimental section. All files you find here are part of the 
+minimum requirements. See also the official Oscar documentation.
+
+## Status
+
+We plan to also provide a function to automatically copy this template 
+for you to start your own package. 
+

experimental/GraphicalModels/docs/doc.main

@@ -0,0 +1,5 @@
+[
+   "GraphicalModels" => [
+      "graphicalmodels.md",
+   ]
+]

experimental/GraphicalModels/docs/src/graphicalmodels.md

@@ -0,0 +1,67 @@
+```@meta
+CurrentModule = Oscar
+DocTestSetup = quote
+  using Oscar
+end
+```
+# Algebraic Phylogenetics
+
+The Algebraic Phylogenetics part of OSCAR provides functionality for
+- specifying phylogenetic models
+- parametrizing such models
+- calculating their algebraic invariants
+
+## Models
+
+The five most common models in algebraic phylogenetics can automatically be specified by calling the below functions, each taking a tree as input and attaching transition matrices to its edges as defined by Jukes Cantor, Kimura, etc., respectively, returning a structure `PhylogeneticModel` or `GroupBasedPhylogeneticModel`.
+
+```@docs
+cavender_farris_neyman_model(graph::Graph{Directed})
+jukes_cantor_model(graph::Graph{Directed})
+kimura2_model(graph::Graph{Directed})
+kimura3_model(graph::Graph{Directed})
+general_markov_model(graph::Graph{Directed})
+```
+
+The models are by default in projective space. For their affine versions, call
+
+```@docs
+affine_phylogenetic_model!(pm::PhylogeneticModel)
+```
+
+## Components of a model
+
+`PhylogeneticModel` specifies any phylogenetic tree model in probability coordinates and `GroupBasedPhylogeneticModel` can specify group-based, e.g. Fourier, coordinates. For any model, we can call its graph, transition matrices attached to the graph's edges, the number of states of each vertex random variable, and the corresponding polynomial rings. For instance for Jukes Cantor on the star with three leaves:
+
+```@docs
+graph(pm::PhylogeneticModel)
+transition_matrices(pm::PhylogeneticModel)
+number_states(pm::PhylogeneticModel)
+probability_ring(pm::PhylogeneticModel)
+fourier_ring(pm::GroupBasedPhylogeneticModel)
+fourier_parameters(pm::GroupBasedPhylogeneticModel)
+group_of_model(pm::GroupBasedPhylogeneticModel)
+```
+
+
+## Parametrization
+
+For each phylogenetic model, we can calculate the parametrization, a map from transition matrices to probabilities, parametrized in probability or Fourier coordinates. For group-based models, we can reparametrize between these and return the transformation matrix, and we can calculate equivalence classes of probabilities with the same parametrization.
+
+```@docs
+probability_map(pm::PhylogeneticModel)
+fourier_map(pm::GroupBasedPhylogeneticModel)
+compute_equivalent_classes(parametrization::Dict{Tuple{Vararg{Int64}}, QQMPolyRingElem})
+sum_equivalent_classes(equivalent_classes::Dict{Vector{Tuple{Vararg{Int64}}}, QQMPolyRingElem})
+specialized_fourier_transform(pm::GroupBasedPhylogeneticModel, p_equivclasses::Dict{Vector{Tuple{Vararg{Int64}}}, QQMPolyRingElem}, f_equivclasses::Dict{Vector{Tuple{Vararg{Int64}}}, QQMPolyRingElem})
+inverse_specialized_fourier_transform(pm::GroupBasedPhylogeneticModel, p_equivclasses::Dict{Vector{Tuple{Vararg{Int64}}}, QQMPolyRingElem}, f_equivclasses::Dict{Vector{Tuple{Vararg{Int64}}}, QQMPolyRingElem})
+```
+
+## Contact
+
+Please direct questions about this part of OSCAR to the following people:
+* [Marina Garrote López](https://sites.google.com/view/marinagarrotelopez),
+* possibly others.
+
+You can ask questions in the [OSCAR Slack](https://www.oscar-system.org/community/#slack).
+Alternatively, you can [raise an issue on github](https://www.oscar-system.org/community/#how-to-report-issues).

experimental/GraphicalModels/src/GraphicalModels.jl

@@ -0,0 +1,38 @@
+include("PhylogeneticModels.jl")
+include("PhylogeneticAuxiliary.jl")
+include("PhylogeneticParametrization.jl")
+
+#export models
+export cavender_farris_neyman_model 
+export jukes_cantor_model
+export kimura2_model
+export kimura3_model
+export general_markov_model
+export affine_phylogenetic_model!
+
+#export phylogenetic models attributes
+export phylogenetic_model
+export graph
+export number_states
+export transition_matrices
+export probability_ring
+export root_distribution
+export fourier_parameters
+export fourier_ring
+export group_of_model
+
+#export probability and fourier map
+export probability_map
+export fourier_map
+
+# export functions to calculate equivalent classes
+export compute_equivalent_classes 
+export sum_equivalent_classes 
+
+# export transformation matrices
+export specialized_fourier_transform
+export inverse_specialized_fourier_transform
+
+# export structs for GroupBasedPhylogeneticModel,PhylogeneticModel
+export PhylogeneticModel
+export GroupBasedPhylogeneticModel

experimental/GraphicalModels/src/PhylogeneticAuxiliary.jl

@@ -0,0 +1,89 @@
+##########################
+#### GROUP OPERATIONS ####
+##########################
+# This function implements the group operation on Z2 x Z2 and is subsequently 
+# used to check the condition given in the computation of the Fourier
+# coordinates and to calculate the sum of the leaves beneath a given edge. 
+
+function group_sum(pm::GroupBasedPhylogeneticModel, states::Vector{Int})
+  group = group_of_model(pm)
+  return sum(group[states]).%2
+end
+
+function is_zero_group_sum(pm::GroupBasedPhylogeneticModel, states::Vector{Int})
+  ng = length(states)
+  return group_sum(pm, [states[1]]) == group_sum(pm, states[2:ng])
+end
+
+function which_group_element(pm::GroupBasedPhylogeneticModel, elem::Vector{Int64})
+  group = group_of_model(pm)
+  return findall([all(group[i].==elem) for i in 1:length(group)])[1]
+end
+
+
+########################################
+#### AUXILIARY FUNCTIONS FOR GRAPHS ####
+########################################
+
+function interior_nodes(graph::Graph)
+  big_graph = Polymake.graph.Graph(ADJACENCY = pm_object(graph))
+  degrees = big_graph.NODE_DEGREES
+  return findall(x -> x > 1, degrees)
+end
+
+function leaves(graph::Graph)
+  big_graph = Polymake.graph.Graph(ADJACENCY = pm_object(graph))
+  degrees = big_graph.NODE_DEGREES
+  return findall(x -> x == 1, degrees)
+end
+
+function vertex_descendants(v::Int, gr::Graph, desc::Vector{Any})
+  lvs = leaves(gr)
+  outn = outneighbors(gr, v)
+
+  if v in lvs
+    return [v]
+  end
+
+  innodes = setdiff(outn, lvs)
+  d = unique(append!(desc, intersect(outn, lvs)))
+
+  if length(innodes) > 0
+    for i in innodes
+      d = vertex_descendants(i, gr, d)
+    end
+    return d 
+  end
+
+  return d
+end
+
+function cherries(graph::Graph)
+  lvs = leaves(graph)
+  cherry = []
+  for l in lvs
+    in_node = inneighbors(graph,l)[1]
+    lvs_cherr = outneighbors(graph, inneighbors(graph,l)[1])
+    if issubset(lvs_cherr, lvs) == 2
+      cherry = append!(cherry, [[Edge(in_node, lvs_cherr[1]), Edge(in_node, lvs_cherr[2])]])
+    end
+  end
+
+  return unique(cherry)
+end
+
+function root(graph::Graph)
+  n_parents = [length(inneighbors(graph, v)) for v in 1:n_vertices(graph)]
+  return findall(x -> x == 0, n_parents)[1]
+end
+
+
+########################
+#### SORT FUNCTIONS ####
+########################
+
+function sort_edges(graph::Graph)
+  edgs = collect(edges(graph))
+  leaves_idx = findall(edge -> dst(edge) in Oscar.leaves(graph), edgs)
+  return edgs[vcat(leaves_idx, setdiff(1:length(edgs), leaves_idx))]
+end