Fixed results comparison and added duration

Author: etpeterson

conftest.py

@@ -2,7 +2,7 @@
 import pathlib
 import json
 import csv
-import datetime
+# import datetime
 
 
 def pytest_addoption(parser):
@@ -49,6 +49,12 @@ def pytest_addoption(parser):
         type=int,
         help="Number of fits to perform on the same parameters",
     )
+    parser.addoption(
+        "--saveDurationFileName",
+        default="",
+        type=str,
+        help="Saved duration results file name",
+    )
 
 
 @pytest.fixture(scope="session")
@@ -60,9 +66,25 @@ def save_file(request):
         # print(filename)
         # filename.unlink(missing_ok=True)
         filename = filename.as_posix()
-        with open(filename, "a") as csv_file:
-            writer = csv.writer(csv_file, delimiter='#')
-            writer.writerow(["", datetime.datetime.now()])
+        with open(filename, "w") as csv_file:
+            writer = csv.writer(csv_file, delimiter=',')
+            writer.writerow(("Algorithm", "Region", "SNR", "f", "Dp", "D", "f_fitted", "Dp_fitted", "D_fitted"))
+            # writer.writerow(["", datetime.datetime.now()])
+    return filename
+
+@pytest.fixture(scope="session")
+def save_duration_file(request):
+    filename = request.config.getoption("--saveDurationFileName")
+    if filename:
+        current_folder = pathlib.Path.cwd()
+        filename = current_folder / filename
+        # print(filename)
+        # filename.unlink(missing_ok=True)
+        filename = filename.as_posix()
+        with open(filename, "w") as csv_file:
+            writer = csv.writer(csv_file, delimiter=',')
+            writer.writerow(("Algorithm", "Region", "SNR", "Duration [us]", "Count"))
+            # writer.writerow(["", datetime.datetime.now()])
     return filename
 
 @pytest.fixture(scope="session")

src/original/ETP_SRI/LinearFitting.py

@@ -66,12 +66,20 @@ def ivim_fit(self, bvalues, signal):
         gt_cutoff = bvalues >= self.linear_cutoff
         linear_signal = np.log(signal)
         D = self.linear_fit(bvalues[gt_cutoff], linear_signal[gt_cutoff])
+        # print(D)
         D[1] = max(D[1], 0)  # constrain to positive values
 
         if lt_cutoff.sum() > 0:
             signal_Dp = linear_signal[lt_cutoff] - gd.linear_signal(D[1], bvalues[lt_cutoff], np.log(D[0]))
-            
-            Dp_prime = self.linear_fit(bvalues[lt_cutoff], np.log(signal_Dp))
+            # print(signal_Dp)
+            signal_valid = signal_Dp > 0
+            lt_cutoff_dual = np.logical_and(lt_cutoff[:len(signal_Dp)], signal_valid)
+            # print(lt_cutoff_dual)
+            Dp_prime = [-1, -1]
+            if lt_cutoff_dual.sum() > 0:
+                # print(np.log(signal_Dp[lt_cutoff_dual]))
+                Dp_prime = self.linear_fit(bvalues[:len(signal_Dp)][lt_cutoff_dual], np.log(signal_Dp[lt_cutoff_dual]))
+                # print(Dp_prime)
 
             if np.any(np.asarray(Dp_prime) < 0) or not np.all(np.isfinite(Dp_prime)):
                 print('Perfusion fit failed')

tests/IVIMmodels/unit_tests/analyze.r

@@ -1,21 +1,25 @@
 library(tidyverse)
