Merge pull request #26 from OSIPI/analysis/improvements

Improved plotting and random number generation

Author: oscarjalnefjord

.github/workflows/analysis.yml

@@ -61,3 +61,5 @@ jobs:
               f_limited.pdf
               Dp_limited.pdf
               durations.pdf
+              curve_plot.pdf
+              fitted_curves.pdf

conftest.py

@@ -78,11 +78,20 @@ def save_file(request):
         # print(filename)
         # filename.unlink(missing_ok=True)
         filename = filename.as_posix()
+
+        data = data_list(request.config.getoption("--dataFile"))  # TODO: clean up this hacky way to get bvalues
+        [_, bvalues, _] = next(data)
+        bvalue_string = ["bval_" + str(bvalue) for bvalue in bvalues]
+        # bvalue_string = ["b_0.0","b_1.0","b_2.0","b_5.0","b_10.0","b_20.0","b_30.0","b_50.0","b_75.0","b_100.0","b_150.0","b_250.0","b_350.0","b_400.0","b_550.0","b_700.0","b_850.0","b_1000.0"]
+
         with open(filename, "w") as csv_file:
             writer = csv.writer(csv_file, delimiter=',')
-            writer.writerow(("Algorithm", "Region", "SNR", "index", "f", "Dp", "D", "f_fitted", "Dp_fitted", "D_fitted"))
+            writer.writerow(("Algorithm", "Region", "SNR", "index", "f", "Dp", "D", "f_fitted", "Dp_fitted", "D_fitted", *bvalue_string))
+            yield writer
             # writer.writerow(["", datetime.datetime.now()])
-    return filename
+    else:
+        yield None
+    # return filename
 
 @pytest.fixture(scope="session")
 def save_duration_file(request):
@@ -96,8 +105,11 @@ def save_duration_file(request):
         with open(filename, "w") as csv_file:
             writer = csv.writer(csv_file, delimiter=',')
             writer.writerow(("Algorithm", "Region", "SNR", "Duration [us]", "Count"))
+            yield writer
             # writer.writerow(["", datetime.datetime.now()])
-    return filename
+    else:
+        yield None
+    # return filename
 
 @pytest.fixture(scope="session")
 def rtol(request):

tests/IVIMmodels/unit_tests/analyze.r

