Commit of work from summer

dataset_size.py

@@ -0,0 +1,37 @@
+# This file is just to tell how many total images there are in the dataset
+
+import threshold_refac as tr
+import torch
+
+config = tr.load_config()
+ENSEMBLE_PATH = f"{config['paths']['model_output']}{config['ensemble']['name']}"
+
+test_dset = torch.load(f'{ENSEMBLE_PATH}/test_dataset.pt')
+val_dset = torch.load(f'{ENSEMBLE_PATH}/val_dataset.pt')
+train_dset = torch.load(f'{ENSEMBLE_PATH}/train_dataset.pt')
+
+
+print(
+    f'Total number of images in dataset: {len(test_dset) + len(val_dset) + len(train_dset)}'
+)
+print(f'Test: {len(test_dset)}, Val: {len(val_dset)}, Train: {len(train_dset)}')
+
+
+def preprocess_data(data, device):
+    mri, xls = data
+    mri = mri.unsqueeze(0).to(device)
+    xls = xls.unsqueeze(0).to(device)
+    return (mri, xls)
+
+
+# Loop through images and determine how many are positive and negative
+positive = 0
+negative = 0
+for _, (_, target) in enumerate(test_dset + train_dset + val_dset):
+    actual = list(target.cpu().numpy())[1].item()
+    if actual == 1:
+        positive += 1
+    else:
+        negative += 1
+
+print(f'Positive: {positive}, Negative: {negative}')

threshold_xarray.py

@@ -12,6 +12,7 @@ import matplotlib.pyplot as plt
 import matplotlib.ticker as mtick
 
 
+
 # The datastructures for this file are as follows
 # models_dict: Dictionary - {model_id: model}
 # predictions: DataArray - (data_id, model_id, prediction_value) - Prediction value has coords ['negative_prediction', 'positive_prediction', 'negative_actual', 'positive_actual']
@@ -78,7 +79,7 @@ def get_ensemble_predictions(models, dataset, device, id_offset=0):
         zeros,
         dims=('data_id', 'model_id', 'prediction_value'),
         coords={
-            'data_id': range(len(dataset)),
+            'data_id': range(id_offset, len(dataset) + id_offset),
             'model_id': list(models.keys()),
             'prediction_value': [
                 'negative_prediction',
@@ -160,7 +161,7 @@ def compute_ensemble_statistics(predictions: xr.DataArray):
 
 # Compute the thresholded predictions given an array of predictions
 def compute_thresholded_predictions(input_stats: xr.DataArray):
-    quantiles = np.linspace(0.05, 0.95, 19) * 100
+    quantiles = np.linspace(0.00, 1.00, 21) * 100
     metrics = ['accuracy', 'f1']
     statistics = ['stdev', 'entropy', 'confidence']
 
@@ -218,6 +219,12 @@ def compute_metric(arr, metric):
         return met.F1(
             arr.loc[{'statistic': 'predicted'}], arr.loc[{'statistic': 'actual'}]
         )
+    elif metric == 'ece':
+        true_labels = arr.loc[{'statistic': 'actual'}].values
+        predicted_labels = arr.loc[{'statistic': 'predicted'}].values
+        confidences = arr.loc[{'statistic': 'confidence'}].values
+
+        return calculate_ece_stats(confidences, predicted_labels, true_labels)
 
     else:
         raise ValueError('Invalid metric: ' + metric)
@@ -373,7 +380,9 @@ def compute_individual_statistics(predictions: xr.DataArray):
         },
     )
 
-    for data_id in predictions.data_id:
+    for data_id in tqdm(
+        predictions.data_id, total=len(predictions.data_id), unit='images'
+    ):
         for model_id in predictions.model_id:
             data = predictions.loc[{'data_id': data_id, 'model_id': model_id}]
             mean = data[0:2]
@@ -416,7 +425,9 @@ def compute_individual_thresholds(input_stats: xr.DataArray):
         },
     )
 
-    for model_id in input_stats.model_id:
+    for model_id in tqdm(
+        input_stats.model_id, total=len(input_stats.model_id), unit='models'
+    ):
         for statistic in statistics:
             # First, we must compute the quantiles for the statistic
             quantile_values = np.percentile(
@@ -550,8 +561,9 @@ def graph_all_individual_thresholded_predictions(
 
 # Calculate statistics of subsets of models for sensitivity analysis
 def calculate_subset_statistics(predictions: xr.DataArray):
-    # Calculate subsets for 1-50 models
-    subsets = range(1, len(predictions.model_id) + 1)
+    # Calculate subsets for 1-49 models
+    subsets = range(1, len(predictions.model_id))
+
     zeros = np.zeros(
         (len(predictions.data_id), len(subsets), 7)
     )  # Include stdev, but for 1 models set to NaN
@@ -574,7 +586,9 @@ def calculate_subset_statistics(predictions: xr.DataArray):
         },
     )
 
