Skip to content

Commit

Permalink
docs: Add explanation on distribution placement.
Browse files Browse the repository at this point in the history
fix and comment distrib_placement.py
  • Loading branch information
drodarie committed Nov 22, 2024
1 parent 45752c6 commit 78be5e7
Show file tree
Hide file tree
Showing 2 changed files with 159 additions and 13 deletions.
103 changes: 102 additions & 1 deletion docs/getting-started/guide_components.rst
Original file line number Diff line number Diff line change
Expand Up @@ -186,7 +186,7 @@ Finally, to import our classes in our configuration file, we will modify the
},
"connectivity": {
"A_to_B": {
"strategy": "connectivity.DistanceConnectivity",
"strategy": "connectome.DistanceConnectivity",
"radius": 100,
"presynaptic": {
"cell_types": ["base_type"]
Expand Down Expand Up @@ -214,6 +214,107 @@ The parameters of `place` includes a dictionary linking each cell type name to i
:class:`PlacementIndicator <.placement.indicator.PlacementIndicator>`, and the `Chunk`
in which to place the cells.

First, let's assign a random position to each cell that need to be within the Chunk.
We will create a for loop to go through each `PlacementIndicator`.
Then, we will use the `guess` function that estimates the number of cells to place
within the `Chunk`.
Finally, we will draw a random position for each cell to place and store them with
the ``place_cells`` function.

.. code-block:: python
def place(self, chunk, indicators):
# For each placement indicator
for name_indic, indicator in indicators.items():
# Prepare an array to store positions
all_positions = np.empty((0, 3))
# Guess the number of cells to place within the chunk.
num_to_place = indicator.guess(chunk=chunk)
if num_to_place > 0:
# Assign a random position to the cells within this Chunk
positions = (
np.random.rand(num_to_place, 3) * chunk.dimensions + chunk.ldc
)
all_positions = np.concatenate([all_positions, positions])
self.place_cells(indicator, all_positions, chunk)
Now, we need to apply our distribution to each Partition of our circuit to see how
the cells are distributed within, along our directed axis.
Let's make a new for loop to loop through Partitions of each indicator.
Then, we extract the number of cells to place within the total Partition, still using
the `guess` function. For that, we convert the partition into a list of voxels.

.. literalinclude:: /../examples/tutorials/distrib_placement.py
:language: python
:lines: 54-62

Now, our function is supposed to place cells within only one Chunk.Fortunately,
Partitions can be decomposed into Chunks. So, we can retrieve from the
distribution the number of cells to place within the ``chunk`` parameter of the
function, according to its position along the directed ``axis``.

To do so, we need to define the interval occupied by the chunk within the partition.
We will leverage the lowest and highest coordinates of the chunk and partition with
respectively the functions ``ldc`` and ``mdc``.

.. literalinclude:: /../examples/tutorials/distrib_placement.py
:language: python
:lines: 63-68

``bounds`` is here the interval of ratios of the space occupied by a Chunk within the Partition
along the chosen axis.

We also need to take into account the case where the direction is negative. In this case, we should
invert the interval. E.g., if the ``bounds`` is [0.2, 0.3] then the inverted interval is [0.7, 0.8].

.. literalinclude:: /../examples/tutorials/distrib_placement.py
:language: python
:lines: 68-73

Great, we have the interval on which we want to place our cells.
Now, we will check how many cells to place within this interval according to our distribution
along the provided axis, knowing the total number of cells to place within the partition.
We will create a separate function for this called ``draw_interval``.
Additionally, we also need to take into account the two other dimensions. We will compute the
ratio of area occupied by the chunk along the two other directions:

.. literalinclude:: /../examples/tutorials/distrib_placement.py
:language: python
:lines: 75-82

So, your final place function should look like this

.. literalinclude:: /../examples/tutorials/distrib_placement.py
:language: python
:lines: 54-89

``draw_interval`` will leverage an
`acceptance-rejection method <https://en.wikipedia.org/wiki/Rejection_sampling>`_ .
In short, this method will draw n random values within [0, 1] and return the number which value
is less than the probability according to the the distribution to fall in the provided interval
boundaries.
We will retrieve the interval of definition of the distribution and within boundaries of our
provided ratio interval.

.. literalinclude:: /../examples/tutorials/distrib_placement.py
:language: python
:lines: 34-41

From the distribution, we can retrieve the probability for a drawn random value be lesser than
the upper bound and the probability for it to be greater than the lower bound.

.. literalinclude:: /../examples/tutorials/distrib_placement.py
:language: python
:lines: 42-45

Finally, we apply the acceptance-rejection algorithm to see how much of the cells to place in
the partition should be placed in the chunk:

.. literalinclude:: /../examples/tutorials/distrib_placement.py
:language: python
:lines: 46-52

We are done with the Placement! Here is how the full strategy looks like:

.. literalinclude:: /../examples/tutorials/distrib_placement.py
Expand Down
69 changes: 57 additions & 12 deletions examples/tutorials/distrib_placement.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,38 +5,83 @@

@config.node
class DistributionPlacement(PlacementStrategy):
"""
Place cells based on a scipy distribution along a specified axis.
"""

distribution = config.attr(type=types.distribution(), required=True)
"""Scipy distribution to apply to the cell coordinate along the axis"""
axis: int = config.attr(type=types.int(min=0, max=2), required=False, default=2)
"""Axis along which to apply the distribution (i.e. x, y or z)"""
direction: str = config.attr(
type=types.in_(["positive", "negative"]), required=False, default="positive"
)
"""Specify if the distribution is applied along positive or negative axis direction"""

def draw_interval(self, n, lower, upper):
"""
Apply an acceptance-rejection method to n points according to the provided
interval of the distribution.
This method draws n random values and returns the number which value is
less than the probability to fall in the provided interval boundaries.
:param int n: Number of points to draw
:param float lower: Lower bound of the interval within [0, 1]
:param float upper: Upper bound of the interval within [0, 1]
:return: Number of points passing the acceptance-rejection test.
:rtype: int
"""
# Extract the interval of values that can be generated by the distribution
# Since some distribution values can reach infinite, we clamp the interval
# so that the probability to be out of this interval is 1e-9
distrib_interval = self.distribution.definition_interval(1e-9)
selected_lt = self.distribution.cdf(
upper * np.diff(distrib_interval) + distrib_interval[0]
)
selected_gt = self.distribution.sf(
lower * np.diff(distrib_interval) + distrib_interval[0]
)
# Retrieve the value at the lower and upper bound ratios of the
# distribution's interval
value_upper = upper * np.diff(distrib_interval) + distrib_interval[0]
value_lower = lower * np.diff(distrib_interval) + distrib_interval[0]
# Retrieve the probability of a value to be lower than the upper value
selected_lt = self.distribution.cdf(value_upper)
# Retrieve the probability of a value to be greater than the lower value
selected_gt = self.distribution.sf(value_lower)
# Apply the acceptance-rejection method: a random point within [0,1] is
# accepted if it is lower than the probability to be less than the upper
# value and the probability to be greater than the lower value
random_numbers = np.random.rand(n)
selected = random_numbers <= selected_lt * selected_gt
# Returns the number of point passing the test
return np.count_nonzero(selected)

def place(self, chunk, indicators):
# For each placement indicator
for name_indic, indicator in indicators.items():
# Prepare an array to store positions
all_positions = np.empty((0, 3))
# For each partitions
for p in indicator.partitions:
# Guess the number of cells to place within the partition.
num_to_place = indicator.guess(voxels=p.to_voxels())
partition_size = (p._data.mdc - p._data.ldc)[self.axis]
chunk_borders = np.array([chunk.ldc[self.axis], chunk.mdc[self.axis]])
ratios = chunk_borders / partition_size
# Extract the size of the partition
partition_size = p._data.mdc - p._data.ldc
# Retrieve the ratio interval occupied by the current Chunk along the axis
chunk_borders = np.array([chunk.ldc, chunk.mdc])
ratios = (chunk_borders - p._data.ldc) / partition_size
bounds = ratios[:, self.axis]
if self.direction == "negative":
ratios = 1 - ratios
ratios = ratios[::-1]
# If the direction on which to apply the distribution is inverted,
# then the ratio interval should be inverted too.
bounds = 1 - bounds
bounds = bounds[::-1]

num_selected = self.draw_interval(num_to_place, *ratios)
# Draw according to the distribution the random number of cells to
# place in the Chunk
num_selected = self.draw_interval(
num_to_place, lower=bounds[0], upper=bounds[1]
)
# ratio of area occupied by the chunk along the two other dimensions
ratio_area = np.diff(np.delete(ratios, self.axis, axis=1), axis=0)
num_selected *= np.prod(ratio_area)
if num_selected > 0:
# Assign a random position to the cells within this Chunk
positions = (
np.random.rand(num_selected, 3) * chunk.dimensions + chunk.ldc
)
Expand Down

0 comments on commit 78be5e7

Please sign in to comment.