@@ -65,6 +65,7 @@ if (runPrediction) {
 ivim_decay <- function(f, D, Dp, bvalues) {
     return(f*exp(-Dp*bvalues) + (1-f)*exp(-D*bvalues))
 }
+#TODO: read bvalues from file
 bvalues <- c(0.0, 1.0, 2.0, 5.0, 10.0, 20.0, 30.0, 50.0, 75.0, 100.0, 150.0, 250.0, 350.0, 400.0, 550.0, 700.0, 850.0, 1000.0)
 
 generate_curves <- function(data, bvalues, appended_string) {
@@ -82,3 +83,13 @@ data_curves_restricted <- cbind(curves_restricted, signal_fitted=curves_restrict
 curve_plot <- ggplot(data_curves_restricted, aes(x=bvalue))  + facet_grid(Region ~ SNR) + geom_line(alpha=0.2, aes(y=signal_fitted, group=interaction(Algorithm, index), color=Algorithm)) + geom_line(aes(y=signal)) + ylim(0, 1) + ylab("Signal (a.u.)")
 print(curve_plot)
 ggsave("curve_plot.pdf", plot=curve_plot, width = 30, height = 30, units = "cm")
+
+
+data_points_restricted <- data_restricted[data_restricted$index==0,]
+data_points_restricted <- pivot_longer(data_points_restricted, starts_with("bval_"), names_to="bvalue", values_to="fitted_data", names_prefix="bval_", names_transform=list(bvalue = as.numeric))
+# data_curves_restricted_idx0 <- data_curves_restricted[data_curves_restricted$index==0,]
+data_points_restricted <- cbind(data_curves_restricted[data_curves_restricted$index==0,], fitted_data=data_points_restricted$fitted_data)
+# fitted_curves <- ggplot(data_points_restricted, aes(x=bvalue))  + facet_grid(Region ~ SNR) + geom_line(alpha=0.5, aes(y=signal_fitted, group=Algorithm, color=Algorithm)) + geom_point(aes(y=fitted_data, group=Algorithm, color=Algorithm)) + ylim(0, 1) + ylab("Signal (a.u.)")
+fitted_curves <- ggplot(data_points_restricted, aes(x=bvalue))  + facet_grid(Region ~ SNR) + geom_line(alpha=0.5, aes(y=signal_fitted, group=Algorithm, color=Algorithm)) + geom_point(aes(y=fitted_data, shape="Data")) + ylim(0, 1) + ylab("Signal (a.u.)") + guides(shape = guide_legend("Fitted Data"))
+print(fitted_curves)
+ggsave("fitted_curves.pdf", plot=fitted_curves, width = 30, height = 30, units = "cm")

tests/IVIMmodels/unit_tests/test_ivim_synthetic.py

@@ -16,17 +16,18 @@
 @pytest.mark.slow
 def test_generated(ivim_algorithm, ivim_data, SNR, rtol, atol, fit_count, rician_noise, save_file, save_duration_file, use_prior):
     # assert save_file == "test"
-    random.seed(42)
+    rng = np.random.RandomState(42)
+    # random.seed(42)
     S0 = 1
-    gd = GenerateData()
+    gd = GenerateData(rng=rng)
     name, bvals, data = ivim_data
     D = data["D"]
     f = data["f"]
     Dp = data["Dp"]
     fit = OsipiBase(algorithm=ivim_algorithm)
     # here is a prior
     if use_prior and hasattr(fit, "accepts_priors") and fit.accepts_priors:
-        prior = [np.random.normal(D, D/3, 10), np.random.normal(f, f/3, 10), np.random.normal(Dp, Dp/3, 10), np.random.normal(1, 1/3, 10)]
+        prior = [rng.normal(D, D/3, 10), rng.normal(f, f/3, 10), rng.normal(Dp, Dp/3, 10), rng.normal(1, 1/3, 10)]
         # prior = [np.repeat(D, 10)+np.random.normal(0,D/3,np.shape(np.repeat(D, 10))), np.repeat(f, 10)+np.random.normal(0,f/3,np.shape(np.repeat(D, 10))), np.repeat(Dp, 10)+np.random.normal(0,Dp/3,np.shape(np.repeat(D, 10))),np.repeat(np.ones_like(Dp), 10)+np.random.normal(0,1/3,np.shape(np.repeat(D, 10)))]  # D, f, D*
         fit.initialize(prior_in=prior)
     time_delta = datetime.timedelta()
@@ -38,11 +39,13 @@ def test_generated(ivim_algorithm, ivim_data, SNR, rtol, atol, fit_count, rician
         start_time = datetime.datetime.now()
         [f_fit, Dp_fit, D_fit] = fit.osipi_fit(signal, bvals) #, prior_in=prior
         time_delta += datetime.datetime.now() - start_time
-        if save_file:
-            save_results(save_file, ivim_algorithm, name, SNR, idx, [f, Dp, D], [f_fit, Dp_fit, D_fit])
+        if save_file is not None:
+            save_file.writerow([ivim_algorithm, name, SNR, idx, f, Dp, D, f_fit, Dp_fit, D_fit, *signal])
+            # 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)
+    if save_duration_file is not None:
+        save_duration_file.writerow([ivim_algorithm, name, SNR, time_delta/datetime.timedelta(microseconds=1), fit_count])
+        # save_duration(save_duration_file, ivim_algorithm, name, SNR, time_delta, fit_count)
 
 
 def save_results(filename, algorithm, name, SNR, idx, truth, fit):

utilities/data_simulation/GenerateData.py

@@ -4,17 +4,23 @@ class GenerateData:
     """
     Generate exponential and linear data
     """
-    def __init__(self, operator=None):
+    def __init__(self, operator=None, rng=None):
         """
         Parameters
         ----------
         operator : numpy or torch
             Must provide mathematical operators
+        rng : random number generator
+            Must provide normal() operator
         """
         if operator is None:
             self._op = np
         else:
             self._op = operator
+        if rng is None:
+            self._rng = self._op.random
+        else:
+            self._rng = rng
 
     def ivim_signal(self, D, Dp, f, S0, bvalues, snr=None, rician_noise=False):
         """
@@ -91,9 +97,9 @@ def add_noise(self, real_signal, snr=None, rician_noise=True, imag_signal=None):
         real_noise = self._op.zeros_like(real_signal)
         imag_noise = self._op.zeros_like(real_signal)
         if snr is not None:
-            real_noise = self._op.random.normal(0, 1 / snr, real_signal.shape)
+            real_noise = self._rng.normal(0, 1 / snr, real_signal.shape)
             if rician_noise:
-                noisy_data = self._op.sqrt(self._op.power(real_signal, 2) + self._op.power(imag_signal, 2))
+                imag_noise = self._rng.normal(0, 1 / snr, imag_signal.shape)
         noisy_data = self._op.sqrt(self._op.power(real_signal + real_noise, 2) + self._op.power(imag_signal + imag_noise, 2))
         return noisy_data