Distribution-based Placement Strategy¶
Note
This example presents in more advanced tutorial to write BSB components. If you are not familiarized with BSB components, check out the getting started section on writing components.
We will start from the following configuration file (corresponds to the first network file from the getting started tutorial):
{
"name": "Starting example",
"storage": {
"engine": "hdf5",
"root": "network.hdf5"
},
"network": {
"x": 200.0,
"y": 200.0,
"z": 200.0
},
"regions": {
"brain_region": {
"type": "stack",
"children": ["base_layer", "top_layer"]
}
},
"partitions": {
"base_layer": {
"type": "layer",
"thickness": 100
},
"top_layer": {
"type": "layer",
"thickness": 100
}
},
"cell_types": {
"base_type": {
"spatial": {
"radius": 2.5,
"density": 3.9e-4
}
},
"top_type": {
"spatial": {
"radius": 7,
"count": 40
}
}
},
"placement": {
"example_placement": {
"strategy": "bsb.placement.RandomPlacement",
"cell_types": ["base_type"],
"partitions": ["base_layer"]
},
"top_placement": {
"strategy": "bsb.placement.RandomPlacement",
"cell_types": ["top_type"],
"partitions": ["top_layer"]
}
},
"connectivity": {
"A_to_B": {
"strategy": "bsb.connectivity.AllToAll",
"presynaptic": {
"cell_types": ["base_type"]
},
"postsynaptic": {
"cell_types": ["top_type"]
}
}
}
}
Let’s save this new configuration in our project folder under the name config_placement.json
Description of the strategy to implement¶
We want here to implement a distribution placement strategy: cells will be placed within
their Partition following a probability distribution along a certain axis and
direction. For instance, let us use the
alpha random distribution:
The distribution should be a density function that produces random numbers, according
to the distance along the axis from a border of the Partition.
Here, we want to control the distribution of the cells within the base_layer
with respect to the border with the top_layer.
Components boiler plate¶
We will write our PlacementStrategy class in a placement.py file,
from bsb import config, PlacementStrategy
@config.node
class DistributionPlacement(PlacementStrategy):
# add your attributes here
def place(self, chunk, indicators):
# write your code here
pass
Here, our class will extend from
PlacementStrategy which is an
abstract class that requires you to implement the
place function.
Note that this strategy leverages the @config.node python decorator.
The configuration node decorator allows you to pass the parameters
defined in the configuration file to the class. It will also handle the
type testing of your configuration attributes (e.g., make sure your
axis parameter is a positive integer). We will see in the following sections how
to create your class configuration attributes.
Add configuration attributes¶
For our strategy, we need to pass a list of parameters to its class through our configuration file:
a density function
distribution, defined by the user, from the list of scipy stats functions.an
axisalong which to apply the distribution. The latter should be in [0, 1, 2]a
directionthat if set will apply the distribution along the"negative"direction, i.e from top to bottom ("positive"direction if not set).
This translates into 3 configuration attributes that you can add to your class:
from bsb import config, PlacementStrategy, types
@config.node
class DistributionPlacement(PlacementStrategy):
distribution = config.attr(type=types.distribution(), required=True)
axis: int = config.attr(type=types.int(min=0, max=2), required=False, default=2)
direction: str = config.attr(type=types.in_(["positive", "negative"]),
required=False, default="positive")
def place(self, chunk, indicators):
# write your code here
pass
distribution is required, and should correspond to a
distribution node which interface scipy distributions.axis here is an optional integer attribute with a default value set to 2.direction is an optional string attribute that can be either the string "positive"
or "negative" (see in_).At this stage, you have created a python class with minimal code implementation, you should now link it to your configuration file. To import our class in our configuration file, we will modify the placement block:
"placement": {
"alpha_placement": {
"strategy": "placement.DistributionPlacement",
"distribution": {
"distribution": "alpha",
"a": 8
},
"axis": 0,
"direction": "negative",
"cell_types": ["base_type"],
"partitions": ["base_layer"]
},
"top_placement": {
"strategy": "bsb.placement.RandomPlacement",
"cell_types": ["top_type"],
"partitions": ["top_layer"]
}
}
Implement the python methods¶
The place function will be used here to produce and store a
PlacementSet for each cell type population
to place in the selected Partition.
BSB is parallelizing placement jobs for each Chunk concerned.
The parameters of place includes a dictionary linking each cell type name to its
PlacementIndicator, and the Chunk
in which to place the cells.
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 two for loops to iterate over the Partitions of each indicator.
Then, we extract the number of cells to place within the total Partition, using
the guess function.
For that, we will convert the partition into a list of voxels.
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())
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 attributes ldc and mdc.
# 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]
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].
bounds = ratios[:, self.axis]
if self.direction == "negative":
# 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]
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:
# 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)
So, your final place function should look like this
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())
# 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":
# 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]
# 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
)
all_positions = np.concatenate([all_positions, positions])
self.place_cells(indicator, all_positions, chunk)
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.
# 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)
# 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]
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.
# 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)
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:
# 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)
We are done with the Placement! Here is how the full strategy looks like:
import numpy as np
from bsb import PlacementStrategy, config, types
@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)
# 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())
# 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":
# 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]
# 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
)
all_positions = np.concatenate([all_positions, positions])
self.place_cells(indicator, all_positions, chunk)
Enjoy¶
You have done the hardest part! Now, you should be able to run the reconstruction once again with your brand new component.
bsb compile config_placement.json --verbosity 3