Version 2 Roadmap

[schlegelp]

This issue is to suggest & discuss (read: spitball) potential breaking changes that might warrant a new major version.

Plotting

Over the past months I have experimented more with the pygfx library and I'm pretty much set on replacing the current vispy-based 3d viewer with octarine. This may well happen already in the next minor version (see #138) but I wonder if a general refactoring of the plotting interface is in order. In particular:

• Do we want to add pygfx as a 2d plotting backend (see also fastplotlib)? It's much faster and perspectively correct (unlike matplotlib).
• In general, I think the plotting should be refactored into classes under the hood instead of the current convoluted functions.
• Further to the above: I would love to play around trying a "Grammar of Graphics"-like interface e.g. via method chaining.

Additional backends

Currently, navis can utilize multiple CPUs (via pathos) but it doesn't scale beyond a single machine. It'd be nice if there was an easy way to divvy up e.g. a big NBLAST across multiple nodes. @aschampion has already made a start in #111 for a Dask backend and the implementation is generic so that it would be easy to tack on more.

Neuron Classes

Two things that have been on my mind with regards to neurons:

1. As discussed with @bdpedigo in Berlin, can we generate a "unified" neuron object that combines multiple representations (skeleton, mesh, image)? @ceesem has done something similar in MeshParty.

2. The current neuron names (TreeNeuron, MeshNeuron, VoxelNeuron) are... ugly

It would be nice if there was only a single Neuron class that could be either a single type or multiple. Something along these lines:

>>> swc = pd.read_table('neuron.swc')  # load some SWC table
>>> n = navis.Neuron(swc)  # automatically interprets that this is a skeleton
>>> # Alternatively: n = navis.Neuron.from_swc('neuron.swc') or n = navis.read_swc('neuron.swc')
>>> n
<navis.Neuron(id=12345, type=skeleton, nodes=4320)>
>>> n.type
("skeleton", )

Combining multiple representations could then look like this:

>>> n = navis.Neuron.from_swc('neuron.swc')
>>> n.add_mesh(mesh)
>>> Alternatively: n = navis.Neuron(skeleton=swc, mesh=mesh)
>>> n 
<navis.Neuron(id=12345, type=(skeleton, mesh), nodes=4320, faces=20432)>
>>> n.type
("skeleton", "mesh")

Individual representations can be accessed like this:

>>> n.skeleton  # this can also be generated on-the-fly from the mesh
<navis.Skeleton(id=12345, faces=20432)>
>>> n.mesh
<navis.Mesh(id=12345, faces=20432)>
>>> navis.plot3d(n)  # plots the "first" representation (what that is tbd)
>>> navis.plot3d(n.skeleton)  # explicitly plot the skeleton

Under-the-hood navis already combines representations when needed. For example, navis.prune_twigs works on MeshNeurons by (a) skeletonizing the mesh, (b) pruning the skeleton and (c) removing faces corresponding to the removed skeleton nodes from the mesh. The big question then is: does making this explicity actually simplify things for the user? I imagine that most users will work with either skeletons or meshes and don't need to know how things happen under-the-hood.

I will add above list as I think of more stuff.

Perhaps some of the above can be tackled in a longer hackathon?

[bdpedigo]

big +1!

this kind of neuron.mesh, neuron.skeleton is exactly what I had in mind.

Some challenges I could see:

• defining what happens for all representations when you do an operation on another. your prune_twigs is a perfect example where you've figured out how to define that operation mapping, but I am wondering if it is easy to generalize to all cases/operations or not
• same as the above, but whether to map operations from one representation to another always, or in some kind of lazy manner
• as you mention, how to hide the details of these backend implementation details when the user doesnt care about the underlying representation, while also allowing the customizability for people who do

[schlegelp]

The way this is currently implemented is via a simple decorator:

@utils.meshneuron_skeleton(method='subset')
def prune_twigs(skeleton, ...):

The decorator intercepts the input to prune_twigs and if it's a mesh it will take its skeleton (generate if not yet present) and pass that through to the function. The function does its thing on the skeleton and when the decorator gets the pruned skeleton back, it applies the changes to the mesh. Bottom line: you write your function to operate on skeletons only and let the decorator deal with meshes.

Your idea for lazy mapping between skeleton and mesh is cool. I can imagine wanting to chain a couple functions that work on the skeleton - in that case, I don't need to re-generate the mesh in between.

>>> # These functions just edit the skeleton and set a flag that the mesh needs to be regenerated from the skeleton 
>>> navis.prune_twigs(n, 5, inplace=True)  
>>> navis.prune_by_strahler(n, 1, inplace=True)
>>> # Only now is the mesh actually regenerated 
>>> n.mesh

