Traumatic brain injury recovery prediction by harmonizing real brain CT and synthetic brain MRI: a pilot study.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches

[1] § Materials and methods › Statistical analysis ↔ Stage_2_CLS/5fold/analysis.ipynb, lines 270 to 397 · score 0.73 · receiver operating characteristic, DeLong, correlated, curves, ROC, AUCs · Report a mismatch
[2] § Materials and methods › GAN as harmonization technique for synthetic imaging rendering ↔ Stage_1_FPGAN/model.py, lines 22 to 59 · score 0.60 · residual blocks, domain information, network, layers, models · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 3,122 lines · 119 KB · no license · 1 match

# %%
%cd
# %% [markdown]
# # Binary classification
# %%
import logging
import os
import sys
from pathlib import Path
import pandas as pd
import monai
import argparse
import numpy as np
import torch
import torch.nn as nn
from torch.utils.data import DataLoader as _TorchDataLoader
from torch.utils.data import Dataset
from torch.utils.tensorboard import SummaryWriter
from monai.data import decollate_batch, CSVSaver, DataLoader
from monai.data.utils import list_data_collate, set_rnd, worker_init_fn
from monai.metrics import ROCAUCMetric
from monai.transforms import ScaleIntensityRanged, NormalizeIntensityd, Activations, AsDiscrete, Compose, LoadImaged, RandRotated, Resized, ScaleIntensityd, RandFlipd
from monai.networks.nets import resnet10, resnet18, resnet34, resnet50
from sklearn.model_selection import StratifiedKFold
from statistics import mean, mode
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, classification_report, accuracy_score, roc_auc_score
# %%
def calculate_scores(y_true, y_pred_act_all):
print("="*20)
print("Test Set Results:\n")
y_pred_act_all = np.array(y_pred_act_all)
sum5 = np.sum(np.array(y_pred_act_all), axis=0)
probs = sum5/5
ypred_soft_votes = np.argmax(probs, axis=1)
acc_score = accuracy_score(y_true, ypred_soft_votes)
probs_auc = y_pred_act_all[:,:,1]
sum5 = np.sum(np.array(probs_auc), axis=0)
probs_auc = sum5/5
auc_score = roc_auc_score(y_true, probs_auc)
cm = confusion_matrix(y_true, ypred_soft_votes)
# Calculate confusion matrix
tn, fp, fn, tp = cm.ravel()
# Calculate Precision
precision = precision_score(y_true, ypred_soft_votes)
# Calculate Recall (Sensitivity)
recall = recall_score(y_true, ypred_soft_votes)
# Calculate Specificity
specificity = tn / (tn + fp)
# Calculate F1 Score
f1 = f1_score(y_true, ypred_soft_votes)
return acc_score, auc_score, precision, recall, specificity, f1, cm, ypred_soft_votes
# %%
import sys
from mri_classification.scripts.sbatch.bi_5fold_solo_eval import main
from scipy.stats import mode
from sklearn.metrics import roc_auc_score
from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score
sys.argv = ['']
parser = argparse.ArgumentParser()
parser.add_argument('--modality', type=str, default='MR')
parser.add_argument('--opt', type=str, default='Adam')
parser.add_argument('--lr', type=float, default=1e-3)
parser.add_argument('--wdecay', type=float, default=1e-4)
config = parser.parse_args()
y_true, y_pred_act_all = main(config)
acc_score, auc_score, precision, recall, specificity, f1, cm, ypred_soft_votes = calculate_scores(y_true, y_pred_act_all)
# Print the results
print(f"Accuracy: {round(acc_score, 4)}")
print(f"AUC Score: {round(auc_score, 4)}")
print(f"Precision: {round(precision,4)}")
print(f"Recall (Sensitivity): {round(recall,4)}")
print(f"Specificity: {round(specificity, 4)}")
print(f"F1 Score: {round(f1, 4)}")
target_names = ['class 0+1', 'class 2']
print(classification_report(y_true, ypred_soft_votes, target_names=target_names))
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot(cmap = 'Blues')
plt.show()
# %%
import sys
from mri_classification.scripts.sbatch.bi_5fold_solo_eval import main
from scipy.stats import mode
from sklearn.metrics import roc_auc_score
from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score
sys.argv = ['']
parser = argparse.ArgumentParser()
parser.add_argument('--modality', type=str, default='CT')
parser.add_argument('--opt', type=str, default='Adam')
parser.add_argument('--lr', type=float, default=1e-5)
parser.add_argument('--wdecay', type=float, default=1e-6)
config = parser.parse_args()
y_true, y_pred_act_all = main(config)
acc_score, auc_score, precision, recall, specificity, f1, cm, ypred_soft_votes = calculate_scores(y_true, y_pred_act_all)
# Print the results
print(f"Accuracy: {round(acc_score, 4)}")
print(f"AUC Score: {round(auc_score, 4)}")
print(f"Precision: {round(precision,4)}")
print(f"Recall (Sensitivity): {round(recall,4)}")
print(f"Specificity: {round(specificity, 4)}")
print(f"F1 Score: {round(f1, 4)}")
target_names = ['class 0+1', 'class 2']
print(classification_report(y_true, ypred_soft_votes, target_names=target_names))
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot(cmap = 'Blues')
plt.show()
# %%
import sys
from mri_classification.scripts.sbatch.bi_5fold_ensemble_eval import main
from scipy.stats import mode
from sklearn.metrics import roc_auc_score
from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score
sys.argv = ['']
parser = argparse.ArgumentParser()
parser.add_argument('--ct_type', type=str, default='Real')
parser.add_argument('--mr_type', type=str, default='Real')
parser.add_argument('--opt', type=str, default='Adam')
parser.add_argument('--epochs', type=int, default=60)
parser.add_argument('--lr', type=float, default=1e-5)
parser.add_argument('--wdecay', type=float, default=1e-3)
config = parser.parse_args()
y_true, y_pred_act_all = main(config)
acc_score, auc_score, precision, recall, specificity, f1, cm, ypred_soft_votes = calculate_scores(y_true, y_pred_act_all)
# Print the results
print(f"Accuracy: {round(acc_score, 4)}")
print(f"AUC Score: {round(auc_score, 4)}")
print(f"Precision: {round(precision,4)}")
print(f"Recall (Sensitivity): {round(recall,4)}")
print(f"Specificity: {round(specificity, 4)}")
print(f"F1 Score: {round(f1, 4)}")
target_names = ['class 0+1', 'class 2']
print(classification_report(y_true, ypred_soft_votes, target_names=target_names))
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot(cmap = 'Blues')
plt.show()
# %%
import sys
from mri_classification.scripts.sbatch.bi_5fold_ensemble_eval import main
from scipy.stats import mode
from sklearn.metrics import roc_auc_score
from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score
sys.argv = ['']
parser = argparse.ArgumentParser()
parser.add_argument('--ct_type', type=str, default='Fake')
parser.add_argument('--mr_type', type=str, default='Real')
parser.add_argument('--opt', type=str, default='Adam')
parser.add_argument('--epochs', type=int, default=60)
parser.add_argument('--lr', type=float, default=1e-5)
parser.add_argument('--wdecay', type=float, default=1e-2)
config = parser.parse_args()
y_true, y_pred_act_all = main(config)
acc_score, auc_score, precision, recall, specificity, f1, cm, ypred_soft_votes = calculate_scores(y_true, y_pred_act_all)
# Print the results
print(f"Accuracy: {round(acc_score, 4)}")
print(f"AUC Score: {round(auc_score, 4)}")
print(f"Precision: {round(precision,4)}")
print(f"Recall (Sensitivity): {round(recall,4)}")
print(f"Specificity: {round(specificity, 4)}")
print(f"F1 Score: {round(f1, 4)}")
target_names = ['class 0+1', 'class 2']
print(classification_report(y_true, ypred_soft_votes, target_names=target_names))
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot(cmap = 'Blues')
plt.show()
# %%
import sys
from mri_classification.scripts.sbatch.bi_5fold_ensemble_eval import main
from scipy.stats import mode
from sklearn.metrics import roc_auc_score
from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score
sys.argv = ['']
parser = argparse.ArgumentParser()
parser.add_argument('--ct_type', type=str, default='Fake')
parser.add_argument('--mr_type', type=str, default='Real')
parser.add_argument('--opt', type=str, default='Adam')
parser.add_argument('--epochs', type=int, default=60)
parser.add_argument('--lr', type=float, default=1e-5)
parser.add_argument('--wdecay', type=float, default=1e-4)
config = parser.parse_args()
y_true, y_pred_act_all = main(config)
acc_score, auc_score, precision, recall, specificity, f1, cm, ypred_soft_votes = calculate_scores(y_true, y_pred_act_all)
# Print the results
print(f"Accuracy: {round(acc_score, 4)}")
print(f"AUC Score: {round(auc_score, 4)}")
print(f"Precision: {round(precision,4)}")
print(f"Recall (Sensitivity): {round(recall,4)}")
print(f"Specificity: {round(specificity, 4)}")
print(f"F1 Score: {round(f1, 4)}")
target_names = ['class 0+1', 'class 2']
print(classification_report(y_true, ypred_soft_votes, target_names=target_names))
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot(cmap = 'Blues')
plt.show()
# %% [markdown]
# ### DeLong Test
# %%
import pandas as pd
import numpy as np
import scipy.stats
# AUC comparison adapted from
# https://github.com/Netflix/vmaf/
def compute_midrank(x):
"""Computes midranks.
Args:
x - a 1D numpy array
Returns:
array of midranks
"""
J = np.argsort(x)
Z = x[J]
N = len(x)
T = np.zeros(N, dtype=float)
i = 0
while i < N:
j = i
while j < N and Z[j] == Z[i]:
j += 1
T[i:j] = 0.5*(i + j - 1)
i = j
T2 = np.empty(N, dtype=float)
# Note(kazeevn) +1 is due to Python using 0-based indexing
# instead of 1-based in the AUC formula in the paper
T2[J] = T + 1
return T2
def fastDeLong(predictions_sorted_transposed, label_1_count):
"""
The fast version of DeLong's method for computing the covariance of
unadjusted AUC.
Args:
predictions_sorted_transposed: a 2D numpy.array[n_classifiers, n_examples]
sorted such as the examples with label "1" are first
Returns:
(AUC value, DeLong covariance)
Reference:
@article{sun2014fast,
title={Fast Implementation of DeLong's Algorithm for
Comparing the Areas Under Correlated Receiver Operating Characteristic Curves},
author={Xu Sun and Weichao Xu},
journal={IEEE Signal Processing Letters},
volume={21},
number={11},
pages={1389--1393},
year={2014},
publisher={IEEE}
}
"""
# Short variables are named as they are in the paper
m = label_1_count
n = predictions_sorted_transposed.shape[1] - m
positive_examples = predictions_sorted_transposed[:, :m]
negative_examples = predictions_sorted_transposed[:, m:]
k = predictions_sorted_transposed.shape[0]
tx = np.empty([k, m], dtype=float)
ty = np.empty([k, n], dtype=float)
tz = np.empty([k, m + n], dtype=np.float)
for r in range(k):
tx[r, :] = compute_midrank(positive_examples[r, :])
ty[r, :] = compute_midrank(negative_examples[r, :])
tz[r, :] = compute_midrank(predictions_sorted_transposed[r, :])
aucs = tz[:, :m].sum(axis=1) / m / n - float(m + 1.0) / 2.0 / n
v01 = (tz[:, :m] - tx[:, :]) / n
v10 = 1.0 - (tz[:, m:] - ty[:, :]) / m
sx = np.cov(v01)
sy = np.cov(v10)
delongcov = sx / m + sy / n
return aucs, delongcov
def calc_pvalue(aucs, sigma):
"""Computes log(10) of p-values.
Args:
aucs: 1D array of AUCs
sigma: AUC DeLong covariances
Returns:
log10(pvalue)
"""
l = np.array([[1, -1]])
z = np.abs(np.diff(aucs)) / np.sqrt(np.dot(np.dot(l, sigma), l.T))
return np.log10(2) + scipy.stats.norm.logsf(z, loc=0, scale=1) / np.log(10)
def compute_ground_truth_statistics(ground_truth):
assert np.array_equal(np.unique(ground_truth), [0, 1])
order = (-ground_truth).argsort()
label_1_count = int(ground_truth.sum())
return order, label_1_count
def delong_roc_variance(ground_truth, predictions):
"""
Computes ROC AUC variance for a single set of predictions
Args:
ground_truth: np.array of 0 and 1
predictions: np.array of floats of the probability of being class 1
"""
order, label_1_count = compute_ground_truth_statistics(ground_truth)
predictions_sorted_transposed = predictions[np.newaxis, order]
aucs, delongcov = fastDeLong(predictions_sorted_transposed, label_1_count)
assert len(aucs) == 1, "There is a bug in the code, please forward this to the developers"
return aucs[0], delongcov
def delong_roc_test(ground_truth, predictions_one, predictions_two):
"""
Computes log(p-value) for hypothesis that two ROC AUCs are different
Args:
ground_truth: np.array of 0 and 1
predictions_one: predictions of the first model,
np.array of floats of the probability of being class 1
predictions_two: predictions of the second model,
np.array of floats of the probability of being class 1
"""
order, label_1_count = compute_ground_truth_statistics(ground_truth)
predictions_sorted_transposed = np.vstack((predictions_one, predictions_two))[:, order]
aucs, delongcov = fastDeLong(predictions_sorted_transposed, label_1_count)
return calc_pvalue(aucs, delongcov)
# %%
import sys
from mri_classification.scripts.sbatch.bi_5fold_ensemble_eval import main
from scipy.stats import mode
from sklearn.metrics import roc_auc_score
from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score
# Fake CT + Real MR
sys.argv = ['']
parser = argparse.ArgumentParser()
parser.add_argument('--ct_type', type=str, default='Fake')
parser.add_argument('--mr_type', type=str, default='Real')
parser.add_argument('--opt', type=str, default='Adam')
parser.add_argument('--epochs', type=int, default=60)
parser.add_argument('--lr', type=float, default=1e-5)
parser.add_argument('--wdecay', type=float, default=1e-4)
config = parser.parse_args()
y_true, y_pred_act_all = main(config)
acc_score, auc_score, precision, recall, specificity, f1, cm, fakect_ypred = calculate_scores(y_true, y_pred_act_all)
# Real MR only
from mri_classification.scripts.sbatch.bi_5fold_solo_eval import main
sys.argv = ['']
parser = argparse.ArgumentParser()
parser.add_argument('--modality', type=str, default='MR')
parser.add_argument('--opt', type=str, default='Adam')
parser.add_argument('--lr', type=float, default=1e-3)
parser.add_argument('--wdecay', type=float, default=1e-4)
config = parser.parse_args()
y_true, y_pred_act_all = main(config)
acc_score, auc_score, precision, recall, specificity, f1, cm, mr_ypred = calculate_scores(y_true, y_pred_act_all)
# %%
import pandas as pd
import numpy as np
import scipy.stats
# AUC comparison adapted from
# https://github.com/Netflix/vmaf/
def compute_midrank(x):
"""Computes midranks.
Args:
x - a 1D numpy array
Returns:
array of midranks
"""
J = np.argsort(x)
Z = x[J]
N = len(x)
T = np.zeros(N, dtype=float)
i = 0
while i < N:
j = i
while j < N and Z[j] == Z[i]:
j += 1
T[i:j] = 0.5*(i + j - 1)
i = j
T2 = np.empty(N, dtype=float)
# Note(kazeevn) +1 is due to Python using 0-based indexing
# instead of 1-based in the AUC formula in the paper
T2[J] = T + 1
return T2
def fastDeLong(predictions_sorted_transposed, label_1_count):
"""
The fast version of DeLong's method for computing the covariance of
unadjusted AUC.
Args:
predictions_sorted_transposed: a 2D numpy.array[n_classifiers, n_examples]
sorted such as the examples with label "1" are first
Returns:
(AUC value, DeLong covariance)
Reference:
@article{sun2014fast,
title={Fast Implementation of DeLong's Algorithm for
Comparing the Areas Under Correlated Receiver Operating Characteristic Curves},
author={Xu Sun and Weichao Xu},
journal={IEEE Signal Processing Letters},
volume={21},
number={11},
pages={1389--1393},
year={2014},
publisher={IEEE}
}
"""
# Short variables are named as they are in the paper
m = label_1_count
n = predictions_sorted_transposed.shape[1] - m
positive_examples = predictions_sorted_transposed[:, :m]
negative_examples = predictions_sorted_transposed[:, m:]
k = predictions_sorted_transposed.shape[0]
tx = np.empty([k, m], dtype=float)
ty = np.empty([k, n], dtype=float)
tz = np.empty([k, m + n], dtype=float)
for r in range(k):
tx[r, :] = compute_midrank(positive_examples[r, :])
ty[r, :] = compute_midrank(negative_examples[r, :])
tz[r, :] = compute_midrank(predictions_sorted_transposed[r, :])
aucs = tz[:, :m].sum(axis=1) / m / n - float(m + 1.0) / 2.0 / n
v01 = (tz[:, :m] - tx[:, :]) / n
v10 = 1.0 - (tz[:, m:] - ty[:, :]) / m
sx = np.cov(v01)
sy = np.cov(v10)
delongcov = sx / m + sy / n
return aucs, delongcov
def calc_pvalue(aucs, sigma):
"""Computes log(10) of p-values.
Args:
aucs: 1D array of AUCs
sigma: AUC DeLong covariances
Returns:
log10(pvalue)
"""
l = np.array([[1, -1]])
z = np.abs(np.diff(aucs)) / np.sqrt(np.dot(np.dot(l, sigma), l.T))
return np.log10(2) + scipy.stats.norm.logsf(z, loc=0, scale=1) / np.log(10)
def compute_ground_truth_statistics(ground_truth):
assert np.array_equal(np.unique(ground_truth), [0, 1])
order = (-ground_truth).argsort()
label_1_count = int(ground_truth.sum())
return order, label_1_count
def delong_roc_variance(ground_truth, predictions):
"""
Computes ROC AUC variance for a single set of predictions
Args:
ground_truth: np.array of 0 and 1
predictions: np.array of floats of the probability of being class 1
"""
order, label_1_count = compute_ground_truth_statistics(ground_truth)
predictions_sorted_transposed = predictions[np.newaxis, order]
aucs, delongcov = fastDeLong(predictions_sorted_transposed, label_1_count)
assert len(aucs) == 1, "There is a bug in the code, please forward this to the developers"
return aucs[0], delongcov
def delong_roc_test(ground_truth, predictions_one, predictions_two):
"""
Computes log(p-value) for hypothesis that two ROC AUCs are different
Args:
ground_truth: np.array of 0 and 1
predictions_one: predictions of the first model,
np.array of floats of the probability of being class 1
predictions_two: predictions of the second model,
np.array of floats of the probability of being class 1
"""
order, label_1_count = compute_ground_truth_statistics(ground_truth)
predictions_sorted_transposed = np.vstack((predictions_one, predictions_two))[:, order]
aucs, delongcov = fastDeLong(predictions_sorted_transposed, label_1_count)
return calc_pvalue(aucs, delongcov)
# %%
pvalue = delong_roc_test(y_true, fakect_ypred, mr_ypred)
pvalue
# %%
p_value_test = np.exp(np.log(10)*pvalue)
p_value_test
# %%
# %%
# %%
import pandas as pd
import numpy as np
import scipy.stats
from scipy import stats
# AUC comparison adapted from
# https://github.com/Netflix/vmaf/
def compute_midrank(x):
"""Computes midranks.
Args:
x - a 1D numpy array
Returns:
array of midranks
"""
J = np.argsort(x)
Z = x[J]
N = len(x)
T = np.zeros(N, dtype=float)
i = 0
while i < N:
j = i
while j < N and Z[j] == Z[i]:
j += 1
T[i:j] = 0.5*(i + j - 1)
i = j
T2 = np.empty(N, dtype=float)
# Note(kazeevn) +1 is due to Python using 0-based indexing
# instead of 1-based in the AUC formula in the paper
T2[J] = T + 1
return T2
def compute_midrank_weight(x, sample_weight):
"""Computes midranks.
Args:
x - a 1D numpy array
Returns:
array of midranks
"""
J = np.argsort(x)
Z = x[J]
cumulative_weight = np.cumsum(sample_weight[J])
N = len(x)
T = np.zeros(N, dtype=float)
i = 0
while i < N:
j = i
while j < N and Z[j] == Z[i]:
j += 1
T[i:j] = cumulative_weight[i:j].mean()
i = j
T2 = np.empty(N, dtype=float)
T2[J] = T
return T2
def fastDeLong(predictions_sorted_transposed, label_1_count, sample_weight=None):
if sample_weight is None:
return fastDeLong_no_weights(predictions_sorted_transposed, label_1_count)
else:
return fastDeLong_weights(predictions_sorted_transposed, label_1_count, sample_weight)
def fastDeLong_weights(predictions_sorted_transposed, label_1_count, sample_weight):
"""
The fast version of DeLong's method for computing the covariance of
unadjusted AUC.
Args:
predictions_sorted_transposed: a 2D numpy.array[n_classifiers, n_examples]
sorted such as the examples with label "1" are first
Returns:
(AUC value, DeLong covariance)
Reference:
@article{sun2014fast,
title={Fast Implementation of DeLong's Algorithm for
Comparing the Areas Under Correlated Receiver Oerating Characteristic Curves},
author={Xu Sun and Weichao Xu},
journal={IEEE Signal Processing Letters},
volume={21},
number={11},
pages={1389--1393},
year={2014},
publisher={IEEE}
}
"""
# Short variables are named as they are in the paper
m = label_1_count
n = predictions_sorted_transposed.shape[1] - m
positive_examples = predictions_sorted_transposed[:, :m]
negative_examples = predictions_sorted_transposed[:, m:]
k = predictions_sorted_transposed.shape[0]
tx = np.empty([k, m], dtype=float)
ty = np.empty([k, n], dtype=float)
tz = np.empty([k, m + n], dtype=float)
for r in range(k):
tx[r, :] = compute_midrank_weight(positive_examples[r, :], sample_weight[:m])
ty[r, :] = compute_midrank_weight(negative_examples[r, :], sample_weight[m:])
tz[r, :] = compute_midrank_weight(predictions_sorted_transposed[r, :], sample_weight)
total_positive_weights = sample_weight[:m].sum()
total_negative_weights = sample_weight[m:].sum()
pair_weights = np.dot(sample_weight[:m, np.newaxis], sample_weight[np.newaxis, m:])
total_pair_weights = pair_weights.sum()
aucs = (sample_weight[:m]*(tz[:, :m] - tx)).sum(axis=1) / total_pair_weights
v01 = (tz[:, :m] - tx[:, :]) / total_negative_weights
v10 = 1. - (tz[:, m:] - ty[:, :]) / total_positive_weights
sx = np.cov(v01)
sy = np.cov(v10)
delongcov = sx / m + sy / n
return aucs, delongcov
def fastDeLong_no_weights(predictions_sorted_transposed, label_1_count):
"""
The fast version of DeLong's method for computing the covariance of
unadjusted AUC.
Args:
predictions_sorted_transposed: a 2D numpy.array[n_classifiers, n_examples]
sorted such as the examples with label "1" are first
Returns:
(AUC value, DeLong covariance)
Reference:
@article{sun2014fast,
title={Fast Implementation of DeLong's Algorithm for
Comparing the Areas Under Correlated Receiver Oerating
Characteristic Curves},
author={Xu Sun and Weichao Xu},
journal={IEEE Signal Processing Letters},
volume={21},
number={11},
pages={1389--1393},
year={2014},
publisher={IEEE}
}
"""
# Short variables are named as they are in the paper
m = label_1_count
n = predictions_sorted_transposed.shape[1] - m
positive_examples = predictions_sorted_transposed[:, :m]
negative_examples = predictions_sorted_transposed[:, m:]
k = predictions_sorted_transposed.shape[0]
tx = np.empty([k, m], dtype=float)
ty = np.empty([k, n], dtype=float)
tz = np.empty([k, m + n], dtype=float)
for r in range(k):
tx[r, :] = compute_midrank(positive_examples[r, :])
ty[r, :] = compute_midrank(negative_examples[r, :])
tz[r, :] = compute_midrank(predictions_sorted_transposed[r, :])
aucs = tz[:, :m].sum(axis=1) / m / n - float(m + 1.0) / 2.0 / n
v01 = (tz[:, :m] - tx[:, :]) / n
v10 = 1.0 - (tz[:, m:] - ty[:, :]) / m
sx = np.cov(v01)
sy = np.cov(v10)
delongcov = sx / m + sy / n
return aucs, delongcov
def calc_pvalue(aucs, sigma):
"""Computes log(10) of p-values.
Args:
aucs: 1D array of AUCs
sigma: AUC DeLong covariances
Returns:
log10(pvalue)
"""
l = np.array([[1, -1]])
z = np.abs(np.diff(aucs)) / (np.sqrt(np.dot(np.dot(l, sigma), l.T)) + 1e-8)
pvalue = 2 * (1 - scipy.stats.norm.cdf(np.abs(z)))
# print(10**(np.log10(2) + scipy.stats.norm.logsf(z, loc=0, scale=1) / np.log(10)))
return pvalue
def compute_ground_truth_statistics(ground_truth, sample_weight=None):
assert np.array_equal(np.unique(ground_truth), [0, 1])
order = (-ground_truth).argsort()
label_1_count = int(ground_truth.sum())
if sample_weight is None:
ordered_sample_weight = None
else:
ordered_sample_weight = sample_weight[order]
return order, label_1_count, ordered_sample_weight
def delong_roc_variance(ground_truth, predictions):
"""
Computes ROC AUC variance for a single set of predictions
Args:
ground_truth: np.array of 0 and 1
predictions: np.array of floats of the probability of being class 1
"""
sample_weight = None
order, label_1_count, ordered_sample_weight = compute_ground_truth_statistics(
ground_truth, sample_weight)
predictions_sorted_transposed = predictions[np.newaxis, order]
aucs, delongcov = fastDeLong(predictions_sorted_transposed, label_1_count)
assert len(aucs) == 1, "There is a bug in the code, please forward this to the developers"
return aucs[0], delongcov
def delong_roc_test(ground_truth, predictions_one, predictions_two):
"""
Computes log(p-value) for hypothesis that two ROC AUCs are different
Args:
ground_truth: np.array of 0 and 1
predictions_one: predictions of the first model,
np.array of floats of the probability of being class 1
predictions_two: predictions of the second model,
np.array of floats of the probability of being class 1
"""
sample_weight = None
order, label_1_count,ordered_sample_weight = compute_ground_truth_statistics(ground_truth)
predictions_sorted_transposed = np.vstack((predictions_one, predictions_two))[:, order]
aucs, delongcov = fastDeLong(predictions_sorted_transposed, label_1_count,sample_weight)
return calc_pvalue(aucs, delongcov)
# def delong_roc_ci(y_true,y_pred):
# aucs, auc_cov = delong_roc_variance(y_true, y_pred)
# auc_std = np.sqrt(auc_cov)
# lower_upper_q = np.abs(np.array([0, 1]) - (1 - alpha) / 2)
# ci = stats.norm.ppf(
# lower_upper_q,
# loc=aucs,
# scale=auc_std)
# ci[ci > 1] = 1
# return aucs,ci
def delong_roc_ci(y_true,y_pred):
aucs, auc_cov = delong_roc_variance(y_true, y_pred)
auc_std = np.sqrt(auc_cov)
lower_upper_q = np.abs(np.array([0, 1]) - (1 - alpha) / 2)
ci = stats.norm.ppf(
lower_upper_q,
loc=aucs,
scale=auc_std)
ci[ci > 1] = 1
return aucs,ci
def get_95CI(y_true,y_pred_1):
"""
Return the 95% CI and AUC of prediction
Args:
labels: array (n,) the ground truth
scores1: array (n,) the predicted probability
"""
alpha = .95
auc_1, auc_cov_1 = delong_roc_variance(y_true, y_pred_1)
auc_std = np.sqrt(auc_cov_1)
lower_upper_q = np.abs(np.array([0, 1]) - (1 - alpha) / 2)
# 95% CI
ci = stats.norm.ppf(
lower_upper_q,
loc=auc_1,
scale=auc_std)
ci[ci > 1] = 1
return ci,auc_1
# threshold
from sklearn import metrics
def get_optimal_threshold(labels,y_pred_1):
'''
get the threshold according to youden index
Args:
labels:<numpy.ndarray> (n,) groundtruth
y_pred_1:<numpy.ndarray> (n,) predicted probabilities
'''
fpr, tpr, thresholds = metrics.roc_curve(labels,y_pred_1)
optimal_index = np.argmax(+tpr-fpr)
optimal_thresholds = thresholds[optimal_index]
# print(thresholds)
# print(optimal_index)
# print(fpr)
return optimal_thresholds
def get_metric(y_true, y_prob,threshold,verbose = True):
'''
Return the commanly used metric value according to the given threshold
Args:
y_true:<numpy.ndarray> (n,) groundtruth
y_prob:<numpy.ndarray> (n,) predicted probabilities
threshold: <float>
Return:
scores: <dict> commanly used metric
'''
scores = {}
y_pred = (y_prob>=(threshold-1E-4)).astype(int)
# print report
target_names = ['class 0', 'class 1']
text = metrics.classification_report(y_true, y_pred, target_names=target_names)
conf_mat=pd.crosstab(y_true, y_pred,rownames=['label'],colnames=['pre'])
if verbose:
print(conf_mat)
print(text)
# accuracy
scores['accuracy'] = metrics.accuracy_score(y_true, y_pred)
# precision
try:
scores['PPV'] = metrics.precision_score(y_true, y_pred)
except:
scores['PPV'] = None
#NPV
try:
scores['NPV'] = conf_mat[0][0]/(conf_mat[0][0]+conf_mat[0][1])
except:
scores['NPV'] = None
# recall
scores['recall (sensitivity)'] = metrics.recall_score(y_true, y_pred)
scores['recall_neg (specificity)'] = metrics.recall_score(y_true==0, y_pred==0)
# F1-score
scores['f1_score'] = metrics.f1_score(y_true, y_pred)
# ROC/AUC
scores['AUC95%CI'], scores['AUC'] = get_95CI(y_true,y_pred)
return scores
# %%
pvalue = delong_roc_test(y_true, fakect_ypred, mr_ypred)
print('p_value:', pvalue)
# %%
# %%
# %%
# %% [markdown]
# # Three-class classification
# %%
import logging
import os
import sys
from pathlib import Path
import pandas as pd
import monai
import numpy as np
import torch
import torch.nn as nn
from torch.utils.data import DataLoader as _TorchDataLoader
from torch.utils.data import Dataset
from torch.utils.tensorboard import SummaryWriter
from monai.data import decollate_batch, CSVSaver, DataLoader
from monai.data.utils import list_data_collate, set_rnd, worker_init_fn
from monai.metrics import ROCAUCMetric
from monai.transforms import ScaleIntensityRanged, NormalizeIntensityd, Activations, AsDiscrete, Compose, LoadImaged, RandRotated, Resized, ScaleIntensityd, RandFlipd
from monai.networks.nets import resnet10, resnet18, resnet34, resnet50
from sklearn.model_selection import StratifiedKFold
from statistics import mean, mode
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, classification_report, accuracy_score, roc_auc_score
class CustomDataset(monai.data.Dataset):
def __init__(self, d1, d2):
self.d1 = d1
self.d2 = d2
def __getitem__(self, idx):
dict1 = self.d1.__getitem__(idx)
image1, label1 = dict1["img"], dict1["label"]
dict2 = self.d2.__getitem__(idx)
image2, label2 = dict2["img"], dict2["label"]
assert label1==label2
dict_1 = dict()
dict_1["img"] = image1
dict_1["label"] = label1
dict_2 = dict()
dict_2["img"] = image2
dict_2["label"] = label2
return dict_1, dict_2
def __len__(self):
return len(self.d1)
class MyEnsemble(nn.Module):
def __init__(self, modelA, modelB):
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
super(MyEnsemble, self).__init__()
self.modelA = modelA
self.modelB = modelB
# Remove last linear layer
self.modelA.fc = nn.Identity()
self.modelB.fc = nn.Identity()
# Create new classifier
self.mlp1 = nn.Linear(1024,256).to(device)
self.mlp2 = nn.Linear(256,32).to(device)
self.classifier = nn.Linear(32,3).to(device)
def forward(self, i1, i2):
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
x1 = self.modelA(i1) # clone to make sure x is not changed by inplace methods
x1 = x1.view(x1.size(0), -1).to(device)
x2 = self.modelB(i2)
x2 = x2.view(x2.size(0), -1).to(device)
x = torch.cat((x1, x2), dim=1).to(device)
x = nn.functional.relu(self.mlp1(x)).to(device)
x = nn.functional.relu(self.mlp2(x)).to(device)
x = self.classifier(x).to(device)
return x
def load_data(train_images, val_images, train_labels, val_labels, test_images, test_labels):
train_labels_array = np.array(train_labels, dtype=np.int64)
val_labels_array = np.array(val_labels, dtype=np.int64)
test_labels_array = np.array(test_labels, dtype=np.int64)
#loading the CT images
ct_path = Path('/scratch/ajoshi83/Data/Reg_CT')
#ct_path = Path('/scratch/ajoshi83/generated_ct')
#ct_path = Path('/data/amciilab/ajoshi83/generated_ct2')
ct_train_images_path = [os.path.join(ct_path, f) for f in train_images]
ct_val_images_path = [os.path.join(ct_path, f) for f in val_images]
ct_test_images_path = [os.path.join(ct_path, f) for f in test_images]
ct_train_files = [{"img": img, "label": label} for img, label in zip(ct_train_images_path, train_labels_array)]
ct_val_files = [{"img": img, "label": label} for img, label in zip(ct_val_images_path, val_labels_array)]
ct_test_files = [{"img": img, "label": label} for img, label in zip(ct_test_images_path, test_labels_array)]
#loading the MR images
mr_path = Path('/scratch/ajoshi83/Data/Reg_MR')
#mr_path = Path('/data/amciilab/ajoshi83/generated_mr2')
mr_train_images_path = [os.path.join(mr_path, f) for f in train_images]
mr_val_images_path = [os.path.join(mr_path, f) for f in val_images]
mr_test_images_path = [os.path.join(mr_path, f) for f in test_images]
mr_train_files = [{"img": img, "label": label} for img, label in zip(mr_train_images_path, train_labels_array)]
mr_val_files = [{"img": img, "label": label} for img, label in zip(mr_val_images_path, val_labels_array)]
mr_test_files = [{"img": img, "label": label} for img, label in zip(mr_test_images_path, test_labels_array)]
return ct_train_files, ct_val_files, mr_train_files, mr_val_files, ct_test_files, mr_test_files
#monai.config.print_config()
# old_stdout = sys.stdout
# log_file = open("mri_classification/logs/5fold_eval_fakemr.log","w")
# sys.stdout = log_file
print("this will be written to message.log")
# file_handler = logging.FileHandler(filename='logs/tmp.log')
# #stdout_handler = logging.StreamHandler(stream=sys.stdout)
# handlers = [file_handler]
# logging.basicConfig(
# level=logging.DEBUG,
# format='[%(asctime)s] {%(filename)s:%(lineno)d} %(levelname)s - %(message)s',
# handlers=handlers
# )
# logger = logging.getLogger('LOGGER_NAME')
# edit the path accordingly
gose = pd.read_csv('mri_classification/Pilot_169_train_val_test_split.csv')
X = []
Y = []
test_subjects = []
test_labels = []
torch.cuda.empty_cache()
for i in range(gose.shape[0]):
subj = gose['Main.GUID'][i]
subj_id = str(subj)[4:]
name = str(subj_id) + '.nii'
if gose['Set'][i]=='Train':
X.append(name)
Y.append(gose['Class'][i])
elif gose['Set'][i]=='Val':
X.append(name)
Y.append(gose['Class'][i])
elif gose['Set'][i]=='Test':
test_subjects.append(name)
test_labels.append(gose['Class'][i])
else:
print("Unknown Set: ", gose['Set'][i])
print("Total Subjects for 5-Fold (Train+Val):", len(X))
print("Test Subjects (kept separate):", len(test_subjects))
print("5-Fold Subjects Information:\n")
print("Subjects of Class 0: ", Y.count(0))
print("Subjects of Class 1: ", Y.count(1))
print("Subjects of Class 2: ", Y.count(2))
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=1)
fold_number = 1
lr =1e-5
decay = 1e-5
X = np.array(X)
Y = np.array(Y)
train_acc, train_auc = [], []
val_acc, val_auc = [], []
test_acc, test_auc = [], []
y_pred_test = []
y_pred_auc = []
for train_index, test_index in skf.split(X, Y):
print("Starting Fold {}..!".format(fold_number))
x_train_fold, x_test_fold = X[train_index], X[test_index]
y_train_fold, y_test_fold = Y[train_index], Y[test_index]
y_train_fold = list(y_train_fold)
y_test_fold = list(y_test_fold)
x_train_fold = list(x_train_fold)
x_test_fold = list(x_test_fold)
# print("\n")
# print("Fold {} statistics:\n".format(fold_number))
# print("Train Subjects: {}".format(len(y_train_fold)))
# print("Subjects of Class 0: ", y_train_fold.count(0))
# print("Subjects of Class 1: ", y_train_fold.count(1))
# print("Subjects of Class 2: {}\n".format(y_train_fold.count(2)))
# print("Val Subjects: {}".format(len(y_test_fold)))
# print("Subjects of Class 0: ", y_test_fold.count(0))
# print("Subjects of Class 1: ", y_test_fold.count(1))
# print("Subjects of Class 2: ", y_test_fold.count(2))
ct_train_files, ct_val_files, mr_train_files, mr_val_files, ct_test_files, mr_test_files = load_data(x_train_fold, x_test_fold, y_train_fold, y_test_fold, test_subjects, test_labels)
# Define transforms for CT and MR respectively
ct_transforms = Compose(
[
LoadImaged(keys=["img"], ensure_channel_first=True),
ScaleIntensityRanged(keys=["img"],a_min=0, a_max=85, b_min=0, b_max=1),
NormalizeIntensityd(keys=["img"]),
#Resized(keys=["img"], spatial_size=(96, 96, 96)),
]
)
mr_transforms = Compose(
[
LoadImaged(keys=["img"], ensure_channel_first=True),
#ScaleIntensityRanged(keys=["img"],a_min=0, a_max=85, b_min=0, b_max=1),
NormalizeIntensityd(keys=["img"]),
# RandFlipd(keys=["img"], prob=1, spatial_axis=2),
# RandRotated(keys=["img"], prob=1, range_x=[0.4,0.4])
#ScaleIntensityRanged(keys=["img"],a_min=0, a_max=85, b_min=0, b_max=1),
]
)
post_pred = Compose([Activations(softmax=True)])
post_label = Compose([AsDiscrete(to_onehot=3)])
# create a training data loader
ct_train_ds = monai.data.Dataset(data=ct_train_files, transform=ct_transforms)
mr_train_ds = monai.data.Dataset(data=mr_train_files, transform=mr_transforms)
combined_train_ds = CustomDataset(d1 = mr_train_ds, d2 = ct_train_ds)
train_loader = DataLoader(combined_train_ds, batch_size=1, shuffle=False, pin_memory=torch.cuda.is_available())
check_data1, check_data2 = monai.utils.misc.first(train_loader)
#print(check_data1["img"].shape, check_data1["label"], check_data2["img"].shape, check_data2["label"])
# create a validation data loader
ct_val_ds = monai.data.Dataset(data=ct_val_files, transform=ct_transforms)
mr_val_ds = monai.data.Dataset(data=mr_val_files, transform=mr_transforms)
combined_val_ds = CustomDataset(d1 = mr_val_ds, d2 = ct_val_ds)
val_loader = DataLoader(combined_val_ds, batch_size=1, shuffle=False, pin_memory=torch.cuda)
# create a test data loader
ct_test_ds = monai.data.Dataset(data=ct_test_files, transform=ct_transforms)
mr_test_ds = monai.data.Dataset(data=mr_test_files, transform=mr_transforms)
combined_test_ds = CustomDataset(d1 = mr_test_ds, d2 = ct_test_ds)
test_loader = DataLoader(combined_test_ds, batch_size=1, shuffle=False, pin_memory=torch.cuda)
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
#print(device)
model_ct = resnet18(pretrained=False, spatial_dims=3, n_input_channels=1, num_classes=3).to(device)
#model_ct.load_state_dict(torch.load("models/best_val_resnet18.pth"))
model_mr = resnet18(pretrained=False, spatial_dims=3, n_input_channels=1, num_classes=3).to(device)
#model_mr.load_state_dict(torch.load("models/best_mr_val_resnet18.pth"))
# Freeze these models
for param in model_mr.parameters():
param.requires_grad_(True)
for param in model_ct.parameters():
param.requires_grad_(True)
# Create ensemble model
model = MyEnsemble(model_ct, model_mr)
#model.load_state_dict(torch.load("/data/amciilab/ajoshi83/models/CT_Fold{}.pth".format(fold_number)))
#model.load_state_dict(torch.load("/scratch/ajoshi83/models_august/Fake_CT_Fold{}_ensemble_adam_slower.pth".format(int(fold_number))))
model.load_state_dict(torch.load("/scratch/ajoshi83/models_sbatch/Both_Real_Fold{}_ensemble_Adam_{}_{}.pth".format(int(fold_number), lr, decay)))
loss_function = torch.nn.CrossEntropyLoss()
optimizer = torch.optim.Adam(model.parameters(),lr= 1e-5, weight_decay=1e-3)
auc_metric = ROCAUCMetric(average="weighted")
# starting evaluation
val_interval = 1
best_metric = -1
best_metric_epoch = -1
best_val_loss = 2
writer = SummaryWriter()
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/train_output")
for batch_data in train_loader:
#step += 1
#print(batch_data["img"])
ct_batch, mr_batch = batch_data[0], batch_data[1]
#print(ct_batch)
#print(mr_batch)
ct_inputs, ct_labels = ct_batch["img"].to(device), ct_batch["label"].to(device)
mr_inputs, mr_labels = mr_batch["img"].to(device), mr_batch["label"].to(device)
train_outputs = model(ct_inputs, mr_inputs).argmax(dim=1)
y_pred = torch.cat([y_pred, model(ct_inputs, mr_inputs)], dim=0)
y = torch.cat([y, mr_labels], dim=0)
value = torch.eq(train_outputs, mr_labels)
metric_count += len(value)
num_correct += value.sum().item()
#saver.save_batch(train_outputs, train_data["img"].meta)
metric = num_correct / metric_count
# print("val evaluation metric:", metric)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred, detach = True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
saver.finalize()
train_acc.append(round(metric,3))
train_auc.append(round(auc_result,3))
print("Train Set: Accuracy: {:.4f}, AUC: {:.4f}".format(metric, auc_result))
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/val_output")
for val_data in val_loader:
ct_val_data, mr_val_data = val_data[0], val_data[1]
ct_val_images, ct_val_labels = ct_val_data["img"].to(device), ct_val_data["label"].to(device)
mr_val_images, mr_val_labels = mr_val_data["img"].to(device), mr_val_data["label"].to(device)
y_pred = torch.cat([y_pred, model(ct_val_images, mr_val_images)], dim=0)
#print(y_pred)
y = torch.cat([y, mr_val_labels], dim=0)
acc_value = torch.eq(y_pred.argmax(dim=1), y)
acc_metric = acc_value.sum().item() / len(acc_value)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred,detach=True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
saver.finalize()
val_acc.append(round(acc_metric,3))
val_auc.append(round(auc_result,3))
print("Val Set: Accuracy: {:.4f}, AUC: {:.4f}".format(acc_metric, auc_result))
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/test_output")
y_pred_mid = []
for test_data in test_loader:
ct_test_data, mr_test_data = test_data[0], test_data[1]
ct_test_images, ct_test_labels = ct_test_data["img"].to(device), ct_test_data["label"].to(device)
mr_test_images, mr_test_labels = mr_test_data["img"].to(device), mr_test_data["label"].to(device)
test_outputs = model(ct_test_images, mr_test_images).argmax(dim=1)
y_pred_mid.append(test_outputs.cpu().numpy()[0])
#softmax_op = torch.nn.functional.softmax(model(ct_test_images, mr_test_images)).argmax(dim=1)
#y_pred_mid.append(softmax_op.cpu().numpy()[0])
y_pred = torch.cat([y_pred, model(ct_test_images, mr_test_images)], dim=0)
y = torch.cat([y, mr_test_labels], dim=0)
value = torch.eq(test_outputs, mr_test_labels)
metric_count += len(value)
num_correct += value.sum().item()
#saver.save_batch(test_outputs, test_data["img"].meta)
y_pred_test.append(y_pred_mid)
metric = num_correct / metric_count
# print("test evaluation metric:", metric)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred, detach=True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
y_pred_auc.append(y_pred_act)
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
# auc_score =roc_auc_score(y_onehot, y_pred_act, multi_class='ovr', average='weighted')
# print(auc_score)
saver.finalize()
test_acc.append(round(metric,3))
test_auc.append(round(auc_result,3))
print("Test Set: Accuracy: {:.4f}, AUC: {:.4f}".format(metric, auc_result))
print("Fold {} completed...Next Fold starting...".format(fold_number))
print("\n")
fold_number += 1
writer.close()
print("\n")
print("Train Acc: {:.4f} +/- {:.4f}, Train AUC: {:.4f} +/- {:.4f}".format(mean(train_acc), np.std(train_acc), mean(train_auc), np.std(train_auc)))
print("Val Acc: {:.4f} +/- {:.4f}, Val AUC: {:.4f} +/- {:.4f}".format(mean(val_acc), np.std(val_acc), mean(val_auc), np.std(val_auc)))
# # print("Val Accuracies of 5 Folds:", val_acc)
# print("Val AUCs of 5 Folds:", val_auc)
print("Test Acc: {:.4f} +/- {:.4f}, Test AUC: {:.4f} +/- {:.4f}".format(mean(test_acc), np.std(test_acc), mean(test_auc), np.std(test_auc)))
# print("Test Accuracies of 5 Folds:", test_acc)
# print("Test AUCs of 5 Folds:", test_auc)
y_pred_tr = np.transpose(y_pred_test)
final = []
for i in range(y_pred_tr.shape[0]):
final.append(mode(y_pred_tr[i]))
y_true = y.cpu().numpy()
final_np = np.array(final)
y_pred_auc = np.array(y_pred_auc)
sum = np.sum(y_pred_auc, axis=0)
sum = sum/5
target_names = ['class 0', 'class 1', 'class 2']
print("\n")
print("Test Accuracy (Hard Voting of 5 models):", round(accuracy_score(y_true, final_np), 4))
print("Test AUC (Hard Voting of 5 models):", round(roc_auc_score(y_onehot, sum, multi_class='ovr', average='weighted'), 4))
cm = confusion_matrix(y_true, final_np)
print(classification_report(y_true, final_np, target_names=target_names))
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot()
plt.show()
# %%
# %%
import logging
import os
import sys
from pathlib import Path
import pandas as pd
import monai
import numpy as np
import torch
import torch.nn as nn
from torch.utils.data import DataLoader as _TorchDataLoader
from torch.utils.data import Dataset
from torch.utils.tensorboard import SummaryWriter
from monai.data import decollate_batch, CSVSaver, DataLoader
from monai.data.utils import list_data_collate, set_rnd, worker_init_fn
from monai.metrics import ROCAUCMetric
from monai.transforms import ScaleIntensityRanged, NormalizeIntensityd, Activations, AsDiscrete, Compose, LoadImaged, RandRotated, Resized, ScaleIntensityd, RandFlipd
from monai.networks.nets import resnet10, resnet18, resnet34, resnet50
from sklearn.model_selection import StratifiedKFold
from statistics import mean, mode
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, classification_report, accuracy_score, roc_auc_score
class CustomDataset(monai.data.Dataset):
def __init__(self, d1, d2):
self.d1 = d1
self.d2 = d2
def __getitem__(self, idx):
dict1 = self.d1.__getitem__(idx)
image1, label1 = dict1["img"], dict1["label"]
dict2 = self.d2.__getitem__(idx)
image2, label2 = dict2["img"], dict2["label"]
assert label1==label2
dict_1 = dict()
dict_1["img"] = image1
dict_1["label"] = label1
dict_2 = dict()
dict_2["img"] = image2
dict_2["label"] = label2
return dict_1, dict_2
def __len__(self):
return len(self.d1)
class MyEnsemble(nn.Module):
def __init__(self, modelA, modelB):
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
super(MyEnsemble, self).__init__()
self.modelA = modelA
self.modelB = modelB
# Remove last linear layer
self.modelA.fc = nn.Identity()
self.modelB.fc = nn.Identity()
# Create new classifier
self.mlp1 = nn.Linear(1024,256).to(device)
self.mlp2 = nn.Linear(256,32).to(device)
self.classifier = nn.Linear(32,3).to(device)
def forward(self, i1, i2):
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
x1 = self.modelA(i1) # clone to make sure x is not changed by inplace methods
x1 = x1.view(x1.size(0), -1).to(device)
x2 = self.modelB(i2)
x2 = x2.view(x2.size(0), -1).to(device)
x = torch.cat((x1, x2), dim=1).to(device)
x = nn.functional.relu(self.mlp1(x)).to(device)
x = nn.functional.relu(self.mlp2(x)).to(device)
x = self.classifier(x).to(device)
return x
def load_data(train_images, val_images, train_labels, val_labels, test_images, test_labels):
train_labels_array = np.array(train_labels, dtype=np.int64)
val_labels_array = np.array(val_labels, dtype=np.int64)
test_labels_array = np.array(test_labels, dtype=np.int64)
#loading the CT images
#ct_path = Path('/scratch/ajoshi83/Data/Reg_CT')
ct_path = Path('/scratch/ajoshi83/generated_ct')
#ct_path = Path('/data/amciilab/ajoshi83/generated_ct2')
ct_train_images_path = [os.path.join(ct_path, f) for f in train_images]
ct_val_images_path = [os.path.join(ct_path, f) for f in val_images]
ct_test_images_path = [os.path.join(ct_path, f) for f in test_images]
ct_train_files = [{"img": img, "label": label} for img, label in zip(ct_train_images_path, train_labels_array)]
ct_val_files = [{"img": img, "label": label} for img, label in zip(ct_val_images_path, val_labels_array)]
ct_test_files = [{"img": img, "label": label} for img, label in zip(ct_test_images_path, test_labels_array)]
#loading the MR images
mr_path = Path('/scratch/ajoshi83/Data/Reg_MR')
#mr_path = Path('/data/amciilab/ajoshi83/generated_mr2')
mr_train_images_path = [os.path.join(mr_path, f) for f in train_images]
mr_val_images_path = [os.path.join(mr_path, f) for f in val_images]
mr_test_images_path = [os.path.join(mr_path, f) for f in test_images]
mr_train_files = [{"img": img, "label": label} for img, label in zip(mr_train_images_path, train_labels_array)]
mr_val_files = [{"img": img, "label": label} for img, label in zip(mr_val_images_path, val_labels_array)]
mr_test_files = [{"img": img, "label": label} for img, label in zip(mr_test_images_path, test_labels_array)]
return ct_train_files, ct_val_files, mr_train_files, mr_val_files, ct_test_files, mr_test_files
#monai.config.print_config()
# old_stdout = sys.stdout
# log_file = open("mri_classification/logs/5fold_eval_fakemr.log","w")
# sys.stdout = log_file
print("this will be written to message.log")
# file_handler = logging.FileHandler(filename='logs/tmp.log')
# #stdout_handler = logging.StreamHandler(stream=sys.stdout)
# handlers = [file_handler]
# logging.basicConfig(
# level=logging.DEBUG,
# format='[%(asctime)s] {%(filename)s:%(lineno)d} %(levelname)s - %(message)s',
# handlers=handlers
# )
# logger = logging.getLogger('LOGGER_NAME')
# edit the path accordingly
gose = pd.read_csv('mri_classification/Pilot_169_train_val_test_split.csv')
X = []
Y = []
test_subjects = []
test_labels = []
torch.cuda.empty_cache()
for i in range(gose.shape[0]):
subj = gose['Main.GUID'][i]
subj_id = str(subj)[4:]
name = str(subj_id) + '.nii'
if gose['Set'][i]=='Train':
X.append(name)
Y.append(gose['Class'][i])
elif gose['Set'][i]=='Val':
X.append(name)
Y.append(gose['Class'][i])
elif gose['Set'][i]=='Test':
test_subjects.append(name)
test_labels.append(gose['Class'][i])
else:
print("Unknown Set: ", gose['Set'][i])
print("Total Subjects for 5-Fold (Train+Val):", len(X))
print("Test Subjects (kept separate):", len(test_subjects))
print("5-Fold Subjects Information:\n")
print("Subjects of Class 0: ", Y.count(0))
print("Subjects of Class 1: ", Y.count(1))
print("Subjects of Class 2: ", Y.count(2))
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=1)
fold_number = 1
lr =1e-5
decay = 1e-5
X = np.array(X)
Y = np.array(Y)
train_acc, train_auc = [], []
val_acc, val_auc = [], []
test_acc, test_auc = [], []
y_pred_test = []
y_pred_auc = []
for train_index, test_index in skf.split(X, Y):
print("Starting Fold {}..!".format(fold_number))
x_train_fold, x_test_fold = X[train_index], X[test_index]
y_train_fold, y_test_fold = Y[train_index], Y[test_index]
y_train_fold = list(y_train_fold)
y_test_fold = list(y_test_fold)
x_train_fold = list(x_train_fold)
x_test_fold = list(x_test_fold)
# print("\n")
# print("Fold {} statistics:\n".format(fold_number))
# print("Train Subjects: {}".format(len(y_train_fold)))
# print("Subjects of Class 0: ", y_train_fold.count(0))
# print("Subjects of Class 1: ", y_train_fold.count(1))
# print("Subjects of Class 2: {}\n".format(y_train_fold.count(2)))
# print("Val Subjects: {}".format(len(y_test_fold)))
# print("Subjects of Class 0: ", y_test_fold.count(0))
# print("Subjects of Class 1: ", y_test_fold.count(1))
# print("Subjects of Class 2: ", y_test_fold.count(2))
ct_train_files, ct_val_files, mr_train_files, mr_val_files, ct_test_files, mr_test_files = load_data(x_train_fold, x_test_fold, y_train_fold, y_test_fold, test_subjects, test_labels)
# Define transforms for CT and MR respectively
ct_transforms = Compose(
[
LoadImaged(keys=["img"], ensure_channel_first=True),
#ScaleIntensityRanged(keys=["img"],a_min=0, a_max=85, b_min=0, b_max=1),
NormalizeIntensityd(keys=["img"]),
#Resized(keys=["img"], spatial_size=(96, 96, 96)),
]
)
mr_transforms = Compose(
[
LoadImaged(keys=["img"], ensure_channel_first=True),
#ScaleIntensityRanged(keys=["img"],a_min=0, a_max=85, b_min=0, b_max=1),
NormalizeIntensityd(keys=["img"]),
# RandFlipd(keys=["img"], prob=1, spatial_axis=2),
# RandRotated(keys=["img"], prob=1, range_x=[0.4,0.4])
#ScaleIntensityRanged(keys=["img"],a_min=0, a_max=85, b_min=0, b_max=1),
]
)
post_pred = Compose([Activations(softmax=True)])
post_label = Compose([AsDiscrete(to_onehot=3)])
# create a training data loader
ct_train_ds = monai.data.Dataset(data=ct_train_files, transform=ct_transforms)
mr_train_ds = monai.data.Dataset(data=mr_train_files, transform=mr_transforms)
combined_train_ds = CustomDataset(d1 = mr_train_ds, d2 = ct_train_ds)
train_loader = DataLoader(combined_train_ds, batch_size=1, shuffle=False, pin_memory=torch.cuda.is_available())
check_data1, check_data2 = monai.utils.misc.first(train_loader)
#print(check_data1["img"].shape, check_data1["label"], check_data2["img"].shape, check_data2["label"])
# create a validation data loader
ct_val_ds = monai.data.Dataset(data=ct_val_files, transform=ct_transforms)
mr_val_ds = monai.data.Dataset(data=mr_val_files, transform=mr_transforms)
combined_val_ds = CustomDataset(d1 = mr_val_ds, d2 = ct_val_ds)
val_loader = DataLoader(combined_val_ds, batch_size=1, shuffle=False, pin_memory=torch.cuda)
# create a test data loader
ct_test_ds = monai.data.Dataset(data=ct_test_files, transform=ct_transforms)
mr_test_ds = monai.data.Dataset(data=mr_test_files, transform=mr_transforms)
combined_test_ds = CustomDataset(d1 = mr_test_ds, d2 = ct_test_ds)
test_loader = DataLoader(combined_test_ds, batch_size=1, shuffle=False, pin_memory=torch.cuda)
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
#print(device)
model_ct = resnet18(pretrained=False, spatial_dims=3, n_input_channels=1, num_classes=3).to(device)
#model_ct.load_state_dict(torch.load("models/best_val_resnet18.pth"))
model_mr = resnet18(pretrained=False, spatial_dims=3, n_input_channels=1, num_classes=3).to(device)
#model_mr.load_state_dict(torch.load("models/best_mr_val_resnet18.pth"))
# Freeze these models
for param in model_mr.parameters():
param.requires_grad_(True)
for param in model_ct.parameters():
param.requires_grad_(True)
# Create ensemble model
model = MyEnsemble(model_ct, model_mr)
#model.load_state_dict(torch.load("/data/amciilab/ajoshi83/models/CT_Fold{}.pth".format(fold_number)))
model.load_state_dict(torch.load("/scratch/ajoshi83/models_sbatch/Fake_CT_Fold{}_ensemble_Adam_{}.pth".format(int(fold_number), lr)))
#model.load_state_dict(torch.load("/scratch/ajoshi83/models_sbatch/Both_Real_Fold{}_ensemble_Adam_{}_{}.pth".format(int(fold_number), lr, decay)))
# loss_function = torch.nn.CrossEntropyLoss()
# optimizer = torch.optim.Adam(model.parameters(),lr= 1e-5, weight_decay=1e-3)
auc_metric = ROCAUCMetric(average="weighted")
# starting evaluation
val_interval = 1
best_metric = -1
best_metric_epoch = -1
best_val_loss = 2
writer = SummaryWriter()
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/train_output")
for batch_data in train_loader:
#step += 1
#print(batch_data["img"])
ct_batch, mr_batch = batch_data[0], batch_data[1]
#print(ct_batch)
#print(mr_batch)
ct_inputs, ct_labels = ct_batch["img"].to(device), ct_batch["label"].to(device)
mr_inputs, mr_labels = mr_batch["img"].to(device), mr_batch["label"].to(device)
train_outputs = model(ct_inputs, mr_inputs).argmax(dim=1)
y_pred = torch.cat([y_pred, model(ct_inputs, mr_inputs)], dim=0)
y = torch.cat([y, mr_labels], dim=0)
value = torch.eq(train_outputs, mr_labels)
metric_count += len(value)
num_correct += value.sum().item()
#saver.save_batch(train_outputs, train_data["img"].meta)
metric = num_correct / metric_count
# print("val evaluation metric:", metric)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred, detach = True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
saver.finalize()
train_acc.append(round(metric,3))
train_auc.append(round(auc_result,3))
print("Train Set: Accuracy: {:.4f}, AUC: {:.4f}".format(metric, auc_result))
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/val_output")
for val_data in val_loader:
ct_val_data, mr_val_data = val_data[0], val_data[1]
ct_val_images, ct_val_labels = ct_val_data["img"].to(device), ct_val_data["label"].to(device)
mr_val_images, mr_val_labels = mr_val_data["img"].to(device), mr_val_data["label"].to(device)
y_pred = torch.cat([y_pred, model(ct_val_images, mr_val_images)], dim=0)
#print(y_pred)
y = torch.cat([y, mr_val_labels], dim=0)
acc_value = torch.eq(y_pred.argmax(dim=1), y)
acc_metric = acc_value.sum().item() / len(acc_value)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred,detach=True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
saver.finalize()
val_acc.append(round(acc_metric,3))
val_auc.append(round(auc_result,3))
print("Val Set: Accuracy: {:.4f}, AUC: {:.4f}".format(acc_metric, auc_result))
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/test_output")
y_pred_mid = []
for test_data in test_loader:
ct_test_data, mr_test_data = test_data[0], test_data[1]
ct_test_images, ct_test_labels = ct_test_data["img"].to(device), ct_test_data["label"].to(device)
mr_test_images, mr_test_labels = mr_test_data["img"].to(device), mr_test_data["label"].to(device)
test_outputs = model(ct_test_images, mr_test_images).argmax(dim=1)
y_pred_mid.append(test_outputs.cpu().numpy()[0])
#softmax_op = torch.nn.functional.softmax(model(ct_test_images, mr_test_images)).argmax(dim=1)
#y_pred_mid.append(softmax_op.cpu().numpy()[0])
y_pred = torch.cat([y_pred, model(ct_test_images, mr_test_images)], dim=0)
y = torch.cat([y, mr_test_labels], dim=0)
value = torch.eq(test_outputs, mr_test_labels)
metric_count += len(value)
num_correct += value.sum().item()
#saver.save_batch(test_outputs, test_data["img"].meta)
y_pred_test.append(y_pred_mid)
metric = num_correct / metric_count
# print("test evaluation metric:", metric)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred, detach=True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
y_pred_auc.append(y_pred_act)
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
# auc_score =roc_auc_score(y_onehot, y_pred_act, multi_class='ovr', average='weighted')
# print(auc_score)
saver.finalize()
test_acc.append(round(metric,3))
test_auc.append(round(auc_result,3))
print("Test Set: Accuracy: {:.4f}, AUC: {:.4f}".format(metric, auc_result))
print("Fold {} completed...Next Fold starting...".format(fold_number))
print("\n")
fold_number += 1
writer.close()
print("\n")
print("Train Acc: {:.4f} +/- {:.4f}, Train AUC: {:.4f} +/- {:.4f}".format(mean(train_acc), np.std(train_acc), mean(train_auc), np.std(train_auc)))
print("Val Acc: {:.4f} +/- {:.4f}, Val AUC: {:.4f} +/- {:.4f}".format(mean(val_acc), np.std(val_acc), mean(val_auc), np.std(val_auc)))
# # print("Val Accuracies of 5 Folds:", val_acc)
# print("Val AUCs of 5 Folds:", val_auc)
print("Test Acc: {:.4f} +/- {:.4f}, Test AUC: {:.4f} +/- {:.4f}".format(mean(test_acc), np.std(test_acc), mean(test_auc), np.std(test_auc)))
# print("Test Accuracies of 5 Folds:", test_acc)
# print("Test AUCs of 5 Folds:", test_auc)
y_pred_tr = np.transpose(y_pred_test)
final = []
for i in range(y_pred_tr.shape[0]):
final.append(mode(y_pred_tr[i]))
y_true = y.cpu().numpy()
final_np = np.array(final)
y_pred_auc = np.array(y_pred_auc)
sum = np.sum(y_pred_auc, axis=0)
sum = sum/5
target_names = ['class 0', 'class 1', 'class 2']
print("\n")
print("Test Accuracy (Hard Voting of 5 models):", round(accuracy_score(y_true, final_np), 4))
print("Test AUC (Hard Voting of 5 models):", round(roc_auc_score(y_onehot, sum, multi_class='ovr', average='weighted'), 4))
cm = confusion_matrix(y_true, final_np)
print(classification_report(y_true, final_np, target_names=target_names))
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot()
plt.show()
# %%
# %%
import logging
import os
import sys
from pathlib import Path
import pandas as pd
import monai
import numpy as np
import torch
import torch.nn as nn
from torch.utils.data import DataLoader as _TorchDataLoader
from torch.utils.data import Dataset
from torch.utils.tensorboard import SummaryWriter
from monai.data import decollate_batch, CSVSaver, DataLoader
from monai.data.utils import list_data_collate, set_rnd, worker_init_fn
from monai.metrics import ROCAUCMetric
from monai.transforms import ScaleIntensityRanged, NormalizeIntensityd, Activations, AsDiscrete, Compose, LoadImaged, RandRotated, Resized, ScaleIntensityd, RandFlipd
from monai.networks.nets import resnet10, resnet18, resnet34, resnet50
from sklearn.model_selection import StratifiedKFold
from statistics import mean, mode
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, classification_report, accuracy_score, roc_auc_score
class CustomDataset(monai.data.Dataset):
def __init__(self, d1, d2):
self.d1 = d1
self.d2 = d2
def __getitem__(self, idx):
dict1 = self.d1.__getitem__(idx)
image1, label1 = dict1["img"], dict1["label"]
dict2 = self.d2.__getitem__(idx)
image2, label2 = dict2["img"], dict2["label"]
assert label1==label2
dict_1 = dict()
dict_1["img"] = image1
dict_1["label"] = label1
dict_2 = dict()
dict_2["img"] = image2
dict_2["label"] = label2
return dict_1, dict_2
def __len__(self):
return len(self.d1)
class MyEnsemble(nn.Module):
def __init__(self, modelA, modelB):
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
super(MyEnsemble, self).__init__()
self.modelA = modelA
self.modelB = modelB
# Remove last linear layer
self.modelA.fc = nn.Identity()
self.modelB.fc = nn.Identity()
# Create new classifier
self.mlp1 = nn.Linear(1024,256).to(device)
self.mlp2 = nn.Linear(256,32).to(device)
self.classifier = nn.Linear(32,3).to(device)
def forward(self, i1, i2):
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
x1 = self.modelA(i1) # clone to make sure x is not changed by inplace methods
x1 = x1.view(x1.size(0), -1).to(device)
x2 = self.modelB(i2)
x2 = x2.view(x2.size(0), -1).to(device)
x = torch.cat((x1, x2), dim=1).to(device)
x = nn.functional.relu(self.mlp1(x)).to(device)
x = nn.functional.relu(self.mlp2(x)).to(device)
x = self.classifier(x).to(device)
return x
def load_data(train_images, val_images, train_labels, val_labels, test_images, test_labels):
train_labels_array = np.array(train_labels, dtype=np.int64)
val_labels_array = np.array(val_labels, dtype=np.int64)
test_labels_array = np.array(test_labels, dtype=np.int64)
#loading the CT images
ct_path = Path('/scratch/ajoshi83/Data/Reg_CT')
#ct_path = Path('/scratch/ajoshi83/generated_ct')
#ct_path = Path('/data/amciilab/ajoshi83/generated_ct2')
ct_train_images_path = [os.path.join(ct_path, f) for f in train_images]
ct_val_images_path = [os.path.join(ct_path, f) for f in val_images]
ct_test_images_path = [os.path.join(ct_path, f) for f in test_images]
ct_train_files = [{"img": img, "label": label} for img, label in zip(ct_train_images_path, train_labels_array)]
ct_val_files = [{"img": img, "label": label} for img, label in zip(ct_val_images_path, val_labels_array)]
ct_test_files = [{"img": img, "label": label} for img, label in zip(ct_test_images_path, test_labels_array)]
#loading the MR images
#mr_path = Path('/scratch/ajoshi83/Data/Reg_MR')
mr_path = Path('/scratch/ajoshi83/generated_mr')
#mr_path = Path('/data/amciilab/ajoshi83/generated_mr2')
mr_train_images_path = [os.path.join(mr_path, f) for f in train_images]
mr_val_images_path = [os.path.join(mr_path, f) for f in val_images]
mr_test_images_path = [os.path.join(mr_path, f) for f in test_images]
mr_train_files = [{"img": img, "label": label} for img, label in zip(mr_train_images_path, train_labels_array)]
mr_val_files = [{"img": img, "label": label} for img, label in zip(mr_val_images_path, val_labels_array)]
mr_test_files = [{"img": img, "label": label} for img, label in zip(mr_test_images_path, test_labels_array)]
return ct_train_files, ct_val_files, mr_train_files, mr_val_files, ct_test_files, mr_test_files
#monai.config.print_config()
# old_stdout = sys.stdout
# log_file = open("mri_classification/logs/5fold_eval_fakemr.log","w")
# sys.stdout = log_file
print("this will be written to message.log")
# file_handler = logging.FileHandler(filename='logs/tmp.log')
# #stdout_handler = logging.StreamHandler(stream=sys.stdout)
# handlers = [file_handler]
# logging.basicConfig(
# level=logging.DEBUG,
# format='[%(asctime)s] {%(filename)s:%(lineno)d} %(levelname)s - %(message)s',
# handlers=handlers
# )
# logger = logging.getLogger('LOGGER_NAME')
# edit the path accordingly
gose = pd.read_csv('mri_classification/Pilot_169_train_val_test_split.csv')
X = []
Y = []
test_subjects = []
test_labels = []
torch.cuda.empty_cache()
for i in range(gose.shape[0]):
subj = gose['Main.GUID'][i]
subj_id = str(subj)[4:]
name = str(subj_id) + '.nii'
if gose['Set'][i]=='Train':
X.append(name)
Y.append(gose['Class'][i])
elif gose['Set'][i]=='Val':
X.append(name)
Y.append(gose['Class'][i])
elif gose['Set'][i]=='Test':
test_subjects.append(name)
test_labels.append(gose['Class'][i])
else:
print("Unknown Set: ", gose['Set'][i])
print("Total Subjects for 5-Fold (Train+Val):", len(X))
print("Test Subjects (kept separate):", len(test_subjects))
print("5-Fold Subjects Information:\n")
print("Subjects of Class 0: ", Y.count(0))
print("Subjects of Class 1: ", Y.count(1))
print("Subjects of Class 2: ", Y.count(2))
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=1)
fold_number = 1
lr =1e-5
decay = 1e-5
X = np.array(X)
Y = np.array(Y)
train_acc, train_auc = [], []
val_acc, val_auc = [], []
test_acc, test_auc = [], []
y_pred_test = []
y_pred_auc = []
for train_index, test_index in skf.split(X, Y):
print("Starting Fold {}..!".format(fold_number))
x_train_fold, x_test_fold = X[train_index], X[test_index]
y_train_fold, y_test_fold = Y[train_index], Y[test_index]
y_train_fold = list(y_train_fold)
y_test_fold = list(y_test_fold)
x_train_fold = list(x_train_fold)
x_test_fold = list(x_test_fold)
# print("\n")
# print("Fold {} statistics:\n".format(fold_number))
# print("Train Subjects: {}".format(len(y_train_fold)))
# print("Subjects of Class 0: ", y_train_fold.count(0))
# print("Subjects of Class 1: ", y_train_fold.count(1))
# print("Subjects of Class 2: {}\n".format(y_train_fold.count(2)))
# print("Val Subjects: {}".format(len(y_test_fold)))
# print("Subjects of Class 0: ", y_test_fold.count(0))
# print("Subjects of Class 1: ", y_test_fold.count(1))
# print("Subjects of Class 2: ", y_test_fold.count(2))
ct_train_files, ct_val_files, mr_train_files, mr_val_files, ct_test_files, mr_test_files = load_data(x_train_fold, x_test_fold, y_train_fold, y_test_fold, test_subjects, test_labels)
# Define transforms for CT and MR respectively
ct_transforms = Compose(
[
LoadImaged(keys=["img"], ensure_channel_first=True),
ScaleIntensityRanged(keys=["img"],a_min=0, a_max=85, b_min=0, b_max=1),
NormalizeIntensityd(keys=["img"]),
#Resized(keys=["img"], spatial_size=(96, 96, 96)),
]
)
mr_transforms = Compose(
[
LoadImaged(keys=["img"], ensure_channel_first=True),
#ScaleIntensityRanged(keys=["img"],a_min=0, a_max=85, b_min=0, b_max=1),
NormalizeIntensityd(keys=["img"]),
# RandFlipd(keys=["img"], prob=1, spatial_axis=2),
# RandRotated(keys=["img"], prob=1, range_x=[0.4,0.4])
#ScaleIntensityRanged(keys=["img"],a_min=0, a_max=85, b_min=0, b_max=1),
]
)
post_pred = Compose([Activations(softmax=True)])
post_label = Compose([AsDiscrete(to_onehot=3)])
# create a training data loader
ct_train_ds = monai.data.Dataset(data=ct_train_files, transform=ct_transforms)
mr_train_ds = monai.data.Dataset(data=mr_train_files, transform=mr_transforms)
combined_train_ds = CustomDataset(d1 = mr_train_ds, d2 = ct_train_ds)
train_loader = DataLoader(combined_train_ds, batch_size=1, shuffle=False, pin_memory=torch.cuda.is_available())
check_data1, check_data2 = monai.utils.misc.first(train_loader)
#print(check_data1["img"].shape, check_data1["label"], check_data2["img"].shape, check_data2["label"])
# create a validation data loader
ct_val_ds = monai.data.Dataset(data=ct_val_files, transform=ct_transforms)
mr_val_ds = monai.data.Dataset(data=mr_val_files, transform=mr_transforms)
combined_val_ds = CustomDataset(d1 = mr_val_ds, d2 = ct_val_ds)
val_loader = DataLoader(combined_val_ds, batch_size=1, shuffle=False, pin_memory=torch.cuda)
# create a test data loader
ct_test_ds = monai.data.Dataset(data=ct_test_files, transform=ct_transforms)
mr_test_ds = monai.data.Dataset(data=mr_test_files, transform=mr_transforms)
combined_test_ds = CustomDataset(d1 = mr_test_ds, d2 = ct_test_ds)
test_loader = DataLoader(combined_test_ds, batch_size=1, shuffle=False, pin_memory=torch.cuda)
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
#print(device)
model_ct = resnet18(pretrained=False, spatial_dims=3, n_input_channels=1, num_classes=3).to(device)
#model_ct.load_state_dict(torch.load("models/best_val_resnet18.pth"))
model_mr = resnet18(pretrained=False, spatial_dims=3, n_input_channels=1, num_classes=3).to(device)
#model_mr.load_state_dict(torch.load("models/best_mr_val_resnet18.pth"))
# Freeze these models
for param in model_mr.parameters():
param.requires_grad_(True)
for param in model_ct.parameters():
param.requires_grad_(True)
# Create ensemble model
model = MyEnsemble(model_ct, model_mr)
#model.load_state_dict(torch.load("/data/amciilab/ajoshi83/models/CT_Fold{}.pth".format(fold_number)))
model.load_state_dict(torch.load("/scratch/ajoshi83/models_sbatch/Fake_MR_Fold{}_ensemble_AdamW_{}.pth".format(int(fold_number), lr)))
#model.load_state_dict(torch.load("/scratch/ajoshi83/models_sbatch/Both_Real_Fold{}_ensemble_Adam_{}_{}.pth".format(int(fold_number), lr, decay)))
# loss_function = torch.nn.CrossEntropyLoss()
# optimizer = torch.optim.Adam(model.parameters(),lr= 1e-5, weight_decay=1e-3)
auc_metric = ROCAUCMetric(average="weighted")
# starting evaluation
val_interval = 1
best_metric = -1
best_metric_epoch = -1
best_val_loss = 2
writer = SummaryWriter()
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/train_output")
for batch_data in train_loader:
#step += 1
#print(batch_data["img"])
ct_batch, mr_batch = batch_data[0], batch_data[1]
#print(ct_batch)
#print(mr_batch)
ct_inputs, ct_labels = ct_batch["img"].to(device), ct_batch["label"].to(device)
mr_inputs, mr_labels = mr_batch["img"].to(device), mr_batch["label"].to(device)
train_outputs = model(ct_inputs, mr_inputs).argmax(dim=1)
y_pred = torch.cat([y_pred, model(ct_inputs, mr_inputs)], dim=0)
y = torch.cat([y, mr_labels], dim=0)
value = torch.eq(train_outputs, mr_labels)
metric_count += len(value)
num_correct += value.sum().item()
#saver.save_batch(train_outputs, train_data["img"].meta)
metric = num_correct / metric_count
# print("val evaluation metric:", metric)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred, detach = True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
saver.finalize()
train_acc.append(round(metric,3))
train_auc.append(round(auc_result,3))
print("Train Set: Accuracy: {:.4f}, AUC: {:.4f}".format(metric, auc_result))
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/val_output")
for val_data in val_loader:
ct_val_data, mr_val_data = val_data[0], val_data[1]
ct_val_images, ct_val_labels = ct_val_data["img"].to(device), ct_val_data["label"].to(device)
mr_val_images, mr_val_labels = mr_val_data["img"].to(device), mr_val_data["label"].to(device)
y_pred = torch.cat([y_pred, model(ct_val_images, mr_val_images)], dim=0)
#print(y_pred)
y = torch.cat([y, mr_val_labels], dim=0)
acc_value = torch.eq(y_pred.argmax(dim=1), y)
acc_metric = acc_value.sum().item() / len(acc_value)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred,detach=True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
saver.finalize()
val_acc.append(round(acc_metric,3))
val_auc.append(round(auc_result,3))
print("Val Set: Accuracy: {:.4f}, AUC: {:.4f}".format(acc_metric, auc_result))
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/test_output")
y_pred_mid = []
for test_data in test_loader:
ct_test_data, mr_test_data = test_data[0], test_data[1]
ct_test_images, ct_test_labels = ct_test_data["img"].to(device), ct_test_data["label"].to(device)
mr_test_images, mr_test_labels = mr_test_data["img"].to(device), mr_test_data["label"].to(device)
test_outputs = model(ct_test_images, mr_test_images).argmax(dim=1)
y_pred_mid.append(test_outputs.cpu().numpy()[0])
#softmax_op = torch.nn.functional.softmax(model(ct_test_images, mr_test_images)).argmax(dim=1)
#y_pred_mid.append(softmax_op.cpu().numpy()[0])
y_pred = torch.cat([y_pred, model(ct_test_images, mr_test_images)], dim=0)
y = torch.cat([y, mr_test_labels], dim=0)
value = torch.eq(test_outputs, mr_test_labels)
metric_count += len(value)
num_correct += value.sum().item()
#saver.save_batch(test_outputs, test_data["img"].meta)
y_pred_test.append(y_pred_mid)
metric = num_correct / metric_count
# print("test evaluation metric:", metric)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred, detach=True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
y_pred_auc.append(y_pred_act)
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
# auc_score =roc_auc_score(y_onehot, y_pred_act, multi_class='ovr', average='weighted')
# print(auc_score)
saver.finalize()
test_acc.append(round(metric,3))
test_auc.append(round(auc_result,3))
print("Test Set: Accuracy: {:.4f}, AUC: {:.4f}".format(metric, auc_result))
print("Fold {} completed...Next Fold starting...".format(fold_number))
print("\n")
fold_number += 1
writer.close()
print("\n")
print("Train Acc: {:.4f} +/- {:.4f}, Train AUC: {:.4f} +/- {:.4f}".format(mean(train_acc), np.std(train_acc), mean(train_auc), np.std(train_auc)))
print("Val Acc: {:.4f} +/- {:.4f}, Val AUC: {:.4f} +/- {:.4f}".format(mean(val_acc), np.std(val_acc), mean(val_auc), np.std(val_auc)))
# # print("Val Accuracies of 5 Folds:", val_acc)
# print("Val AUCs of 5 Folds:", val_auc)
print("Test Acc: {:.4f} +/- {:.4f}, Test AUC: {:.4f} +/- {:.4f}".format(mean(test_acc), np.std(test_acc), mean(test_auc), np.std(test_auc)))
# print("Test Accuracies of 5 Folds:", test_acc)
# print("Test AUCs of 5 Folds:", test_auc)
y_pred_tr = np.transpose(y_pred_test)
final = []
for i in range(y_pred_tr.shape[0]):
final.append(mode(y_pred_tr[i]))
y_true = y.cpu().numpy()
final_np = np.array(final)
y_pred_auc = np.array(y_pred_auc)
sum = np.sum(y_pred_auc, axis=0)
sum = sum/5
target_names = ['class 0', 'class 1', 'class 2']
print("\n")
print("Test Accuracy (Hard Voting of 5 models):", round(accuracy_score(y_true, final_np), 4))
print("Test AUC (Hard Voting of 5 models):", round(roc_auc_score(y_onehot, sum, multi_class='ovr', average='weighted'), 4))
cm = confusion_matrix(y_true, final_np)
print(classification_report(y_true, final_np, target_names=target_names))
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot()
plt.show()
# %%
# %%
import logging
import os
import sys
from pathlib import Path
import pandas as pd
import monai
import numpy as np
import torch
import torch.nn as nn
from torch.utils.data import DataLoader as _TorchDataLoader
from torch.utils.data import Dataset
from torch.utils.tensorboard import SummaryWriter
from monai.data import decollate_batch, CSVSaver, DataLoader
from monai.data.utils import list_data_collate, set_rnd, worker_init_fn
from monai.metrics import ROCAUCMetric
from monai.transforms import ScaleIntensityRanged, NormalizeIntensityd, Activations, AsDiscrete, Compose, LoadImaged, RandRotated, Resized, ScaleIntensityd, RandFlipd
from monai.networks.nets import resnet10, resnet18, resnet34, resnet50
from sklearn.model_selection import StratifiedKFold
from statistics import mean, mode
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, classification_report, accuracy_score, roc_auc_score
class CustomDataset(monai.data.Dataset):
def __init__(self, d1, d2):
self.d1 = d1
self.d2 = d2
def __getitem__(self, idx):
dict1 = self.d1.__getitem__(idx)
image1, label1 = dict1["img"], dict1["label"]
dict2 = self.d2.__getitem__(idx)
image2, label2 = dict2["img"], dict2["label"]
assert label1==label2
dict_1 = dict()
dict_1["img"] = image1
dict_1["label"] = label1
dict_2 = dict()
dict_2["img"] = image2
dict_2["label"] = label2
return dict_1, dict_2
def __len__(self):
return len(self.d1)
class MyEnsemble(nn.Module):
def __init__(self, modelA, modelB):
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
super(MyEnsemble, self).__init__()
self.modelA = modelA
self.modelB = modelB
# Remove last linear layer
self.modelA.fc = nn.Identity()
self.modelB.fc = nn.Identity()
# Create new classifier
self.mlp1 = nn.Linear(1024,256).to(device)
self.mlp2 = nn.Linear(256,32).to(device)
self.classifier = nn.Linear(32,3).to(device)
def forward(self, i1, i2):
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
x1 = self.modelA(i1) # clone to make sure x is not changed by inplace methods
x1 = x1.view(x1.size(0), -1).to(device)
x2 = self.modelB(i2)
x2 = x2.view(x2.size(0), -1).to(device)
x = torch.cat((x1, x2), dim=1).to(device)
x = nn.functional.relu(self.mlp1(x)).to(device)
x = nn.functional.relu(self.mlp2(x)).to(device)
x = self.classifier(x).to(device)
return x
def load_data(train_images, val_images, train_labels, val_labels, test_images, test_labels):
train_labels_array = np.array(train_labels, dtype=np.int64)
val_labels_array = np.array(val_labels, dtype=np.int64)
test_labels_array = np.array(test_labels, dtype=np.int64)
#loading the CT images
ct_path = Path('Data/Reg_MR')
#ct_path = Path('/data/amciilab/ajoshi83/generated_ct2')
ct_train_images_path = [os.path.join(ct_path, f) for f in train_images]
ct_val_images_path = [os.path.join(ct_path, f) for f in val_images]
ct_test_images_path = [os.path.join(ct_path, f) for f in test_images]
ct_train_files = [{"img": img, "label": label} for img, label in zip(ct_train_images_path, train_labels_array)]
ct_val_files = [{"img": img, "label": label} for img, label in zip(ct_val_images_path, val_labels_array)]
ct_test_files = [{"img": img, "label": label} for img, label in zip(ct_test_images_path, test_labels_array)]
#loading the MR images
mr_path = Path('Data/Reg_MR')
#mr_path = Path('/data/amciilab/ajoshi83/generated_mr2')
mr_train_images_path = [os.path.join(mr_path, f) for f in train_images]
mr_val_images_path = [os.path.join(mr_path, f) for f in val_images]
mr_test_images_path = [os.path.join(mr_path, f) for f in test_images]
mr_train_files = [{"img": img, "label": label} for img, label in zip(mr_train_images_path, train_labels_array)]
mr_val_files = [{"img": img, "label": label} for img, label in zip(mr_val_images_path, val_labels_array)]
mr_test_files = [{"img": img, "label": label} for img, label in zip(mr_test_images_path, test_labels_array)]
return ct_train_files, ct_val_files, mr_train_files, mr_val_files, ct_test_files, mr_test_files
#monai.config.print_config()
# old_stdout = sys.stdout
# log_file = open("mri_classification/logs/5fold_eval_fakemr.log","w")
# sys.stdout = log_file
print("this will be written to message.log")
# file_handler = logging.FileHandler(filename='logs/tmp.log')
# #stdout_handler = logging.StreamHandler(stream=sys.stdout)
# handlers = [file_handler]
# logging.basicConfig(
# level=logging.DEBUG,
# format='[%(asctime)s] {%(filename)s:%(lineno)d} %(levelname)s - %(message)s',
# handlers=handlers
# )
# logger = logging.getLogger('LOGGER_NAME')
# edit the path accordingly
gose = pd.read_csv('mri_classification/Pilot_169_train_val_test_split.csv')
X = []
Y = []
test_subjects = []
test_labels = []
torch.cuda.empty_cache()
for i in range(gose.shape[0]):
subj = gose['Main.GUID'][i]
subj_id = str(subj)[4:]
name = str(subj_id) + '.nii'
if gose['Set'][i]=='Train':
X.append(name)
Y.append(gose['Class'][i])
elif gose['Set'][i]=='Val':
X.append(name)
Y.append(gose['Class'][i])
elif gose['Set'][i]=='Test':
test_subjects.append(name)
test_labels.append(gose['Class'][i])
else:
print("Unknown Set: ", gose['Set'][i])
print("Total Subjects for 5-Fold (Train+Val):", len(X))
print("Test Subjects (kept separate):", len(test_subjects))
print("5-Fold Subjects Information:\n")
print("Subjects of Class 0: ", Y.count(0))
print("Subjects of Class 1: ", Y.count(1))
print("Subjects of Class 2: ", Y.count(2))
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=1)
fold_number = 1
X = np.array(X)
Y = np.array(Y)
train_acc, train_auc = [], []
val_acc, val_auc = [], []
test_acc, test_auc = [], []
y_pred_test = []
y_pred_auc = []
for train_index, test_index in skf.split(X, Y):
print("Starting Fold {}..!".format(fold_number))
x_train_fold, x_test_fold = X[train_index], X[test_index]
y_train_fold, y_test_fold = Y[train_index], Y[test_index]
y_train_fold = list(y_train_fold)
y_test_fold = list(y_test_fold)
x_train_fold = list(x_train_fold)
x_test_fold = list(x_test_fold)
# print("\n")
# print("Fold {} statistics:\n".format(fold_number))
# print("Train Subjects: {}".format(len(y_train_fold)))
# print("Subjects of Class 0: ", y_train_fold.count(0))
# print("Subjects of Class 1: ", y_train_fold.count(1))
# print("Subjects of Class 2: {}\n".format(y_train_fold.count(2)))
# print("Val Subjects: {}".format(len(y_test_fold)))
# print("Subjects of Class 0: ", y_test_fold.count(0))
# print("Subjects of Class 1: ", y_test_fold.count(1))
# print("Subjects of Class 2: ", y_test_fold.count(2))
ct_train_files, ct_val_files, mr_train_files, mr_val_files, ct_test_files, mr_test_files = load_data(x_train_fold, x_test_fold, y_train_fold, y_test_fold, test_subjects, test_labels)
# Define transforms for CT and MR respectively
ct_transforms = Compose(
[
LoadImaged(keys=["img"], ensure_channel_first=True),
#ScaleIntensityRanged(keys=["img"],a_min=0, a_max=85, b_min=0, b_max=1),
NormalizeIntensityd(keys=["img"]),
#Resized(keys=["img"], spatial_size=(96, 96, 96)),
]
)
mr_transforms = Compose(
[
LoadImaged(keys=["img"], ensure_channel_first=True),
NormalizeIntensityd(keys=["img"]),
RandFlipd(keys=["img"], prob=1, spatial_axis=2),
RandRotated(keys=["img"], prob=1, range_x=[0.4,0.4])
#ScaleIntensityRanged(keys=["img"],a_min=0, a_max=85, b_min=0, b_max=1),
]
)
post_pred = Compose([Activations(softmax=True)])
post_label = Compose([AsDiscrete(to_onehot=3)])
# create a training data loader
ct_train_ds = monai.data.Dataset(data=ct_train_files, transform=ct_transforms)
mr_train_ds = monai.data.Dataset(data=mr_train_files, transform=mr_transforms)
combined_train_ds = CustomDataset(d1 = mr_train_ds, d2 = ct_train_ds)
train_loader = DataLoader(combined_train_ds, batch_size=1, shuffle=False, pin_memory=torch.cuda.is_available())
check_data1, check_data2 = monai.utils.misc.first(train_loader)
#print(check_data1["img"].shape, check_data1["label"], check_data2["img"].shape, check_data2["label"])
# create a validation data loader
ct_val_ds = monai.data.Dataset(data=ct_val_files, transform=ct_transforms)
mr_val_ds = monai.data.Dataset(data=mr_val_files, transform=mr_transforms)
combined_val_ds = CustomDataset(d1 = mr_val_ds, d2 = ct_val_ds)
val_loader = DataLoader(combined_val_ds, batch_size=1, shuffle=False, pin_memory=torch.cuda)
# create a test data loader
ct_test_ds = monai.data.Dataset(data=ct_test_files, transform=ct_transforms)
mr_test_ds = monai.data.Dataset(data=mr_test_files, transform=mr_transforms)
combined_test_ds = CustomDataset(d1 = mr_test_ds, d2 = ct_test_ds)
test_loader = DataLoader(combined_test_ds, batch_size=1, shuffle=False, pin_memory=torch.cuda)
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
#print(device)
model_ct = resnet18(pretrained=False, spatial_dims=3, n_input_channels=1, num_classes=3).to(device)
#model_ct.load_state_dict(torch.load("models/best_val_resnet18.pth"))
model_mr = resnet18(pretrained=False, spatial_dims=3, n_input_channels=1, num_classes=3).to(device)
#model_mr.load_state_dict(torch.load("models/best_mr_val_resnet18.pth"))
# Freeze these models
for param in model_mr.parameters():
param.requires_grad_(True)
for param in model_ct.parameters():
param.requires_grad_(True)
# Create ensemble model
model = MyEnsemble(model_ct, model_mr)
#model.load_state_dict(torch.load("mri_classification/models/5-Fold/Final_Fold{}_val_resnet18.pth".format(int(fold_number))))
#model.load_state_dict(torch.load("mri_classification/models_gen/5-Fold/Fake_MR_Fold{}_ensemble.pth".format(int(fold_number))))
model.load_state_dict(torch.load("/data/amciilab/ajoshi83/models/MR_Fold{}.pth".format(fold_number)))
#print("Ensemble Model arch check:")
#print(model)
loss_function = torch.nn.CrossEntropyLoss()
optimizer = torch.optim.Adam(model.parameters(),lr= 1e-5, weight_decay=1e-3)
auc_metric = ROCAUCMetric(average="weighted")
# starting evaluation
val_interval = 1
best_metric = -1
best_metric_epoch = -1
best_val_loss = 2
writer = SummaryWriter()
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/train_output")
for batch_data in train_loader:
#step += 1
#print(batch_data["img"])
ct_batch, mr_batch = batch_data[0], batch_data[1]
#print(ct_batch)
#print(mr_batch)
ct_inputs, ct_labels = ct_batch["img"].to(device), ct_batch["label"].to(device)
mr_inputs, mr_labels = mr_batch["img"].to(device), mr_batch["label"].to(device)
train_outputs = model(ct_inputs, mr_inputs).argmax(dim=1)
y_pred = torch.cat([y_pred, model(ct_inputs, mr_inputs)], dim=0)
y = torch.cat([y, mr_labels], dim=0)
value = torch.eq(train_outputs, mr_labels)
metric_count += len(value)
num_correct += value.sum().item()
#saver.save_batch(train_outputs, train_data["img"].meta)
metric = num_correct / metric_count
# print("val evaluation metric:", metric)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred, detach = True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
saver.finalize()
train_acc.append(round(metric,3))
train_auc.append(round(auc_result,3))
print("Train Set: Accuracy: {:.4f}, AUC: {:.4f}".format(metric, auc_result))
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/val_output")
for val_data in val_loader:
ct_val_data, mr_val_data = val_data[0], val_data[1]
ct_val_images, ct_val_labels = ct_val_data["img"].to(device), ct_val_data["label"].to(device)
mr_val_images, mr_val_labels = mr_val_data["img"].to(device), mr_val_data["label"].to(device)
y_pred = torch.cat([y_pred, model(ct_val_images, mr_val_images)], dim=0)
#print(y_pred)
y = torch.cat([y, mr_val_labels], dim=0)
acc_value = torch.eq(y_pred.argmax(dim=1), y)
acc_metric = acc_value.sum().item() / len(acc_value)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred,detach=True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
saver.finalize()
val_acc.append(round(acc_metric,3))
val_auc.append(round(auc_result,3))
print("Val Set: Accuracy: {:.4f}, AUC: {:.4f}".format(acc_metric, auc_result))
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/test_output")
y_pred_mid = []
for test_data in test_loader:
ct_test_data, mr_test_data = test_data[0], test_data[1]
ct_test_images, ct_test_labels = ct_test_data["img"].to(device), ct_test_data["label"].to(device)
mr_test_images, mr_test_labels = mr_test_data["img"].to(device), mr_test_data["label"].to(device)
test_outputs = model(ct_test_images, mr_test_images).argmax(dim=1)
y_pred_mid.append(test_outputs.cpu().numpy()[0])
#softmax_op = torch.nn.functional.softmax(model(ct_test_images, mr_test_images)).argmax(dim=1)
#y_pred_mid.append(softmax_op.cpu().numpy()[0])
y_pred = torch.cat([y_pred, model(ct_test_images, mr_test_images)], dim=0)
y = torch.cat([y, mr_test_labels], dim=0)
value = torch.eq(test_outputs, mr_test_labels)
metric_count += len(value)
num_correct += value.sum().item()
#saver.save_batch(test_outputs, test_data["img"].meta)
y_pred_test.append(y_pred_mid)
metric = num_correct / metric_count
# print("test evaluation metric:", metric)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred, detach=True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
y_pred_auc.append(y_pred_act)
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
# auc_score =roc_auc_score(y_onehot, y_pred_act, multi_class='ovr', average='weighted')
# print(auc_score)
saver.finalize()
test_acc.append(round(metric,3))
test_auc.append(round(auc_result,3))
print("Test Set: Accuracy: {:.4f}, AUC: {:.4f}".format(metric, auc_result))
print("Fold {} completed...Next Fold starting...".format(fold_number))
print("\n")
fold_number += 1
writer.close()
print("\n")
print("Train Acc: {:.4f} +/- {:.4f}, Train AUC: {:.4f} +/- {:.4f}".format(mean(train_acc), np.std(train_acc), mean(train_auc), np.std(train_auc)))
print("Val Acc: {:.4f} +/- {:.4f}, Val AUC: {:.4f} +/- {:.4f}".format(mean(val_acc), np.std(val_acc), mean(val_auc), np.std(val_auc)))
# # print("Val Accuracies of 5 Folds:", val_acc)
# print("Val AUCs of 5 Folds:", val_auc)
print("Test Acc: {:.4f} +/- {:.4f}, Test AUC: {:.4f} +/- {:.4f}".format(mean(test_acc), np.std(test_acc), mean(test_auc), np.std(test_auc)))
# print("Test Accuracies of 5 Folds:", test_acc)
# print("Test AUCs of 5 Folds:", test_auc)
y_pred_tr = np.transpose(y_pred_test)
final = []
for i in range(y_pred_tr.shape[0]):
final.append(mode(y_pred_tr[i]))
y_true = y.cpu().numpy()
final_np = np.array(final)
y_pred_auc = np.array(y_pred_auc)
sum = np.sum(y_pred_auc, axis=0)
sum = sum/5
target_names = ['class 0', 'class 1', 'class 2']
print("\n")
print("Test Accuracy (Hard Voting of 5 models):", round(accuracy_score(y_true, final_np), 4))
print("Test AUC (Hard Voting of 5 models):", round(roc_auc_score(y_onehot, sum, multi_class='ovr', average='weighted'), 4))
cm = confusion_matrix(y_true, final_np)
print(classification_report(y_true, final_np, target_names=target_names))
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot()
plt.show()
# %%
# %%
# %%
import logging
import os
import sys
from pathlib import Path
import pandas as pd
import numpy as np
import torch
from torch.utils.tensorboard import SummaryWriter
import monai
from monai.data import decollate_batch, CSVSaver, DataLoader
from monai.metrics import ROCAUCMetric
from monai.transforms import ScaleIntensityRanged, NormalizeIntensityd, Activations, AsDiscrete, Compose, LoadImaged, RandRotate90d, Resized, ScaleIntensityd
from monai.networks.nets import resnet10, resnet18, resnet34, resnet50
import torch
import torch.nn as nn
from torch.utils.data import DataLoader as _TorchDataLoader
from torch.utils.data import Dataset
from monai.data.utils import list_data_collate, set_rnd, worker_init_fn
from sklearn.model_selection import StratifiedKFold
from statistics import mean, mode
from sklearn.metrics import accuracy_score
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, classification_report, accuracy_score, roc_auc_score
class CustomDataset(monai.data.Dataset):
def __init__(self, d1, d2):
self.d1 = d1
self.d2 = d2
def __getitem__(self, idx):
dict1 = self.d1.__getitem__(idx)
image1, label1 = dict1["img"], dict1["label"]
dict2 = self.d2.__getitem__(idx)
image2, label2 = dict2["img"], dict2["label"]
assert label1==label2
dict_1 = dict()
dict_1["img"] = image1
dict_1["label"] = label1
dict_2 = dict()
dict_2["img"] = image2
dict_2["label"] = label2
return dict_1, dict_2
def __len__(self):
return len(self.d1)
class MyEnsemble(nn.Module):
def __init__(self, modelA, modelB):
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
super(MyEnsemble, self).__init__()
self.modelA = modelA
self.modelB = modelB
# Remove last linear layer
self.modelA.fc = nn.Identity()
self.modelB.fc = nn.Identity()
# Create new classifier
self.mlp1 = nn.Linear(1024,256).to(device)
self.mlp2 = nn.Linear(256,32).to(device)
self.classifier = nn.Linear(32,3).to(device)
def forward(self, i1, i2):
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
x1 = self.modelA(i1) # clone to make sure x is not changed by inplace methods
x1 = x1.view(x1.size(0), -1).to(device)
x2 = self.modelB(i2)
x2 = x2.view(x2.size(0), -1).to(device)
x = torch.cat((x1, x2), dim=1).to(device)
x = nn.functional.relu(self.mlp1(x)).to(device)
x = nn.functional.relu(self.mlp2(x)).to(device)
x = self.classifier(x).to(device)
return x
def load_data(train_images, val_images, train_labels, val_labels, test_images, test_labels):
train_labels_array = np.array(train_labels, dtype=np.int64)
val_labels_array = np.array(val_labels, dtype=np.int64)
test_labels_array = np.array(test_labels, dtype=np.int64)
#loading the MR images
mr_path = Path('Data/Reg_MR')
mr_train_images_path = [os.path.join(mr_path, f) for f in train_images]
mr_val_images_path = [os.path.join(mr_path, f) for f in val_images]
mr_test_images_path = [os.path.join(mr_path, f) for f in test_images]
mr_train_files = [{"img": img, "label": label} for img, label in zip(mr_train_images_path, train_labels_array)]
mr_val_files = [{"img": img, "label": label} for img, label in zip(mr_val_images_path, val_labels_array)]
mr_test_files = [{"img": img, "label": label} for img, label in zip(mr_test_images_path, test_labels_array)]
return mr_train_files, mr_val_files, mr_test_files
#monai.config.print_config()
# old_stdout = sys.stdout
# log_file = open("mri_classification/logs/MR_5Fold_eval.log","w")
# sys.stdout = log_file
# edit the path accordingly
gose = pd.read_csv('mri_classification/Pilot_169_train_val_test_split.csv')
X = []
Y = []
test_subjects = []
test_labels = []
torch.cuda.empty_cache()
for i in range(gose.shape[0]):
subj = gose['Main.GUID'][i]
subj_id = str(subj)[4:]
name = str(subj_id) + '.nii'
if gose['Set'][i]=='Train':
X.append(name)
Y.append(gose['Class'][i])
elif gose['Set'][i]=='Val':
X.append(name)
Y.append(gose['Class'][i])
elif gose['Set'][i]=='Test':
test_subjects.append(name)
test_labels.append(gose['Class'][i])
else:
print("Unknown Set: ", gose['Set'][i])
print("Total Subjects for 5-Fold (Train+Val):", len(X))
print("Test Subjects (kept separate):", len(test_subjects))
print("5-Fold Subjects Information:\n")
print("Subjects of Class 0: ", Y.count(0))
print("Subjects of Class 1: ", Y.count(1))
print("Subjects of Class 2: ", Y.count(2))
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=1)
fold_number = 1
X = np.array(X)
Y = np.array(Y)
train_acc, train_auc = [], []
val_acc, val_auc = [], []
test_acc, test_auc = [], []
y_pred_test = []
y_pred_auc = []
for train_index, val_index in skf.split(X, Y):
print("Starting Fold {}..!".format(fold_number))
# if fold_number <5:
# fold_number += 1
# continue
x_train_fold, x_val_fold = X[train_index], X[val_index]
y_train_fold, y_val_fold = Y[train_index], Y[val_index]
y_train_fold = list(y_train_fold)
y_val_fold = list(y_val_fold)
x_train_fold = list(x_train_fold)
x_val_fold = list(x_val_fold)
# test_subjects = list(test_subjects)
# test_labels = list(test_labels)
print("\n")
print("Fold {} statistics:\n".format(fold_number))
print("Train Subjects: {}".format(len(y_train_fold)))
print("Subjects of Class 0: ", y_train_fold.count(0))
print("Subjects of Class 1: ", y_train_fold.count(1))
print("Subjects of Class 2: {}\n".format(y_train_fold.count(2)))
print("val Subjects: {}".format(len(y_val_fold)))
print("Subjects of Class 0: ", y_val_fold.count(0))
print("Subjects of Class 1: ", y_val_fold.count(1))
print("Subjects of Class 2: ", y_val_fold.count(2))
mr_train_files, mr_val_files, mr_test_files = load_data(x_train_fold, x_val_fold, y_train_fold, y_val_fold, test_subjects, test_labels)
mr_transforms = Compose(
[
LoadImaged(keys=["img"], ensure_channel_first=True),
NormalizeIntensityd(keys=["img"]),
]
)
post_pred = Compose([Activations(softmax=True)])
post_label = Compose([AsDiscrete(to_onehot=3)])
# create a training data loader
mr_train_ds = monai.data.Dataset(data=mr_train_files, transform=mr_transforms)
mr_train_loader = DataLoader(mr_train_ds, batch_size=1, shuffle=False, num_workers=4, pin_memory=torch.cuda.is_available())
check_data1 = monai.utils.misc.first(mr_train_loader)
print(check_data1["img"].shape, check_data1["label"])
# create a validation data loader
mr_val_ds = monai.data.Dataset(data=mr_val_files, transform=mr_transforms)
mr_val_loader = DataLoader(mr_val_ds, batch_size=1, shuffle=False, num_workers=4, pin_memory=torch.cuda.is_available())
# create a test data loader
mr_test_ds = monai.data.Dataset(data= mr_test_files, transform=mr_transforms)
mr_test_loader = DataLoader(mr_test_ds, batch_size=1, shuffle=False, num_workers=4, pin_memory=torch.cuda.is_available())
# Create Model, CrossEntropy Loss and Adam optimizer
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print(device)
model = resnet18(pretrained=False, spatial_dims=3, n_input_channels=1, num_classes=3).to(device)
#model.load_state_dict(torch.load("mri_classification/models/5-Fold/MR_Fold{}.pth".format(int(fold_number))))
model.load_state_dict(torch.load("/data/amciilab/ajoshi83/models/MR_Fold{}.pth".format(fold_number)))
loss_function = torch.nn.CrossEntropyLoss()
optimizer = torch.optim.Adam(model.parameters(),lr= 1e-5, weight_decay=1e-3)
auc_metric = ROCAUCMetric(average="weighted")
# starting evaluation
val_interval = 1
best_metric = -1
best_metric_epoch = -1
best_val_loss = 2
writer = SummaryWriter()
# with torch.no_grad():
# num_correct = 0.0
# metric_count = 0
# y_pred = torch.tensor([], dtype=torch.float32, device=device)
# y = torch.tensor([], dtype=torch.long, device=device)
# saver = CSVSaver(output_dir="./train_output")
# for batch_data in mr_train_loader:
# inputs, labels = batch_data["img"].to(device), batch_data["label"].to(device)
# train_outputs = model(inputs).argmax(dim=1)
# y_pred = torch.cat([y_pred, model(inputs)], dim=0)
# y = torch.cat([y, labels], dim=0)
# value = torch.eq(train_outputs, labels)
# metric_count += len(value)
# num_correct += value.sum().item()
# #saver.save_batch(train_outputs, train_data["img"].meta)
# metric = num_correct / metric_count
# # print("val evaluation metric:", metric)
# y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
# y_pred_act = [post_pred(i) for i in decollate_batch(y_pred, detach = True)]
# y_onehot = torch.stack(y_onehot, dim=0)
# y_pred_act = torch.stack(y_pred_act, dim=0)
# y_onehot = y_onehot.to(device="cpu")
# y_pred_act = y_pred_act.to(device="cpu")
# auc_metric(y_pred_act, y_onehot)
# auc_result = auc_metric.aggregate()
# auc_metric.reset()
# saver.finalize()
# train_acc.append(round(metric,3))
# train_auc.append(round(auc_result,3))
# print("Train Set: Accuracy: {:.4f}, AUC: {:.4f}".format(metric, auc_result))
# with torch.no_grad():
# num_correct = 0.0
# metric_count = 0
# y_pred = torch.tensor([], dtype=torch.float32, device=device)
# y = torch.tensor([], dtype=torch.long, device=device)
# saver = CSVSaver(output_dir="./val_output")
# for val_data in mr_val_loader:
# val_images, val_labels = val_data["img"].to(device), val_data["label"].to(device)
# y_pred = torch.cat([y_pred, model(val_images)], dim=0)
# #print(y_pred)
# y = torch.cat([y, val_labels], dim=0)
# acc_value = torch.eq(y_pred.argmax(dim=1), y)
# acc_metric = acc_value.sum().item() / len(acc_value)
# y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
# y_pred_act = [post_pred(i) for i in decollate_batch(y_pred,detach=True)]
# y_onehot = torch.stack(y_onehot, dim=0)
# y_pred_act = torch.stack(y_pred_act, dim=0)
# y_onehot = y_onehot.to(device="cpu")
# y_pred_act = y_pred_act.to(device="cpu")
# auc_metric(y_pred_act, y_onehot)
# auc_result = auc_metric.aggregate()
# auc_metric.reset()
# saver.finalize()
# val_acc.append(round(acc_metric,3))
# val_auc.append(round(auc_result,3))
# print("Val Set: Accuracy: {:.4f}, AUC: {:.4f}".format(acc_metric, auc_result))
# with torch.no_grad():
# num_correct = 0.0
# metric_count = 0
# y_pred = torch.tensor([], dtype=torch.float32, device=device)
# y = torch.tensor([], dtype=torch.long, device=device)
# #saver = CSVSaver(output_dir="./test_output")
# for test_data in mr_test_loader:
# mr_test_images, mr_test_labels = test_data["img"].to(device), test_data["label"].to(device)
# test_outputs = model(mr_test_images).argmax(dim=1)
# y_pred = torch.cat([y_pred, model(mr_test_images)], dim=0)
# y = torch.cat([y, mr_test_labels], dim=0)
# value = torch.eq(test_outputs, mr_test_labels)
# metric_count += len(value)
# num_correct += value.sum().item()
# #saver.save_batch(test_outputs, test_data["img"].meta)
# metric = num_correct / metric_count
# # print("test evaluation metric:", metric)
# y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
# y_pred_act = [post_pred(i) for i in decollate_batch(y_pred, detach=True)]
# y_onehot = torch.stack(y_onehot, dim=0)
# y_pred_act = torch.stack(y_pred_act, dim=0)
# y_onehot = y_onehot.to(device="cpu")
# y_pred_act = y_pred_act.to(device="cpu")
# auc_metric(y_pred_act, y_onehot)
# auc_result = auc_metric.aggregate()
# auc_metric.reset()
# saver.finalize()
# test_acc.append(round(metric,3))
# test_auc.append(round(auc_result,3))
# print("Test Set: Accuracy: {:.4f}, AUC: {:.4f}".format(metric, auc_result))
# print("Fold {} completed...Next Fold starting...".format(fold_number))
# fold_number += 1
# writer.close()
with torch.no_grad():
num_correct = 0.0
metric_count = 0
y_pred = torch.tensor([], dtype=torch.float32, device=device)
y = torch.tensor([], dtype=torch.long, device=device)
saver = CSVSaver(output_dir="mri_classification/test_output")
y_pred_mid = []
for test_data in mr_test_loader:
mr_test_images, mr_test_labels = test_data["img"].to(device), test_data["label"].to(device)
test_outputs = model(mr_test_images).argmax(dim=1)
#y_pred_mid.append(torch.nn.functional.softmax(model(mr_test_images), dim=1).cpu().numpy()[0])
#y_pred_mid.append(test_outputs.cpu().numpy()[0])
softmax_op = torch.nn.functional.softmax(model(mr_test_images)).argmax(dim=1)
y_pred_mid.append(softmax_op.cpu().numpy()[0])
y_pred = torch.cat([y_pred, model(mr_test_images)], dim=0)
y = torch.cat([y, mr_test_labels], dim=0)
value = torch.eq(test_outputs, mr_test_labels)
metric_count += len(value)
num_correct += value.sum().item()
#saver.save_batch(test_outputs, test_data["img"].meta)
y_pred_test.append(y_pred_mid)
# print(num_correct)
# print(metric_count)
metric = num_correct / metric_count
# print("test evaluation metric:", metric)
y_onehot = [post_label(i) for i in decollate_batch(y, detach=True)]
y_pred_act = [post_pred(i) for i in decollate_batch(y_pred, detach=True)]
y_onehot = torch.stack(y_onehot, dim=0)
y_pred_act = torch.stack(y_pred_act, dim=0)
y_onehot = y_onehot.to(device="cpu")
y_pred_act = y_pred_act.to(device="cpu")
auc_metric(y_pred_act, y_onehot)
auc_result = auc_metric.aggregate()
auc_metric.reset()
saver.finalize()
test_acc.append(round(metric,3))
test_auc.append(round(auc_result,3))
print("Test Set: Accuracy: {:.4f}, AUC: {:.4f}".format(metric, auc_result))
print("Fold {} completed...Next Fold starting...".format(fold_number))
fold_number += 1
#writer.close()
print("\n")
# print("Train Acc: {:.4f} +/- {:.4f}, Train AUC: {:.4f} +/- {:.4f}".format(mean(train_acc), np.std(train_acc), mean(train_auc), np.std(train_auc)))
# print("Val Acc: {:.4f} +/- {:.4f}, Val AUC: {:.4f} +/- {:.4f}".format(mean(val_acc), np.std(val_acc), mean(val_auc), np.std(val_auc)))
# # # print("Val Accuracies of 5 Folds:", val_acc)
# print("Val AUCs of 5 Folds:", val_auc)
print("Test Acc: {:.4f} +/- {:.4f}, Test AUC: {:.4f} +/- {:.4f}".format(mean(test_acc), np.std(test_acc), mean(test_auc), np.std(test_auc)))
# print("Test Accuracies of 5 Folds:", test_acc)
# print("Test AUCs of 5 Folds:", test_auc)
y_pred_tr = np.transpose(y_pred_test)
final = []
for i in range(y_pred_tr.shape[0]):
final.append(mode(y_pred_tr[i]))
y_true = y.cpu().numpy()
final_np = np.array(final)
y_pred_auc = np.array(y_pred_auc)
sum = np.sum(y_pred_auc, axis=0)
sum = sum/5
target_names = ['class 0', 'class 1', 'class 2']
print("\n")
print("Test Accuracy (Hard Voting of 5 models):", round(accuracy_score(y_true, final_np), 4))
#print("Test AUC (Hard Voting of 5 models):", round(roc_auc_score(y_onehot, sum, multi_class='ovr', average='weighted'), 4))
cm = confusion_matrix(y_true, final_np)
print(classification_report(y_true, final_np, target_names=target_names))
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot()
plt.show()
# %%
torch.nn.functional.softmax(model(mr_test_images), dim=1).detach().cpu().numpy()[0]
# %%
y_pred_tr.shape
# %%
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot()
plt.show()
# %%
roc_auc_score(y_onehot, sum, multi_class='ovr', average='weighted')
# %%
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, classification_report, accuracy_score
target_names = ['class 0', 'class 1', 'class 2']
print("Test Accuracy (Hard Voting of 5 models):", accuracy_score(y.cpu().numpy(), final_np))
cm = confusion_matrix(y.cpu().numpy(), final_np)
print(classification_report(y.cpu().numpy(), final_np, target_names=target_names))
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot()
plt.show()
# %%