-    for data_id in predictions.data_id:
+    for data_id in tqdm(
+        predictions.data_id, total=len(predictions.data_id), unit='images'
+    ):
         for subset in subsets:
             data = predictions.sel(
                 data_id=data_id, model_id=predictions.model_id[:subset]
@@ -603,22 +617,24 @@ def calculate_subset_statistics(predictions: xr.DataArray):
 # Calculate Accuracy, F1 and ECE for subset stats - sensityvity analysis
 def calculate_sensitivity_analysis(subset_stats: xr.DataArray):
     subsets = subset_stats.model_count
-    stats = ['accuracy', 'f1']
+    stats = ['accuracy', 'f1', 'ece']
 
     zeros = np.zeros((len(subsets), len(stats)))
 
     sens_analysis = xr.DataArray(
         zeros,
         dims=('model_count', 'statistic'),
-        coords={'model_count': subsets, 'statistic': ['accuracy', 'f1']},
+        coords={'model_count': subsets, 'statistic': stats},
     )
 
-    for subset in subsets:
+    for subset in tqdm(subsets, total=len(subsets), unit='model subsets'):
+
         data = subset_stats.sel(model_count=subset)
-        acc = compute_metric(data, 'accuracy')
-        f1 = compute_metric(data, 'f1')
+        acc = compute_metric(data, 'accuracy').item()
+        f1 = compute_metric(data, 'f1').item()
+        ece = compute_metric(data, 'ece').item()
 
-        sens_analysis.loc[{'model_count': subset}] = [acc, f1]
+        sens_analysis.loc[{'model_count': subset.item()}] = [acc, f1, ece]
 
     return sens_analysis
 
@@ -648,18 +664,12 @@ def calculate_overall_stats(ensemble_statistics: xr.DataArray):
 
 
 # https://towardsdatascience.com/expected-calibration-error-ece-a-step-by-step-visual-explanation-with-python-code-c3e9aa12937d
-def calculate_ece_stats(statistics, bins=10):
+def calculate_ece_stats(confidences, predicted_labels, true_labels, bins=10):
     bin_boundaries = np.linspace(0, 1, bins + 1)
     bin_lowers = bin_boundaries[:-1]
     bin_uppers = bin_boundaries[1:]
 
-    confidences = ((statistics.sel(statistic='mean').values) - 0.5) * 2
-    accuracies = statistics.sel(statistic='correct').values
-
     ece = np.zeros(1)
-    bin_accuracies = xr.DataArray(
-        np.zeros(bins), dims=('lower_bound'), coords={'lower_bound': bin_lowers}
-    )
 
     for bin_lower, bin_upper in zip(bin_lowers, bin_uppers):
         in_bin = np.logical_and(
@@ -668,16 +678,13 @@ def calculate_ece_stats(statistics, bins=10):
         prob_in_bin = in_bin.mean()
 
         if prob_in_bin.item() > 0:
-            accuracy_in_bin = accuracies[in_bin].mean()
+            accuracy_in_bin = true_labels[in_bin].mean()
 
-            bin_accuracies.loc[{'lower_bound': bin_lower}]
             avg_confidence_in_bin = confidences[in_bin].mean()
-            ece += np.abs(avg_confidence_in_bin - accuracy_in_bin) * prob_in_bin
 
-    bin_accuracies.attrs['ece'] = ece
-    bin_accuracies.attrs['bin_number'] = bins
+            ece += np.abs(avg_confidence_in_bin - accuracy_in_bin) * prob_in_bin
 
-    return bin_accuracies
+    return ece
 
 
 def plot_ece_graph(ece_stats, title, xlabel, ylabel, save_path):
@@ -788,6 +795,7 @@ def main():
     print('Individual Thresholded Predictions Graphed')
 
     # Compute subset statistics and graph
+    print('Computing Sensitivity Analysis...')
     subset_stats = calculate_subset_statistics(predictions)
     sens_analysis = calculate_sensitivity_analysis(subset_stats)
     graph_sensitivity_analysis(
@@ -798,6 +806,14 @@ def main():
         '# of Models',
         'Accuracy',
     )
+    graph_sensitivity_analysis(
+        sens_analysis,
+        'ece',
+        f'{V4_PATH}/sens_analysis_ece.png',
+        'Sensitivity Analysis of ECE vs. # of Models',
+        '# of Models',
+        'ECE',
+    )
     print(sens_analysis.sel(statistic='accuracy'))
     print(calculate_overall_stats(ensemble_statistics))
 

xarray_images.py

@@ -0,0 +1,21 @@
+# I need a couple of images of positive and negative class members for my poster
+
+import xarray as xr
+import numpy as np
+import matplotlib.pyplot as plt
+import torch as th
+
+
+# Dataset is a torch dataset
+def get_image(dataset, idx):
+    img = dataset[idx][0].numpy()
+    return img
+
+
+def plot_image(img, path):
+    plt.imshow(img)
+    plt.savefig(path)
+
+
+def get_random_positive_image(dataset):
+    idx =

xarray_sensitivity.py

@@ -0,0 +1,98 @@
+# For this project, we want to calculate the accuracy for many different numbers of models and model selections
+# We cannot calculate every possible permutation and combination of models, so we will use a random selection
+
+
+import math
+import itertools as it
+
+import xarray as xr
+import numpy as np
+import matplotlib.pyplot as plt
+
+import torch as th
+import random as rand
+
+import threshold_xarray as txr
+import os
+
+
+# Generate a selection of combinations given an iterable, combination length, and number of combinations
+# This checks the number of possible combinations and compares it to the requested number of combinations
+# If the number of requested combinations is more than the possible number of combinations, it wil error
+# If it is less, it will generate all possible combinations and select the requested number of combinations
+# If it is much less, it will randomly generate the requested number of combinations and check that they are unique
+# If it is equal, it will generate and return all possible combinations\
+def get_combinations(iterable, r, n_combinations):
+    possible_combinations = math.comb(len(iterable), r)
+
+    if n_combinations < possible_combinations:
+        raise ValueError(
+            f'Number of requested combinations {n_combinations} of length {r} on set of length {len(iterable)} is less than the possible number of combinations {possible_combinations}'
+        )
+    elif n_combinations == possible_combinations:
+        return list(it.combinations(iterable, r))
+    else:
+        if n_combinations < possible_combinations / 5:
+            combinations = []
+            while len(combinations) < n_combinations:
+                combination = rand.sample(iterable, r)
+
+                if combination not in combinations:
+                    combinations.append(combination)
+            return combinations
+        else:
+            combinations = list(it.combinations(iterable, r))  # All possible combinations
+            return rand.sample(
+                combinations, n_combinations
+            )  # Randomly select n_combinations
+        
+
+# Now that we have a function to generate combinations, we can generate a list of 49 * 50 + 1 combinations
+# This will be a list of 2451 combinations of 50 models
+
+models = list(range(50))
+combos = {}
+for i in range(49):
+    combos[i] = get_combinations(models, i + 1, 50)
+
+combos[50] = [models]
+
+
+# Now that we have the list of combinations, we need the predictions
+print('Loading Config...')
+config = txr.load_config()
+ENSEMBLE_PATH = f"{config['paths']['model_output']}{config['ensemble']['name']}"
+V4_PATH = ENSEMBLE_PATH + '/v4'
+
+if not os.path.exists(V4_PATH):
+    os.makedirs(V4_PATH)
+print('Config Loaded')
+
+test_predictions = xr.open_dataarray(f'{V4_PATH}/test_predictions.nc')
+val_predictions = xr.open_dataarray(f'{V4_PATH}/val_predictions.nc')
+
+# Prune Data
+print('Pruning Data...')
+if config['operation']['exclude_blank_ids']:
+    excluded_data_ids = config['ensemble']['excluded_ids']
+    test_predictions = txr.prune_data(test_predictions, excluded_data_ids)
+    val_predictions = txr.prune_data(val_predictions, excluded_data_ids)
+
+# Concatenate Predictions
+predictions = xr.concat([test_predictions, val_predictions], dim='data_id')
+
+# Now that we have the list of predictions, we can calculate the accuracy and other stats for each combination
+# We will calculate the accuracy for each combination of models and save the results
+
+# Calculate the accuracy for each combination of models
+for num_models, model_combinations in combos.items():
+    print(f'Calculating Accuracy for {num_models} Models')
+    for i, model_combination in enumerate(model_combinations):
+        print(f'Calculating Accuracy for Combination {i} of {len(model_combinations)}')
+        model_predictions = predictions.sel(model=model_combination)
+        
+        # Calculate the accuracy
+        num_correct = 
+
+
+