The lazy evaluation would be super easy to implement:

class Neuron:
...
@property 
def mesh(self):
   if self._mesh_stale:
      self._regenerate_mesh_from_skeleton()
   return self._mesh

[bdpedigo]

yeah that kind of pattern sounds wonderful to use, i worried it would be hard to keep track of things if multiple steps are chained together but maybe if there is some robust code for keeping track of the mapping between representations it would not be so bad

[bdpedigo]

re: visualization, i've been using https://github.com/pyvista/pyvista a lot lately - i wonder what the pros/cons are relative to pygfx. one thing that I appreciate about pyvista is that it seems to have a big user base, which helps with googling for solutions.

[schlegelp]

I remember looking at PyVista way back when I added the 3D viewer to navis. A bit fuzzy on why I decided to go for vispy instead of it but it probably had something to with the clunky VTK dependency/interface. To be fair: that was a while ago and I probably didn't try very hard.

In any event: pygfx is super light-weight, blazingly fast (based on modern WGPU instead of OpenGL) and has a beautiful API that is just the right spot between low- and high-level interface to be super customisable without too much effort. VTK on the other hand has the advantage of coming with "batteries included" (i.e. it has things like mesh-simplification built in) and having a bunch of bells and whistles that pygfx does not (yet).

[bdpedigo]

re: representations:

• non-tree spatial graphs: for example the level2 graph used in CAVE, which isn't a DAG, but also isn't exactly a mesh as we normally think of it. however meshparty actually just treats this spatial graph as a mesh in terms of the code, though, so it might be easy to support this representation too
• condensed trees / segment graphs: taking the skeleton but converting to where each node is an entire segment (branch to branch or branch to tip) in the skeleton. ive used these for visualizations and can be helpful for some computations too

[schlegelp]

Re non-tree spatial graphs: minimally, the spatial graph could be used to fix breaks the skeleton representation. That should be fairly easy to implement but I can see how having the spatial graph would be cool to have elsewhere too.

RE condensed trees / segments graphs: there is already navis.TreeNeuron.simple that will give you that. It's used under-the-hood for some graph traversal stuff but using it for visualizations hadn't crossed my mind. Curious to see what you did with it!

[bdpedigo]

But yes, would be a big fan of a morphology representations hackathon! Perhaps as a way to get started, it could be worth writing down a kind of desired API, and working backwards from there in terms of how to implement things?

[schlegelp]

💯 I think writing down some example API and small workflows would be super helpful.

I might sit down, take some of the existing navis tutorials and mock up how they could look like in v2. The Neuron Basics tutorial would be an interesting exercise to walk through as it supposed to explain to the user how neurons work.

[ceesem]

Unsurprisingly, Ben and I are on the same page here with our goals.

A few things that I have come to believe are useful generalizations/distinctions for such a project:

• Having a generic category of "SpatialGraph" is useful. That is to say, a non-skeleton structure with vertices in space and edges between them. This could be the mesh, this could be the level 2 graph, or any other wireframe of a neuron. The key feature is that it's a single connected component and that any annotation you might want to is can be associated with a vertex of this graph, but it's not necessarily tree-like.
• It is useful to distinguish point-like annotations (e.g. synapses, mitochondria) that can have a any-number-to-one relationship with vertices from labels, which are necessarily in 1:1 correspondence with vertices. Labels could be anything from neuronal compartment to radius to SegCLR embeddings.

[schlegelp]

I'm on-board with both of these and I understand the power of a spatial graph from the L2 graphs in FlyWire - although on a side note: I did a bit of experimenting trying to fix breaks in full-res skeletons using the L2 graph but that turned out to be a mixed bag because users don't always click on the neighbouring supervoxels when merging.

For annotations vs labels: would that be like connectors vs tags in CATMAID?

[ceesem]

In fact, to not just totally repeat what @bdpedigo said while I was already typing, a key distinction between the mesh and a spatial graph is that the mesh is built of faces and is typically not a single connected component. We always have to remove disconnected faces and graft in connecting edges in order to build a "nice" single component graph out of it, which already makes it not a true triangle mesh.

[schlegelp]

I totally see how the spatial graph can be very powerful. In my head, I'm already thinking about having a SpatialGraph class that does the heavy lifting of generating, fixing and syncing between mesh and skeleton representation.

[bdpedigo]

Having a generic category of "SpatialGraph" is useful.

Exactly, and now a TreeNeuron/skeleton can inherit from spatial graph, but could have extra methods specific to trees (e.g. a notion of tips and a root)