analysis.ipynb at commit ae8ad46, no license · at the source

Overview

Authors: Yiming Che^1,2, Amogh Manoj Joshi^1,2, Jay Shah^1,2, Md Mahfuzur Rahman Siddiquee^1,2, Catherine D Chong^2,3, Simona Nikolova³, Gina Dumkrieger³, Baoxin Li^1,2, Teresa Wu^1,2, Todd J Schwedt^2,3

ORCID iDs: Jay Shah, Md Mahfuzur Rahman Siddiquee

Institutions: Arizona State University (United States); Mayo Clinic (United States); Mayo Clinic in Arizona (United States); Mayo Clinic Hospital (United States)

Journal: Brain communications, volume 8, issue 2, article fcag123

Dates: received 1 May 2025; accepted 2 April 2026; published online 6 April 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1093/braincomms/fcag123 · PMID 42004011 · PMCID PMC13084558 · OpenAlex W7152176423

Open access: gold, a free copy (OpenAlex)

Status: code verified (how the code is found and verified)

Categories: structural MRI / diffusion (modality), other (modality), human (organism), traumatic brain injury (population)

Methods: Connectivity, Machine learning, Statistics

Keywords: deep learning, traumatic brain injury, concussion, neuroimaging, data harmonization

Topic: Traumatic Brain Injury Research (Clinical Neurology, Medicine), according to OpenAlex

