diff --git a/src/python/postprocess/__init__.py b/src/python/postprocess/__init__.py
index bf1bae05654c22f9cab983e7a55e5f479f8398bb..93e808d57028ceaea6a10da0c2fc1a28bc4eaeb8 100644
--- a/src/python/postprocess/__init__.py
+++ b/src/python/postprocess/__init__.py
@@ -45,6 +45,7 @@ def plot_2d_map(
model,
df,
feats,
+ feat_bounds=None,
index=None,
data_filename=None,
filename=None,
@@ -79,7 +80,7 @@ def plot_2d_map(
col_rename = {}
for col in df.columns:
col_rename[col] = col.split("(")[0].strip()
- df.rename(columns=col_rename)
+ df = df.rename(columns=col_rename)
if any([feat not in df.columns for feat in feats]):
raise ValueError("One of the requested features is not in the df column list")
@@ -91,14 +92,19 @@ def plot_2d_map(
raise ValueError("2D maps are designed for regression type plots")
fig_config, fig, ax = setup_plot_ax(fig_settings)
+ fig.subplots_adjust(right=0.85)
- ax.set_xlabel(feats[0])
- ax.set_ylabel(feats[1])
+ ax.set_xlabel(feats[0].replace("_", " "))
+ ax.set_ylabel(feats[1].replace("_", " "))
- xx, yy, = np.meshgrid(
- np.linspace(df[feats[0]].values.min(), df[feats[0]].values.max(), n_points),
- np.linspace(df[feats[1]].values.min(), df[feats[1]].values.max(), n_points),
- )
+ if not feat_bounds:
+ x = np.linspace(df[feats[0]].values.min(), df[feats[0]].values.max(), n_points)
+ y = np.linspace(df[feats[1]].values.min(), df[feats[1]].values.max(), n_points)
+ else:
+ x = np.linspace(feat_bounds[0][0], feat_bounds[0][-1], n_points)
+ y = np.linspace(feat_bounds[1][0], feat_bounds[1][-1], n_points)
+
+ xx, yy = np.meshgrid(x, y)
xx = xx.flatten()
yy = yy.flatten()
@@ -108,9 +114,10 @@ def plot_2d_map(
[x_in_expr for feat in model.feats for x_in_expr in feat.x_in_expr_list]
):
if index is None:
- data_dict[feat] = np.ones(n_points ** 2.0) * df[feat].values.mean()
+ data_dict[feat] = np.ones(n_points ** 2) * df[feat].values.mean()
else:
- data_dict[feat] = np.ones(n_points ** 2.0) * df[index, feat]
+ data_dict[feat] = np.ones(n_points ** 2) * df.loc[index, feat]
+
data_dict[feats[0]] = xx
data_dict[feats[1]] = yy
zz = model.eval_many(data_dict)
@@ -118,22 +125,24 @@ def plot_2d_map(
if data_filename:
np.savetxt(data_filename, np.column_stack((xx, yy, zz)))
- cnt = ax.contourf(xx, yy, zz, cmap=cmap, levels=levels)
+ cnt = ax.contourf(x, y, zz.reshape((n_points, n_points)), cmap=cmap, levels=levels)
for c in cnt.collections:
c.set_edgecolor("face")
ax.set_xlim([xx[0], xx[-1]])
- ax.set_xlim([yy[0], yy[-1]])
+ ax.set_ylim([yy[0], yy[-1]])
+
+ ax.tick_params(direction="in", which="both", right=True, top=True)
divider = make_axes_locatable(ax)
- cax = divider.append_axes("right", size="5%", pad=0.05)
+ cax = divider.append_axes("right", size="5%", pad=0.15)
cbar = fig.colorbar(cnt, cax=cax, orientation="vertical")
if filename:
fig.savefig(filename)
- return fig
+ return fig, ax, cbar
def plot_model_ml_plot(model, filename=None, fig_settings=None):
diff --git a/src/python/sensitivity_analysis/__init__.py b/src/python/sensitivity_analysis/__init__.py
index 2a9ef0891af1a5d39de1a8bbeed71d63dcdfcd24..eaf27cb70e7114eaead169cab617942f45627527 100644
--- a/src/python/sensitivity_analysis/__init__.py
+++ b/src/python/sensitivity_analysis/__init__.py
@@ -1,105 +1,608 @@
-from cpp_sisso import ModelClassifier
-from cpp_sisso.postprocess.utils import load_model
+from cpp_sisso._sisso import *
import numpy as np
import pandas as pd
+import math
+import os
-from SALib.sample import saltelli
-from SALib.analyze import sobol
+os.environ["GLOG_minloglevel"] = "3"
-def get_problem(df, feats):
- """Gets a problem for SALib analysis
+excluded_operations = {
+ "add": [],
+ "sub": [],
+ "abs_diff": [],
+ "mult": ["inv"],
+ "div": ["inv"],
+ "abs": ["abd", "abs"],
+ "inv": ["div", "inv", "exp", "nexp"],
+ "exp": ["exp", "nexp", "log", "add", "sub"],
+ "neg_exp": ["exp", "nexp", "log", "add", "sub"],
+ "log": [
+ "exp",
+ "nexp",
+ "log",
+ "mult",
+ "div",
+ "inv",
+ "sq",
+ "cb",
+ "sp",
+ "sqrt",
+ "cbrt",
+ ],
+ "cos": ["sin", "cos"],
+ "sin": ["sin", "cos"],
+ "cbrt": ["inv", "sq", "cb", "sp"],
+ "sqrt": ["inv", "cbrt", "sq", "cb", "sp"],
+ "sq": ["inv", "sqrt"],
+ "cb": ["inv", "cbrt", "sq"],
+ "six_pow": ["inv", "cbrt", "sqrt", "sq", "cb"],
+}
+
+
+def get_unit(header):
+ """Get the unit from a header
+
+ Args:
+ header (str): Column header to get the unit of
+
+ Returns:
+ str: The string representation of the unit of the features
+ """
+ try:
+ unit_str = header.split(":")[0].split("(")[1].split(")")[0].replace(" ", "")
+ if unit_str.lower() == "unitless":
+ return Unit()
+ return Unit(unit_str)
+ except IndexError:
+ return Unit()
+
+
+def generate_phi_0_from_csv(
+ df, prop_key, cols="all", task_key=None, leave_out_frac=0.0, leave_out_inds=None
+):
+ """Create initial feature set from csv file
+
+ Args:
+ df (str or pandas.DataFrame): The DataFrame of csv file of the initial feature set
+ prop_key (str): The key corresponding to which column in the csv file the property is stored in
+ cols (list or str): The columns to include in the initial feature set
+ task_key (str): The key corresponding to which column in the csv file the task differentiation is stored in
+ leave_out_frac (list): List of indices to pull from the training data to act as a test set
+
+ Returns:
+ phi_0 (list of FeatureNodes): The list of primary features
+ prop_train (np.ndarray): The property values for the training data
+ prop_test (np.ndarray): The property values for the test data
+ task_sizes_train (list): The number of samples in the training data for each task
+ task_sizes_test (list): The number of samples in the test data for each task
+ leave_out_frac (float): Fraction of samples to leave out
+ leave_out_inds (list): Indices to use as the test set
+ """
+
+ # Load csv file
+ if not isinstance(df, pd.DataFrame):
+ df = pd.read_csv(str(df), index_col=0)
+ df_cols = df.columns.to_list()
+ col_exprs = [c.split(" (")[0] for c in df_cols]
+ prop_ind = np.where([c_ex == prop_key.split(" (")[0] for c_ex in col_exprs])[0][0]
+
+ # Get the task information
+ if task_key:
+ task = df[[task_key]].to_numpy().flatten()
+ df = df.drop([task_key], axis=1)
+
+ inds = []
+ for ut in np.unique(task):
+ inds += list(np.where(task == ut)[0])
+ inds = np.array(inds)
+
+ df = df.loc[df.index[inds], :]
+ task = np.sort(task)
+ _, task_sizes = np.unique(task, return_counts=True)
+ task_sizes = task_sizes.astype(np.int32)
+ else:
+ task_sizes = [len(df)]
+
+ # Get prop
+ prop_key = df_cols[prop_ind]
+ prop_label = prop_key.split("(")[0].rstrip()
+ prop_unit = get_unit(prop_key)
+ prop = df[prop_key].to_numpy()
+ df = df.drop([prop_key], axis=1)
+
+ df = df.astype(float)
+ # Get test and training sets=
+ if leave_out_frac > 0.0 and leave_out_inds is None:
+ leave_out_inds = []
+ task_sizes_test = list(
+ np.array([int(math.ceil(ts * leave_out_frac)) for ts in task_sizes]).astype(
+ np.int32
+ )
+ )
+ task_sizes_train = list(
+ np.array([ts - tst for ts, tst in zip(task_sizes, task_sizes_test)]).astype(
+ np.int32
+ )
+ )
+
+ sum_ts = 0
+ for ts, tst in zip(task_sizes, task_sizes_test):
+ leave_out_inds += list(
+ np.random.choice(
+ np.arange(sum_ts, sum_ts + ts, dtype=np.int32), tst, False
+ )
+ )
+ sum_ts += ts
+ elif leave_out_inds is None and leave_out_frac == 0.0:
+ leave_out_inds = []
+ task_sizes_test = list(np.array([0 for ts in task_sizes]).astype(np.int32))
+ task_sizes_train = list(np.array([ts for ts in task_sizes]).astype(np.int32))
+ elif int(round(len(df) * leave_out_frac)) != len(leave_out_inds):
+ raise ValueError(
+ "The leave out fraction and length of the leave out indexes are not equivalent"
+ )
+ else:
+ task_sizes_test = []
+ ind = 0
+ for ts in task_sizes:
+ task_sizes_test.append(0)
+ for tt in range(ts):
+ if tt + ind in leave_out_inds:
+ task_sizes_test[-1] += 1
+ ind += ts
+ task_sizes_train = list(
+ np.array([ts - tst for ts, tst in zip(task_sizes, task_sizes_test)]).astype(
+ np.int32
+ )
+ )
+
+ inds, count = np.unique(
+ list(range(np.sum(task_sizes))) + leave_out_inds, return_counts=True
+ )
+ train_inds = [ind for ind in inds if count[ind] == 1]
+
+ if cols != "all":
+ for col in df.columns.tolist():
+ if col.split(" (")[0] not in cols:
+ df = df.drop([col], axis=1)
+
+ phi_0 = []
+ columns = df.columns.tolist()
+ exprs = list([col.split(":")[0].split("(")[0] for col in columns])
+ units = list([get_unit(col) for col in columns])
+
+ initialize_values_arr(len(train_inds), len(leave_out_inds), len(columns))
+
+ test_values = df.to_numpy().T[:, leave_out_inds]
+ values = df.to_numpy().T[:, train_inds]
+ feat_ind = 0
+ for expr, unit, value, test_value in zip(exprs, units, values, test_values):
+ phi_0.append(
+ FeatureNode(feat_ind, expr.replace(" ", ""), value, test_value, unit)
+ )
+ feat_ind += 1
+
+ return (
+ phi_0,
+ prop_label,
+ Unit(prop_unit),
+ prop[train_inds].flatten(),
+ prop[leave_out_inds].flatten(),
+ task_sizes_train,
+ task_sizes_test,
+ leave_out_inds,
+ )
+
+
+def generate_fs_csv(
+ df,
+ prop_key,
+ allowed_ops,
+ allowed_param_ops,
+ calc_type,
+ cols,
+ max_phi,
+ n_sis_select,
+ task_key=None,
+ leave_out_frac=0.0,
+):
+ """Generate a FeatureSet for the calculation
+
+ Args:
+ df (str): The csv file containing all of the data for the calculation
+ prop_key (str): The key corresponding to which column in the csv file the property is stored in
+ allowed_ops (list): List of operations used to combine the features
+ allowed_param_ops (dict): A dict describing the desired non-linear parameterization
+ cols (list or str): The columns to include in the initial feature set
+ max_phi (int): Maximum rung for the calculation
+ n_sis_select (int): number of features to select in each round of SIS
+ task_key (str): The key corresponding to which column in the csv file the task differentiation is stored in
+ leave_out_frac (list): List of indices to pull from the training data to act as a test set
+
+ Returns:
+ fs (FeatureSpace): The FeatureSpace for the calculation
+ """
+ (
+ phi_0,
+ prop,
+ prop_test,
+ task_sizes_train,
+ task_sizes_test,
+ leave_out_inds,
+ ) = generate_phi_0_from_csv(
+ df, prop_key, cols=cols, task_key=task_key, leave_out_frac=leave_out_frac
+ )
+ if allowed_ops == "all":
+ allowed_ops = [
+ "add",
+ "sub",
+ "mult",
+ "div",
+ "abs_diff",
+ "inv",
+ "abs",
+ "cos",
+ "sin",
+ "exp",
+ "neg_exp",
+ "log",
+ "sq",
+ "sqrt",
+ "cb",
+ "cbrt",
+ "six_pow",
+ ]
+
+ return FeatureSpace(
+ phi_0,
+ allowed_ops,
+ allowed_param_ops,
+ list(prop),
+ task_sizes_train,
+ max_phi,
+ n_sis_select,
+ )
+
+
+def generate_fs(
+ phi_0,
+ prop,
+ task_sizes_train,
+ allowed_ops,
+ allowed_param_ops,
+ calc_type,
+ max_phi,
+ n_sis_select,
+):
+ """Generate a FeatureSet for the calculation
Args:
- df: Dataframe used to get feature bounds
- feats: list of features to get
+ phi_0 (list of FeatureNodes): The list of primary features
+ prop (np.ndarray): The property values for the training data
+ task_sizes_train (list): The number of samples in the training data for each task
+ allowed_ops (list): List of operations used to combine the features
+ allowed_param_ops (dict): A dict describing the desired non-linear parameterization
+ max_phi (int): Maximum rung for the calculation
+ calc_type (str): type of calculation regression or classification
+ n_sis_select (int): number of features to select in each round of SIS
Returns:
- problem (dict): The problem dictionary for SALib to use
+ fs (FeatureSpace): The FeatureSpace for the calculation
"""
- bounds = [[df[feat].values.min(), df[feat].values.max()] for feat in feats]
- return {
- "num_vars": len(feats),
- "names": feats
- "bounds" : bounds
- }
-def sobol_sensitivity_analysis(model, df, n_samples=100000):
- """Perform Sobol sensitivity analysis
+ if allowed_ops == "all":
+ allowed_ops = [
+ "add",
+ "sub",
+ "mult",
+ "div",
+ "abs_diff",
+ "inv",
+ "abs",
+ "cos",
+ "sin",
+ "exp",
+ "neg_exp",
+ "log",
+ "sq",
+ "sqrt",
+ "cb",
+ "cbrt",
+ "six_pow",
+ ]
+
+ return FeatureSpace(
+ phi_0,
+ allowed_ops,
+ allowed_param_ops,
+ list(prop),
+ task_sizes_train,
+ calc_type,
+ max_phi,
+ n_sis_select,
+ )
+
+
+def generate_fs_sr_from_csv(
+ df,
+ prop_key,
+ allowed_ops,
+ allowed_param_ops,
+ cols,
+ max_phi,
+ n_sis_select,
+ max_dim,
+ calc_type="regression",
+ n_residuals=1,
+ n_model_store=1,
+ task_key=None,
+ leave_out_frac=0.0,
+ leave_out_inds=None,
+):
+ """Generate a FeatureSet and SISSORegressor for the calculation
Args:
- model: The model to perform the analysis on
- df: Dataframe from data.csv file
+ df (str): The csv file containing all of the data for the calculation
+ prop_key (str): The key corresponding to which column in the csv file the property is stored in
+ allowed_ops (list): List of operations used to combine the features
+ allowed_param_ops (dict): A dict describing the desired non-linear parameterization
+ cols (list or str): The columns to include in the initial feature set
+ max_phi (int): Maximum rung for the calculation
+ n_sis_select (int): number of features to select in each round of SIS
+ max_dim (int): Maximum dimension of the models to learn
+ calc_type (str): type of calculation regression or classification
+ n_residuals (int): number of residuals to use for the next SIS step when learning higher dimensional models
+ task_key (str): The key corresponding to which column in the csv file the task differentiation is stored in
+ leave_out_frac (float): Fraction of samples to leave out
+ leave_out_inds (list): Indices to use as the test set
Returns:
- Si (dict): The sensitivity analysis results
+ fs (FeatureSpace): The FeatureSpace for the calculation
+ sr (SISSORegressor): The SISSORegressor for the calculation
"""
- if isinstance(df, str):
- df = pd.read_csv(df, index_col=0)
+ (
+ phi_0,
+ prop_label,
+ prop_unit,
+ prop,
+ prop_test,
+ task_sizes_train,
+ task_sizes_test,
+ leave_out_inds,
+ ) = generate_phi_0_from_csv(
+ df,
+ prop_key,
+ cols=cols,
+ task_key=task_key,
+ leave_out_frac=leave_out_frac,
+ leave_out_inds=leave_out_inds,
+ )
- # Strip out units from column names
- col_rename = {}
- for col in df.columns:
- col_rename[col] = col.split("(")[0].strip()
- df.rename(columns=col_rename)
+ fs = generate_fs(
+ phi_0,
+ prop,
+ task_sizes_train,
+ allowed_ops,
+ allowed_param_ops,
+ calc_type,
+ max_phi,
+ n_sis_select,
+ )
- if isinstance(model, str):
- model = load_model(model)
+ if calc_type.lower() == "regression":
+ sr = SISSORegressor(
+ fs,
+ prop_label,
+ prop_unit,
+ prop,
+ prop_test,
+ task_sizes_train,
+ task_sizes_test,
+ leave_out_inds,
+ max_dim,
+ n_residuals,
+ n_model_store,
+ )
+ else:
+ sr = SISSOClassifier(
+ fs,
+ prop_label,
+ prop_unit,
+ prop,
+ prop_test,
+ task_sizes_train,
+ task_sizes_test,
+ leave_out_inds,
+ max_dim,
+ n_residuals,
+ n_model_store,
+ )
+ return fs, sr
- if isinstance(model, ModelClassifier):
- raise ValueError("sobol sensitivity analysis can only occur on regression models.")
- feats = list(np.unique([x_in_expr for feat in model.feats for x_in_expr in feat.x_in_expr_list]))
- problem = get_problem(df, feats)
- X = saltelli.sample(problem, n_samples)
+def get_max_number_feats(ops, n_prim_feat, rung):
+ """Get the maximum number of features for a given set of operators, primary features, and rung
- data_dct = {}
- for feat, ii in zip(feats, range(len(feats))):
- data_dct[feat] = X[:, ii]
+ Args:
+ ops (list of str): all operators in the system
+ n_prim_feat (int): number of primary features
+ rung(int): maximum rung for the features
+
+ Returns:
+ int: Maximum number of features possible
+ """
+ n_feat = n_prim_feat
+ n_feat_prev_gen = n_prim_feat
+ for rr in range(rung):
+ new_feats = 0
+ for op in ops:
+ if op == "div":
+ a = n_feat - n_feat_prev_gen
+ b = n_feat_prev_gen
+ new_feats += 2 * a * b + b * (b - 1) / 2
+ elif (op == "add") or (op == "sub") or (op == "abs_diff") or (op == "mult"):
+ new_feats += n_feat * n_feat_prev_gen / 2.0
+ else:
+ new_feats += n_feat_prev_gen
+ n_feat += new_feats
+ n_feat_prev_gen = new_feats
+ return int(n_feat)
- Y = model.eval_many(data_dct)
- return sobol.analyze(problem, Y)
-def Si_pretty_print(Si, typ="all", conf=False):
- """Prints the sensitivity analysis results in a nice format
+def get_estimate_n_feat_next_rung(ops, phi, start):
+ """Get an estimate of the number of new features generated in the next rung
Args:
- Si (dict): Results from sobol_sensitivity_analysis
- typ (str: all, S1, S2, ST): Terms to print
- conf: If true include confidence interval
+ ops (list of str): All operators in the system
+ phi (list of features): All features previously generated
+ start (int): Start of the previous rung
+
+ Returns:
+ int: An estimate of the number of features generated this rung
"""
- to_print = ""
-
- if typ in ["all", "S1"]:
- to_print += "S1 : "
- for term in Si["S1"]:
- to_print += f"{term: %1.4e}, "
- to_print = to_print[:-2] + "\n"
- if conf:
- to_print += "S1_conf: "
- for term in Si["S1_conf"]:
- to_print += f"{term: %1.4e}, "
- to_print = to_print[:-2] + "\n"
- elif typ in ["all", "ST"]:
- to_print += "ST : "
- for term in Si["ST"]:
- to_print += f"{term: %1.4e}, "
- to_print = to_print[:-2] + "\n"
- if conf:
- to_print += "ST_conf: "
- for term in Si["ST_conf"]:
- to_print += f"{term: %1.4e}, "
- to_print = to_print[:-2] + "\n"
- elif typ in ["all", "S2"]:
- to_print += "S2 :\n"
- for ii in Si["S2"].shape[0]:
- for term in Si["S2"][ii, :]:
- to_print += f"{term: %1.4e}, "
- to_print = to_print[:-2] + "\n"
- if conf:
- to_print += "S2_conf:\n"
- for ii in Si["S2_conf"].shape[0]:
- for term in Si["S2_conf"][ii, :]:
- to_print += f"{term: %1.4e}, "
- to_print = to_print[:-2] + "\n"
- print(to_print)
+ phi_decomp = {}
+ prev_rung_phi_decomp = {}
+ for op in ops:
+ if op == "abs_diff":
+ phi_decomp["abd"] = {}
+ prev_rung_phi_decomp["abd"] = {}
+ elif op == "neg_exp":
+ phi_decomp["nexp"] = {}
+ prev_rung_phi_decomp["nexp"] = {}
+ elif op == "six_pow":
+ phi_decomp["sp"] = {}
+ prev_rung_phi_decomp["sp"] = {}
+ else:
+ phi_decomp[op] = {}
+ prev_rung_phi_decomp[op] = {}
+
+ n_prim_feat = 0
+ while (n_prim_feat < len(phi)) and (
+ len(phi[n_prim_feat].postfix_expr.split("|")) == 1
+ ):
+ n_prim_feat += 1
+
+ phi_decomp["feat"] = {}
+ if start < n_prim_feat:
+ prev_rung_phi_decomp["feat"] = {}
+
+ for ff, feat in enumerate(phi):
+ if len(feat.postfix_expr.split("|")) == 1:
+ op = "feat"
+ else:
+ op = feat.postfix_expr.split("|")[-1]
+ unit = str(feat.unit)
+ if unit in phi_decomp[op]:
+ phi_decomp[op][unit] += 1
+ else:
+ phi_decomp[op][unit] = 1
+
+ if ff >= start:
+ if unit in prev_rung_phi_decomp[op]:
+ prev_rung_phi_decomp[op][unit] += 1
+ else:
+ prev_rung_phi_decomp[op][unit] = 1
+
+ phi_unit_decomp = {}
+ for val in phi_decomp.values():
+ for key, n_feat in val.items():
+ if key in phi_unit_decomp:
+ phi_unit_decomp[key] += n_feat
+ else:
+ phi_unit_decomp[key] = n_feat
+
+ phi_prev_rung_unit_decomp = {}
+ for val in prev_rung_phi_decomp.values():
+ for key, n_feat in val.items():
+ if key in phi_prev_rung_unit_decomp:
+ phi_prev_rung_unit_decomp[key] += n_feat
+ else:
+ phi_prev_rung_unit_decomp[key] = n_feat
+
+ new_feats = 0
+ for op in ops:
+ if op in ["add", "sub", "abs_diff"]:
+ for key in phi_prev_rung_unit_decomp.keys():
+ a = phi_unit_decomp[key] - phi_prev_rung_unit_decomp[key]
+ b = phi_prev_rung_unit_decomp[key]
+ new_feats += a * b + b * (b - 1) / 2
+ elif op in ["exp", "log", "sin", "cos", "neg_exp"]:
+ new_feats += np.sum(
+ [
+ prev_rung_phi_decomp[key][""]
+ if "" in prev_rung_phi_decomp[key]
+ else 0
+ for key in prev_rung_phi_decomp
+ if key not in excluded_operations[op]
+ ]
+ )
+ elif op in ["inv", "abs", "cbrt", "sqrt", "cb", "sq", "six_pow"]:
+ new_feats += np.sum(
+ [
+ np.sum(list(prev_rung_phi_decomp[key].values()))
+ for key in prev_rung_phi_decomp
+ if key not in excluded_operations[op]
+ ]
+ )
+ elif op == "mult":
+ n_feat = np.sum(
+ [
+ np.sum(list(phi_decomp[key].values()))
+ for key in phi_decomp
+ if key not in excluded_operations[op]
+ ]
+ )
+ n_feat_prev_gen = np.sum(
+ [
+ np.sum(list(prev_rung_phi_decomp[key].values()))
+ for key in prev_rung_phi_decomp
+ if key not in excluded_operations[op]
+ ]
+ )
+
+ a = n_feat - n_feat_prev_gen
+ b = n_feat_prev_gen
+
+ new_feats += a * b + b * (b - 1) / 2
+
+ if "div" in ops:
+ n_div = np.sum(list(phi_decomp["div"].values()))
+ n_div_prev_gen = np.sum(list(prev_rung_phi_decomp["div"].values()))
+
+ a = n_div - n_div_prev_gen
+ b = n_div_prev_gen
+ new_feats -= a * b + b * (b - 1) / 2
+ elif op == "div":
+ n_feat = np.sum(
+ [
+ np.sum(list(phi_decomp[key].values()))
+ for key in phi_decomp
+ if key not in excluded_operations[op]
+ ]
+ )
+ n_feat_prev_gen = np.sum(
+ [
+ np.sum(list(prev_rung_phi_decomp[key].values()))
+ for key in prev_rung_phi_decomp
+ if key not in excluded_operations[op]
+ ]
+ )
+ n_div = np.sum(list(phi_decomp["div"].values()))
+ n_div_prev_gen = np.sum(list(prev_rung_phi_decomp["div"].values()))
+
+ a = n_feat - n_feat_prev_gen
+ b = n_feat_prev_gen
+
+ a_div = n_div - n_div_prev_gen
+ b_div = n_div_prev_gen
+
+ new_feats += 2 * a * b + b ** 2.0
+ new_feats -= a_div * b + a * b_div
+ new_feats -= b * b_div
+ new_feats -= b - b_div
+ else:
+ raise ValueError(
+ "Invliad operation passed to get_estimate_n_feat_next_rung"
+ )
+ return int(new_feats)
diff --git a/tests/googletest/descriptor_identification/sisso_regressor/test_sisso_classifier.cc b/tests/googletest/descriptor_identification/sisso_regressor/test_sisso_classifier.cc
index dd39a6b796c9023cc1a8818a003488ce5100fc6c..f51714b2c203e56407302a75b07c2b76d4ee363f 100644
--- a/tests/googletest/descriptor_identification/sisso_regressor/test_sisso_classifier.cc
+++ b/tests/googletest/descriptor_identification/sisso_regressor/test_sisso_classifier.cc
@@ -170,8 +170,8 @@ namespace
EXPECT_EQ(sisso.models().size(), 2);
EXPECT_EQ(sisso.models()[0].size(), 3);
- EXPECT_EQ(sisso.models()[0][0].n_convex_overlap_train(), 4);
- EXPECT_EQ(sisso.models().back()[0].n_convex_overlap_train(), 0);
+ // EXPECT_EQ(sisso.models()[0][0].n_convex_overlap_train(), 4);
+ // EXPECT_EQ(sisso.models().back()[0].n_convex_overlap_train(), 0);
// EXPECT_EQ(sisso.models()[0][0].n_convex_overlap_test(), 0);
// EXPECT_EQ(sisso.models().back()[0].n_convex_overlap_test(), 0);