[ceesem]

And because there’s a mapping between graph and tree vertices, the tree concept of “tip” or “root” maps onto one or more graph vertices.

[schlegelp]

FYI: I converted this to a discussion since it allows directly responding to a post.

[schlegelp]

General question for spatial graphs: is there a reason why the spatial graph can't just be the mesh's skeleton, fixed up and without breaks? Or the other way around: is there an advantage of having the spatial graph consist of the mesh's vertices and edges plus a set of extra edges that make it a single connected component?

[bdpedigo]

but to your more general point I go back and forth about whether the mesh should be in this category or not. im torn between leveraging existing code to operate on meshes like trimesh or VTK/pyvista or {name your mesh package} vs writing something that is a subset of this spatial graph, or maybe some hybrid of the two

[ceesem]

The main reason is that you want to have around the most complete representation of the neuron available while retaining provenance of the various data sources. For example, if you want to analyze the local morphology of the structure around synapses or even the nanometer-scale distances between inputs, synapses live most naturally with either the surface mesh or even a voxel graph. However, once you start getting to branch level statistics, you want to be using the skeleton. We have definitely seen opportunity in our work to have both. Maybe it’s over-engineering to have a method that lets you do both, but I think most of the algorithms and data structures are similar. If you’re building the skeleton from the spatial graph, for example, you already need the graph to skeleton mapping to put e.g. synapses in the correct spot, so why not keep that around for bidirectional analyses?

[bdpedigo]

that all makes sense but im still missing how that implies one connected component

[schlegelp]

I guess we expect a real neuron to be a single connected component so the spatial graph should (ideally) reflect that?

For what it's worth: I didn't think we would enforce connectedness.

[aschampion]

Even without "noisy" multi-component neurons from meshes, there's always been a case for multi-component neurons due to truncated volumes or reconstruction. It would be good to factor explicit interfaces for what works on multi-component graphs (label selection, morphometrics, etc.) vs single-components (global synaptic flows, etc.). For your magic default Neuron type, then, just implement both, and have any interface methods requiring single-component first check (with some cached dirty flag) that it's single-component, with any mutations resetting the flag. Users with non-arborescence spatial graphs can then just skip implementing that interface if it doesn't suit the data, and other parts of navis can clarify with hints/docs what capabilities they actually require.

[clbarnes]

My utility for and ability to contribute to navis has severely diminished, but my two cents about the "single neuron representation" conversation is that there is plenty of functionality which only applies to certain representations. Given that, it is nice to be able to rely on the type system to enforce what kind of representation you're working with rather than having a mixed bag and try to keep track yourself. That would mean separate types for different representations, and any given function could Union[Rep1, Rep2] for all the types it works with.

You could, for convenience, additionally have a joint representation which composes over them all and has some basic functionality which can call on any of the representations (e.g. a bounding box which prefers meshes, then trees, then point clouds) but I'd be quite cautious about over-using it when it comes to type annotations.

[schlegelp]

The sanity gains from a codebase which can be typed end-to-end are huge.

I remember trying that when type hinting and mypy came out and concluding that it wasn't worth the effort for complex data types. Back then, there weren't even any stubs for numpy or pandas - I'm sure they exist now - and making mypy 100% happy was just impossible. I'm a big fan of type safety in Rust but for navis I always thought that a decent number of doc tests and assertions is the way to go.

Being able to distinguish between neuron representations and avoid automagic conversions between them is, IMO, very important for scaling - skeletons are a very efficient representation for certain analyses and people (/CPU clusters) dealing with a large number of neurons may get upset if their 100,000 skeletons all got secretly turned into voxels.

Very true. I had already been thinking about ways to make sure a function only gets the data it actually needs. For example, we wouldn't want to pass more than the dotprops into an NBLAST. The same applies for other functions where we'd need to avoid passing more data than necessary to the child processes when parallel=True.

I imagine that in practice automagic conversion will only ever be simplification. At least with the current set of (non-plotting) functions, we'd always end up with a skeleton. Still, you never know what the users might try to do...

Re your UnionNeuron suggestion: that would be the explicit option, similar to @ceesem's meshwork. Perhaps it'd be worth implementing something like this while still in version 1.x and see how it goes?

Might even be cool to play with it in downstream libraries: for example, I was thinking that we could have some on-demand-data fetching for FlyWire neurons in fafbseg (similar to CatmaidNeurons in pymaid):

