import uuid
import numpy as np
import pandas as pd
from scipy import stats
from scipy.stats import beta
from scipy.stats import lognorm
from scipy.stats import norm
from scipy.stats import truncnorm
from scipy.stats import uniform
from nuctransportdb.generate_id import get_entry_str
from nuctransportdb.generate_id import ntd_namespace
from nuctransportdb.property2dataframe import preserve_value_type
from nuctransportdb.recommended_CV import rock_CV
# --- Create masks based on different inputs ---#
[docs]
def value_empty_mask(input_pd_df):
"""Filter data in the input DataFrame that is empty.
Args:
input_pd_df (pd.DataFrame): Input DataFrame.
Returns:
pd.Series: A Boolean series indicating which data is empty.
"""
return (
input_pd_df["value"].isna()
& input_pd_df["value_std"].isna()
& input_pd_df["value_min"].isna()
& input_pd_df["value_max"].isna()
& input_pd_df["sampled_data"].isna()
)
[docs]
def value_invalid_mask(input_pd_df):
"""Filter some data in the input DataFrame that contains invalid entries.
Args:
input_pd_df (pd.DataFrame): Input DataFrame.
Returns:
pd.Series: A Boolean Series indicating which data is invalid.
"""
# The value is not provided, but the standard deviation is given.
invalid_mask_1 = (
input_pd_df["value"].isna()
& input_pd_df["value_std"].notna()
& input_pd_df["sampled_data"].isna()
)
# Only the minimum value is provided
invalid_mask_2 = (
input_pd_df["value"].isna()
& input_pd_df["value_min"].notna()
& input_pd_df["value_max"].isna()
& input_pd_df["value_std"].isna()
& input_pd_df["sampled_data"].isna()
)
# Only the maximum value is provided
invalid_mask_3 = (
input_pd_df["value"].isna()
& input_pd_df["value_min"].isna()
& input_pd_df["value_max"].notna()
& input_pd_df["value_std"].isna()
& input_pd_df["sampled_data"].isna()
)
# Rows that satisfy value_min < value < value_max
is_scalar = input_pd_df["type"].eq("scalar")
scalar_df = input_pd_df.loc[is_scalar]
has_value = scalar_df["value"].notna()
has_min = scalar_df["value_min"].notna()
has_max = scalar_df["value_max"].notna()
# convert None to np.nan for comparison
value = scalar_df["value"].where(scalar_df["value"].notnull(), np.nan)
value_min = scalar_df["value_min"].where(scalar_df["value_min"].notnull(), np.nan)
value_max = scalar_df["value_max"].where(scalar_df["value_max"].notnull(), np.nan)
mask = pd.Series(False, index=input_pd_df.index)
mask.loc[is_scalar] = (
(has_value & has_min & (value < value_min)) # value < value_min
| (has_value & has_max & (value > value_max)) # value > value_max
| (has_min & has_max & (value_min > value_max)) # value_min > value_max
)
return invalid_mask_1 | invalid_mask_2 | invalid_mask_3 | mask
[docs]
def value_pdf_mask(input_pd_df):
"""Filter some data in the input DataFrame given by a probability distribution.
Args:
input_pd_df (pd.DataFrame): Input DataFrame.
Returns:
pd.Series: A Boolean series indicating which data entries have a corresponding probability distribution.
"""
is_pdf_mask = input_pd_df["sampled_data"].notna() & (
input_pd_df["type"] == "scalar"
)
invalid = value_invalid_mask(input_pd_df)
return is_pdf_mask & (~invalid)
[docs]
def value_truncnorm_mask(input_pd_df):
"""Filter some data in the input DataFrame that can be assumed with a truncated normal distribution.
Args:
input_pd_df (pd.DataFrame): Input DataFrame.
Returns:
pd.Series: A Boolean series indicating which data is truncated normal.
"""
is_truncnorm_mask = (
input_pd_df["value"].notna()
& input_pd_df["value_min"].notna()
& input_pd_df["value_max"].notna()
& input_pd_df["value_std"].notna()
& input_pd_df["sampled_data"].isna()
& (input_pd_df["type"] == "scalar")
)
invalid = value_invalid_mask(input_pd_df)
return is_truncnorm_mask & (~invalid)
[docs]
def value_lognorm_mask(input_pd_df):
"""Filter some data in the input DataFrame that can be assumed with a log-normal distribution.
Args:
input_pd_df (pd.DataFrame): Input DataFrame.
Returns:
pd.Series: A Boolean series indicating which data is log-normal.
"""
is_lognorm_mask1 = (
input_pd_df["value"].notna()
& input_pd_df["value_min"].isna()
& input_pd_df["value_max"].isna()
& input_pd_df["value_std"].notna()
& input_pd_df["sampled_data"].isna()
& (input_pd_df["type"] == "scalar")
)
is_lognorm_mask2 = (
input_pd_df["value"].notna()
& input_pd_df["value_min"].isna()
& input_pd_df["value_max"].notna()
& input_pd_df["value_std"].notna()
& input_pd_df["sampled_data"].isna()
& (input_pd_df["type"] == "scalar")
)
is_lognorm_mask3 = (
input_pd_df["value"].notna()
& input_pd_df["value_min"].notna()
& input_pd_df["value_max"].isna()
& input_pd_df["value_std"].notna()
& input_pd_df["sampled_data"].isna()
& (input_pd_df["type"] == "scalar")
)
is_lognorm_mask4 = (
input_pd_df["value"].notna()
& input_pd_df["value_min"].notna()
& input_pd_df["value_max"].isna()
& input_pd_df["value_std"].isna()
& input_pd_df["sampled_data"].isna()
& (input_pd_df["type"] == "scalar")
)
is_lognorm_mask5 = (
input_pd_df["value"].notna()
& input_pd_df["value_min"].isna()
& input_pd_df["value_max"].notna()
& input_pd_df["value_std"].isna()
& input_pd_df["sampled_data"].isna()
& (input_pd_df["type"] == "scalar")
)
is_lognorm_mask_single = (
input_pd_df["value"].notna()
& input_pd_df["value_min"].isna()
& input_pd_df["value_max"].isna()
& input_pd_df["value_std"].isna()
& input_pd_df["sampled_data"].isna()
& (input_pd_df["type"] == "scalar")
)
invalid = value_invalid_mask(input_pd_df)
return (
is_lognorm_mask1
| is_lognorm_mask2
| is_lognorm_mask3
| is_lognorm_mask4
| is_lognorm_mask5
| is_lognorm_mask_single
) & (~invalid)
[docs]
def value_PERT_mask(input_pd_df):
"""Filter some data in the input DataFrame that can be assumed with a PERT distribution.
Args:
input_pd_df (pd.DataFrame): Input DataFrame.
Returns:
pd.Series: A Boolean series indicating which data is PERT.
"""
is_PERT_mask = (
input_pd_df["value"].notna()
& input_pd_df["value_min"].notna()
& input_pd_df["value_max"].notna()
& input_pd_df["value_std"].isna()
& input_pd_df["sampled_data"].isna()
& (input_pd_df["type"] == "scalar")
)
invalid = value_invalid_mask(input_pd_df)
return is_PERT_mask & (~invalid)
# Present numbers with adaptive precision
# --- Generate different distributions based on the input parameters---#
[docs]
def generate_truncnorm(value, value_std, value_min, value_max, as_string=False, **kwargs):
"""Generate a truncated normal distribution based on the given parameters.
Args:
value (float): The mean of the original dataset.
value_std (float): The standard deviation of the original dataset.
value_min (float): The minimum value of the original dataset.
value_max (float): The maximum value of the original dataset.
as_string (bool): If True, return a string representation of the truncated normal distribution. If False, return a frozen distribution object. Default is False.
**kwargs: Extra keyword arguments accepted for compatibility and ignored.
Returns:
scipy.stats._distn_infrastructure.rv_continuous_frozen | str: A truncated normal distribution object if `as_string` is False, otherwise a string representation.
"""
a = (value_min - value) / value_std
b = (value_max - value) / value_std
if as_string:
return f"scipy.stats.truncnorm({format_number_adaptive(a)}, {format_number_adaptive(b)}, loc={format_number_adaptive(value)}, scale={format_number_adaptive(value_std)})"
return truncnorm(a, b, loc=value, scale=value_std)
[docs]
def generate_PERT(value, value_min, value_max, as_string=False, **kwargs):
"""Generate a PERT distribution based on the given parameters.
Args:
value (float): The mean of the original dataset.
value_min (float): The minimum value of the original dataset.
value_max (float): The maximum value of the original dataset.
as_string (bool): If True, return a string representation of the PERT distribution. If False, return a beta distribution object. Default is False.
**kwargs: Extra keyword arguments accepted for compatibility and ignored.
Returns:
scipy.stats._distn_infrastructure.rv_continuous_frozen | str: A beta distribution object if `as_string` is False, otherwise a string representation.
"""
# PERT distribution parameters calculated based on the book Risk Analysis: A Quantitative Guide by David Vose
value_mu = 4
value_a = 1 + value_mu * ((value - value_min) / (value_max - value_min))
value_b = 1 + value_mu * ((value_max - value) / (value_max - value_min))
if as_string:
return f"scipy.stats.beta({format_number_adaptive(value_a)}, {format_number_adaptive(value_b)}, loc={format_number_adaptive(value_min)}, scale={format_number_adaptive(value_max - value_min)})"
return beta(value_a, value_b, loc=value_min, scale=value_max - value_min)
[docs]
def generate_lognorm(value, value_std, value_min, as_string=False, **kwargs):
"""Generate a log-normal distribution based on the given parameters.
Args:
value (float): The mean of the original dataset.
value_std (float): The standard deviation of the original dataset.
value_min (float): Shift the distribution to start from the minimum value.
as_string (bool): If True, return a string representation of the log-normal distribution. If False, return a log-normal distribution object. Default is False.
**kwargs: Extra keyword arguments accepted for compatibility and ignored.
Returns:
scipy.stats._distn_infrastructure.rv_continuous_frozen | str: A log-normal distribution object if `as_string` is False, otherwise a string representation.
"""
if pd.isna(value_min): # check both None and NaN
value_min = 0.0
adjusted_mean = value - value_min
# lognorm parameters calculated based on https://doi.org/10.5281/ZENODO.4305949
log_std = np.sqrt(np.log(1 + (value_std / adjusted_mean) ** 2))
log_mean = np.log(adjusted_mean) - 0.5 * log_std**2
if as_string:
return f"scipy.stats.lognorm(s={format_number_adaptive(log_std)}, scale={format_number_adaptive(np.exp(log_mean))}, loc={format_number_adaptive(value_min)})"
return lognorm(s=log_std, scale=np.exp(log_mean), loc=value_min)
[docs]
def generate_norm(value, value_std, as_string=False, **kwargs):
"""Generate a normal distribution based on the given parameters.
Args:
value (float): The mean of the original dataset.
value_std (float): The standard deviation of the original dataset.
as_string (bool): If True, return a string representation of the normal distribution. If False, return a normal distribution object. Default is False.
**kwargs: Extra keyword arguments accepted for compatibility and ignored.
Returns:
scipy.stats._distn_infrastructure.rv_continuous_frozen | str: A normal distribution object if `as_string` is False, otherwise a string representation.
"""
norm_std = value_std
norm_mean = value
if as_string:
return f"scipy.stats.norm(loc={format_number_adaptive(norm_mean)}, scale={format_number_adaptive(norm_std)})"
return norm(loc=norm_mean, scale=norm_std)
[docs]
dict_sorption_rock_cv = rock_CV()
# --- Generate samples from different distributions ---#
[docs]
def generate_samples(
input_df_group,
df_pdf_type,
random_state=21,
):
"""Generate samples for a group of input data according to its specified distribution type, including uniform, truncated normal, log-normal, and PERT.
Args:
input_df_group (pd.DataFrame): The input DataFrame containing the data for one specific property and one type of pdf.
df_pdf_type (str): The type of probability distribution function to use. It can be "is_pdf_df", "is_uniform_df", "is_truncnorm_df", "is_lognorm_df", or "is_PERT_df".
random_state (int, optional): The random state for reproducibility. Defaults to 21.
Raises:
ValueError: If the distribution type is not recognized.
Returns:
np.ndarray: The generated samples falls within three standard deviations of the mean.
"""
samples_combined = []
# Get the nuclide name
rock_type = input_df_group["rock_type"].iloc[0]
nuclide_rock_property = input_df_group["nuclide_property"].iloc[0]
for (
sample_size,
sampled_data,
value,
value_std,
value_min,
value_max,
) in zip(
input_df_group["sample_size"],
input_df_group["sampled_data"],
input_df_group["value"],
input_df_group["value_std"],
input_df_group["value_min"],
input_df_group["value_max"], strict=False,
):
if sample_size is None or np.isnan(sample_size):
sample_size = 1000000
if df_pdf_type == "is_pdf_df":
# if the scipy.stats.rvs or numpy array is provided
samples = eval(sampled_data) if isinstance(sampled_data, str) else np.array(sampled_data)
elif df_pdf_type == "is_uniform_df":
uniform_samples = generate_uniform(value_min, value_max).rvs(
size=sample_size, random_state=random_state,
)
samples = uniform_samples
elif df_pdf_type == "is_truncnorm_df":
truncnorm_samples = generate_truncnorm(
value, value_std, value_min, value_max,
).rvs(size=sample_size, random_state=random_state)
samples = truncnorm_samples
elif df_pdf_type == "is_lognorm_df":
if pd.isna(value_std): # check both None and NaN
if nuclide_rock_property == "sorption_coefficient":
value_std = value * dict_sorption_rock_cv[rock_type]
if value == 0.0 and value_std == 0.0:
lognorm_samples = np.zeros(sample_size)
else:
lognorm_samples = generate_lognorm(value, value_std, value_min).rvs(
size=sample_size, random_state=random_state,
)
samples = lognorm_samples
elif df_pdf_type == "is_PERT_df":
PERT_samples = generate_PERT(value, value_min, value_max).rvs(
size=sample_size, random_state=random_state,
)
samples = PERT_samples
else:
msg = f"The df_pdf_type {df_pdf_type} is not recognized!"
raise ValueError(msg)
# Remove negative samples
samples = samples[samples >= 0]
if np.all(samples == 0):
samples_combined.append(samples)
else:
# Filter samples within three standard deviation of the mean
samples_zscore = stats.zscore(samples)
samples_filtered = samples[np.abs(samples_zscore) < 3]
samples_combined.append(samples_filtered)
return np.concatenate(samples_combined)
[docs]
def get_sample_statistics(input_property_group):
"""This function concatenates samples from different distributions for a group of input data and computes statistical properties (mean, std, min, max).
Args:
input_property_group (pd.DataFrame): DataFrame containing values for one specific property.
Returns:
tuple: A tuple containing the mean, standard deviation, minimum value, maximum value of the combined samples.
The tuple is in the following order:
(samples_mean, samples_std, samples_min, samples_max)
"""
# Generate samples for each distribution type
sample_combined = []
# Get samples from the existing pdf
is_pdf_mask = value_pdf_mask(input_property_group)
is_pdf_df = input_property_group[is_pdf_mask].copy()
if not is_pdf_df.empty:
sample_combined.append(generate_samples(is_pdf_df, "is_pdf_df"))
# Generate uniform samples
is_uniform_mask = value_uniform_mask(input_property_group)
is_uniform_df = input_property_group[is_uniform_mask].copy()
if not is_uniform_df.empty:
sample_combined.append(generate_samples(is_uniform_df, "is_uniform_df"))
# Generate truncated normal samples
is_truncnorm_mask = value_truncnorm_mask(input_property_group)
is_truncnorm_df = input_property_group[is_truncnorm_mask].copy()
if not is_truncnorm_df.empty:
sample_combined.append(generate_samples(is_truncnorm_df, "is_truncnorm_df"))
# Generate log-normal samples
is_lognorm_mask = value_lognorm_mask(input_property_group)
is_lognorm_df = input_property_group[is_lognorm_mask].copy()
if not is_lognorm_df.empty:
sample_combined.append(generate_samples(is_lognorm_df, "is_lognorm_df"))
# Generate PERT samples
is_PERT_mask = value_PERT_mask(input_property_group)
is_PERT_df = input_property_group[is_PERT_mask].copy()
if not is_PERT_df.empty:
sample_combined.append(generate_samples(is_PERT_df, "is_PERT_df"))
samples_all = np.concatenate(sample_combined)
# compute statistical properties
samples_mean = np.mean(samples_all)
samples_std = np.std(samples_all)
samples_min = np.min(samples_all)
samples_max = np.max(samples_all)
return (
format_number_adaptive(samples_mean),
format_number_adaptive(samples_std),
format_number_adaptive(samples_min),
format_number_adaptive(samples_max),
)
# --- Merge property values from multiple datasets into a single DataFrame with combined statistics ---#
[docs]
def generate_unit(property_name):
"""Generate a unit string and base for a given property name.
Args:
property_name (str): The name of the property for which to generate the unit.
Raises:
TypeError: If the property name is not recognized.
Returns:
str: The unit string for the property.
list: The unit base for the property.
"""
if property_name == "sorption_coefficient":
unit_str = "m^3/kg"
unit_base = [-1, 3, 0, 0, 0, 0, 0]
else:
msg = f"The property name {property_name} is not recognized!"
raise TypeError(msg)
return unit_str, unit_base
# Distribution descriptions for merged dataset
[docs]
distribution_descriptions = {
generate_norm: "A normal distribution",
generate_lognorm: "A lognormal distribution",
generate_PERT: "A PERT distribution",
generate_truncnorm: "A truncated normal distribution",
generate_uniform: "A uniform distribution",
}
[docs]
def merge_property_value(input_property, sample_size=1000000, source_type="merged", sampling_functions_by_property=None):
"""Merge property values from multiple datasets into a single DataFrame with combined statistics. The merged results include summary statistics and truncated normal distribution parameters for each property, along with metadata such as description and combined IDs.
Args:
input_property (pd.DataFrame): DataFrame containing property values.
sample_size (int): The number of samples to generate for each merged dataset. Default is 1000000.
source_type (str): "merged" or "default". Default is "merged".
sampling_functions_by_property (dict): A dictionary mapping property names to desired sampling functions. If provided, the function will use the specified sampling functions for each property instead of the truncnorm.
Returns:
pd.DataFrame: DataFrame containing merged property values, summary statistics, distribution parameters, and metadata for each property.
"""
# Replace np.nan with None
input_property = input_property.replace({np.nan: None, "": None})
# Preserve float for scalar and string for expression and dictionary
input_property[["value", "value_min", "value_max", "value_std"]] = input_property[
["value", "value_min", "value_max", "value_std"]
].map(preserve_value_type)
# Filter out empty and unreasonable values
is_empty_mask = value_empty_mask(input_property)
is_unreasonable_mask = value_invalid_mask(input_property)
input_property = input_property[~(is_empty_mask | is_unreasonable_mask)].copy()
# Filter out scalar type and non-scalar type data
input_prop_scalar = input_property.loc[input_property["type"] == "scalar"]
input_property_nonscalar = input_property.loc[input_property["type"] != "scalar"]
# Create an empty DataFrame with selected columns
selected_columns = [
"nuclide",
"source",
"type",
"value",
"value_min",
"value_max",
"value_std",
"sample_size",
"sampled_data",
"unit_str",
"unit_base",
"variable_name",
"variable_unit_str",
"variable_unit_base",
"description",
"agency",
"location",
"simplified_lithology",
"ID",
]
merged_property = pd.DataFrame(columns=selected_columns)
# Group the input property by 'nuclide'
input_nuclide_grouped = input_prop_scalar.groupby(["nuclide"])
# Get SDH namespace for generating unique IDs
NTD_NAMESPACE = ntd_namespace()
# Merge property values for each group
for nuclide_group_keys, input_nuclide_group in input_nuclide_grouped:
is_pdf_mask = value_pdf_mask(input_nuclide_group)
is_pdf_df = input_nuclide_group[is_pdf_mask].copy()
# If only one dataset is provided by a probability distribution, then the merge process is unnecessary
if (input_nuclide_group.shape[0] == 1) and (not is_pdf_df.empty):
merged_property = pd.concat(
[
property_df
for property_df in [merged_property, input_nuclide_group]
if not property_df.empty
],
ignore_index=True,
)
else:
# Get statistics from the combined samples
(
samples_mean,
samples_std,
samples_min,
samples_max,
) = get_sample_statistics(input_nuclide_group)
sample_statistics_config = {"value": samples_mean, "value_std": samples_std, "value_min": samples_min, "value_max": samples_max}
nuclide_name = nuclide_group_keys[0]
nuclide_rock_property = input_nuclide_group["nuclide_property"].iloc[0]
# if all values are identical, then no sampling is needed
if samples_std == 0.0:
samples_min, samples_max, sampling_rvs = None, None, None
description_dist = "All values are identical. All samples are"
else:
if (sampling_functions_by_property) and (nuclide_name in sampling_functions_by_property):
# use desired sampling functions
sampling_function = sampling_functions_by_property[nuclide_name]
sampling_rvs = sampling_function(as_string=True,**sample_statistics_config) + f".rvs(size={sample_size}, random_state=21)"
# get distribution description according to the sampling function
description_dist = distribution_descriptions.get(
sampling_function,
"A distribution",
)
else:
# use log normal distribution as default sampling function
sampling_rvs = generate_lognorm(as_string=True,**sample_statistics_config) + f".rvs(size={sample_size}, random_state=21)"
# get distribution description for log normal distribution
description_dist = distribution_descriptions[generate_lognorm]
# Get meta information
ids_list = list(input_nuclide_group["ID"])
ids_combined = ",".join(ids_list)
print_before_id = "combined datasets with ids:" if len(ids_list) > 1 else "dataset with id:"
number_of_datasets = input_nuclide_group.shape[0]
property_dict = {
"nuclide": nuclide_name,
"type": "scalar",
"source": source_type,
"value": samples_mean,
"value_min": samples_min,
"value_max": samples_max,
"value_std": samples_std,
"sample_size": sample_size,
"sampled_data": sampling_rvs,
"unit_str": generate_unit(nuclide_rock_property)[0],
"unit_base": generate_unit(nuclide_rock_property)[1],
"variable_name": None,
"variable_unit_str": None,
"variable_unit_base": None,
"description": None,
"agency": None,
"location": None,
"simplified_lithology": None,
"ID": None,
}
if source_type in {"default", "merged"}:
property_dict["description"] = (
f"{description_dist} fitted from {number_of_datasets} {print_before_id} {ids_combined}."
)
property_dict["ID"] = str(
uuid.uuid5(
NTD_NAMESPACE, get_entry_str(property_dict, nuclide_name),
),
)
else:
msg = f"The allowed source_type entries are 'default' and 'merged'. The source_type '{source_type}' is not recognized!"
raise ValueError(
msg,
)
# Create a new row with the merged property values
input_property_group_merged = pd.DataFrame([property_dict])
# Append the merged data of each property group to the merged DataFrame
merged_property = pd.concat(
[
property_df
for property_df in [merged_property, input_property_group_merged]
if not property_df.empty
],
ignore_index=True,
)
if (not input_property_nonscalar.empty) and (source_type == "merged"):
merged_property = pd.concat(
[
property_df
for property_df in [
merged_property,
input_property_nonscalar[selected_columns],
]
if not property_df.empty
],
ignore_index=True,
)
# Replace np.nan with None
merged_property = merged_property.replace({np.nan: None, "": None})
# Drop unnecessary columns where all values are None
minimum_required_columns = {
"nuclide", "source", "type",
"value", "value_min", "value_max", "value_std", "sample_size",
"sampled_data", "unit_str", "unit_base", "description", "ID"}
# Columns where every value is None/NaN
empty_columns = merged_property.columns[merged_property.isna().all()]
# Drop only the empty columns that are not required
cols_to_drop = [col for col in empty_columns if col not in minimum_required_columns]
return merged_property.drop(columns=cols_to_drop)