Skip to content

Commit

Permalink
Revert format changes
Browse files Browse the repository at this point in the history
  • Loading branch information
@eleftherioszisis
eleftherioszisis committed Mar 25, 2024
1 parent 85e76fb commit dc78193
Showing 1 changed file with 24 additions and 73 deletions.
97 changes: 24 additions & 73 deletions atlas_densities/densities/fitting.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,9 +92,7 @@ def create_dataframe_from_known_densities(
for region_name in result.index:
density_mask = average_densities["brain_region"] == region_name
for cell_type in cell_types:
cell_type_mask = density_mask & (
average_densities["cell_type"] == cell_type
)
cell_type_mask = density_mask & (average_densities["cell_type"] == cell_type)
cell_type = cell_type.lower().replace(" ", "_")
result.at[region_name, cell_type] = np.mean(
average_densities[cell_type_mask]["measurement"]
Expand Down Expand Up @@ -160,9 +158,7 @@ def fill_in_homogenous_regions(
density = np.mean(neuron_density[region_mask])
densities.at[child_region_name, "inhibitory_neuron"] = density
if cell_density_stddevs is None:
densities.at[
child_region_name, "inhibitory_neuron_standard_deviation"
] = density
densities.at[child_region_name, "inhibitory_neuron_standard_deviation"] = density
else:
densities.at[
child_region_name, "inhibitory_neuron_standard_deviation"
Expand All @@ -175,9 +171,7 @@ def fill_in_homogenous_regions(
id_mask = [id_ in desc_id_set for id_ in hierarchy_info["id"]]
for child_region_name in hierarchy_info[id_mask].index:
densities.at[child_region_name, "inhibitory_neuron"] = 0.0
densities.at[
child_region_name, "inhibitory_neuron_standard_deviation"
] = 0.0
densities.at[child_region_name, "inhibitory_neuron_standard_deviation"] = 0.0


def compute_average_intensity(
Expand Down Expand Up @@ -226,9 +220,7 @@ def _add_depths(region_map_df):
]
depth = 1
while parent_ids:
children_ids = region_map_df[
np.isin(region_map_df["parent_id"], list(parent_ids))
].index
children_ids = region_map_df[np.isin(region_map_df["parent_id"], list(parent_ids))].index
region_map_df.loc[children_ids, "depth"] = depth
parent_ids = set(children_ids)
depth += 1
Expand Down Expand Up @@ -257,9 +249,7 @@ def _fill_densities(region_map, region_map_df, df):
count = voxel_count.sum()
if count:
df.loc[id_, "voxel_count"] = count
df.loc[id_, "density"] = np.average(
df.loc[children, "density"], weights=voxel_count
)
df.loc[id_, "density"] = np.average(df.loc[children, "density"], weights=voxel_count)


def _apply_density_slices(gene_marker_volumes):
Expand Down Expand Up @@ -288,7 +278,6 @@ def _compute_average_intensities_helper(index, gene_marker_volumes, id_):
for marker, intensity in gene_marker_volumes.items():
mask_voxels = index.ravel(intensity["mask"])[voxel_ids]

# restricted_mask = np.logical_and(intensity["mask"], mask)
count = mask_voxels.sum()

if count <= 0:
Expand Down Expand Up @@ -405,20 +394,16 @@ def compute_average_intensities(

index = ValueToIndexVoxels(annotation)

# _helper = delayed(_compute_average_intensities_helper)
_helper = _compute_average_intensities_helper
work = []
for id_ in index.values:
if id_ == 0:
continue
work.append(_helper(index, gene_marker_volumes, id_))

# res = Parallel()(work)
res = work
densities = (
pd.DataFrame(
list(it.chain(*res)), columns=["marker", "id", "voxel_count", "density"]
)
pd.DataFrame(list(it.chain(*res)), columns=["marker", "id", "voxel_count", "density"])
.set_index("id")
.pivot(columns="marker")
.swaplevel(axis=1)
Expand All @@ -435,9 +420,7 @@ def compute_average_intensities(
result["brain_region"] = region_map_df.loc[result.index].name

for marker in gene_marker_volumes:
df = pd.DataFrame(
data=0.0, index=result.index, columns=["voxel_count", "density"]
)
df = pd.DataFrame(data=0.0, index=result.index, columns=["voxel_count", "density"])
df.update(densities[marker.lower()])
_fill_densities(region_map, region_map_df, df)
result[marker.lower()] = df["density"]
Expand Down Expand Up @@ -501,10 +484,7 @@ def _optimize_func(x, coefficient):
absolute_sigma=True,
)
ss_reg = np.sum((_optimize_func(xdata, parameters[0][0]) - np.mean(ydata)) ** 2)
ss_tot = (
np.sum((np.array(ydata) - _optimize_func(xdata, parameters[0][0])) ** 2)
+ ss_reg
)
ss_tot = np.sum((np.array(ydata) - _optimize_func(xdata, parameters[0][0])) ** 2) + ss_reg
# if total sum of square is null, no variance can be explained by the fitting
return {
"coefficient": parameters[0][0],
Expand All @@ -517,9 +497,7 @@ def _optimize_func(x, coefficient):


def compute_fitting_coefficients(
groups: dict[str, set[str]],
average_intensities: pd.DataFrame,
densities: pd.DataFrame,
groups: dict[str, set[str]], average_intensities: pd.DataFrame, densities: pd.DataFrame
) -> FittingData:
"""
Compute the linear fitting coefficient of the cloud of 2D points (average marker intensity,
Expand Down Expand Up @@ -581,21 +559,14 @@ def compute_fitting_coefficients(

result = {
group_name: {
cell_type: {
"coefficient": np.nan,
"standard_deviation": np.nan,
"r_square": np.nan,
}
cell_type: {"coefficient": np.nan, "standard_deviation": np.nan, "r_square": np.nan}
for cell_type in cell_types
}
for group_name in groups
}

clouds: dict[str, dict[str, dict[str, list[float]]]] = {
group_name: {
cell_type: {"xdata": [], "ydata": [], "sigma": []}
for cell_type in cell_types
}
group_name: {cell_type: {"xdata": [], "ydata": [], "sigma": []} for cell_type in cell_types}
for group_name in groups
}
L.info("Computing fitting coefficients in %d groups ...", len(groups))
Expand All @@ -611,31 +582,22 @@ def compute_fitting_coefficients(
intensity = average_intensities.at[region_name, cell_type[:-1]]
density = densities.at[region_name, cell_type]
if not (
np.isnan(density)
or np.isnan(intensity)
or intensity == 0.0
or density == 0.0
np.isnan(density) or np.isnan(intensity) or intensity == 0.0 or density == 0.0
):
assert_msg = f"in region {region_name} for cell type {cell_type}"
assert intensity >= 0.0, "Negative average intensity " + assert_msg
assert density >= 0.0, "Negative density " + assert_msg
standard_deviation = densities.at[
region_name, cell_type + "_standard_deviation"
]
assert not np.isnan(standard_deviation), (
"NaN standard deviation " + assert_msg
)
assert standard_deviation >= 0.0, (
"Negative standard deviation " + assert_msg
)
assert not np.isnan(standard_deviation), "NaN standard deviation " + assert_msg
assert standard_deviation >= 0.0, "Negative standard deviation " + assert_msg

clouds[group_name][cell_type]["xdata"].append(intensity)
clouds[group_name][cell_type]["ydata"].append(density)
clouds[group_name][cell_type]["sigma"].append(standard_deviation)

L.info(
"Computing regression coefficients for %d cell types ...", len(cell_types)
)
L.info("Computing regression coefficients for %d cell types ...", len(cell_types))
for cell_type in tqdm(cell_types):
cloud = clouds[group_name][cell_type]
result[group_name][cell_type] = linear_fitting_xy(
Expand Down Expand Up @@ -679,9 +641,7 @@ def fit_unknown_densities(
for marker in average_intensities.columns:
intensity = average_intensities.at[region_name, marker]
cell_type = marker + "+"
if np.isnan(densities.at[region_name, cell_type]) and not np.isnan(
intensity
):
if np.isnan(densities.at[region_name, cell_type]) and not np.isnan(intensity):
fitting = fitting_coefficients[group_name][cell_type]
if np.isnan(fitting["coefficient"]):
warnings.warn(
Expand All @@ -707,9 +667,7 @@ def _check_homogenous_regions_sanity(homogenous_regions: pd.DataFrame) -> None:
supported_cell_types = {"inhibitory", "excitatory"}
diff = actual_cell_types - supported_cell_types
if len(diff) > 0:
raise AtlasDensitiesError(
f"`homogenous_regions` has unexpected cell types: {diff}"
)
raise AtlasDensitiesError(f"`homogenous_regions` has unexpected cell types: {diff}")


def _check_average_densities_sanity(average_densities: pd.DataFrame) -> None:
Expand All @@ -719,13 +677,10 @@ def _check_average_densities_sanity(average_densities: pd.DataFrame) -> None:
actual_measurement_types = set(average_densities["measurement_type"])
diff = actual_measurement_types - {"cell density"}
if len(diff) > 0:
raise AtlasDensitiesError(
f"`average_densities` has unexpected measurement types: {diff}"
)
raise AtlasDensitiesError(f"`average_densities` has unexpected measurement types: {diff}")

nan_mask = (
average_densities["measurement"].isna()
| average_densities["standard_deviation"].isna()
average_densities["measurement"].isna() | average_densities["standard_deviation"].isna()
)

if np.any(nan_mask):
Expand All @@ -745,9 +700,7 @@ def _check_average_densities_sanity(average_densities: pd.DataFrame) -> None:
)


def _get_group_names(
region_map: "RegionMap", group_ids_config: dict
) -> dict[str, set[str]]:
def _get_group_names(region_map: "RegionMap", group_ids_config: dict) -> dict[str, set[str]]:
"""
Get AIBS names for regions in several region groups of interest.
Expand All @@ -769,9 +722,9 @@ def _get_group_names(
# the set of names of the Rest group is closed under taking descendants.
isocortex_id = region_map.find("Isocortex", attr="name").pop()
cerebellum_id = region_map.find("Cerebellum", attr="name").pop()
ascendant_ids = set(
region_map.get(isocortex_id, attr="id", with_ascendants=True)
) | set(region_map.get(cerebellum_id, attr="id", with_ascendants=True))
ascendant_ids = set(region_map.get(isocortex_id, attr="id", with_ascendants=True)) | set(
region_map.get(cerebellum_id, attr="id", with_ascendants=True)
)
group_ids["Rest"] -= ascendant_ids

return {
Expand Down Expand Up @@ -883,9 +836,7 @@ def linear_fitting( # pylint: disable=too-many-arguments
groups = _get_group_names(region_map, group_ids_config)

L.info("Computing fitting coefficients ...")
fitting_coefficients = compute_fitting_coefficients(
groups, average_intensities, densities
)
fitting_coefficients = compute_fitting_coefficients(groups, average_intensities, densities)
L.info("Fitting unknown average densities ...")
fit_unknown_densities(groups, average_intensities, densities, fitting_coefficients)

Expand Down

0 comments on commit dc78193

Please sign in to comment.