[WIP] Branch length estimation #50
Conversation
…ngths to conform to API
| import copy | ||
| from typing import List, Tuple | ||
|
|
||
| import cvxpy as cp |
There was a problem hiding this comment.
We'll need to add this to the setup.py, my tests break because this wasn't installed during make install
There was a problem hiding this comment.
Agreed!
QQ: I noticed many of the requirements have a minimum version requirement, e.g. numpy > 1.17, is there a reason for that? I was thinking to just add cvxpy without a version requirement.
There was a problem hiding this comment.
Having version requirements is helpful for a couple of reasons -
(1) Sometimes, APIs are different between versions (e.g., networkx changes the way it calls some functionality in versions <2 and >2, I believe)
(2) Sometimes, earlier versions of packages throw warnings (this is especially relevant if, for example, a particular package relies on another package that has updated and is deprecating some feature). Requiring a minimal version (or maximal version, for that matter) may help you get rid of warnings being thrown around.
In general, it's likely that having an outdated version of a particular package will be fine and the worst thing that will happen is that a warning might be thrown. However, I think it's good for software reproducibility so that when a new user comes along they can create a very similar environment to the developer. Hope this makes sense!
| log_likelihood: The log-likelihood of the training data under the | ||
| estimated model. | ||
| log_loss: The log-loss of the training data under the estimated model. | ||
| This is the log likhelihood plus the regularization terms. |
There was a problem hiding this comment.
typo -> "likelihood"
|
|
||
| def estimate_branch_lengths(self, tree: Tree) -> None: | ||
| r""" | ||
| See base class. The only caveat is that this method raises if it fails |
There was a problem hiding this comment.
incomplete sentence: "this method raises [an error?]..."
There was a problem hiding this comment.
I've always used "raises" for short to mean that it "raises an error", but it might not be proper english.
There was a problem hiding this comment.
Ah - you might have more software engineering lingo than I do! If it's not too much of a problem, let's be more verbose here in the docstring.
Can you also add a Raises section to the docstring (ala Args and Returns)?
| def dfs(node: int, tree: Tree): | ||
| node_state = tree.get_state(node) | ||
| for child in tree.children(node): | ||
| # Compute the state of the child | ||
| child_state = "" | ||
| edge_length = tree.get_age(node) - tree.get_age(child) | ||
| # print(f"{node} -> {child}, length {edge_length}") | ||
| assert edge_length >= 0 | ||
| for i in range(num_characters): | ||
| # See what happens to character i | ||
| if node_state[i] != "0": | ||
| # The character has already mutated; there in nothing | ||
| # to do | ||
| child_state += node_state[i] | ||
| continue | ||
| else: | ||
| # Determine if the character will mutate. | ||
| mutates = ( | ||
| np.random.exponential(1.0 / mutation_rate) | ||
| < edge_length | ||
| ) | ||
| if mutates: | ||
| child_state += "1" | ||
| else: | ||
| child_state += "0" | ||
| tree.set_state(child, child_state) | ||
| dfs(child, tree) |
There was a problem hiding this comment.
A general comment - a depth-first traversal function already exists in networkx.
But, if you want to do this a special way, I think this belongs as a method for the Tree class, not as a nested function here. Thoughts?
There was a problem hiding this comment.
I should look into the function from networkx that you mention. What exactly does it do? Can you please point me to it?
We need to do two things in overlay_lineage_tracing_data: (1) iterate over all edges down the tree (the tree traversal) (2) do some custom computation on each edge as we visit it. If we were to any of this behavior to the tree class I guess is would be (1)? IMO (2) should rest in the LineageTracingSimulator since it is something that varies.
There was a problem hiding this comment.
Yeah, that's a good point. The action items at each internal node in the DF traversal will change so it's not necessarily a general method for the tree class.
However, having a DFS method in the class that just returns nodes in post-order might be worthwhile and could cut down on some duplicated code in the other classes that operate on the CassiopeiaTree object. I'll add this a todo in the class.
| np.random.exponential(1.0 / mutation_rate) | ||
| < edge_length |
There was a problem hiding this comment.
I could be misremembering, but I believe np.random.exponential already inverts the scale parameter (i.e. np.random.exponential(b) draws from 1/b * exp(x/b)). [You can read the docs here: https://numpy.org/doc/stable/reference/random/generated/numpy.random.exponential.html -- indeed, it looks like the function inverts the scale parameter for you]
If this is the case, do you want to pass instead use np.random.exponential(mutation_rate)?
There was a problem hiding this comment.
np.random.exponential uses the scale parameterization, so we should invert the rate as in np.random.exponential(1.0 / mutation_rate), right?
There was a problem hiding this comment.
Sorry - I'm a bit confused. If the density function you want is lambda * exp(x*lambda) (where lambda is the mutation rate) then you are good to go. O.w. you should just pass in lambda to the exponential density function.
cassiopeia/tools/tree.py
Outdated
| def reconstruct_ancestral_states(self): | ||
| r""" | ||
| Reconstructs ancestral states with maximum parsimony. | ||
| """ | ||
| root = self.root() | ||
|
|
||
| def dfs(v: int) -> None: | ||
| children = self.children(v) | ||
| n_children = len(children) | ||
| if n_children == 0: | ||
| return | ||
| for child in children: | ||
| dfs(child) | ||
| children_states = [self.get_state(child) for child in children] | ||
| n_characters = len(children_states[0]) | ||
| state = "" | ||
| for character_id in range(n_characters): | ||
| states_for_this_character = set( | ||
| [ | ||
| children_states[i][character_id] | ||
| for i in range(n_children) | ||
| ] | ||
| ) | ||
| if len(states_for_this_character) == 1: | ||
| state += states_for_this_character.pop() | ||
| else: | ||
| state += "0" | ||
| self.set_state(v, state) | ||
| if v == root: | ||
| # Reset state to all zeros! | ||
| self.set_state(v, "0" * n_characters) | ||
|
|
||
| dfs(root) |
There was a problem hiding this comment.
Richard wrote a function that does this in cassiopeia.solver.solver_utilities. You might just want to take that because it takes into account missing data, etc.
There was a problem hiding this comment.
Yes, you are right! I just found it easier to code this up myself then to figure out how to get the right inputs to annotate_ancestral_characters, because e.g. the Tree class I am using holds character vectors as a str, not as a List[int] :( Probably a bad design choice of mine? This is one thing I wanted to discuss! I left the character vector as a str because it met all my needs, and changing it to List[int] meant changing the code prematurely. (We definitely want to have only one implementation of parsimony reconstruction when we merge, the question is what's the API/where it lives.)
There was a problem hiding this comment.
Makes sense! This is definitely something out of scope of your branch length estimator.
But, yes I feel that using lists to represent the character states of a sample is more desirable. One major reason for this is the separation of characters in a string is not always intuitive (e.g., if you have double-digit or triple-digit state representations). One can get around this by setting a delimiter (e.g., the | symbol) but if you always end up splitting the string to change something and then rejoining, why not just ditch the string data structure altogether for this?
… cuts_p <= x <= cuts_v
I migrated my branch length estimation code here, under
cassiopeia/tools. This PR should be see as a big draft to get discussions started around APIs. For example, note that my APIs currently depend on aTreeclass, rather than onnetworkx.DiGraph.Note also that my code is totally self-contained right now, segregated from any utilities that Cassiopeia might already provide, such as ancestral state reconstruction via maximum parsimony (I have a method for this in the
Treeclass which I just implemented de-novo).