Merge pull request #109 from OSIPI/DL_wrapper

Dl wrapper + DL testing

Author: oliverchampion

.gitignore

@@ -23,6 +23,7 @@ tests/IVIMmodels/unit_tests/*.log
 junit/*
 ivim_simulation.bval
 ivim_simulation.bvec
+*.pt
 
 # Unit test / coverage reports
 .tox/
@@ -32,4 +33,5 @@ nosetests.xml
 coverage.xml
 *.pyc
 phantoms/MR_XCAT_qMRI/*.json
-phantoms/MR_XCAT_qMRI/*.txt
+phantoms/MR_XCAT_qMRI/*.txt
+tests/IVIMmodels/unit_tests/models

conftest.py

@@ -37,6 +37,12 @@ def pytest_addoption(parser):
         type=str,
         help="Default data file name",
     )
+    parser.addoption(
+        "--dataFileDL",
+        default="tests/IVIMmodels/unit_tests/generic_DL.json",
+        type=str,
+        help="Default data file name",
+    )
     parser.addoption(
         "--saveFileName",
         default="",
@@ -179,6 +185,10 @@ def pytest_generate_tests(metafunc):
     if "bound_input" in metafunc.fixturenames:
         args = bound_input(metafunc.config.getoption("dataFile"),metafunc.config.getoption("algorithmFile"))
         metafunc.parametrize("bound_input", args)
+    if "deep_learning_algorithms" in metafunc.fixturenames:
+        args = deep_learning_algorithms(metafunc.config.getoption("dataFileDL"),metafunc.config.getoption("algorithmFile"))
+        metafunc.parametrize("deep_learning_algorithms", args)
+
 
 
 def data_list(filename):
@@ -210,17 +220,18 @@ def data_ivim_fit_saved(datafile, algorithmFile):
         first = True
         for name, data in all_data.items():
             algorithm_dict = algorithm_information.get(algorithm, {})
-            xfail = {"xfail": name in algorithm_dict.get("xfail_names", {}),
-                "strict": algorithm_dict.get("xfail_names", {}).get(name, True)}
-            kwargs = algorithm_dict.get("options", {})
-            tolerances = algorithm_dict.get("tolerances", {})
-            skiptime=False
-            if first:
-                if algorithm_dict.get("fail_first_time", False):
-                    skiptime = True
-                    first = False
-            requires_matlab = algorithm_dict.get("requires_matlab", False)
-            yield name, bvals, data, algorithm, xfail, kwargs, tolerances, skiptime, requires_matlab
+            if not algorithm_dict.get('deep_learning',False):
+                xfail = {"xfail": name in algorithm_dict.get("xfail_names", {}),
+                    "strict": algorithm_dict.get("xfail_names", {}).get(name, True)}
+                kwargs = algorithm_dict.get("options", {})
+                tolerances = algorithm_dict.get("tolerances", {})
+                skiptime=False
+                if first:
+                    if algorithm_dict.get("fail_first_time", False):
+                        skiptime = True
+                        first = False
+                requires_matlab = algorithm_dict.get("requires_matlab", False)
+                yield name, bvals, data, algorithm, xfail, kwargs, tolerances, skiptime, requires_matlab
 
 def algorithmlist(algorithmFile):
     # Find the algorithms from algorithms.json
@@ -233,7 +244,7 @@ def algorithmlist(algorithmFile):
     for algorithm in algorithms:
         algorithm_dict = algorithm_information.get(algorithm, {})
         requires_matlab = algorithm_dict.get("requires_matlab", False)
-        yield algorithm, requires_matlab
+        yield algorithm, requires_matlab, algorithm_dict.get('deep_learning', False)
 
 def bound_input(datafile,algorithmFile):
     # Find the algorithms from algorithms.json
@@ -251,9 +262,31 @@ def bound_input(datafile,algorithmFile):
     for name, data in all_data.items():
         for algorithm in algorithms:
             algorithm_dict = algorithm_information.get(algorithm, {})
-            xfail = {"xfail": name in algorithm_dict.get("xfail_names", {}),
-                "strict": algorithm_dict.get("xfail_names", {}).get(name, True)}
+            if not algorithm_dict.get('deep_learning',False):
+                xfail = {"xfail": name in algorithm_dict.get("xfail_names", {}),
+                    "strict": algorithm_dict.get("xfail_names", {}).get(name, True)}
+                kwargs = algorithm_dict.get("options", {})
+                tolerances = algorithm_dict.get("tolerances", {})
+                requires_matlab = algorithm_dict.get("requires_matlab", False)
+                yield name, bvals, data, algorithm, xfail, kwargs, tolerances, requires_matlab
+
+def deep_learning_algorithms(datafile,algorithmFile):
+    # Find the algorithms from algorithms.json
+    current_folder = pathlib.Path.cwd()
+    algorithm_path = current_folder / algorithmFile
+    with algorithm_path.open() as f:
+        algorithm_information = json.load(f)
+    # Load generic test data generated from the included phantom: phantoms/MR_XCAT_qMRI
+    generic = current_folder / datafile
+    with generic.open() as f:
+        all_data = json.load(f)
+    algorithms = algorithm_information["algorithms"]
+    bvals = all_data.pop('config')
+    bvals = bvals['bvalues']
+    for algorithm in algorithms:
+        algorithm_dict = algorithm_information.get(algorithm, {})
+        if algorithm_dict.get('deep_learning',False):
             kwargs = algorithm_dict.get("options", {})
-            tolerances = algorithm_dict.get("tolerances", {})
             requires_matlab = algorithm_dict.get("requires_matlab", False)
-            yield name, bvals, data, algorithm, xfail, kwargs, tolerances, requires_matlab
+            tolerances = algorithm_dict.get("tolerances", {"atol":{"f": 2e-1, "D": 8e-4, "Dp": 8e-2},"rtol":{"f": 0.2, "D": 0.3, "Dp": 0.4}})
+            yield algorithm, all_data, bvals, kwargs, requires_matlab, tolerances

src/standardized/IVIM_NEToptim.py

@@ -0,0 +1,185 @@
+from src.wrappers.OsipiBase import OsipiBase
+import numpy as np
+import IVIMNET.deep as deep
+import torch
+import warnings
+from utilities.data_simulation.GenerateData import GenerateData
+
+class IVIM_NEToptim(OsipiBase):
+    """
+    Bi-exponential fitting algorithm by Oliver Gurney-Champion, Amsterdam UMC
+    """
+
+    # I'm thinking that we define default attributes for each submission like this
+    # And in __init__, we can call the OsipiBase control functions to check whether
+    # the user inputs fulfil the requirements
+
+    # Some basic stuff that identifies the algorithm
+    id_author = "Oliver Gurney Champion, Amsterdam UMC"
+    id_algorithm_type = "Deep learnt bi-exponential fit"
+    id_return_parameters = "f, D*, D, S0"
+    id_units = "seconds per milli metre squared or milliseconds per micro metre squared"
+
+    # Algorithm requirements
+    required_bvalues = 4
+    required_thresholds = [0,
+                           0]  # Interval from "at least" to "at most", in case submissions allow a custom number of thresholds
+    required_bounds = False
+    required_bounds_optional = True  # Bounds may not be required but are optional
+    required_initial_guess = False
+    required_initial_guess_optional = False
+    accepted_dimensions = 1  # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most?
+
+
+    # Supported inputs in the standardized class
+    supported_bounds = True
+    supported_initial_guess = False
+    supported_thresholds = False
+
+    def __init__(self, SNR=None, bvalues=None, thresholds=None, bounds=None, initial_guess=None, fitS0=True, traindata=None):
+        """
+            Everything this algorithm requires should be implemented here.
+            Number of segmentation thresholds, bounds, etc.
+
+            Our OsipiBase object could contain functions that compare the inputs with
+            the requirements.
+        """
+        if bvalues == None:
+            raise ValueError("for deep learning models, bvalues need defining at initiaition")
+        #super(OGC_AmsterdamUMC_biexp, self).__init__(bvalues, bounds, initial_guess, fitS0)
+        super(IVIM_NEToptim, self).__init__(bvalues=bvalues, bounds=bounds, initial_guess=initial_guess)
+        self.fitS0=fitS0
+        self.bvalues=np.array(bvalues)
+        self.initialize(bounds, initial_guess, fitS0, traindata, SNR)
+
+    def initialize(self, bounds, initial_guess, fitS0, traindata, SNR):
+        self.fitS0=fitS0
+        self.deep_learning = True
+        self.supervised = False
+        if traindata is None:
+            warnings.warn('no training data provided (traindata = None). Training data will be simulated')
+            if SNR is None:
+                warnings.warn('No SNR indicated. Data simulated with SNR = (5-1000)')
+                SNR = (5, 1000)
+            self.training_data(self.bvalues,n=1000000,SNR=SNR)
+        self.arg=Arg()
+        if bounds is not None:
+            self.arg.net_pars.cons_min = bounds[0]  # Dt, Fp, Ds, S0
+            self.arg.net_pars.cons_max = bounds[1]  # Dt, Fp, Ds, S0
+        if traindata is None:
+            self.net = deep.learn_IVIM(self.train_data['data'], self.bvalues, self.arg)
+        else:
+            self.net = deep.learn_IVIM(traindata, self.bvalues, self.arg)
+        self.algorithm =lambda data: deep.predict_IVIM(data, self.bvalues, self.net, self.arg)
+
+
+    def ivim_fit(self, signals, bvalues, **kwargs):
+        """Perform the IVIM fit
+
+        Args:
+            signals (array-like)
+            bvalues (array-like): b-values for the signals. If None, self.bvalues will be used. Default is None.
+
+        Returns:
+            _type_: _description_
+        """
+        if not np.array_equal(bvalues, self.bvalues):
+            raise ValueError("bvalue list at fitting must be identical as the one at initiation, otherwise it will not run")
+
+        paramsNN = deep.predict_IVIM(signals, self.bvalues, self.net, self.arg)
+
+        results = {}
+        results["D"] = paramsNN[0]
+        results["f"] = paramsNN[1]
+        results["Dp"] = paramsNN[2]
+
+        return results
+
+
+    def ivim_fit_full_volume(self, signals, bvalues, retrain_on_input_data=False, **kwargs):
+        """Perform the IVIM fit
+
+        Args:
+            signals (array-like)
+            bvalues (array-like): b-values for the signals. If None, self.bvalues will be used. Default is None.
+
+        Returns:
+            _type_: _description_
+        """
+        if not np.array_equal(bvalues, self.bvalues):
+            raise ValueError("bvalue list at fitting must be identical as the one at initiation, otherwise it will not run")
+
+        signals = self.reshape_to_voxelwise(signals)
+        if retrain_on_input_data:
+            self.net = deep.learn_IVIM(signals, self.bvalues, self.arg, net=self.net)
+        paramsNN = deep.predict_IVIM(signals, self.bvalues, self.net, self.arg)
+
+        results = {}
+        results["D"] = paramsNN[0]
+        results["f"] = paramsNN[1]
+        results["Dp"] = paramsNN[2]
+
+        return results
+
+    def reshape_to_voxelwise(self, data):
+        """
+        reshapes multi-D input (spatial dims, bvvalue) data to 2D voxel-wise array
+        Args:
+            data (array): mulit-D array (data x b-values)
+        Returns:
+            out (array): 2D array (voxel x b-value)
+        """
+        B = data.shape[-1]
+        voxels = int(np.prod(data.shape[:-1]))  # e.g., X*Y*Z
+        return data.reshape(voxels, B)
+
+
+    def training_data(self, bvalues, data=None, SNR=(5,1000), n=1000000,Drange=(0.0005,0.0034),frange=(0,1),Dprange=(0.005,0.1),rician_noise=False):
+        rng = np.random.RandomState(42)
+        if data is None:
+            gen = GenerateData(rng=rng)
+            data, D, f, Dp = gen.simulate_training_data(bvalues, SNR=SNR, n=n,Drange=Drange,frange=frange,Dprange=Dprange,rician_noise=rician_noise)
+            if self.supervised:
+                self.train_data = {'data':data,'D':D,'f':f,'Dp':Dp}
+            else:
+                self.train_data = {'data': data}
+
+class NetArgs:
+    def __init__(self):
+        self.optim = 'adam'  # these are the optimisers implementd. Choices are: 'sgd'; 'sgdr'; 'adagrad' adam
+        self.lr = 0.00003  # this is the learning rate.
+        self.patience = 10  # this is the number of epochs without improvement that the network waits untill determining it found its optimum
+        self.batch_size = 128  # number of datasets taken along per iteration
+        self.maxit = 500  # max iterations per epoch
+        self.split = 0.9  # split of test and validation data
+        self.load_nn = False  # load the neural network instead of retraining
+        self.loss_fun = 'rms'  # what is the loss used for the model. rms is root mean square (linear regression-like); L1 is L1 normalisation (less focus on outliers)
+        self.skip_net = False  # skip the network training and evaluation
+        self.scheduler = False  # as discussed in the article, LR is important. This approach allows to reduce the LR itteratively when there is no improvement throughout an 5 consecutive epochs
+        # use GPU if available
+        self.use_cuda = torch.cuda.is_available()
+        self.device = torch.device("cuda:0" if self.use_cuda else "cpu")
+        self.select_best = False
+        # the optimized network settings
+
+class NetPars:
+    def __init__(self):
+        self.dropout = 0.1  # 0.0/0.1 chose how much dropout one likes. 0=no dropout; internet says roughly 20% (0.20) is good, although it also states that smaller networks might desire smaller amount of dropout
+        self.batch_norm = True  # False/True turns on batch normalistion
+        self.parallel = 'parallel'  # defines whether the network exstimates each parameter seperately (each parameter has its own network) or whether 1 shared network is used instead
+        self.con = 'sigmoid'  # defines the constraint function; 'sigmoid' gives a sigmoid function giving the max/min; 'abs' gives the absolute of the output, 'none' does not constrain the output
+        self.tri_exp = False
+        #### only if sigmoid constraint is used!
+        self.cons_min = [0, 0, 0.005, 0]  # Dt, Fp, Ds, S0
+        self.cons_max = [0.005, 0.8, 0.2, 2.0]  # Dt, Fp, Ds, S0
+        ####
+        self.fitS0 = True  # indicates whether to fit S0 (True) or fix it to 1 (for normalised signals); I prefer fitting S0 as it takes along the potential error is S0.
+        self.depth = 2  # number of layers
+        self.width = 0  # new option that determines network width. Putting to 0 makes it as wide as the number of b-values
+        boundsrange = 0.3 * (np.array(self.cons_max)-np.array(self.cons_min)) # ensure that we are on the most lineair bit of the sigmoid function
+        self.cons_min = np.array(self.cons_min) - boundsrange
+        self.cons_max = np.array(self.cons_max) + boundsrange
+class Arg:
+    def __init__(self):
+        self.train_pars = NetArgs()
+        self.net_pars = NetPars()