-data <- read.csv("test_output.csv", skip=1, col.names=c("Algorithm", "Region", "SNR", "Truth", "Fitted"))
-data<-separate_wider_delim(data, cols = Truth, delim = ",", names = c("f", "D", "Dp"))
-data$f <- substring(data$f, 2)
-data$Dp <- substring(data$Dp, 1, nchar(data$Dp) - 1)
-data<-separate_wider_delim(data, cols = Fitted, delim = ",", names = c("f_fitted", "D_fitted", "Dp_fitted"))
-data$f_fitted <- substring(data$f_fitted, 2)
-data$Dp_fitted <- substring(data$Dp_fitted, 1, nchar(data$Dp_fitted) - 1)
-data[c("Algorithm", "Region")] <- sapply(data[c("Algorithm", "Region")], as.factor)
-data[c("SNR", "f", "D", "Dp", "f_fitted", "D_fitted", "Dp_fitted")] <- sapply(data[c("SNR", "f", "D", "Dp", "f_fitted", "D_fitted", "Dp_fitted")], as.numeric)
-ggplot(data, aes(x=Algorithm)) + geom_boxplot(aes(y=f_fitted)) + geom_boxplot(color="red", aes(y=f)) + facet_grid(SNR ~ Region) + scale_x_discrete(guide = guide_axis(angle = 90)) + ylim(0, 1) + ggtitle("Perfusion fraction grid") + xlab("Perfusion fraction")
-ggsave("f.pdf", width = 50, height = 50, units = "cm")
-ggplot(data, aes(x=Algorithm)) + geom_boxplot(aes(y=D_fitted)) + geom_boxplot(color="red", aes(y=D)) + facet_grid(SNR ~ Region) + scale_x_discrete(guide = guide_axis(angle = 90)) + ggtitle("Diffusion grid") + xlab("Diffusion")
-ggsave("D.pdf", width = 50, height = 50, units = "cm")
-ggplot(data, aes(x=Algorithm)) + geom_boxplot(aes(y=Dp_fitted)) + geom_boxplot(color="red", aes(y=Dp)) + facet_grid(SNR ~ Region) + scale_x_discrete(guide = guide_axis(angle = 90)) + ggtitle("Perfusion grid") + xlab("Perfusion")
-ggsave("Dp.pdf", width = 50, height = 50, units = "cm")
-#why?
-ggplot(data, aes(x=Algorithm)) + geom_boxplot(aes(y=D_fitted)) + geom_boxplot(color="red", aes(y=Dp)) + facet_grid(SNR ~ Region) + scale_x_discrete(guide = guide_axis(angle = 90)) + ggtitle("Diffusion grid") + xlab("Diffusion")
-ggsave("D_tweak.pdf", width = 50, height = 50, units = "cm")
-ggplot(data, aes(x=Algorithm)) + geom_boxplot(aes(y=Dp_fitted)) + geom_boxplot(color="red", aes(y=D)) + facet_grid(SNR ~ Region) + scale_x_discrete(guide = guide_axis(angle = 90)) + ggtitle("Perfusion grid") + xlab("Perfusion")
-ggsave("Dp_tweak.pdf", width = 50, height = 50, units = "cm")
+library(plyr)
+
+plot_ivim <- function(data, fileExtension) {
+    ggplot(data, aes(x=Algorithm)) + geom_boxplot(aes(y=f_fitted)) + geom_boxplot(color="red", aes(y=f)) + facet_grid(SNR ~ Region) + scale_x_discrete(guide = guide_axis(angle = 90)) + ylim(0, 1) + ggtitle("Perfusion fraction grid") + ylab("Perfusion fraction")
+    ggsave(paste("f", fileExtension, sep=""), width = 50, height = 50, units = "cm")
+    ggplot(data, aes(x=Algorithm)) + geom_boxplot(aes(y=D_fitted)) + geom_boxplot(color="red", aes(y=D)) + facet_grid(SNR ~ Region) + scale_x_discrete(guide = guide_axis(angle = 90)) + ggtitle("Diffusion grid") + ylab("Diffusion")
+    ggsave(paste("D", fileExtension, sep=""), width = 50, height = 50, units = "cm")
+    ggplot(data, aes(x=Algorithm)) + geom_boxplot(aes(y=Dp_fitted)) + geom_boxplot(color="red", aes(y=Dp)) + facet_grid(SNR ~ Region) + scale_x_discrete(guide = guide_axis(angle = 90)) + ylim(0, 0.25) + ggtitle("Perfusion grid") + ylab("Perfusion")
+    ggsave(paste("Dp", fileExtension, sep=""), width = 50, height = 50, units = "cm")
+}
+
+data <- read.csv("test_output.csv")
+data <- data %>% mutate_if(is.character, as.factor)
+plot_ivim(data, ".pdf")
+
+data_restricted <- data[data$Region %in% c("Liver", "spleen", "Right kydney cortex", "right kidney medulla"),]
+plot_ivim(data_restricted, "_limited.pdf")
+
+data_duration <- read.csv("test_duration.csv")
+data_duration <- data_duration %>% mutate_if(is.character, as.factor)
+data_duration$ms <- data_duration$Duration..us./data_duration$Count/1000
+ggplot(data_duration, aes(x=Algorithm, y=ms)) + geom_boxplot() + scale_x_discrete(guide = guide_axis(angle = 90)) + ggtitle("Fit Duration") + ylab("Time (ms)")
+ggsave("durations.pdf", width = 20, height = 20, units = "cm")
+