class FlywireNeuron(navis.UnionNeuron):
    def __init__(self, root_id):
        self.root_id = root_id
        self._mesh = None 
        self._skeleton = None 
    
    @property
    def mesh(self) -> Mesh:
        if self._mesh is None:
            self._mesh = fafbseg.flywire.get_mesh_neuron(self.root_id)
        return self._mesh

    @property 
    def skeleton(self) -> Skeleton:
        if self._skeleton is None:
            self._skeleton = fafbseg.flywire.get_skeleton(self.root_id)
        return self._skeleton

[aschampion]

While I agree with Chris re: having some interface to make decision points regarding representation explicit, and Philipp re: having an ergonomic wrapper implicitly providing the most common conversions, the trick is as always what happens once there's one layer of indirection. That is, when these representations are being used in a library-provided analysis algorithm or dynamic resampling for vis., etc., how do you expose policy about representation mixing without a mess of either being stuck with the default wrapper, verbose interface, or hidden mixed-implementation bugs. Either the implementations have to be very careful to only use caller-supplied types and type constructors (which even first-party implementations struggle with, look at Pandas' awful extended dataframe attributes or subclass support), or you have to do some kind of dependency injection which is too much Engineer for a Biologist library like navis.

Some nits in lieu of a detailed but important implementation discussion to be had: a UnionNeuron is actually a ProductNeuron in the sense of type theory (not a union or sum), and I'm not sure if it really extends Mesh and Skeleton rather than just containing them. There's two concerns here: the representations the type contains, and how it decides to mix/match/convert those. The reason you might want to decouple those concerns is motivated by Philipp's example of FlywireNeuron. I might want an easy interface to retrieve multiple representations from a backend, but when it comes time to wrangle them like Philipp's cable_length example, make a different decision about which to use for provenance or fidelity. That decision might not be backend-specific, so being able to mix and match a easy backend representation provider with a multi-representation wrangler is valuable. This is all possible in python it just takes some care to avoid stuffing too much into a single type/hierarchy or knotting yourself into inheritance diamonds.

[bdpedigo]

being able to mix and match a easy backend representation provider with a multi-representation wrangler is valuable

I like the mix and match idea. It feels like there is already a lot of code for representing meshes in particular, and I might want to use one or another depending on what I'm doing. But I feel like the key problem that I havent seen already solved besides in (navis or meshparty) is this tracking of the mapping between representations and storing annotations, so that would still be the value added here.

I guess the downside is it's harder to know what you can expect your neuron can do? If neuron.skeleton can be anything that implements a tree (lets say it just has .nodes and .edges or .parents), then you might not know if neuron.skeleton.cable_length() is going to exist, depending on the backend representation.

if im misunderstanding what you mean though, let me know

[bdpedigo]

Going to sketch out some thoughts — perhaps repeating parts of what others have said above. Feedback welcome:

I have a use case where I have both high and low resolution versions of meshes for the same neuron, in addition to skeletons, synapses, etc. So I think the most general form of a ProductNeuron (this and all other names up for debate) is just a container for some arbitrary set of morphology representations that are references to the same object, e.g.

neuron = ProductNeuron(mesh1, mesh2, skeleton, ...)

and something like neuron['mesh1'] or neuron.mesh1 just gives me back that underlying representation.

I'm wondering if because of all of the issues above, it is easier to just punt on convenience methods like ProductNeuron.cable_length, at least for this most general form. ProductNeuron would then only be responsible for

1. keeping track of correspondences between representations, where applicable. For instance, if a skeleton was computed from a mesh, storing that vertex correspondence
2. storing annotations, e.g. synapses or edits. I think in terms of implementation, this could even be a special case of 1). For instance, the synapses could be an additional PointCloud layer, and ProductNeuron knows the mapping from that point cloud to the other representations. Would be up for debate whether to expose annotations to the user in this way, but I think implementation-wise that could work.
3. keeping track of how to change the other representations when one is manipulated

A ProductNeuron is then just a container for a finite set of representations that we could define clear interfaces for. I think the obvious ones are Voxels, Mesh, SpatialGraph/Skeleton, and Points. I think each of those has very clear requirements, and we could define interfaces for each, e.g. a Mesh is anything that implements .vertices and.faces. As a start this could just come in the form of a single MeshNeuron class, but that could even evolve over time. For instance I've been using PyVista a lot, and if I wanted a PyVistaMeshNeuron class someday, I could roll my own as long as it has .vertices and .faces and whatever else we define for that interface.

Now, for each of those underlying representations, we'd be free to define whatever methods or computations make sense. .cable_length() for a SkeletonNeuron as in the discussion above. I like this approach because as others have pointed out above, I think certain representations (mesh/skeleton/pointcloud) or implementations (trimesh, pyvista mesh, something else) will always have advantages and disadvantages for particular computations. So it feels natural to leave those up to the specific implementation to say what it is good at.