Funding: U.S. Department of Defense

Citations: cited by 1 paper (Europe PMC); 52 references in the paper

Abstract

To enhance the prediction of traumatic brain injury (mTBI) outcomes, we propose a deep learning approach that integrates brain computed tomography (CT) scans with corresponding synthetic T1-weighted magnetic resonance imaging (T1-MRI). Our method significantly outperforms the prediction using CT scans alone. TRACK-TBI Pilot dataset, which includes imaging and clinical outcome data from patients with TBI, is studied. The hypothesis is brain CT and T1-MRI complement each other and together will improve TBI prognosis compared to using either CT or T1-MRI alone. Since CT and T1-MRI may not be available for the same individual, we employed a specialized version of a generative adversarial network (GAN), known as fixed-point GAN (FP-GAN). FP-GAN was trained using unpaired CT and T1-MRI scans to generate synthetic T1-MRIs from real CT scans. This process produced pseudo-paired CT-MRI data, which was then used to train a deep learning classifier for outcome prediction. The classifier consists of dual parallel 3D ResNet-18 models, each independently processing T1-MRI and CT scans. We used Glasgow Outcome Scale-Extended (GOSE) scores at 3 months post-TBI as the measure of patient outcomes. To avoid data leakage, the subjects used in FP-GAN and ResNet-18 model have no overlap. We further divided the paired data, allocating 69 samples for 5-fold cross-validation and 17 samples for testing. Prognostic performance was evaluated using the area under the receiver operating characteristic curve (AUC), F1-score (the harmonic mean of precision and recall), sensitivity (true positive rate) and specificity (true negative rate). For binary classification, we defined good recovery as GOSE ≥ 7 (positive) and poor recovery (negative) as 3 ≤ GOSE ≤ 6. Accordingly, our training set consists of 24 subjects with poor recovery and 45 subjects with good recovery, while the testing set includes 5 subjects with poor recovery and 12 subjects with good recovery. A DeLong test on AUC confirms that the improvement from incorporating synthetic T1-MRI (AUC = 0.76 ± 0.10) is statistically significant (P < 0.05) compared to using CT alone (AUC = 0.68 ± 0.13). The significant improvement from using the combination of real CT and synthetic T1-MRI in sensitivity (SEN = 0.95 ± 0.07) and overall performance metrics, such as F1-score (F1 = 0.84 ± 0.03), suggests that the proposed approach provides a robust and effective prognostic approach compared to using CT alone (SEN = 0.83 ± 0.18 and F1 = 0.76 ± 0.07). This pilot research demonstrates the potential of a deep learning-based harmonization model to bridge the gap between CT and T1-MRI in TBI assessment. By integrating synthetic T1-MRI with CT, prediction performance is substantially enhanced.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