tests/IVIMmodels/unit_tests/test_ivim_synthetic.py

@@ -6,14 +6,15 @@
 import os
 import csv
 import random
+import datetime
 
 from src.wrappers.OsipiBase import OsipiBase
 from utilities.data_simulation.GenerateData import GenerateData
 
 #run using pytest <path_to_this_file> --saveFileName test_output.txt --SNR 50 100 200
 #e.g. pytest -m slow tests/IVIMmodels/unit_tests/test_ivim_synthetic.py  --saveFileName test_output.csv --SNR 10 50 100 200 --fitCount 20
 @pytest.mark.slow
-def test_generated(ivim_algorithm, ivim_data, SNR, rtol, atol, fit_count, save_file):
+def test_generated(ivim_algorithm, ivim_data, SNR, rtol, atol, fit_count, save_file, save_duration_file):
     # assert save_file == "test"
     random.seed(42)
     S0 = 1
@@ -23,19 +24,30 @@ def test_generated(ivim_algorithm, ivim_data, SNR, rtol, atol, fit_count, save_f
     f = data["f"]
     Dp = data["Dp"]
     fit = OsipiBase(algorithm=ivim_algorithm)
+    time_delta = datetime.timedelta()
     for idx in range(fit_count):
         # if "data" not in data:
         signal = gd.ivim_signal(D, Dp, f, S0, bvals, SNR)
         # else:
         #     signal = data["data"]
+        start_time = datetime.datetime.now()
         [f_fit, Dp_fit, D_fit] = fit.osipi_fit(signal, bvals)
+        time_delta += datetime.datetime.now() - start_time
         if save_file:
-            save_results(save_file, ivim_algorithm, name, SNR, [f, D, Dp], [f_fit, Dp_fit, D_fit])
-        npt.assert_allclose([f, D, Dp], [f_fit, D_fit, Dp_fit], rtol, atol)
+            save_results(save_file, ivim_algorithm, name, SNR, [f, Dp, D], [f_fit, Dp_fit, D_fit])
+        npt.assert_allclose([f, Dp, D], [f_fit, Dp_fit, D_fit], rtol, atol)
+    if save_duration_file:
+        save_duration(save_duration_file, ivim_algorithm, name, SNR, time_delta, fit_count)
 
 
 def save_results(filename, algorithm, name, SNR, truth, fit):
     with open(filename, "a") as csv_file:
         writer = csv.writer(csv_file, delimiter=',')
-        data = [algorithm, name, SNR, truth, fit]
+        data = [algorithm, name, SNR, *truth, *fit]
+        writer.writerow(data)
+
+def save_duration(filename, algorithm, name, SNR, duration, count):
+    with open(filename, "a") as csv_file:
+        writer = csv.writer(csv_file, delimiter=',')
+        data = [algorithm, name, SNR, duration/datetime.timedelta(microseconds=1), count]
         writer.writerow(data)