Further, if particular combinations of representations get commonly used, that would also be easy to write down as a special case of ProductNeuron. Say, the mesh/skeleton/synapses combo that comes up a lot in meshparty. If someone wants to roll this particular combination into a class, they are free to do that given these tools. ProductNeuron handles the mapping logic, and the individual layers know what they are good at. So if I have thought about these underlying representations, I know what they are good at, and I want to write a PhilippNeuron class where .cable_length always uses the skeleton layer representation, I can do that. In fact, this is probably what many users would still end up using most of the time.

TL;DR

• A very general ProductNeuron class that just serves as a container for 1 or more views/layers referring to the same object, knows how those layers correspond to each other, and how to manipulate the other representations when one changes. It has no other computational logic. This would probably be for "super users" but also enable development of specific cases for many users (see point 3)
• A set of interfaces defined for each representation type (mesh/skeleton/voxels/...), and some specific implementations for each. As a start I am guessing these could be the current TreeNeuron,MeshNeuron,etc implementations basically dropped in. These know what computations they are good for.
• If desired, specific special cases of ProductNeuron x {choose your representations} which can add additional logic for what representations to use for what computations. The designer could choose to promote said computations to this class level.
• Could also have functions in navis.algorithms, and they could be as strict as they want about what functions apply to what representations/implementations

I like this division of labor because I constantly find myself using different representations for different tasks, and I think those representations will always be in flux depending on what people are doing with them. But at the end of the day the task I always find challenging is storing the mapping between those layers. This would be a way of abstracting that logic for the mapping from what meshparty already does, and allowing people to drop in things like TreeNeuron as one of the underlying representations if they wanted.

[bdpedigo]

https://bdpedigo.github.io/morphlink/examples/microns_example/ I mocked up a notebook with essentially version of what I described above. The code itself is very fragile but I hope illustrates what this kind of approach could look like.

Feedback very much welcome - I know @ceesem already has some thoughts on further separating things like annotations/labels from the morphology representation itself, using the same mapping logic that is in there.

This multiple representations + mappings is something that is coming up a lot for me right now, so I would love to steadily push on this somehow.

[clbarnes]

Something else to consider: a good chunk of navis' data is stored using pandas, and pandas supports the arrow backend. With increased interest in speeding up some hot loops with rust, as well as interoperating with nat in some places, v2 might be an opportunity to start insisting on the arrow backend internally.

To my understanding, pandas' arrow implementation isn't complete - pandas dataframes are basically a bag of arrow arrays rather than using an arrow table directly, so that interop may not be quite such an easy win. Additionally, while there are arrow implementations in rust, I have a sneaking suspicion that the arrow support in rust is not as good as the numpy support (somewhat conjecture, as I haven't used the arrow crates, just had a look over the docs).

A much larger change and probably too pie-in-the-sky would be to look into polars, which uses a more complete arrow representation and has better rust interop too.

[aschampion]

Rust's arrow support is great, largely due to the datafusion/ballista work.

Pandas' arrow support has come far and is now, like most of pandas, 90% great and the rest a minefield of deep, un-workaroundable bugs, especially with their new nullable types.

My opinion regarding polars vs pandas for navis is that pandas is a necessary evil to be useful to most folks (myself included) because it is just what you're going to be using everywhere else.

[aschampion]

Currently, navis can utilize multiple CPUs (via pathos) but it doesn't scale beyond a single machine. It'd be nice if there was an easy way to divvy up e.g. a big NBLAST across multiple nodes. @aschampion has already made a start in #111 for a Dask backend and the implementation is generic so that it would be easy to tack on more.

The right thing to do here to cover the most use cases is probably just use joblib, which allows for an ok-ish Dask and Ray backend as well. This might not scale to the largest things, but those will always require custom plumbing.

[schlegelp]

Not sure if that's would have to be a separate things but it'd be cool if we had different backends for NBLAST. Currently, NBLAST has it's own implementation for parallelisation but perhaps that can be unified?

Markus recently reminded me of https://github.com/navis-org/navis-nblast-gpu where I tested running NBLAST on GPUs. Back then my laptop had a pretty measly GPU but I should probably test this again with the new ARM chips.

[schlegelp]

One more thing to consider is loosing some weight. Candidates for deprecation and subsequent removal in navis 2.0:

• Cytoscape interface: this is a stow-away from back when navis was build from pymaid
• R interface: not sure if anyone uses this actively
• the vispy-based Viewer (replaced by Octarine)