SoloChe/TBI-Recovery-Prediction-Harmonization

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: ae8ad4680577f45f5178d327ad6ae18f6202d7ca, 8 March 2026

Languages: Python (30), Shell (4), Jupyter (3)

Size: 39 files, 37 scripts

Software Heritage: not archived

Found in: “Data availability”

Holds: README, 2 notebooks

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: PyTorch (32 files), NumPy (30 files), scikit-learn (29 files), MONAI (28 files), pandas (28 files), Matplotlib (4 files), Pillow (3 files), SciPy (2 files), scikit-image (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

38 files

Stage_1_FPGAN/data_loader.py, Python, 152 lines
Stage_1_FPGAN/main.py, Python, 171 lines
Stage_1_FPGAN/model.py, Python, 85 lines, 1 match
Stage_1_FPGAN/solver.py, Python, 937 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/5fold_CT_double_arch-checkpoint.py, Python, 345 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/5fold_CT_double_arch_eval-checkpoint.py, Python, 389 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/5fold_MR_double_arch-checkpoint.py, Python, 343 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/5fold_MR_double_arch_eval-checkpoint.py, Python, 388 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/5fold_MR_eval-checkpoint.py, Python, 369 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/5fold_MR_train-checkpoint.py, Python, 289 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/5fold_eval-checkpoint.py, Python, 386 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/5fold_eval_fakemr-checkpoint.py, Python, 385 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/5fold_training-checkpoint.py, Python, 308 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/5fold_training_with_weights-checkpoint.py, Python, 336 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/5fold_training_with_weights_fakemr-checkpoint.py, Python, 344 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/analysis-checkpoint.ipynb, Jupyter, 2,110 lines
Stage_2_CLS/5fold/.ipynb_checkpoints/eval_cm-checkpoint.py, Python, 1 line
Stage_2_CLS/5fold/5fold_CT_double_arch.py, Python, 345 lines
Stage_2_CLS/5fold/5fold_CT_double_arch_eval.py, Python, 389 lines
Stage_2_CLS/5fold/5fold_CT_eval.py, Python, 317 lines
Stage_2_CLS/5fold/5fold_CT_train.py, Python, 289 lines
Stage_2_CLS/5fold/5fold_MR_double_arch.py, Python, 344 lines
Stage_2_CLS/5fold/5fold_MR_double_arch_eval.py, Python, 388 lines
Stage_2_CLS/5fold/5fold_MR_eval.py, Python, 369 lines
Stage_2_CLS/5fold/5fold_MR_train.py, Python, 289 lines
Stage_2_CLS/5fold/5fold_eval.py, Python, 386 lines
Stage_2_CLS/5fold/5fold_eval_fakemr.py, Python, 386 lines
Stage_2_CLS/5fold/5fold_training.py, Python, 308 lines
Stage_2_CLS/5fold/5fold_training_with_weights.py, Python, 336 lines
Stage_2_CLS/5fold/5fold_training_with_weights_fakemr.py, Python, 344 lines
Stage_2_CLS/5fold/analysis.ipynb, Jupyter, 3,122 lines, 1 match
Stage_2_CLS/5fold/analysis2.ipynb, Jupyter, 147 lines
Stage_2_CLS/5fold/eval.sh, Shell, 34 lines
Stage_2_CLS/5fold/eval2.sh, Shell, 84 lines
Stage_2_CLS/5fold/eval_actual.sh, Shell, 42 lines
Stage_2_CLS/5fold/eval_cm.py, Python, 377 lines
Stage_2_CLS/5fold/eval_tri.sh, Shell, 56 lines
README.md, Text, 5 lines

At the source: github.com/SoloChe/TBI-Recovery-Prediction-Harmonization

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
37 scripts, each with its path and the digest of its content;
2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

All data used in this study were obtained from the Federal Interagency Traumatic Brain Injury Research (FITBIR) Informatics System. Access to the FITBIR datasets requires proper authorization and adherence to their data use agreements. Researchers interested in accessing the data can apply through the FITBIR Data Access Request process. Data sharing is not applicable to this article as no new data were created or analysed in this study. Model code, training scripts and pre/postprocessing pipelines are available at: https://github.com/SoloChe/TBI-Recovery-Prediction-Harmonization.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 10 authors, 5 keywords, 1 funder, 42 references.

Cite

This paper

Che, Y., Joshi, A. M., Shah, J., Rahman Siddiquee, M. M., Chong, C. D., Nikolova, S., Dumkrieger, G., Li, B., Wu, T., & Schwedt, T. J. (2026). Traumatic brain injury recovery prediction by harmonizing real brain CT and synthetic brain MRI: a pilot study. Brain communications, 8(2), fcag123. https://doi.org/10.1093/braincomms/fcag123

BibTeX

@article{che2026traumatic, author = {Che, Yiming and Joshi, Amogh Manoj and Shah, Jay and Rahman Siddiquee, Md Mahfuzur and Chong, Catherine D and Nikolova, Simona and Dumkrieger, Gina and Li, Baoxin and Wu, Teresa and Schwedt, Todd J}, title = {{Traumatic brain injury recovery prediction by harmonizing real brain CT and synthetic brain MRI: a pilot study}}, journal = {Brain communications}, year = {2026}, month = apr, volume = {8}, number = {2}, pages = {fcag123}, publisher = {Oxford University Press}, issn = {2632-1297}, doi = {10.1093/braincomms/fcag123}, url = {https://doi.org/10.1093/braincomms/fcag123}, pmid = {42004011}, pmcid = {PMC13084558} }

RIS

TY - JOUR AU - Che, Yiming AU - Joshi, Amogh Manoj AU - Shah, Jay AU - Rahman Siddiquee, Md Mahfuzur AU - Chong, Catherine D AU - Nikolova, Simona AU - Dumkrieger, Gina AU - Li, Baoxin AU - Wu, Teresa AU - Schwedt, Todd J TI - Traumatic brain injury recovery prediction by harmonizing real brain CT and synthetic brain MRI: a pilot study T2 - Brain communications J2 - Brain Commun PY - 2026 DA - 2026/04/06 VL - 8 IS - 2 SP - fcag123 SN - 2632-1297 PB - Oxford University Press DO - 10.1093/braincomms/fcag123 UR - https://doi.org/10.1093/braincomms/fcag123 LA - en ER -

CSL-JSON

{ "id": "10.1093/braincomms/fcag123", "type": "article-journal", "title": "Traumatic brain injury recovery prediction by harmonizing real brain CT and synthetic brain MRI: a pilot study", "container-title": "Brain communications", "author": [ { "family": "Che", "given": "Yiming" }, { "family": "Joshi", "given": "Amogh Manoj" }, { "family": "Shah", "given": "Jay" }, { "family": "Rahman Siddiquee", "given": "Md Mahfuzur" }, { "family": "Chong", "given": "Catherine D" }, { "family": "Nikolova", "given": "Simona" }, { "family": "Dumkrieger", "given": "Gina" }, { "family": "Li", "given": "Baoxin" }, { "family": "Wu", "given": "Teresa" }, { "family": "Schwedt", "given": "Todd J" } ], "container-title-short": "Brain Commun", "volume": "8", "issue": "2", "page": "fcag123", "DOI": "10.1093/braincomms/fcag123", "PMID": "42004011", "PMCID": "PMC13084558", "ISSN": "2632-1297", "publisher": "Oxford University Press", "URL": "https://doi.org/10.1093/braincomms/fcag123", "language": "en", "issued": { "date-parts": [ [ 2026, 4, 6 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown

<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1093/braincomms/fcag123.svg)](https://openscicode.org/badge/doi/10.1093/braincomms/fcag123)
<!-- oscr-badge:end -->

reStructuredText

.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1093/braincomms/fcag123.svg
   :target: https://openscicode.org/badge/doi/10.1093/braincomms/fcag123
   :alt: OSCR

.. oscr-badge:end

HTML

<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1093/braincomms/fcag123"><img src="https://openscicode.org/badge/doi/10.1093/braincomms/fcag123.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

SoloChe/TBI-Recovery-Prediction-Harmonization: README.md on GitHub

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.