def tanimoto_pairwise(fps: np.ndarray) -> np.ndarray:
    """
    Compute the pairwise Tanimoto similarity for a set of fingerprints.

    Parameters:
    fps (np.ndarray)
        A 2D numpy array where each row represents a binary fingerprint.

    Returns:
    np.ndarray
        A 2D numpy array containing the Tanimoto similarity coefficients 
        for each pair of fingerprints. The element at position (i, j) in
        the output array represents the Tanimoto similarity coefficient
        between the i-th and j-th fingerprints in the input array.

    See:
        https://doi.org/10.1021/ci300261r
        https://doi.org/10.1111/j.1469-8137.1912.tb05611.x
    """
    c = np.matmul(fps, fps.T)
    ab = np.diag(c)
    ab = np.add.outer(ab, ab)
    return c / (ab - c)

# get pairwise distances based on Tanimoto similarity
dist_matrix = 1 - tanimoto_pairwise(fps)
# flatten lower triangle of symmetric matrix
dist_data = dist_matrix[np.tril_indices(len(dist_matrix), -1)]
# butina clustering
clusters = Butina.ClusterData(
    dist_data, distThresh=0.65, isDistData=True, nPts=len(dist_matrix)
)

# map molecules to clusters in an array
butina_groups = np.zeros(len(fps))
for i in range(len(clusters)):
    for j in clusters[i]:
        butina_groups[j] = i

f'Number of clusters: {len(clusters)}'

# set number of repeats and cross-validation folds
reps = 100
num_folds = 5
total_splits = reps*num_folds

# output results storage
results = {
    'ecfps': {
        'mae': np.zeros(total_splits),
        'r2': np.zeros(total_splits),
    },
    'func_fps': {
        'mae': np.zeros(total_splits),
        'r2': np.zeros(total_splits),
    },
    'pdv': {
        'mae': np.zeros(total_splits),
        'r2': np.zeros(total_splits),
    }
}

counter = 0
pbar = tqdm(total=total_splits)
for i in range(reps):
    # loop over repetitions
    # initiate splitter with new random state
    splitter = GroupKFold(n_splits=num_folds, shuffle=True, random_state=i)
    # split into folds based on butina clusters
    for fold in splitter.split(fps, groups=butina_groups):
        # loop over folds
        train_idx, test_idx = fold
        # get train and test y values
        y_train = y_true[train_idx]
        y_test = y_true[test_idx]

        # fit LGBM model and predict using Morgan Fingerprints
        X_train = fps[train_idx]
        X_test = fps[test_idx]
        model = lgb.LGBMRegressor(verbose=-1)
        model.fit(X_train, y_train)
        preds = model.predict(X_test)
        results['ecfps']['mae'][counter] = mean_absolute_error(y_test, preds)
        results['ecfps']['r2'][counter] = r2_score(y_test, preds)

        # fit LGBM model and predict using functional Morgan Fingerprints
        X_train = func_fps[train_idx]
        X_test = func_fps[test_idx]
        model = lgb.LGBMRegressor(verbose=-1)
        model.fit(X_train, y_train)
        preds = model.predict(X_test)
        results['func_fps']['mae'][counter] = mean_absolute_error(y_test, preds)
        results['func_fps']['r2'][counter] = r2_score(y_test, preds)

        # fit LGBM model and predict using physicochemical descriptor vectors
        X_train = pdvs[train_idx]
        X_test = pdvs[test_idx]
        model = lgb.LGBMRegressor(verbose=-1)
        model.fit(X_train, y_train)
        preds = model.predict(X_test)
        results['pdv']['mae'][counter] = mean_absolute_error(y_test, preds)
        results['pdv']['r2'][counter] = r2_score(y_test, preds)

        counter += 1
        pbar.update(1)
pbar.close()

plt.figure(figsize=(10, 6))
labels = {'ecfps': 'ECFP$4_{1024}$','pdv': 'PDV', 'func_fps': 'FCFP$4_{1024}$'}
for method, scores in results.items():
    x = scores['mae'].mean()
    # mae confidence intervals
    x_ci =stats.t.interval(
        0.95,
        total_splits-1,
        loc=x,
        scale=stats.sem(scores['mae'])
    )
    y = scores['r2'].mean()
    # r2 confidence intervals
    y_ci =stats.t.interval(
        0.95,
        total_splits-1,
        loc=y,
        scale=stats.sem(scores['r2'])
    )
        # plot means
    plt.scatter(x, y, alpha=0.5, label=labels[method])
    # plot confidence intervals
    plt.errorbar(x, y, alpha=0.5, xerr=(x_ci[1]-x_ci[0])/2, yerr=(y_ci[1]-y_ci[0])/2, capsize=5)

plt.title('Coefficent of determination versus mean absolute error on hPPB task')
plt.xlabel('MAE\u2193 [log$_{10}$(%)]')
plt.ylabel('R$^{2}\u2191$')
plt.legend()
plt.grid(True)
plt.show()