|
| 1 | +""" |
| 2 | +September 2020 by Oliver Gurney-Champion |
| 3 | +oliver.gurney.champion@gmail.com / o.j.gurney-champion@amsterdamumc.nl |
| 4 | +https://www.github.com/ochampion |
| 5 | +
|
| 6 | +Code is uploaded as part of our publication in MRM (Kaandorp et al. Improved unsupervised physics-informed deep learning for intravoxel-incoherent motion modeling and evaluation in pancreatic cancer patients. MRM 2021) |
| 7 | +
|
| 8 | +requirements: |
| 9 | +numpy |
| 10 | +torch |
| 11 | +tqdm |
| 12 | +matplotlib |
| 13 | +scipy |
| 14 | +joblib |
| 15 | +""" |
| 16 | +import torch |
| 17 | +import numpy as np |
| 18 | + |
| 19 | + |
| 20 | +#most of these are options from the article and explained in the M&M. |
| 21 | +class train_pars: |
| 22 | + def __init__(self,nets): |
| 23 | + self.optim='adam' #these are the optimisers implementd. Choices are: 'sgd'; 'sgdr'; 'adagrad' adam |
| 24 | + if nets == 'optim': |
| 25 | + self.lr = 0.00003 # this is the learning rate. |
| 26 | + elif nets == 'orig': |
| 27 | + self.lr = 0.001 # this is the learning rate. |
| 28 | + else: |
| 29 | + self.lr = 0.00003 # this is the learning rate. |
| 30 | + self.patience= 10 # this is the number of epochs without improvement that the network waits untill determining it found its optimum |
| 31 | + self.batch_size= 128 # number of datasets taken along per iteration |
| 32 | + self.maxit = 500 # max iterations per epoch |
| 33 | + self.split = 0.9 # split of test and validation data |
| 34 | + self.load_nn= False # load the neural network instead of retraining |
| 35 | + 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) |
| 36 | + self.skip_net = False # skip the network training and evaluation |
| 37 | + 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 |
| 38 | + # use GPU if available |
| 39 | + self.use_cuda = torch.cuda.is_available() |
| 40 | + self.device = torch.device("cuda:0" if self.use_cuda else "cpu") |
| 41 | + self.select_best = False |
| 42 | + |
| 43 | + |
| 44 | +class net_pars: |
| 45 | + def __init__(self,nets): |
| 46 | + # select a network setting |
| 47 | + if (nets == 'optim'): |
| 48 | + # the optimized network settings |
| 49 | + 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 |
| 50 | + self.batch_norm = True # False/True turns on batch normalistion |
| 51 | + 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 |
| 52 | + 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 |
| 53 | + self.tri_exp = False |
| 54 | + #### only if sigmoid constraint is used! |
| 55 | + if self.tri_exp: |
| 56 | + self.cons_min = [0., 0.0003, 0.0, 0.003, 0.0, 0.08] # F0', D0, F1', D1, F2', D2 |
| 57 | + self.cons_max = [2.5, 0.003, 1, 0.08, 1, 5] # F0', D0, F1', D1, F2', D2 |
| 58 | + else: |
| 59 | + self.cons_min = [0, 0, 0.005, 0] # Dt, Fp, Ds, S0 |
| 60 | + self.cons_max = [0.005, 0.7, 0.2, 2.0] # Dt, Fp, Ds, S0 |
| 61 | + #### |
| 62 | + 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. |
| 63 | + self.depth = 2 # number of layers |
| 64 | + self.width = 0 # new option that determines network width. Putting to 0 makes it as wide as the number of b-values |
| 65 | + elif nets == 'orig': |
| 66 | + # as summarized in Table 1 from the main article for the original network |
| 67 | + self.dropout = 0.0 #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 |
| 68 | + self.batch_norm = False # False/True turns on batch normalistion |
| 69 | + self.parallel = 'single' # defines whether the network exstimates each parameter seperately (each parameter has its own network) or whether 1 shared network is used instead |
| 70 | + self.con = 'abs' # 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 |
| 71 | + self.tri_exp = False |
| 72 | + #### only if sigmoid constraint is used! |
| 73 | + if self.tri_exp: |
| 74 | + self.cons_min = [0., 0.0003, 0.0, 0.003, 0.0, 0.08] # F0', D0, F1', D1, F2', D2 |
| 75 | + self.cons_max = [2.5, 0.003, 1, 0.08, 1, 5] # F0', D0, F1', D1, F2', D2 |
| 76 | + else: |
| 77 | + self.cons_min = [0, 0, 0.005, 0] # Dt, Fp, Ds, f0 |
| 78 | + self.cons_max = [0.005, 0.7, 0.2, 2.0] # Dt, Fp, Ds, f0 |
| 79 | + #### |
| 80 | + self.fitS0 = False # 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. |
| 81 | + self.depth = 3 # number of layers |
| 82 | + self.width = 0 # new option that determines network width. Putting to 0 makes it as wide as the number of b-values |
| 83 | + else: |
| 84 | + # chose wisely :) |
| 85 | + self.dropout = 0.3 #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 |
| 86 | + self.batch_norm = True # False/True turns on batch normalistion |
| 87 | + 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 |
| 88 | + 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 |
| 89 | + self.tri_exp = True |
| 90 | + #### only if sigmoid constraint is used! |
| 91 | + if self.tri_exp: |
| 92 | + self.cons_min = [0., 0.0003, 0.0, 0.003, 0.0, 0.08] # F0', D0, F1', D1, F2', D2 |
| 93 | + self.cons_max = [2.5, 0.003, 1, 0.08, 1, 5] # F0', D0, F1', D1, F2', D2 |
| 94 | + else: |
| 95 | + self.cons_min = [0, 0, 0.005, 0] # Dt, Fp, Ds, f0 |
| 96 | + self.cons_max = [0.005, 0.7, 0.2, 2.0] # Dt, Fp, Ds, f0 |
| 97 | + #### |
| 98 | + self.fitS0 = False # 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. |
| 99 | + self.depth = 4 # number of layers |
| 100 | + self.width = 500 # new option that determines network width. Putting to 0 makes it as wide as the number of b-values |
| 101 | + 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 |
| 102 | + self.cons_min = np.array(self.cons_min) - boundsrange |
| 103 | + self.cons_max = np.array(self.cons_max) + boundsrange |
| 104 | + |
| 105 | + |
| 106 | + |
| 107 | +class lsqfit: |
| 108 | + def __init__(self): |
| 109 | + self.method = 'lsq' #seg, bayes or lsq |
| 110 | + self.model = 'bi-exp' #bi-exp or tri-exp |
| 111 | + self.do_fit = False # skip lsq fitting |
| 112 | + self.load_lsq = False # load the last results for lsq fit |
| 113 | + self.fitS0 = True # indicates whether to fit S0 (True) or fix it to 1 in the least squares fit. |
| 114 | + self.jobs = 4 # number of parallel jobs. If set to 1, no parallel computing is used |
| 115 | + self.bvalues_included=[50, 700] # for ADC fit |
| 116 | + if self.model == 'tri-exp': |
| 117 | + self.bounds = ([0, 0, 0, 0.008, 0, 0.06], [2.5, 0.008, 1, 0.08, 1, 5]) # F0', D0, F1', D1, F2', D2 |
| 118 | + elif self.model == 'ADC': |
| 119 | + self.bounds = ([0, 0], [0.005, 2.5]) # Dt, Fp, Ds, S0 |
| 120 | + else: |
| 121 | + self.bounds = ([0.0003, 0, 0.005, 0.5], [0.005, 0.7, 0.3, 2.5]) # Dt, Fp, Ds, S0 |
| 122 | + |
| 123 | +class sim: |
| 124 | + def __init__(self): |
| 125 | + self.bvalues = np.array([0, 5, 10, 20, 30, 40, 60, 150, 300, 500, 700]) # array of b-values |
| 126 | + self.SNR = [15, 20, 30, 50] # the SNRs to simulate at |
| 127 | + self.sims = 1000000 # number of simulations to run |
| 128 | + self.num_samples_eval = 10000 # number of simualtiosn te evaluate. This can be lower than the number run. Particularly to save time when fitting. More simulations help with generating sufficient data for the neural network |
| 129 | + self.repeats = 1 # this is the number of repeats for simulations |
| 130 | + self.rician = False # add rician noise to simulations; if false, gaussian noise is added instead |
| 131 | + self.model = 'bi-exp' |
| 132 | + if self.model == 'bi-exp': |
| 133 | + self.range = ([0.0005, 0.05, 0.01], [0.003, 0.55, 0.1]) |
| 134 | + else: |
| 135 | + self.range = ([0.0005, 0.05, 0.001, 0.05, 0.08], [0.003, 0.5, 0.05, 0.5, 2]) # D0, F1', D1, F2', D2 |
| 136 | + |
| 137 | + |
| 138 | +class hyperparams: |
| 139 | + def __init__(self): |
| 140 | + self.fig = False # plot results and intermediate steps |
| 141 | + self.save_name = 'optim' # orig or optim (or optim_adsig for in vivo) |
| 142 | + self.net_pars = net_pars(self.save_name) |
| 143 | + self.train_pars = train_pars(self.save_name) |
| 144 | + self.fit = lsqfit() |
| 145 | + self.sim = sim() |
0 commit comments