From 78be5e716ba0badae62ae0e412d931e3b1f7f9a5 Mon Sep 17 00:00:00 2001 From: drodarie Date: Fri, 22 Nov 2024 18:01:10 +0100 Subject: [PATCH] docs: Add explanation on distribution placement. fix and comment distrib_placement.py --- docs/getting-started/guide_components.rst | 103 +++++++++++++++++++++- examples/tutorials/distrib_placement.py | 69 ++++++++++++--- 2 files changed, 159 insertions(+), 13 deletions(-) diff --git a/docs/getting-started/guide_components.rst b/docs/getting-started/guide_components.rst index 3bda8c16..66fd9bd3 100644 --- a/docs/getting-started/guide_components.rst +++ b/docs/getting-started/guide_components.rst @@ -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"] @@ -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 `_ . +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 diff --git a/examples/tutorials/distrib_placement.py b/examples/tutorials/distrib_placement.py index cdf8cebd..a4b04f43 100644 --- a/examples/tutorials/distrib_placement.py +++ b/examples/tutorials/distrib_placement.py @@ -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 )