diff --git a/.github/workflows/analysis.yml b/.github/workflows/analysis.yml index ba1d26f..d429c31 100644 --- a/.github/workflows/analysis.yml +++ b/.github/workflows/analysis.yml @@ -61,3 +61,5 @@ jobs: f_limited.pdf Dp_limited.pdf durations.pdf + curve_plot.pdf + fitted_curves.pdf diff --git a/conftest.py b/conftest.py index fd6e0ea..a29d409 100644 --- a/conftest.py +++ b/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): diff --git a/tests/IVIMmodels/unit_tests/analyze.r b/tests/IVIMmodels/unit_tests/analyze.r index 27025a6..d656f9b 100644 --- a/tests/IVIMmodels/unit_tests/analyze.r +++ b/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") diff --git a/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py b/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py index cf4b09a..42017ff 100644 --- a/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py +++ b/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py @@ -16,9 +16,10 @@ @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"] @@ -26,7 +27,7 @@ def test_generated(ivim_algorithm, ivim_data, SNR, rtol, atol, fit_count, rician 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): diff --git a/utilities/data_simulation/GenerateData.py b/utilities/data_simulation/GenerateData.py index 0f91ce8..5d46e57 100644 --- a/utilities/data_simulation/GenerateData.py +++ b/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 + 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