Added OGC to testing

Author: etpeterson

conftest.py

@@ -68,7 +68,7 @@ def save_file(request):
         filename = filename.as_posix()
         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(("Algorithm", "Region", "SNR", "index", "f", "Dp", "D", "f_fitted", "Dp_fitted", "D_fitted"))
             # writer.writerow(["", datetime.datetime.now()])
     return filename
 

src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py

@@ -26,6 +26,7 @@ class OGC_AmsterdamUMC_Bayesian_biexp(OsipiBase):
     required_initial_guess = False
     required_initial_guess_optional = True
     accepted_dimensions = 1  # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most?
+    accepts_priors = True
 
     def __init__(self, bvalues=None, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]), initial_guess=None, fitS0=True, thresholds=None, prior_in=None):
         """
@@ -38,7 +39,11 @@ def __init__(self, bvalues=None, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.
             Args:
                  datain is a 2D array with values of D, f, D* (and S0) that will form the prior.
         """
-        super(OGC_AmsterdamUMC_Bayesian_biexp, self).__init__(bvalues, bounds,initial_guess,fitS0,prior_in)
+        super(OGC_AmsterdamUMC_Bayesian_biexp, self).__init__(bvalues, thresholds, bounds, initial_guess) #, fitS0, prior_in)
+        self.OGC_algorithm = fit_bayesian
+        self.initialize(bounds, initial_guess, fitS0, thresholds, prior_in)
+
+    def initialize(self, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]), initial_guess=None, fitS0=True, thresholds=None, prior_in=None):
         if bounds is None:
             self.bounds=([0, 0, 0.005, 0.7],[0.005, 1, 0.2, 1.3])
         else:
@@ -50,14 +55,13 @@ def __init__(self, bvalues=None, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.
                 self.neg_log_prior = empirical_neg_log_prior(prior_in[0], prior_in[1], prior_in[2],prior_in[3])
             else:
                 self.neg_log_prior = empirical_neg_log_prior(prior_in[0], prior_in[1], prior_in[2])
-        self.OGC_algorithm = fit_bayesian
         if initial_guess is None:
             self.initial_guess = [0.001, 0.5, 0.1, 1]
         else:
             self.initial_guess = initial_guess
         self.fitS0=fitS0
 
-    def ivim_fit(self, signals, bvalues=None):
+    def ivim_fit(self, signals, bounds, initial_guess, bvalues=None):
         """Perform the IVIM fit
 
         Args:

src/standardized/OGC_AmsterdamUMC_biexp.py

@@ -26,6 +26,7 @@ class OGC_AmsterdamUMC_biexp(OsipiBase):
     required_initial_guess = False
     required_initial_guess_optional = True
     accepted_dimensions = 1  # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most?
+    accepts_priors = True
 
     def __init__(self, bvalues=None, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]), initial_guess=None, fitS0=True, thresholds=None):
         """
@@ -37,6 +38,9 @@ def __init__(self, bvalues=None, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.
         """
         super(OGC_AmsterdamUMC_biexp, self).__init__(bvalues, bounds,initial_guess,fitS0)
         self.OGC_algorithm = fit_least_squares
+        self.initialize(bounds, initial_guess, fitS0)
+
+    def initialize(self, bounds, initial_guess, fitS0):
         if bounds is None:
             self.bounds=([0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3])
         else:
@@ -47,7 +51,7 @@ def __init__(self, bvalues=None, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.
             self.initial_guess = initial_guess
         self.fitS0=fitS0
 
-    def ivim_fit(self, signals, bvalues=None):
+    def ivim_fit(self, signals, bounds, initial_guess, bvalues=None):
         """Perform the IVIM fit
 
         Args:

src/standardized/OGC_AmsterdamUMC_biexp_segmented.py

@@ -26,6 +26,7 @@ class OGC_AmsterdamUMC_biexp_segmented(OsipiBase):
     required_initial_guess = False
     required_initial_guess_optional = True
     accepted_dimensions = 1  # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most?
+    accepts_priors = True
 
     def __init__(self, bvalues=None, thresholds=150, bounds=([0, 0, 0.005],[0.005, 0.7, 0.2]), initial_guess=[0.001, 0.01, 0.01,1]):
         """
@@ -37,6 +38,10 @@ def __init__(self, bvalues=None, thresholds=150, bounds=([0, 0, 0.005],[0.005, 0
         """
         super(OGC_AmsterdamUMC_biexp_segmented, self).__init__(bvalues,thresholds, bounds, initial_guess)
         self.OGC_algorithm = fit_segmented
+        self.initialize(bounds, initial_guess, thresholds)
+        
+
+    def initialize(self, bounds, initial_guess, thresholds):
         if bounds is None:
             self.bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3])
         else:
@@ -47,7 +52,7 @@ def __init__(self, bvalues=None, thresholds=150, bounds=([0, 0, 0.005],[0.005, 0
             self.initial_guess = initial_guess
         self.thresholds=thresholds
 
-    def ivim_fit(self, signals, bvalues=None):
+    def ivim_fit(self, signals, bounds, initial_guess, bvalues=None):
         """Perform the IVIM fit
 
         Args:

src/wrappers/OsipiBase.py

@@ -7,7 +7,7 @@
 class OsipiBase:
     """The base class for OSIPI IVIM fitting"""
 
-    def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, algorithm=None):
+    def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, algorithm=None, **kwargs):
         # Define the attributes as numpy arrays only if they are not None
         self.bvalues = np.asarray(bvalues) if bvalues is not None else None
         self.thresholds = np.asarray(thresholds) if thresholds is not None else None
@@ -17,7 +17,7 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
         # If the user inputs an algorithm to OsipiBase, it is intereprete as initiating
         # an algorithm object with that name.
         if algorithm:
-            self.osipi_initiate_algorithm(algorithm, bvalues=bvalues, thresholds=thresholds, bounds=bounds, initial_guess=initial_guess)
+            self.osipi_initiate_algorithm(algorithm, bvalues=bvalues, thresholds=thresholds, bounds=bounds, initial_guess=initial_guess, **kwargs)
 
     def osipi_initiate_algorithm(self, algorithm, **kwargs):
         """Turns the class into a specified one by the input.
@@ -41,6 +41,10 @@ def osipi_initiate_algorithm(self, algorithm, **kwargs):
         # Change the class from OsipiBase to the specified algorithm
         self.__class__ = getattr(importlib.import_module(import_path), algorithm)
         self.__init__(**kwargs)
+    
+    def initialize(**kwargs):
+        """Placeholder for subclass initialization"""
+        pass
 
     def osipi_fit(self, data=None, bvalues=None, thresholds=None, bounds=None, initial_guess=None, **kwargs):
         """Fits the data with the bvalues
@@ -51,17 +55,17 @@ def osipi_fit(self, data=None, bvalues=None, thresholds=None, bounds=None, initi
         # Then check if they are input here, if they are, these should overwrite the attributes
         use_bvalues = bvalues if bvalues is not None else self.bvalues
         kwargs["bvalues"] = use_bvalues
-        use_thresholds = thresholds if self.bvalues is None else self.thresholds
-        use_bounds = bounds if self.bounds is None else self.bounds
-        use_initial_guess = initial_guess if self.initial_guess is None else self.initial_guess
+        use_thresholds = thresholds if thresholds is not None else self.thresholds
+        use_bounds = bounds if bounds is not None else self.bounds
+        use_initial_guess = initial_guess if initial_guess is not None else self.initial_guess
 
         # Make sure we don't make arrays of None's
         if use_bvalues is not None: use_bvalues = np.asarray(use_bvalues) 
         if use_thresholds is not None: use_thresholds = np.asarray(use_thresholds) 
         if use_bounds is not None: use_bounds = np.asarray(use_bounds) 
         if use_initial_guess is not None: use_initial_guess = np.asarray(use_initial_guess) 
 
-        #args = [data, use_bvalues, use_thresholds]
+        # args = [data, use_bvalues, use_thresholds]
         args = [data, use_thresholds]
         if self.required_bounds or self.required_bounds_optional:
             args.append(use_bounds)
@@ -77,6 +81,8 @@ def osipi_fit(self, data=None, bvalues=None, thresholds=None, bounds=None, initi
 
         #args = [data, use_bvalues, use_initial_guess, use_bounds, use_thresholds]
         args = [arg for arg in args if arg is not None]
+        # print(args)
+        # print(kwargs)
         results = self.ivim_fit(*args, **kwargs)
 
         #self.parameter_estimates = self.ivim_fit(data, bvalues)

tests/IVIMmodels/unit_tests/algorithms.json

@@ -6,7 +6,10 @@
         "IAR_LU_modified_topopro",
         "IAR_LU_segmented_2step",
         "IAR_LU_segmented_3step",
-        "IAR_LU_subtracted"
+        "IAR_LU_subtracted",
+        "OGC_AmsterdamUMC_Bayesian_biexp",
+        "OGC_AmsterdamUMC_biexp_segmented",
+        "OGC_AmsterdamUMC_biexp"
     ],
     "IAR_LU_biexp": {
         "tolerances": {

tests/IVIMmodels/unit_tests/analyze.r

@@ -1,13 +1,18 @@
+library(plyr)  # order matters, we need dplyr to come later for grouping
+library(dplyr)
 library(tidyverse)
-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")
+    f_plot <- 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")
+    print(f_plot)
+    ggsave(paste("f", fileExtension, sep=""), plot=f_plot, width = 50, height = 50, units = "cm")
+    D_plot <- 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")
+    print(D_plot)
+    ggsave(paste("D", fileExtension, sep=""), plot=D_plot, width = 50, height = 50, units = "cm")
+    Dp_plot <- 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")
+    print(Dp_plot)
+    ggsave(paste("Dp", fileExtension, sep=""), plot=Dp_plot, width = 50, height = 50, units = "cm")
 }
 
 data <- read.csv("test_output.csv")
@@ -23,3 +28,17 @@ 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")
 
+
+# some kind of black voodoo magic to number unique groups
+# data_new <- data %>% group_by(Algorithm, Region, SNR) %>% mutate(index = row_number()) %>% ungroup()
+# Then this widens it so we can lm()
+data_wide <- data %>% pivot_wider(names_from=Algorithm, values_from=c(f_fitted, Dp_fitted, D_fitted), id_cols=c(Region, SNR, index, f, D, Dp))
+# linear fit for f
+f_model <- lm(f ~ f_fitted_IAR_LU_biexp + f_fitted_IAR_LU_modified_mix + f_fitted_IAR_LU_modified_topopro + f_fitted_IAR_LU_segmented_2step + f_fitted_IAR_LU_segmented_3step + f_fitted_IAR_LU_subtracted, data=data_wide)
+#f_model <- lm(f ~ f_fitted_ETP_SRI_LinearFitting + f_fitted_IAR_LU_biexp + f_fitted_IAR_LU_modified_mix + f_fitted_IAR_LU_modified_topopro + f_fitted_IAR_LU_segmented_2step + f_fitted_IAR_LU_segmented_3step + f_fitted_IAR_LU_subtracted, data=data_new_wide)
+D_model <- lm(D ~ D_fitted_ETP_SRI_LinearFitting + D_fitted_IAR_LU_biexp + D_fitted_IAR_LU_modified_topopro + D_fitted_IAR_LU_segmented_2step + D_fitted_IAR_LU_segmented_3step + D_fitted_IAR_LU_subtracted, data=data_wide)
+# D_model <- lm(D ~ D_fitted_ETP_SRI_LinearFitting + D_fitted_IAR_LU_biexp + D_fitted_IAR_LU_modified_mix + D_fitted_IAR_LU_modified_topopro + D_fitted_IAR_LU_segmented_2step + D_fitted_IAR_LU_segmented_3step + D_fitted_IAR_LU_subtracted, data=data_new_wide)
+Dp_model <- lm(Dp ~ Dp_fitted_IAR_LU_biexp + Dp_fitted_IAR_LU_modified_mix + Dp_fitted_IAR_LU_modified_topopro + Dp_fitted_IAR_LU_segmented_2step + Dp_fitted_IAR_LU_segmented_3step + Dp_fitted_IAR_LU_subtracted, data=data_wide)
+# Dp_model <- lm(Dp ~ Dp_fitted_ETP_SRI_LinearFitting + Dp_fitted_IAR_LU_biexp + Dp_fitted_IAR_LU_modified_mix + Dp_fitted_IAR_LU_modified_topopro + Dp_fitted_IAR_LU_segmented_2step + Dp_fitted_IAR_LU_segmented_3step + Dp_fitted_IAR_LU_subtracted, data=data_new_wide)
+# predict new data from existing model
+predict(object = f_model, newdata = data_wide)

tests/IVIMmodels/unit_tests/test_ivim_synthetic.py

@@ -24,6 +24,10 @@ 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)
+    # here try a prior, but it's not seeing it, why?
+    # if hasattr(fit, "accepts_priors") and fit.accepts_priors:
+    #     prior = [np.repeat(D, 10), np.repeat(f, 10), np.repeat(Dp, 10)]  # D, f, D*
+    #     fit.initialize(prior_in=prior)
     time_delta = datetime.timedelta()
     for idx in range(fit_count):
         # if "data" not in data:
@@ -34,16 +38,16 @@ def test_generated(ivim_algorithm, ivim_data, SNR, rtol, atol, fit_count, save_f
         [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, Dp, D], [f_fit, Dp_fit, D_fit])
+            save_results(save_file, ivim_algorithm, name, SNR, idx, [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):
+def save_results(filename, algorithm, name, SNR, idx, 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, idx, *truth, *fit]
         writer.writerow(data)
 
 def save_duration(filename, algorithm, name, SNR, duration, count):