diff --git a/.gitignore b/.gitignore index 2035425..f537e58 100644 --- a/.gitignore +++ b/.gitignore @@ -63,6 +63,7 @@ cover/ # OS .DS_Store +*.pdf # Django stuff: *.log diff --git a/MANUAL.md b/MANUAL.md index 1d4161e..372765f 100644 --- a/MANUAL.md +++ b/MANUAL.md @@ -2,7 +2,9 @@ ## Working tree -One should perform calculations in a working folder e.g. `adaptfx/work`. There the instruction files can be specified by the user and `adaptfx` will automatically produce log files if the user wishes to. +The package ships a CLI with which calculations can be performed which will be explained here. One should perform calculations in a working folder e.g. `adaptfx/work`. There the instruction files can be specified by the user and `adaptfx` will automatically produce log files if the user wishes to. + +The package also also provides a class that can be used in scripting mode. When using the package in scripts, a class creates an object with according attributes. To understand the attributes and behaviour of optimisation, read the doc-string of the `aft.py` script [here](src/adaptfx/aft.py) ## Format of the instruction file The user specifies following elements of the dictionary for the main entries: @@ -94,7 +96,7 @@ min_dose : float max_dose : float maximal physical doses to be delivered in one fraction. The doses are aimed at PTV 95. If -1 the dose is adapted to the - remaining dose tumor dose to be delivered or the remaining OAR dose + remaining tumor BED to be delivered or the remaining OAR dose allowed to be prescribed. default: -1 ``` @@ -136,9 +138,11 @@ sf_high : float sf_stepsize: float stepsize of the sparing factor stepsize. sf_prob_threshold': float - probability threshold of the sparing factor occuring. + probability threshold of the sparing factor occuring + which defines the range of sparing factors. inf_penalty : float - infinite penalty for certain undesired states. + define infinite penalty for undesired states + choose arbitrarily large compared to highest occuring reward. plot_policy : int starting from which fraction policy should be plotted. plot_values : int @@ -146,14 +150,74 @@ plot_values : int plot_remains : int starting from which fraction expected remaining number of fractions should be plotted. +plot_probability : int + flag if the probability distribution should be plotted ``` -# Example +> :warning: Note:\ +> It is dangerous to use an upper and lower bound in `sf_low` and `sf_high`, as a truncated probability distribution may result that not accurately represents the environment model. Best is to set `sf_prob_threshold` to `1e-3` or lower or leave it at the default which is set at `1e-4`. + +## Note on Plots +Policy, Value and Remaining number of fractions plots are calculated with the probability distribution in the fraction from which the plots should start. That is the value function that is known when iterating backwards through the fractions. E.g. the plotted policy starting to plot in the first fraction i.e `plot_policy = 1` and `prob_update = 0` is the policy which is known throughout the treatment, when observing only the first sparing factor. In the case of probability updating e.g `prob_update = 1` the plotted policy is the optimal policy for the probability distribution known when observing the first sparing factor. As the probability distribution changes also future optimal policies change and one has to keep in mind only policy with the constant probability distribution from fraction `1` is plotted. + +Different is the probability plot: the probability is set for each fraction and updated with additional oberved sparing factors. + +# AFT CLI -Outlined is an example instruction file for fraction minimisation. It simply is a `.json` that is translated into a python dictionary. An example can be found [here](work/oar_example.json) +Outlined is an example instruction file for fraction minimisation. It simply is a `.json` that is translated into a python dictionary. An example can be found [here](work/example_0.json) This `.json` file can be called in with the CLI as: ``` $ aft -f work/oar_example.json +``` + +# AST CLI +There is also a second CLI that allows to plot sparing factors, policy functions, temporal Adaptive Fractionation Therapy etc. +``` +$ ast [options] -f +``` + +The entry for algorithm simulation type is + +``` +algorithm_simulation: histogram, fraction, single_state, all_state +single_distance, single_patient, grid_distance, grid_fraction + allowed plots options +keys_simulation: dict + simulation keys +``` + +``` +keys_simulation +---------------- +# Histogram of applied AFT for sampled patients +n_patients : float + stepsize of the actionspace. +fixed_mean_sample : float + mean of sampled patient sparing factor +fixed_std_sample : float + standard deviation of sampled patient sparing factor + +# Data related plots +c_list : list + list of c parameters +plot_index : int + which index to plot +data_filepath : string + path to sparing factor file +data_selection : list + two elements of header in sparing factor file +data_row_hue : string + seaborn hue or row + +# Settings of plots +figsize : list + two elements of size figure +fontsize : float + fontsize of plots +save : bool + boolean to instruct saving plot +usetex : bool + boolean if TeX font should be used ``` \ No newline at end of file diff --git a/README.md b/README.md index b0d5f76..46c8940 100644 --- a/README.md +++ b/README.md @@ -124,6 +124,10 @@ A last addition is made with graphical user interfaces that facilitate the use o > :warning: Note:\ > The interfaces are not optimized, and thus it is not recommended using them to further develop extensions. +### Reducing Number of Fractions + +For the 2D algorithms there exist the possibility to reduce number of fractions. A constant $c$ can be chosen which introduces a reward (or rather a cost) linear to the number of fractions used to finish the treatment. The cost is added to the immediate reward returned by the environment in the current fraction. There exist a simulative model helping to estimate what the constant $c$ should be chosen in order for the treatment to finish on some target number of fractions $n_{\text{targ}}$. The function can be found [here](src/adaptfx/radiobiology.py) in `c_calc`. + ### Probability Updating The DP algorithm relies on a description of the environment to compute an optimal policy, in this case the probability distribution of the sparing factor $P(\delta)$, which we assume to be a Gaussian distribution truncated at $0$, with patient-specific parameters for mean and standard deviation. At the start of a treatment, only two sparing factors are available for that patient, from the planning scan and the first fraction. In each fraction, an additional sparing factor is measured, which can be used to calculate updated estimates $\mu_t$ and $\sigma_t$ for mean and standard deviation, respectively. @@ -173,12 +177,12 @@ ImportError: No module named '_ctypes' **Solution:** with the specific package manager of the Linux distribution install `libffi-dev` development tool. E.g. in Fedora Linux and derivatives install this tool ``` -sudo dnf install libffi-devel +$ sudo dnf install libffi-devel ``` On Ubuntu: ``` -sudo apt install libffi-dev +$ sudo apt install libffi-dev ``` ### No GUI backend for `matplotlib` @@ -195,19 +199,19 @@ No matching distribution found for tkinter **Solution:** on Fedora Linux and derivative distributions one could solve this by either installing python tkinter ``` -sudo dnf install python3-tkinter +$ sudo dnf install python3-tkinter ``` on Ubuntu ``` -sudo apt-get install python3-tk +$ sudo apt-get install python3-tk ``` **Solution:** on MacOS and Linux one could instead use `pip` to install `pyqt` ``` -pip install pyqt5 +$ pip install pyqt5 ``` diff --git a/pyproject.toml b/pyproject.toml index e93587e..b0fa9a9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,7 +26,9 @@ dependencies = [ "numpy == 1.24.1", "scipy == 1.10.0", "matplotlib == 3.6.3", - "pyqt6 == 6.4.0" + "pyqt6 == 6.4.0", + "pandas==1.5.2", + "seaborn==0.12.2" ] dynamic = ["version"] @@ -47,3 +49,4 @@ include = ["*"] [project.scripts] aft = "adaptfx.aft:main" +ast = "adaptsim.ast:main" diff --git a/src/adaptfx/__init__.py b/src/adaptfx/__init__.py index 40f15df..93a3931 100644 --- a/src/adaptfx/__init__.py +++ b/src/adaptfx/__init__.py @@ -14,7 +14,7 @@ from .radiobiology import * from .visualiser import * from .reinforce import * -from .reinforce_old import * +from .reinforce_old import max_tumor_bed_old, min_oar_bed_old, min_oar_max_tumor_old __all__ = ['bed_calc_matrix', 'convert_to_physical', diff --git a/src/adaptfx/aft.py b/src/adaptfx/aft.py index e88add7..a0f3717 100644 --- a/src/adaptfx/aft.py +++ b/src/adaptfx/aft.py @@ -7,11 +7,80 @@ class RL_object(): """ - Reinforcement Learning class to check instructions - of calculation, invoke keys and define - calculation settings from file - """ - def __init__(self, instruction_filename): + Invokes a class with one dictionary object and multiple functions + to operate on the dictionary. Dictionary is created from instruction + file and contains instructions, keys, settings. Once optimisation + is performed it also provides results. + + Parameters + ---------- + instruction_filename_in : string + string of instruction filename + + Returns + ------- + returns : ``plan`` class + The optimisation result represented as a ``plan`` object. + Important attributes are: ``keys`` and ``settings`` the keys and + settings from the instruction file. To start the optimisation one + can run ``plan.optimise`` which will store the results in ``output``. + Important attributes are: ``output.physical_doses``, + ``output.tumor_doses``, ``output.oar_doses``. + If one is looking for summed doses they are found in: + ``output.tumor_sum``, ``output.oar_sum``. + In case where number of fractions are minimised one can check + utilised number of fractions with ``output.fractions_used``. + + A full list of attributes: + ``plan.filename`` + ``plan.basename`` + ``plan.algorithm`` + ``plan.log`` + ``plan.log_level`` + ``plan.keys`` + ``plan.settings`` + ``plan.output.physical_doses`` + ``plan.output.tumor_doses`` + ``plan.output.oar_doses`` + ``plan.output.tumor_sum`` + ``plan.output.oar_sum`` + ``plan.output.fractions_used`` + + A full list of optional attributes (dependent if the user specifies plot): + ``plan.output.policy`` + ``plan.output.value`` + ``plan.output.remains`` + ``plan.output.probability`` + + A full list of available functions: + ``plan.optimise()`` + ``plan.plot()`` + ``plan.fraction_counter()`` + + + Examples + -------- + Consider the following demonstration: + + >>> import adaptfx as afx + >>> plan = afx.RL_object('path/to/instruction_file') + + >>> plan.keys.sparing_factors + [0.98, 0.97, 0.8, 0.83, 0.8, 0.85, 0.94] + + + >>> plan.optimise() + process duration: 0.0219 s: + >>> plan.output.oar_doses + array([ 0. , 3.2, 3.3, 11.6, 107.1]) + >>> plan.output.tumor_sum + 72.0 + >>> plan.output.fractions_used + 5 + +""" + def __init__(self, instruction_filename_in): + instruction_filename, basename = afx.get_abs_path(instruction_filename_in, nme) try: # check if file can be opened with open(instruction_filename, 'r') as f: read_in = f.read() @@ -48,7 +117,7 @@ def __init__(self, instruction_filename): if not log_level in afx.LOG_LEVEL_LIST: afx.aft_error('invalid "debug" flag was set', nme) - afx.logging_init(instruction_filename, log_bool, log_level) + afx.logging_init(basename, log_bool, log_level) afx.aft_message_info('log level:', log_level, nme) afx.aft_message_info('log to file:', log_bool, nme) @@ -86,6 +155,8 @@ def __init__(self, instruction_filename): settings = afx.setting_reader(afx.SETTING_DICT, user_settings) afx.aft_message_dict('settings', settings, nme, 1) + self.filename = instruction_filename + self.basename = basename self.algorithm = algorithm self.log = log_bool self.log_level = log_level @@ -93,7 +164,10 @@ def __init__(self, instruction_filename): self.settings = afx.DotDict(settings) def optimise(self): + afx.aft_message('optimising...', nme, 1) + start = afx.timing() self.output = afx.multiple(self.algorithm, self.keys, self.settings) + afx.timing(start) def fraction_counter(self): if self.algorithm == 'frac' and self.keys.fraction == 0: @@ -104,14 +178,27 @@ def fraction_counter(self): def plot(self): out = self.output sets = self.settings + figures = [] if self.settings.plot_policy: - afx.plot_val(out.policy.sf, out.policy.states, out.policy.val, out.policy.fractions) + policy_fig = afx.plot_val(out.policy.sf, out.policy.states, out.policy.val, + out.policy.fractions, 'turbo') + figures.append(policy_fig) if self.settings.plot_values: - afx.plot_val(out.value.sf, out.value.states, out.value.val, out.value.fractions) + values_fig = afx.plot_val(out.value.sf, out.value.states, + out.value.val, out.value.fractions, 'viridis') + figures.append(values_fig) if self.settings.plot_remains: - afx.plot_val(out.remains.sf, out.remains.states, out.remains.val, out.remains.fractions) - - if sets.plot_policy or sets.plot_values or sets.plot_remains: + remains_fig = afx.plot_val(out.remains.sf, out.remains.states, + out.remains.val, out.remains.fractions, 'plasma') + figures.append(remains_fig) + if self.settings.plot_probability: + prob_fig = afx.plot_probability(out.probability.sf, out.probability.pdf, + out.probability.fractions) + figures.append(prob_fig) + + if sets.save_plot: + afx.save_plot(self.basename, *figures) + elif sets.plot_policy or sets.plot_values or sets.plot_remains or sets.plot_probability: afx.show_plot() else: afx.aft_message('nothing to plot', nme, 1) @@ -120,7 +207,6 @@ def main(): """ CLI interface to invoke the RL class """ - start = afx.timing() parser = argparse.ArgumentParser( description='Calculate optimal dose per fraction dependent on algorithm type' ) @@ -143,9 +229,7 @@ def main(): args = parser.parse_args(args=None if sys.argv[1:] else ['--help']) plan = RL_object(args.filename) - afx.aft_message('start session...', nme, 1) plan.optimise() - afx.timing(start) # show retrospective dose prescribtion afx.aft_message_list('physical tumor dose:', plan.output.physical_doses, nme, 1) @@ -165,4 +249,4 @@ def main(): if __name__ == '__main__': - main() + main() \ No newline at end of file diff --git a/src/adaptfx/aft_prompt.py b/src/adaptfx/aft_prompt.py index d3fc83c..4240845 100644 --- a/src/adaptfx/aft_prompt.py +++ b/src/adaptfx/aft_prompt.py @@ -1,19 +1,19 @@ # -*- coding: utf-8 -*- import sys import logging -import os +import adaptfx as afx prompt = 'AFT> ' empty = '' -def logging_init(filename, log, level): +def logging_init(basename, log, level): """ log initialisation to write to filename Parameters ---------- - filename : string - filename of log file + basename : string + filename of log file without suffix log : bool if true store log to filename debug: bool @@ -38,24 +38,7 @@ def logging_init(filename, log, level): log_level = logging.ERROR if log: - logfile_extension = "log" - # create logfile name - # get the basename before .json extension - basename = filename.rsplit('.')[0] - # search for existing filename ... - i = 1 - while os.path.exists(f'{basename}_{i}.{logfile_extension}'): - # exponential search if many files exist - i *= 2 - a, b = (i // 2, i) - while a+1 < b: - c = (a + b) // 2 - if os.path.exists(f'{basename}_{c}.{logfile_extension}'): - a, b = (c, b) - else: - a, b = (a, c) - # ... end of search - log_filename = f'{basename}_{b}.{logfile_extension}' + log_filename = afx.create_name(basename, "log") logging.basicConfig( format=format_file, level=log_level, diff --git a/src/adaptfx/aft_utils.py b/src/adaptfx/aft_utils.py index 33d85cb..c3fdad8 100644 --- a/src/adaptfx/aft_utils.py +++ b/src/adaptfx/aft_utils.py @@ -1,6 +1,7 @@ # -*- coding: utf-8 -*- import time -from adaptfx import aft_message, aft_error, aft_warning, aft_message +import os +import adaptfx as afx nme = __name__ def stat_rounding(number, decimal): @@ -30,7 +31,74 @@ def timing(start=None): else: stop = time.perf_counter() time_elapsed = stat_rounding((stop - start), 4) - aft_message(f'process duration: {time_elapsed} s:', nme, 1) + afx.aft_message(f'process duration: {time_elapsed} s:', nme, 1) + +def get_abs_path(filename, name): + """ + from filename create absolute path and basename + without ".appendix" + + Parameters + ---------- + filename : string + filename of instructions + name : string + logger name + + Returns + ------- + norm_path : string + normalised absolute path + basename : string + absolute path to file without suffix + + + """ + if os.path.isfile(filename): # check if filepath exists + if not os.path.isabs(filename): # check if absolute + filename = os.path.abspath(filename) + # collapse redundant separators get base + norm_path = os.path.normpath(filename) + abs_path, name = os.path.split(norm_path) + base, _ = os.path.splitext(name) + basename = os.path.join(abs_path, base) + else: + afx.aft_error(f'did not find "{filename}"', name) + return norm_path, basename + +def create_name(basename, suffix): + """ + from basename create name for a file + that does not yet exist. Increments numbers + + Parameters + ---------- + basename: string + absolute basename of file without suffix + suffix : string + suffix of the output + + Returns + ------- + new_name : string + new name nummerated + + """ + # create logfile name and + # search for existing filename ... + i = 1 + while os.path.exists(f'{basename}_{i}.{suffix}'): + # exponential search if many files exist + i *= 2 + a, b = (i // 2, i) + while a+1 < b: + c = (a + b) // 2 + if os.path.exists(f'{basename}_{c}.{suffix}'): + a, b = (c, b) + else: + a, b = (a, c) + # ... end of search + return f'{basename}_{b}.{suffix}' def key_reader(all_keys, full_dict, user_keys, algorithm): """ @@ -54,7 +122,7 @@ def key_reader(all_keys, full_dict, user_keys, algorithm): all keys copied from parameters """ - aft_message('loading keys...', nme, 1) + afx.aft_message('loading keys...', nme, 1) whole_dict = full_dict.copy() key_dict = all_keys[algorithm] @@ -63,17 +131,17 @@ def key_reader(all_keys, full_dict, user_keys, algorithm): whole_dict[key] = user_keys[key] elif key not in user_keys: if full_dict[key] == None: - aft_error(f'missing mandatory key: "{key}"', nme) + afx.aft_error(f'missing mandatory key: "{key}"', nme) else: whole_dict[key] = full_dict[key] for key in user_keys: if key not in key_dict and key in full_dict: - aft_warning( + afx.aft_warning( f'key: "{key}" is not allowed for "{algorithm}"', nme, 0 ) elif key not in full_dict: - aft_warning(f'unexpected key: "{key}"', nme, 0) + afx.aft_warning(f'unexpected key: "{key}"', nme, 0) return dict((k, whole_dict[k]) for k in key_dict) @@ -95,7 +163,7 @@ def setting_reader(all_settings, user_settings): all settings """ - aft_message('loading settings...', nme, 1) + afx.aft_message('loading settings...', nme, 1) whole_settings = all_settings.copy() for skey in all_settings: @@ -104,7 +172,7 @@ def setting_reader(all_settings, user_settings): for skey in user_settings: if skey not in all_settings: - aft_warning(f'unexpected setting: "{skey}"', nme, 0) + afx.aft_warning(f'unexpected setting: "{skey}"', nme, 0) return whole_settings diff --git a/src/adaptfx/constants.py b/src/adaptfx/constants.py index 974b9e9..a0687ff 100644 --- a/src/adaptfx/constants.py +++ b/src/adaptfx/constants.py @@ -1,9 +1,9 @@ # -*- coding: utf-8 -*- # log settings -LOG_BOOL = 0 -LOG_BOOL_LIST = [0,1] -LOG_LEVEL = 1 -LOG_LEVEL_LIST = [0,1,2] +LOG_BOOL = 0 # no logging to file +LOG_BOOL_LIST = [0,1] # options for logging +LOG_LEVEL = 1 # logging level middle +LOG_LEVEL_LIST = [0,1,2]# options for logging level # dose, sf DOSE_STEPSIZE = 0.1 @@ -11,17 +11,18 @@ SF_LOW = 0 SF_HIGH = 1.7 SF_STEPSIZE = 0.01 -SF_PROB_THRESHOLD = 1e-5 +SF_PROB_THRESHOLD = 1e-4 INF_PENALTY = 1e4 # keys ALPHA_BETA_TUMOR = 10 ALPHA_BETA_OAR = 3 -FULL_DICT = {'number_of_fractions':None, +# "None" are mandatory keys +FULL_DICT = {'number_of_fractions': None, 'fraction': 0, 'sparing_factors': None, - 'prob_update': 0, + 'prob_update': None, 'fixed_mean': None, 'fixed_std': None, 'shape': None, @@ -50,6 +51,8 @@ 'plot_policy': 0, 'plot_values': 0, 'plot_remains': 0, + 'plot_probability': 0, + 'save_plot': 0, } STANDARD_LIST = [ @@ -71,9 +74,10 @@ 'max_dose', ] +# define the list for each optimisation method OAR_LIST = STANDARD_LIST + ['tumor_goal'] -TUMOR_LIST = STANDARD_LIST + ['oar_limit'] +TUMOR_LIST = STANDARD_LIST + ['oar_limit', 'c'] FRAC_LIST = STANDARD_LIST + ['tumor_goal', 'c'] @@ -82,6 +86,6 @@ KEY_DICT = { 'oar': OAR_LIST, 'oar_old': OAR_LIST, 'tumor': TUMOR_LIST, 'tumor_old': TUMOR_LIST, - 'frac': FRAC_LIST, 'frac_old': FRAC_LIST, + 'frac': FRAC_LIST, 'tumor_oar': TUMOR_OAR_LIST, 'tumor_oar_old': TUMOR_OAR_LIST } \ No newline at end of file diff --git a/src/adaptfx/maths.py b/src/adaptfx/maths.py index eb63184..4a2f4b7 100644 --- a/src/adaptfx/maths.py +++ b/src/adaptfx/maths.py @@ -28,7 +28,7 @@ def truncated_normal(mean, std, low, upp): """ normal = truncnorm((low - mean) / std, (upp - mean) / std, - loc=mean, scale=std) + loc=mean, scale=std) return normal @@ -143,8 +143,8 @@ def sf_probdist(X, sf_low, sf_high, sf_stepsize, probability_threshold): # sum probability density from cumulative density function in interval # to get probability half_interval = sf_stepsize/2 - upp_bound = (half_interval - sf_stepsize*1e-6) * np.ones(n_samples) - low_bound = half_interval * np.ones(n_samples) + upp_bound = (half_interval - sf_stepsize*1e-6) + low_bound = half_interval prob = X.cdf(sample_sf + upp_bound) - X.cdf(sample_sf - low_bound) # get rid of all probabilities below given threshold diff --git a/src/adaptfx/planning.py b/src/adaptfx/planning.py index c400841..e4bde3d 100644 --- a/src/adaptfx/planning.py +++ b/src/adaptfx/planning.py @@ -22,18 +22,27 @@ def multiple(algorithm, keys, sets=afx.SETTING_DICT): sets = afx.DotDict(sets) if keys.fraction != 0: - # if only a specific fraction should be calculated - fractions_list = np.array([keys.fraction]) - physical_doses = np.zeros(1) - tumor_doses = np.zeros(1) - oar_doses = np.zeros(1) + # if only up to a specific fraction should be calculated + fractions_list = np.arange(1, keys.fraction + 1, 1) + physical_doses = np.zeros(keys.fraction) + tumor_doses = np.zeros(keys.fraction) + oar_doses = np.zeros(keys.fraction) else: # for calculation whole treatment in retrospect fractions_list = np.arange(1, keys.number_of_fractions + 1, 1) physical_doses = np.zeros(keys.number_of_fractions) tumor_doses = np.zeros(keys.number_of_fractions) oar_doses = np.zeros(keys.number_of_fractions) + + if sets.plot_probability: + # create lists for storing probability distribution + # note that the shape of probability distribution + # is 1) unknown and 2) non-constant --> use list + sf = [] + pdf = [] + # create array for keeping track of fractions + plot_numbering = np.arange(1, keys.number_of_fractions + 1, 1) first_tumor_dose = keys.accumulated_tumor_dose first_oar_dose = keys.accumulated_oar_dose @@ -45,6 +54,8 @@ def multiple(algorithm, keys, sets=afx.SETTING_DICT): output = afx.min_oar_bed(keys, sets) elif algorithm == 'frac': output = afx.min_n_frac(keys, sets) + elif algorithm == 'tumor': + output = afx.max_tumor_bed(keys, sets) elif algorithm == 'oar_old': output = afx.min_oar_bed_old(keys, sets) @@ -66,15 +77,18 @@ def multiple(algorithm, keys, sets=afx.SETTING_DICT): # user specifies to plot policy number, if equal to fraction plot # if both zero than the user doesn't want to plot policy + if sets.plot_probability: + sf.append(list(output.probability.sf)) + pdf.append(list(output.probability.pdf)) if sets.plot_policy == keys.fraction: output_whole.policy = output.policy - output_whole.policy.fractions = fractions_list[sets.plot_policy - 1:] + output_whole.policy.fractions = plot_numbering[sets.plot_policy - 1:] if sets.plot_values == keys.fraction: output_whole.value = output.value - output_whole.value.fractions = fractions_list[sets.plot_values - 1:] + output_whole.value.fractions = plot_numbering[sets.plot_values - 1:] if sets.plot_remains == keys.fraction: output_whole.remains = output.remains - output_whole.remains.fractions = fractions_list[sets.plot_remains - 1:] + output_whole.remains.fractions = plot_numbering[sets.plot_remains - 1:] # store doses exponent = afx.find_exponent(sets.dose_stepsize) - 1 @@ -83,5 +97,11 @@ def multiple(algorithm, keys, sets=afx.SETTING_DICT): output_whole.oar_sum, output_whole.tumor_sum = np.around( [np.nansum(oar_doses), np.nansum(tumor_doses)], -exponent) output_whole.fractions_used = np.count_nonzero(~np.isnan(output_whole.physical_doses)) + if sets.plot_probability: + # store probability distribution + output_whole.probability = {} + output_whole.probability.sf = sf + output_whole.probability.pdf = pdf + output_whole.probability.fractions = fractions_list return output_whole \ No newline at end of file diff --git a/src/adaptfx/radiobiology.py b/src/adaptfx/radiobiology.py index 10e820a..081fddd 100644 --- a/src/adaptfx/radiobiology.py +++ b/src/adaptfx/radiobiology.py @@ -1,5 +1,7 @@ # -*- coding: utf-8 -*- import numpy as np +import adaptfx as afx +import scipy.optimize as opt def bed_calc0(dose, ab, sf=1): """ @@ -61,7 +63,7 @@ def convert_to_physical(bed, ab, sf=1): ab : float alpha-beta ratio. sf : float/array - sparing factor, only specify when OAR BED + sparing factor, only specify when OAR BED. Returns ------- @@ -73,3 +75,116 @@ def convert_to_physical(bed, ab, sf=1): physical_dose = (-sf + np.sqrt(sf**2 + 4 * sf**2 * bed_array / ab)) / ( 2 * sf**2 / ab) return physical_dose + +def cost_func(keys, n_list, n_samples): + """ + For a specified list of maximum number of fractions + simulates average cumulative OAR BED for uniform-, + adaptive- and optimal-fractionation (theoretical optimum) + + Parameters + ---------- + keys : dict + algorithm instructions. + n_list : array + array of maximum number of fractions + n_samples : int + number of patients to sample + + Returns + ------- + uft : array + uniform fractionated average cumulative BED + aft : array + adaptive fractionated average cumulative BED + opt : array + optimal fractionated average cumulative BED + """ + para = keys + BED_matrix = np.zeros((3,len(n_list), n_samples)) + for i, n_max in enumerate(n_list): + dose = keys.abt/(2*n_max) * (np.sqrt(n_max ** 2 + 4*n_max*keys.tumor_goal/keys.abt) - n_max) + para.number_of_fractions = n_max + for j in range(n_samples): + sf_list = np.random.normal(para.fixed_mean, para.fixed_std, n_max+1) + para.sparing_factors = sf_list + # calculate uniform fractionation + BED_matrix[0][i][j] = np.sum(bed_calc0(dose, para.abn, sf_list[1:])) + # calculate adaptive fractionation + BED_matrix[1][i][j] = afx.multiple('oar', para).oar_sum + # calculate optimal fractionation if all sparing factors are known at the beginning + def bed_calc_d(dose): + cum_bed = afx.bed_calc0(dose, para.abn, sf_list[1:]) + return np.sum(cum_bed) + d_in = dose * np.ones(n_max) + # define tumor BED constraint + cons = [{'type': 'eq', 'fun': lambda x: np.sum(afx.bed_calc0(x, para.abt)) - para.tumor_goal}] + BED_matrix[2][i][j] = opt.minimize(bed_calc_d, x0=d_in, constraints=cons).fun + + BED_means = np.mean(BED_matrix, axis=2) + BED_uft, BED_aft, BED_opt = BED_means[0], BED_means[1], BED_means[2] + return BED_uft, BED_aft, BED_opt + + +def c_calc(keys, n_target, n_samples, plot=False): + """ + For a specified targeted number of fractions + gives the optimal C, when minimising OAR BED + + Parameters + ---------- + keys : dict + algorithm instructions. + + n_target : int + targeted number of fractions + + Returns + ------- + c : positive float + optimal parameter for achieving n_pres fractions. + """ + n_upper = keys.number_of_fractions + n_list = np.arange(1, n_upper + 3) + + def cost_fit_func(n_list, a, b): + """ + Fit function for total OAR BED cost, + derived from optimal fraction decision-making + + Parameters + ---------- + a, b : float + fit parameter. + + Returns + ------- + cost : positive float + total OAR BED + """ + curvature = 4 * keys.tumor_goal / keys.abt + cost = a * (n_list - np.sqrt(n_list**2 + curvature * n_list)) + b + return cost + + def d_cost_fit_func(a, n_target): + """ + Negative derivative of cost_fit_func at n_target + """ + curvature = 4 * keys.tumor_goal / keys.abt + d_cost_fit_func = a - (a*curvature + 2 * a * n_target)/(2*np.sqrt(n_target * (curvature + n_target))) + return -d_cost_fit_func + + if n_upper <= n_target: + c_opt = 0 + else: + y_no, y_aft, y_opt = cost_func(keys, n_list, n_samples) + [a_opt, b_opt], _ = opt.curve_fit(cost_fit_func, n_list, y_aft) + c_opt = d_cost_fit_func(a_opt, n_target) + + if plot: + bed_results = {'uniform':y_no, 'adaptive':y_aft, 'optimal':y_opt, + 'aft_fit':cost_fit_func(n_list, a_opt, b_opt)} + fig = afx.plot_accumulated_bed(n_list, bed_results) + + return c_opt + diff --git a/src/adaptfx/reinforce.py b/src/adaptfx/reinforce.py index 646b213..1fdf8ae 100644 --- a/src/adaptfx/reinforce.py +++ b/src/adaptfx/reinforce.py @@ -39,6 +39,7 @@ def min_n_frac(keys, sets=afx.SETTING_DICT): max_dose = keys.max_dose # ---------------------------------------------------------------------- # # check in which fraction data should be returned for plotting + probability_plot = 1 if sets.plot_probability else 0 policy_plot = 1 if sets.plot_policy == fraction else 0 values_plot = 1 if sets.plot_values == fraction else 0 remains_plot = 1 if sets.plot_remains == fraction else 0 @@ -133,11 +134,9 @@ def min_n_frac(keys, sets=afx.SETTING_DICT): future_bedt = accumulated_tumor_dose + bedt_space future_bedt = np.where(future_bedt > tumor_goal, tumor_limit, future_bedt) c_penalties = np.where(np.round(future_bedt, -exp) < tumor_goal, -c, 0) - # for discrete matching - # future_values = future_values_discrete[np.where(np.isclose(bedt_states, future_bedt))] future_values = afx.interpolate(future_bedt, bedt_states, future_values_discrete) vs = -bedn_space + future_values + c_penalties - # argmax of vs along axis 0 to find best action fot the actual sf + # argmax of vs along axis 0 to find best action for the actual sf action_index = vs.argmax(axis=0) if policy_plot or values_plot or remains_plot: @@ -187,13 +186,13 @@ def min_n_frac(keys, sets=afx.SETTING_DICT): # penalties need to be reshaped for broadcasting vs = -last_bedn + penalties.reshape(n_bedt_states, 1) values[fraction_index] = vs - # ensure that for the goal reached the value/poliy is zero (min_dose) + # ensure that for the goal reached the value/policy is zero (min_dose) values[fraction_index][bedt_states==tumor_goal] = 0 if policy_plot: # policy calculation for each bedt, but sf is not considered policy[fraction_index] += last_bed_actions.reshape(n_bedt_states, 1) - # ensure that for the goal reached the value/poliy is zero (min_dose) + # ensure that for the goal reached the value/policy is zero (min_dose) policy[fraction_index][bedt_states==tumor_goal] = 0 elif fraction_state != number_of_fractions: @@ -210,7 +209,7 @@ def min_n_frac(keys, sets=afx.SETTING_DICT): vs = -bedn_sf_space + future_values.reshape(n_bedt_states, n_action, 1) + c_penalties # check vs along the sf axis values[fraction_index] = vs.max(axis=1) - # ensure that for the goal reached the value/poliy is zero (min_dose) + # ensure that for the goal reached the value/policy is zero (min_dose) values[fraction_index][bedt_states==tumor_goal] = 0 if policy_plot or remains_plot: @@ -231,20 +230,266 @@ def min_n_frac(keys, sets=afx.SETTING_DICT): output.tumor_dose = bedt_space[action_index] if not finished else np.nan output.oar_dose = bedn_space[action_index] if not finished else np.nan + if probability_plot: + output.probability = {} + output.probability.sf = sf + output.probability.pdf = rv.pdf(sf) + if policy_plot: + output.policy = {} + output.policy.val = policy + output.policy.sf = sf + output.policy.states = bedt_states + if values_plot: + output.value = {} + output.value.val = values + output.value.sf = sf + output.value.states = bedt_states + if remains_plot: + output.remains = {} + output.remains.val = remains + output.remains.sf = sf + output.remains.states = bedt_states + + return output + + +def max_tumor_bed(keys, sets=afx.SETTING_DICT): + # check if keys is a dictionary from manual user + if isinstance(keys, dict): + keys = afx.DotDict(keys) + + if isinstance(sets, dict): + sets = afx.DotDict(sets) + + fraction = keys.fraction + number_of_fractions = keys.number_of_fractions + accumulated_oar_dose = keys.accumulated_oar_dose + sparing_factors_public = keys.sparing_factors_public + prob_update = keys.prob_update + fixed_mean = keys.fixed_mean + fixed_std = keys.fixed_std + shape = keys.shape + scale = keys.scale + shape_inv = keys.shape_inv + scale_inv = keys.scale_inv + oar_limit = keys.oar_limit + c = keys.c + abt = keys.abt + abn = keys.abn + min_dose = keys.min_dose + max_dose = keys.max_dose + # ---------------------------------------------------------------------- # + # check in which fraction data should be returned for plotting + probability_plot = 1 if sets.plot_probability else 0 + policy_plot = 1 if sets.plot_policy == fraction else 0 + values_plot = 1 if sets.plot_values == fraction else 0 + remains_plot = 1 if sets.plot_remains == fraction else 0 + + # prepare distribution + actual_sf = sparing_factors_public[fraction] + if prob_update == 0: + # fixed normal distribution for sparing factor + mean = fixed_mean + std = fixed_std + rv = afx.truncated_normal(mean, std, sets.sf_low, sets.sf_high) + elif prob_update == 1: + # update normal distribution with gamma prior + mean = np.mean(sparing_factors_public) + std = afx.std_posterior(sparing_factors_public, shape, scale) + rv = afx.truncated_normal(mean, std, sets.sf_low, sets.sf_high) + elif prob_update == 2: + # update distribution with inverse-gamma prior + # posterior predictive distribution is student-t distribution + rv = afx.student_t(sparing_factors_public, shape_inv, scale_inv) + else: + afx.aft_error('invalid "prob_update" key was set', nme) + # initialise distribution from random variable (rv) + [sf, prob] = afx.sf_probdist(rv, sets.sf_low, sets.sf_high, + sets.sf_stepsize, sets.sf_prob_threshold) + n_sf = len(sf) + + # actionspace + exp = afx.find_exponent(sets.dose_stepsize) + accumulated_oar_dose = np.round(accumulated_oar_dose, -exp) + remaining_bed = np.round(oar_limit - accumulated_oar_dose, -exp) + n_bedsteps = int(np.round(remaining_bed / sets.dose_stepsize, 0)) + n_statesteps = int(np.round(remaining_bed / sets.state_stepsize, 0)) + + # automatic max_dose calculation + if max_dose == -1: + max_dose = remaining_bed + # Reduce max_dose to prohibit tumor_goal overshoot (efficiency) + max_dose = min(max_dose, remaining_bed) + + # actionspace in bed dose + # pre calculation + bedn_space_pre = np.linspace(0, remaining_bed, n_bedsteps + 1) + actionspace_pre = afx.convert_to_physical(bedn_space_pre, abn, actual_sf) + bedt_space_pre = afx.bed_calc0(actionspace_pre, abt) + range_action = (bedn_space_pre >= min_dose) & (bedn_space_pre <= max_dose) + if range_action.any(): + # check if actionspace is empty + bedn_space = bedn_space_pre[range_action] + actionspace = actionspace_pre[range_action] + bedt_space = bedt_space_pre[range_action] + n_action = len(bedn_space) + bedn_space_sf = np.zeros((1, n_sf)) + bedn_space.reshape(n_action, 1) + # bed_space to relate actionspace to oar- and tumor-dose + actionspace_sf = afx.convert_to_physical(bedn_space.reshape(n_action, 1), abn, sf) + else: + bedn_space_sf = np.zeros((1, n_sf)) + np.array([min_dose]) + actionspace_sf = afx.convert_to_physical(bedn_space_sf, abn, sf) + n_action = 1 + bedt_space_sf = afx.bed_calc0(actionspace_sf, abt) + + # tumor bed states for tracking dose + oar_upper_bound = oar_limit + sets.state_stepsize + # include at least one more step for bedt + # define number of bed_dose steps to fulfill stepsize + bedn_states = np.linspace(accumulated_oar_dose, + oar_limit, n_statesteps + 1) + remaining_states = bedn_states - accumulated_oar_dose + n_bedn_states = len(bedn_states) + + # values matrix + # dim(values) = dim(policy) = fractions_remaining * bedt * sf + n_remaining_fractions = number_of_fractions - fraction + values = np.zeros((n_remaining_fractions + 1, n_bedn_states, n_sf)) + if policy_plot or values_plot or remains_plot: + policy = np.zeros((n_remaining_fractions + 1, n_bedn_states, n_sf)) + remains = np.zeros((n_remaining_fractions + 1, n_bedn_states, n_sf)) + + finished = False + # ------------------------------------------------------------------------------------- # + remaining_fractions = np.arange(number_of_fractions, fraction - 1, -1) + remaining_index = remaining_fractions - fraction + # note that lowest fraction_state is one not zero + # and remaining_index counts in python indices + for fraction_index, fraction_state in zip(remaining_index, remaining_fractions): + if remaining_bed <= 0: + finished = True + break + + elif fraction_state == fraction and fraction != number_of_fractions: + # state is the actual fraction to calculate + # e.g. in the first fraction_state there is no prior dose delivered + # and future_bedt is equal to bedt_space + future_values_discrete = (values[fraction_index + 1] * prob).sum(axis=1) + future_bedn = accumulated_oar_dose + bedn_space + future_bedn = np.where(future_bedn > oar_limit, oar_upper_bound, future_bedn) + c_penalties = np.where(np.round(future_bedn, -exp) < oar_limit, -c, 0) + future_values = afx.interpolate(future_bedn, bedn_states, future_values_discrete) + vs = bedt_space + future_values + c_penalties + # argmax of vs along axis 0 to find best action for the actual sf + action_index = vs.argmax(axis=0) + + if policy_plot or values_plot or remains_plot: + # for the policy plot + future_bedn_full = bedn_states.reshape(n_bedn_states, 1, 1) + bedn_space_sf.reshape(1, n_action, n_sf) + future_bedn_full = np.where(future_bedn_full > oar_limit, oar_upper_bound, future_bedn_full) + c_penalties_full = np.where(future_bedn_full < oar_limit, -c, 0) + future_values_full = afx.interpolate(future_bedn_full, bedn_states, future_values_discrete) + vs_full = bedt_space_sf.reshape(1, n_action, n_sf) + future_values_full + c_penalties_full + # check vs along the sf axis + current_policy = bedn_space[vs_full.argmax(axis=1)] + # ensure that for the goal reached the value/policy is zero (min_dose) + current_policy[bedn_states == oar_limit] = 0 + future_remains_discrete = (remains[fraction_index + 1] * prob).sum(axis=1) + future_bedn_opt = current_policy + bedn_states.reshape(n_bedn_states, 1) + future_remains = afx.interpolate(future_bedn_opt, bedn_states, future_remains_discrete) + current_remains = np.where((current_policy - remaining_states[::-1].reshape(n_bedn_states, 1)) >= 0, 0, 1) + # write to arrays + policy[fraction_index] = current_policy + values[fraction_index] = vs_full.max(axis=1) + remains[fraction_index] = current_remains + future_remains + + elif fraction == number_of_fractions: + # in the last fraction value is not relevant + # max_dose is the already calculated remaining dose + # this should compensate the discretisation of the state space + action_index = np.abs(remaining_bed - bedn_space).argmin() + + if policy_plot: + policy[fraction_index][0] = bedn_space[action_index] + + elif fraction_state == number_of_fractions: + # final state to initialise terminal reward + # dose remaining to be delivered, this is the actionspace in bedt + last_bed_actions = np.round(oar_limit - bedn_states, -exp) + # cut the actionspace to min and max dose constraints + last_bed_actions = np.where(last_bed_actions < min_dose, min_dose, last_bed_actions) + last_bed_actions = np.where(last_bed_actions > max_dose, max_dose, last_bed_actions) + last_bed_actions_reshaped = last_bed_actions.reshape(n_bedn_states, 1) + last_actions = afx.convert_to_physical(last_bed_actions_reshaped, abn, sf) + last_bedt = afx.bed_calc0(last_actions, abt) + # this smooths out the penalties in underdose and overdose regions + bedn_diff = np.round(bedn_states.reshape(n_bedn_states, 1) + last_bed_actions_reshaped - oar_limit, -exp) + penalties = np.where(bedn_diff == 0, 0, -np.abs(bedn_diff) * sets.inf_penalty) + # to each best action add the according penalties + # penalties need to be reshaped for broadcasting + vs = last_bedt + penalties + values[fraction_index] = vs + # ensure that for the goal reached the value/policy is zero (min_dose) + values[fraction_index][bedn_states==oar_limit] = 0 + + if policy_plot: + # policy calculation for each bedt, but sf is not considered + policy[fraction_index] += last_bed_actions_reshaped + # ensure that for the goal reached the value/policy is zero (min_dose) + policy[fraction_index][bedn_states==oar_limit] = 0 + + elif fraction_state != number_of_fractions: + # every other state but the last + # this calculates the value function in the future fractions + future_values_discrete = (values[fraction_index + 1] * prob).sum(axis=1) + # bedt_states is reshaped such that numpy broadcast leads to 2D array + future_bedn = bedn_states.reshape(n_bedn_states, 1, 1) + bedn_space_sf.reshape(1, n_action, n_sf) + future_bedn = np.where(future_bedn > oar_limit, oar_upper_bound, future_bedn) + future_values = afx.interpolate(future_bedn, bedn_states, future_values_discrete) + c_penalties = np.where(future_bedn < oar_limit, -c, 0) + # dim(bedn_sf_space)=(1,n_action,n_sf),dim(future_values)=(n_states,n_action) + # every row of values_penalties is transposed and copied n_sf times + vs = bedt_space_sf.reshape(1, n_action, n_sf) + future_values + c_penalties + # check vs along the sf axis + values[fraction_index] = vs.max(axis=1) + # ensure that for the goal reached the value/policy is zero (min_dose) + values[fraction_index][bedn_states==oar_limit] = 0 + + if policy_plot or remains_plot: + current_policy = bedn_space[vs.argmax(axis=1)] + # ensure that for the goal reached the value/policy is zero (min_dose) + current_policy[bedn_states == oar_limit] = 0 + future_remains_discrete = (remains[fraction_index + 1] * prob).sum(axis=1) + future_bedn_opt = current_policy + bedn_states.reshape(n_bedn_states, 1) + future_remains = afx.interpolate(future_bedn_opt, bedn_states, future_remains_discrete) + current_remains = np.where((current_policy - remaining_states[::-1].reshape(n_bedn_states, 1)) >= 0, 0, 1) + # write to arrays + policy[fraction_index] = current_policy + remains[fraction_index] = current_remains + future_remains + + output = afx.DotDict({}) + output.physical_dose = actionspace[action_index] if not finished else np.nan + output.tumor_dose = bedt_space[action_index] if not finished else np.nan + output.oar_dose = bedn_space[action_index] if not finished else np.nan + + if probability_plot: + output.probability = {} + output.probability.sf = sf + output.probability.pdf = rv.pdf(sf) if policy_plot: output.policy = {} output.policy.val = policy output.policy.sf = sf - output.policy.states = remaining_states + output.policy.states = bedn_states if values_plot: output.value = {} output.value.val = values output.value.sf = sf - output.value.states = remaining_states + output.value.states = bedn_states if remains_plot: output.remains = {} output.remains.val = remains output.remains.sf = sf - output.remains.states = remaining_states + output.remains.states = bedn_states return output \ No newline at end of file diff --git a/src/adaptfx/reinforce_old.py b/src/adaptfx/reinforce_old.py index 1d3ae0d..775cbef 100644 --- a/src/adaptfx/reinforce_old.py +++ b/src/adaptfx/reinforce_old.py @@ -65,26 +65,26 @@ def min_oar_bed_old(keys, sets=SETTING_DICT): number_of_fractions = keys.number_of_fractions accumulated_tumor_dose = keys.accumulated_tumor_dose sparing_factors_public = keys.sparing_factors_public - alpha = keys.alpha - beta = keys.beta + shape = keys.shape + scale = keys.scale tumor_goal = keys.tumor_goal abt = keys.abt abn = keys.abn min_dose = keys.min_dose max_dose = keys.max_dose - fixed_prob = keys.fixed_prob + prob_update = keys.prob_update fixed_mean = keys.fixed_mean fixed_std = keys.fixed_std # ---------------------------------------------------------------------- # # check in which fraction policy should be returned policy_plot = 1 if sets.plot_policy == fraction else 0 - if fixed_prob != 1: + if prob_update == 1: mean = np.mean( sparing_factors_public ) # extract the mean and std to setup the sparingfactor distribution - standard_deviation = std_posterior(sparing_factors_public, alpha, beta) - if fixed_prob == 1: + standard_deviation = std_posterior(sparing_factors_public, shape, scale) + elif prob_update == 0: mean = fixed_mean standard_deviation = fixed_std X = truncated_normal(mean, standard_deviation, sets.sf_low, sets.sf_high) @@ -224,207 +224,20 @@ def min_oar_bed_old(keys, sets=SETTING_DICT): -def min_n_frac_old(keys, sets=SETTING_DICT): - fraction = keys.fraction - number_of_fractions = keys.number_of_fractions - accumulated_tumor_dose = keys.accumulated_tumor_dose - sparing_factors_public = keys.sparing_factors_public - alpha = keys.alpha - beta = keys.beta - tumor_goal = keys.tumor_goal - c = keys.c - abt = keys.abt - abn = keys.abn - min_dose = keys.min_dose - max_dose = keys.max_dose - fixed_prob = keys.fixed_prob - fixed_mean = keys.fixed_mean - fixed_std = keys.fixed_std - # ---------------------------------------------------------------------- # - # check in which fraction policy should be returned - policy_plot = 1 if sets.plot_policy == fraction else 0 - - if fixed_prob != 1: - mean = np.mean( - sparing_factors_public - ) # extract the mean and std to setup the sparingfactor distribution - standard_deviation = std_posterior(sparing_factors_public, alpha, beta) - if fixed_prob == 1: - mean = fixed_mean - standard_deviation = fixed_std - X = truncated_normal(mean, standard_deviation, sets.sf_low, sets.sf_high) - bedt_states = np.arange(accumulated_tumor_dose, tumor_goal, 1) - bedt_states = np.concatenate( - (bedt_states, [tumor_goal, tumor_goal + 1]) - ) # add an extra step outside of our prescribed tumor dose which will be penalised to make sure that we aim at the prescribe tumor dose - [sf, prob] = sf_probdist(X, sets.sf_low, sets.sf_high, - sets.sf_stepsize, sets.sf_prob_threshold) - # we prepare an empty values list and open an action space which is equal to all the dose numbers that can be given in one fraction - values = np.zeros( - (number_of_fractions - fraction, len(bedt_states), len(sf)) - ) # 2d values list with first indice being the BED and second being the sf - max_physical_dose = convert_to_physical(tumor_goal, abt) - if max_dose == -1: - max_dose = max_physical_dose - elif max_dose > max_physical_dose: - # if the max dose is too large we lower it, so we dont needlessly check too many actions - max_dose = max_physical_dose - if min_dose > max_dose: - min_dose = max_dose - 0.1 - actionspace = np.arange(min_dose, max_dose + 0.1, 0.1) - # now we set up the policy array which has len(BEDT)*len(sf)*len(actionspace) entries. We give each action the same probability to start with - policy = np.zeros((number_of_fractions - fraction, len(bedt_states), len(sf))) - - for state, fraction_state in enumerate( - np.arange(number_of_fractions, fraction -1, -1) - ): - if ( - state == number_of_fractions - 1 #fraction_state == 1 - ): # first state with no prior dose delivered so we dont loop through BEDT - bedn = bed_calc_matrix( - sparing_factors_public[-1], abn, actionspace - ) # calculate all delivered doses to the Normal tissues (the penalty) - future_bedt = bed_calc0(actionspace, abt) - future_values_func = interp1d(bedt_states, (values[state - 1] * prob).sum(axis=1)) - future_values = future_values_func( - future_bedt - ) # for each action and sparing factor calculate the penalty of the action and add the future value we will only have as many future values as we have actions (not sparing dependent) - c_penalties = np.zeros(future_bedt.shape) - c_penalties[future_bedt < tumor_goal] = -c - vs = -bedn + c_penalties + future_values - policy4 = vs.argmax(axis=1) - elif ( - accumulated_tumor_dose >= tumor_goal - ): - best_action = 0 - last_BEDN = bed_calc0(best_action, abn, sparing_factors_public[-1]) - policy4 = 0 - break - elif ( - fraction_state == fraction and fraction != number_of_fractions - ): # actual fraction is first state to calculate - actionspace_clipped = actionspace[ - 0 : max_action(accumulated_tumor_dose, actionspace, tumor_goal) + 1 - ] - bedn = bed_calc_matrix(sparing_factors_public[-1], abn, actionspace_clipped) - future_bedt = accumulated_tumor_dose + bed_calc0(actionspace_clipped, abt) - future_bedt[future_bedt > tumor_goal] = tumor_goal + 1 - penalties = np.zeros(future_bedt.shape) - c_penalties = np.zeros(future_bedt.shape) - penalties[future_bedt > tumor_goal] = 0 #overdosing is indirectly penalised with BEDN - c_penalties[future_bedt < tumor_goal] = -c - future_values_func = interp1d(bedt_states, (values[state - 1] * prob).sum(axis=1)) - future_values = future_values_func( - future_bedt - ) # for each action and sparing factor calculate the penalty of the action and add the future value we will only have as many future values as we have actions (not sparing dependent) - vs = -bedn + c_penalties + future_values + penalties - policy4 = vs.argmax(axis=1) - elif ( - fraction == number_of_fractions - ): # in this state no penalty has to be defined as the value is not relevant - best_action = convert_to_physical(tumor_goal-accumulated_tumor_dose, abt) - if accumulated_tumor_dose > tumor_goal: - best_action = 0 - if best_action < min_dose: - best_action = min_dose - if best_action > max_dose: - best_action = max_dose - last_BEDN = bed_calc0(best_action, abn, sparing_factors_public[-1]) - policy4 = best_action * 10 - else: - future_value_prob = (values[state - 1] * prob).sum(axis=1) - future_values_func = interp1d(bedt_states, future_value_prob) - for tumor_index, tumor_value in enumerate( - bedt_states - ): # this and the next for loop allow us to loop through all states - actionspace_clipped = actionspace[ - 0 : max_action(tumor_value, actionspace, tumor_goal) + 1 - ] # we only allow the actions that do not overshoot - bedn = bed_calc_matrix( - sf, abn, actionspace_clipped - ) # this one could be done outside of the loop and only the clipping would happen inside the loop. - bed = bed_calc_matrix(np.ones(len(sf)), abt, actionspace_clipped) - if state != 0 and tumor_value < tumor_goal: - future_bedt = tumor_value + bed - future_bedt[future_bedt > tumor_goal] = tumor_goal + 1 - future_values = future_values_func( - future_bedt - ) # for each action and sparing factor calculate the penalty of the action and add the future value we will only have as many future values as we have actions (not sparing dependent) - c_penalties = np.zeros(future_bedt.shape) - c_penalties[future_bedt < tumor_goal] = -c - vs = -bedn + c_penalties + future_values - if vs.size == 0: - best_action = np.zeros(len(sf)) - valer = np.zeros(len(sf)) - else: - best_action = actionspace_clipped[vs.argmax(axis=1)] - valer = vs.max(axis=1) - elif tumor_value > tumor_goal: #calculate value for terminal case - best_action = 0 - future_bedt = tumor_value - valer = ( - + np.zeros(sf.shape) - ) - else: # last state no more further values to add - best_action = convert_to_physical(tumor_goal-tumor_value, abt) - if best_action > max_dose: - best_action = max_dose - if best_action < min_dose: - best_action = min_dose - last_BEDN = bed_calc0(best_action, abn, sf) - future_bedt = tumor_value + bed_calc0(best_action, abt) - underdose_penalty = 0 - overdose_penalty = 0 - if future_bedt < tumor_goal: - underdose_penalty = -10000 - if future_bedt > tumor_goal: - overdose_penalty = 0 - valer = ( - -last_BEDN - + underdose_penalty * np.ones(sf.shape) - + overdose_penalty * np.ones(sf.shape) - ) # gives the value of each action for all sparing factors. elements 0-len(sparingfactors) are the Values for - - policy[state][tumor_index] = best_action - values[state][tumor_index] = valer - - if sets.plot_policy or sets.plot_value: - policy = np.concatenate([np.zeros((1, len(bedt_states), len(sf))), policy[::-1]]) - values = np.concatenate([np.zeros((1, len(bedt_states), len(sf))), values[::-1]]) - if fraction != number_of_fractions: - optimal_action = actionspace[policy4] - if fraction == number_of_fractions: - optimal_action = policy4 / 10 - tumor_dose = bed_calc0(optimal_action, abt) - oar_dose = bed_calc0(optimal_action, abn, sparing_factors_public[-1]) - - output = DotDict({}) - - output.physical_dose = optimal_action - output.tumor_dose = tumor_dose - output.oar_dose = oar_dose - - return output - - - - - - def max_tumor_bed_old(keys, sets=SETTING_DICT): fraction = keys.fraction number_of_fractions = keys.number_of_fractions - accumulated_oar_dose = keys.accumulated_tumor_dose + accumulated_oar_dose = keys.accumulated_oar_dose sparing_factors_public = keys.sparing_factors_public - alpha = keys.alpha - beta = keys.beta + shape = keys.shape + scale = keys.scale oar_limit = keys.oar_limit abt = keys.abt abn = keys.abn min_dose = keys.min_dose max_dose = keys.max_dose - fixed_prob = keys.fixed_prob + prob_update = keys.prob_update fixed_mean = keys.fixed_mean fixed_std = keys.fixed_std # ---------------------------------------------------------------------- # @@ -432,12 +245,12 @@ def max_tumor_bed_old(keys, sets=SETTING_DICT): policy_plot = 1 if sets.plot_policy == fraction else 0 actual_sparing = sparing_factors_public[-1] - if fixed_prob != 1: + if prob_update == 1: mean = np.mean( sparing_factors_public ) # extract the mean and std to setup the sparingfactor distribution - standard_deviation = std_posterior(sparing_factors_public, alpha, beta) - if fixed_prob == 1: + standard_deviation = std_posterior(sparing_factors_public, shape, scale) + elif prob_update == 0: mean = fixed_mean standard_deviation = fixed_std X = truncated_normal(mean, standard_deviation, sets.sf_low, sets.sf_high) @@ -634,27 +447,27 @@ def min_oar_max_tumor_old(keys, sets=SETTING_DICT): accumulated_tumor_dose = keys.accumulated_tumor_dose accumulated_oar_dose = keys.accumulated_oar_dose sparing_factors_public = keys.sparing_factors_public - alpha = keys.alpha - beta = keys.beta + shape = keys.shape + scale = keys.scale tumor_goal = keys.tumor_goal oar_limit = keys.oar_limit abt = keys.abt abn = keys.abn min_dose = keys.min_dose max_dose = keys.max_dose - fixed_prob = keys.fixed_prob + prob_update = keys.prob_update fixed_mean = keys.fixed_mean fixed_std = keys.fixed_std # ---------------------------------------------------------------------- # # check in which fraction policy should be returned policy_plot = 1 if sets.plot_policy == fraction else 0 - if fixed_prob != 1: + if prob_update == 1: mean = np.mean( sparing_factors_public ) # extract the mean and std to setup the sparingfactor distribution - standard_deviation = std_posterior(sparing_factors_public, alpha, beta) - if fixed_prob == 1: + standard_deviation = std_posterior(sparing_factors_public, shape, scale) + elif prob_update == 0: mean = fixed_mean standard_deviation = fixed_std X = truncated_normal(mean, standard_deviation, sets.sf_low, sets.sf_high) diff --git a/src/adaptfx/visualiser.py b/src/adaptfx/visualiser.py index c0bc0ec..13b8303 100644 --- a/src/adaptfx/visualiser.py +++ b/src/adaptfx/visualiser.py @@ -3,8 +3,9 @@ import matplotlib.pyplot as plt from matplotlib.colors import Normalize as normalise import matplotlib.cm as cm +import adaptfx as afx -def plot_val(sfs, states, policies, fractions): +def plot_val(sfs, states, data, fractions, colmap='turbo'): """ creates a subplot grid for the policy in each fraction that is given @@ -14,8 +15,12 @@ def plot_val(sfs, states, policies, fractions): 1d array with dimension n states : array 1d array with dimension m - policies : array + data : array 2d array with dimension n_policies:n:m + fractions : array + 1d array with fractions numbers + colmap : string + heat colourmap for matplotlib Returns ------- @@ -23,22 +28,30 @@ def plot_val(sfs, states, policies, fractions): matplotlib pyplot figure """ - [n_policies, _, _] = policies.shape + if colmap == 'turbo': + label = r'Policy $\pi$ in BED$_{10}$ [Gy]' + elif colmap == 'viridis': + label = r'Value $v$' + elif colmap == 'plasma': + label = r'Expected Remaining Number $\varepsilon$' + else: + label = 'empty' + [n_grids, _, _] = data.shape # search for optimal rectangular size of subplot grid - n_rows = n_columns = int(np.sqrt(n_policies)) - while n_rows * n_columns < n_policies: + n_rows = n_columns = int(np.sqrt(n_grids)) + while n_rows * n_columns < n_grids: if n_rows >= n_columns: n_columns += 1 else: n_rows += 1 # initiate plot and parameters fig, ax = plt.subplots(n_rows, n_columns) - x_min, x_max, y_min, y_max = sfs[0], sfs[-1], states[0], states[-1] + x_min, x_max, y_min, y_max = np.min(sfs), np.max(sfs), np.min(states), np.max(states) # create shared colorbar - colmin, colmax = np.min(policies), np.max(policies) + colmin, colmax = np.min(data), np.max(data) normaliser = normalise(colmin, colmax) - colormap = cm.get_cmap('turbo') + colormap = cm.get_cmap(colmap) im = cm.ScalarMappable(cmap=colormap, norm=normaliser) # loop through the axes @@ -47,19 +60,95 @@ def plot_val(sfs, states, policies, fractions): except: # in case ax is a 1x1 subplot axs = np.array([ax]) - - for i, pol in enumerate(policies): + # turn off axes + for a in axs: + a.axis(False) + + for i, pol in enumerate(data): + axs[i].axis(True) axs[i].imshow(pol, interpolation=None, origin='upper', norm=normaliser, cmap=colormap, aspect='auto', - extent=[x_min, x_max, y_min, y_max]) - axs[i].set_title(fractions[i]) - - fig.supxlabel('sparing factor') - fig.supylabel('remaining BED') + extent=[x_min, x_max, y_max, y_min]) + axs[i].set_title(rf'$t = {fractions[i]}$', loc='left') + try: # get rid of inner axes values + axs[i].label_outer() + except: + pass + + T_text = r'{\mathrm{T}}' + fig.supxlabel(r'$\delta$') + fig.supylabel(rf'$B^{T_text}$ [Gy]') fig.tight_layout() - fig.colorbar(mappable=im, ax=axs.tolist()) + fig.colorbar(mappable=im, ax=axs.tolist(), label=label) + + return fig + +def plot_accumulated_bed(n_list, bed_dict): + """ + creates plot for simulated adaptive fractionation therapies + + Parameters + ---------- + n_list : array + 1d array for number of fractions + accumulated_bed : list of arrays + list of 1d array for average accumulated bed + + Returns + ------- + fig : matplotlib.pyplot.figure + matplotlib pyplot figure + """ + fig, ax = plt.subplots(1, 1) + for key in bed_dict: + ax.plot(n_list, bed_dict[key], label=key) + + fig.supylabel('accumulated BED') + fig.supxlabel('total number of fractions $n$') + plt.legend() + + return fig + +def plot_probability(sf_list, pdf_list, fractions_list): + """ + creates plot for multiple probability density functions + + Parameters + ---------- + sf_list : list + list with sf lists + pdf_list : list + list with probability density values list + fractions_list : list + list with specified fraction + Returns + ------- + fig : matplotlib.pyplot.figure + matplotlib pyplot figure + """ + fig, ax = plt.subplots(1, 1) + for i, fraction in enumerate(fractions_list): + ax.plot(sf_list[i], pdf_list[i], label=rf'$t={fraction}$') + + fig.supxlabel(r'$\delta$') + fig.supylabel(r'$PDF(\delta)$') + fig.tight_layout() + plt.legend() + return fig + def show_plot(): - plt.show() \ No newline at end of file + plt.show() + +def save_plot(basename, *figures): + if len(figures)==1: + figures[0].savefig(f'{basename}.pdf', format='pdf') + plt.clf() + plt.close() + else: + for fig in figures: + # create name from base, search for untaken name + fig_name = afx.create_name(basename, 'pdf') + fig.savefig(fig_name, format='pdf') \ No newline at end of file diff --git a/src/adaptsim/__init__.py b/src/adaptsim/__init__.py new file mode 100644 index 0000000..1e58860 --- /dev/null +++ b/src/adaptsim/__init__.py @@ -0,0 +1,13 @@ +# -*- coding: utf-8 -*- +""" +this folder holds the adaptsim package +Created on Wed Jan 27 12:28:31 2023 +@author: janicweber +""" +from .constants import ALL_SIM_DICT, KEY_DICT_SIM, RCPARAMS +from .ast import MC_object +from .visualiser import * +from .visualiser_data import * + +__all__ = [''] + diff --git a/src/adaptsim/ast.py b/src/adaptsim/ast.py new file mode 100644 index 0000000..ad8b6a8 --- /dev/null +++ b/src/adaptsim/ast.py @@ -0,0 +1,209 @@ +# -*- coding: utf-8 -*- +import argparse +import json +import adaptfx as afx +import adaptsim as afs +import numpy as np +import sys +nme = __name__ + +class MC_object(): + """ + Reinforcement Learning class to check instructions + of calculation, invoke keys and define + calculation settings from file + """ + def __init__(self, model_filename): + model = afx.RL_object(model_filename) + with open(model.filename, 'r') as f: + read_in = f.read() + raw_simulation_dict = json.loads(read_in) + + try: + algorithm_simulation = raw_simulation_dict['algorithm_simulation'] + except KeyError as algo_err: + afx.aft_error(f'{algo_err} key missing in: "{model.filename}"', nme) + else: + afx.aft_message_info('algorithm simulation:', algorithm_simulation, nme, 1) + + try: # check if simulation_keys exists and is a dictionnary + raw_keys = raw_simulation_dict['keys_simulation'] + except KeyError: + afx.aft_error(f'"keys_simulation" is missing in : "{model.filename}"', nme) + else: + simulation_dict = afx.key_reader(afs.KEY_DICT_SIM, afs.ALL_SIM_DICT, raw_keys, 'sim') + afx.aft_message_dict('simulation', simulation_dict, nme) + + self.algorithm = model.algorithm + self.filename = model.filename + self.basename = model.basename + self.log = model.log + self.log_level = model.log_level + self.keys_model = model.keys + self.settings = model.settings + self.algorithm_simulation = algorithm_simulation + self.keys_simulation = afx.DotDict(simulation_dict) + + + def simulate(self): + plot_sets = afs.RCPARAMS + if self.settings.usetex: + # if CLI specifies global use of latex + plot_sets["text.usetex"] = True + elif self.keys_simulation.usetex == True: + plot_sets["text.usetex"] = True + plot_sets["figure.figsize"] = self.keys_simulation.figsize + plot_sets["font.size"] = self.keys_simulation.fontsize + n_frac = self.keys_model.number_of_fractions + # case listing for all options + if self.algorithm_simulation == 'histogram': + n_patients = self.keys_simulation.n_patients + mu = self.keys_simulation.fixed_mean_sample + std = self.keys_simulation.fixed_std_sample + plans = np.zeros(n_patients) + for i in range(n_patients): + self.keys_model.sparing_factors = list(np.random.normal(mu , std, n_frac + 1)) + output = afx.multiple(self.algorithm, self.keys_model, self.settings) + # output.oar_sum output.tumor_sum + n_frac_used = np.count_nonzero(~np.isnan(output.physical_doses)) + plans[i] = n_frac_used + end_plot = afs.plot_hist(plans, n_frac, plot_sets) + + if self.algorithm_simulation == 'NEW': + # Set parameters + u = 0 # mean of normal distribution + sigma = 1 # standard deviation of normal distribution + shape = 2 # shape parameter of gamma distribution + scale = 3 # scale parameter of gamma distribution + n = 6 # number of samples in each array + m = 10 # number of arrays to generate + + # Vectorized implementation + u_mean = np.random.normal(u, sigma, size=(m, 1)) + std = np.random.gamma(shape, scale, size=(m, 1)) + samples = np.random.normal(u_mean, std, size=(m, n)) + + for i in range(m): + self.keys_model.sparing_factors = samples[i] + output = afx.multiple(self.algorithm, self.keys_model, self.settings) + + elif self.algorithm_simulation == 'fraction': + # plot applied dose, sparing factor dependent on fraction + c_list = self.keys_simulation.c_list + n_c = len(c_list) + mu = self.keys_model.fixed_mean + std = self.keys_model.fixed_std + sf_list = self.keys_model.sparing_factors + c_dose_array = np.zeros((n_c, n_frac)) + oar_sum_array = np.zeros((n_c)) + for i, c in enumerate(self.keys_simulation.c_list): + self.keys_model.c = c + output = afx.multiple(self.algorithm, self.keys_model, self.settings) + c_dose_array[i] = output.tumor_doses + oar_sum_array[i] = output.oar_sum + self.c_dose_list = c_dose_array + end_plot = afs.plot_dose(self.c_dose_list, sf_list, n_frac, c_list, oar_sum_array, + mu, std, plot_sets) + + elif self.algorithm_simulation == 'single_state': + # plot policy, values or remaining number of fraction for one state + out = afx.multiple(self.algorithm, self.keys_model, self.settings) + if self.settings.plot_policy: + end_plot = afs.plot_val_single(out.policy.sf, out.policy.states, out.policy.val, + out.policy.fractions, self.keys_simulation.plot_index, r'Policy $\pi$ in BED$_{10}$', 'turbo', plot_sets) + if self.settings.plot_values: + end_plot = afs.plot_val_single(out.value.sf, out.value.states, out.value.val, + out.value.fractions, self.keys_simulation.plot_index, r'Value $v$', 'viridis', plot_sets) + if self.settings.plot_remains: + end_plot = afs.plot_val_single(out.remains.sf, out.remains.states, out.remains.val, + out.remains.fractions, self.keys_simulation.plot_index, r'Expected Remaining Number $\varepsilon$', 'plasma', plot_sets) + + elif self.algorithm_simulation == 'all_state': + # plot all policy, values or remaining number of fraction except for the last + out = afx.multiple(self.algorithm, self.keys_model, self.settings) + if self.settings.plot_policy: + end_plot = afs.plot_val_all(out.policy.sf, out.policy.states, out.policy.val, + out.policy.fractions, r'Policy $\pi$ in BED$_{10}$ [Gy]', 'turbo', plot_sets) + if self.settings.plot_values: + end_plot = afs.plot_val_all(out.value.sf, out.value.states, out.value.val, + out.value.fractions, r'Value $v$', 'viridis', plot_sets) + if self.settings.plot_remains: + end_plot = afs.plot_val_all(out.remains.sf, out.remains.states, out.remains.val, + out.remains.fractions, r'Expected Remaining Number $\varepsilon$', 'plasma', plot_sets) + + elif self.algorithm_simulation == 'single_distance': + selec = self.keys_simulation.data_selection + row_hue = self.keys_simulation.data_row_hue + plot_sets["axes.linewidth"] = 1.3 + data_sf = afs.data_reader(self.keys_simulation.data_filepath, selec[0], selec[1]) + end_plot = afs.plot_single(data_sf, 'Distance', 'sparing_factor', + row_hue, r'$w$ [cm]', r'$\delta$', True, 'Set2', plot_sets) + + elif self.algorithm_simulation == 'single_patient': + selec = self.keys_simulation.data_selection + row_hue = self.keys_simulation.data_row_hue + plot_sets["axes.linewidth"] = 1.3 + data_sf = afs.data_reader(self.keys_simulation.data_filepath, selec[0], selec[1], selec[2], selec[3]) + end_plot = afs.plot_single(data_sf, 'Patient', 'sparing_factor', + row_hue, 'Patient', r'$\delta$', False, 'Set2', plot_sets) + + elif self.algorithm_simulation == 'grid_distance': + selec = self.keys_simulation.data_selection + row_hue = self.keys_simulation.data_row_hue + plot_sets["axes.linewidth"] = 1.3 + data_sf = afs.data_reader(self.keys_simulation.data_filepath, selec[0], selec[1], selec[2], selec[3]) + end_plot = afs.plot_grid(data_sf, 'Distance', 'sparing_factor', + 'Patient', row_hue, r'$w$ [cm]', r'$\delta$', 'colorblind', plot_sets) + + elif self.algorithm_simulation == 'grid_fraction': + selec = self.keys_simulation.data_selection + row_hue = self.keys_simulation.data_row_hue + plot_sets["axes.linewidth"] = 1.2 + data_sf = afs.data_reader(self.keys_simulation.data_filepath, selec[0], selec[1], selec[2], selec[3]) + end_plot = afs.plot_twin_grid(data_sf, 'Fraction', 'sparing_factor', + 'Distance', 'Patient', row_hue, r'Fraction $t$', r'$\delta$', r'$w$ [cm]', 'colorblind', plot_sets) + elif end_plot == None: + # catch case where no algorithm simulation is specified + afx.aft_error(f'No such simulation: "{self.algorithm_simulation}"', nme) + + if self.keys_simulation.save: + afx.save_plot(self.basename, end_plot) + else: + afx.show_plot() + +def main(): + """ + CLI interface to invoke the RL class + """ + start = afx.timing() + parser = argparse.ArgumentParser( + description='Patient Monte Carlo Simulation to test adaptive fractionation' + ) + parser.add_argument( + '-f', + '--filenames', + metavar='', + help='input adaptive fractionation instruction filename(s)', + type=str, + nargs='*' + ) + parser.add_argument( + '-t', + '--usetex', + help='matplotlib usetex parameter flag', + action='store_true', + default=False + ) + # In case there is no input show help + args = parser.parse_args(args=None if sys.argv[1:] else ['--help']) + for single_file in args.filenames: + sim = MC_object(single_file) + sim.settings.usetex = args.usetex + afx.aft_message(f'start session for {single_file}', nme, 1) + sim.simulate() + afx.timing(start) + afx.aft_message(f'finish session for {single_file}', nme, 1) + + +if __name__ == '__main__': + main() diff --git a/src/adaptsim/constants.py b/src/adaptsim/constants.py new file mode 100644 index 0000000..faf1597 --- /dev/null +++ b/src/adaptsim/constants.py @@ -0,0 +1,48 @@ +from matplotlib import cycler + +ALL_SIM_DICT = {'n_patients': 0, # histogram + 'fixed_mean_sample': 0, # histogram + 'fixed_std_sample': 0, # histogram + 'c_list': 0, # fraction + 'plot_index': 1, # single state + 'data_filepath': './', # data + 'data_selection': [], # data + 'data_row_hue': 0, # data + 'figsize': [6,4], # settings + 'fontsize': 14, # settings + 'save': 0, # settings + 'usetex': 0 # settings + } + +KEY_DICT_SIM = {'sim': list(ALL_SIM_DICT)} + +plot_line_cycle = cycler('linestyle', ['-', '--', ':', '-.']) + +RCPARAMS = {'figure.figsize' : [6,4], + 'font.size' : 14, + 'font.family' : 'serif', + 'legend.handlelength' : 1.4, + 'legend.fontsize' : 11, + 'legend.title_fontsize' : 11, + 'text.usetex' : False, + 'axes.labelpad' : 10, + 'axes.linewidth' : 1.1, + 'axes.xmargin' : 0.05, + 'axes.prop_cycle' : plot_line_cycle, + 'axes.autolimit_mode' : 'data', + 'axes.titlelocation' : 'left', + 'xtick.major.size' : 7, + 'xtick.minor.size' : 3.5, + 'xtick.major.width' : 1.1, + 'xtick.minor.width' : 1.1, + 'xtick.major.pad' : 5, + 'xtick.minor.visible' : True, + 'ytick.major.size' : 7, + 'ytick.minor.size' : 3.5, + 'ytick.major.width' : 1.1, + 'ytick.minor.width' : 1.1, + 'ytick.major.pad' : 5, + 'ytick.minor.visible' : True, + 'lines.markersize' : 7, + 'lines.markerfacecolor' : 'none', + 'lines.markeredgewidth' : 0.8} \ No newline at end of file diff --git a/src/adaptsim/visualiser.py b/src/adaptsim/visualiser.py new file mode 100644 index 0000000..16d7318 --- /dev/null +++ b/src/adaptsim/visualiser.py @@ -0,0 +1,178 @@ +# -*- coding: utf-8 -*- +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import rcParams +from matplotlib.colors import Normalize as normalise +import matplotlib.cm as cm + +def plot_dose(data, sf_list, n_frac, c_list, oar_array, mean, std, plot_sets=None): + """ + creates a plot of applied dose and corresponding sparing factor + + Parameters + ---------- + data : array + 2d array + sf_list : array + 1d array with dimension n + n_frac : int + + Returns + ------- + ax : matplotlib.pyplot.axes + + """ + if plot_sets: + rcParams.update(plot_sets) + + x = np.arange(1, n_frac+1) + fig, ax = plt.subplots(1,1) + # plot the applied policy + N_text = r'{\mathrm{N}}' + for i, c_raw in enumerate(c_list): + oar_bed = np.round(oar_array[i], 1) + c = np.round(c_raw, 1) + ax.plot(x, data[i], label=rf'$C={c:.1f}$, $B^{N_text}_{n_frac}={oar_bed}$ Gy', + alpha=0.5, color='black') + # plot the sparing factors + ax2 = ax.twinx() + ax2.scatter(x, sf_list[1:], label=rf'$\delta_t \sim \mu={mean}, \sigma={std}$', + marker='^', color='black') + ax2.invert_yaxis() + ax2.set_ylabel(r'$\delta$') + ax.set_ylabel(r'BED$_{10}$ [Gy]') + ax.set_xlabel(r'Fraction $t$') + ax.set_xticks(range(min(x), max(x)+1)) + ax.tick_params(axis='x', which='minor', bottom=False) + lines, labels = ax.get_legend_handles_labels() + cross, clabels = ax2.get_legend_handles_labels() + ax2.legend(lines + cross, labels + clabels, loc=0) + fig.tight_layout() + + return fig + +def plot_hist(data, n_frac, plot_sets=None): + """ + creates a histogram plot of numbers of fractions used + + Parameters + ---------- + data : array + 1d array + n_frac : int + + Returns + ------- + ax : matplotlib.pyplot.axes + + """ + if plot_sets: + rcParams.update(plot_sets) + + x = np.arange(1, n_frac+1) + fig, ax = plt.subplots(1,1) + ax.hist(data, bins=x, alpha=0.4, align= 'left', + histtype= 'stepfilled', color='red') + ax.set_xticks(range(min(x), max(x)+2)) + ax.tick_params(axis='x', which='minor', bottom=False) + ax.set_ylabel(r'Number of Patients') + ax.set_xlabel(r'Fraction $t$') + fig.tight_layout() + + return fig + +def plot_val_single(sfs, states, data, fractions, index, label, colmap='turbo', plot_sets=None): + if plot_sets: + rcParams.update(plot_sets) + [n_grids, _, _] = data.shape + # search for optimal rectangular size of subplot grid + n_rows = n_columns = int(np.sqrt(n_grids)) + while n_rows * n_columns < n_grids: + if n_rows >= n_columns: + n_columns += 1 + else: + n_rows += 1 + # initiate plot and parameters + fig, ax = plt.subplots(1, 1) + x_min, x_max, y_min, y_max = sfs[0], sfs[-1], states[0], states[-1] + + # create shared colorbar + colmin, colmax = np.min(data), np.max(data) + normaliser = normalise(colmin, colmax) + colormap = cm.get_cmap(colmap) + im = cm.ScalarMappable(cmap=colormap, norm=normaliser) + + # loop through the axes + try: + axs = ax.ravel() + except: + # in case ax is a 1x1 subplot + axs = np.array([ax]) + + i = np.where(fractions==index)[0][0] + axs[0].imshow(data[i], interpolation=None, origin='upper', + norm=normaliser, cmap=colormap, aspect='auto', + extent=[x_min, x_max, y_min, y_max]) + T_text = r'{\mathrm{T}}' + axs[0].set_title(rf'$t = {fractions[i]}$') + axs[0].set_xlabel(r'$\delta$') + axs[0].set_ylabel(rf'$B^{T_text}$ [Gy]') + + fig.subplots_adjust(wspace=0.3, hspace=0.3, bottom=0.2, left=0.25, right=0.75) + fig.colorbar(mappable=im, ax=axs.tolist(), label=label) + + return fig + +def plot_val_all(sfs, states, data_full, fractions, label, colmap='turbo', plot_sets=None): + if plot_sets: + rcParams.update(plot_sets) + data = data_full[:-1] + [n_grids, _, _] = data.shape + # search for optimal rectangular size of subplot grid + n_rows = n_columns = int(np.sqrt(n_grids)) + while n_rows * n_columns < n_grids: + if n_rows >= n_columns: + n_columns += 1 + else: + n_rows += 1 + # initiate plot and parameters + fig, ax = plt.subplots(n_rows, n_columns) + x_min, x_max, y_min, y_max = sfs[0], sfs[-1], states[0], states[-1] + + # create shared colorbar + colmin, colmax = np.min(data_full), np.max(data_full) + normaliser = normalise(colmin, colmax) + colormap = cm.get_cmap(colmap) + im = cm.ScalarMappable(cmap=colormap, norm=normaliser) + + # loop through the axes + try: + axs = ax.ravel() + except: + # in case ax is a 1x1 subplot + axs = np.array([ax]) + + # turn off axes + for a in axs: + a.axis(False) + + for i, pol in enumerate(data): + axs[i].axis(True) + axs[i].imshow(pol, interpolation=None, origin='upper', + norm=normaliser, cmap=colormap, aspect='auto', + extent=[x_min, x_max, y_min, y_max]) + axs[i].set_title(rf'$t = {fractions[i]}$') + try: # get rid of inner axes values + axs[i].label_outer() + except: + pass + + T_text = r'{\mathrm{T}}' + fig.supxlabel(r'$\delta$') + fig.supylabel(rf'$B^{T_text}$ [Gy]') + + # fig.tight_layout() + fig.subplots_adjust(wspace=0.3, hspace=0.3, bottom=0.13, left=0.15,right=0.92) + fig.colorbar(mappable=im, ax=axs.tolist(), label=label) + + return fig \ No newline at end of file diff --git a/src/adaptsim/visualiser_data.py b/src/adaptsim/visualiser_data.py new file mode 100644 index 0000000..aed0d8a --- /dev/null +++ b/src/adaptsim/visualiser_data.py @@ -0,0 +1,86 @@ +import pandas as pd +from matplotlib import rcParams +import matplotlib.pyplot as plt +import seaborn as sns + +def data_reader(filename, key_1, entry_1, key_2=False, entry_2=False, sep=','): + df_load = pd.read_csv(filename, sep=sep) + df = df_load.loc[(df_load['Patient'] < 31)] + if key_2 and entry_2: + return df.loc[(df[key_1]==entry_1) & (df[key_2]==entry_2)] + else: + return df.loc[(df[key_1]==entry_1)] + +def plot_single(data, x, y, hue, x_label=None, y_label=None, minor_ticks=True, palette='Set2', plot_sets=None): + if plot_sets: + rcParams.update(plot_sets) + ax = sns.scatterplot(data=data, x=x, y=y, hue=hue, palette=palette, alpha=0.85) + if not minor_ticks: + ax.tick_params(axis='x', which='minor', bottom=False) + fig = ax.get_figure() + if x_label != None and y_label != None: + ax.set_xlabel(x_label) + ax.set_ylabel(y_label) + ax.legend_.set_title(hue) + fig.subplots_adjust(left=0.18, bottom=0.2, right=0.93, top=0.95) + return fig + +def plot_grid(data, x, y, hue, row, x_label=None, y_label=None, palette='Set2', plot_sets=None): + if plot_sets: + rcParams.update(plot_sets) + [length, height] = plot_sets["figure.figsize"] + aspect = length/height + + fig = sns.FacetGrid(data=data, hue=hue, row=row, height=height, aspect=aspect, palette=palette) + fig.map(sns.scatterplot, x, y, alpha=0.85) + fig.despine(top=False, right=False) + fig.add_legend() + fig.legend.set_title(hue) + if x_label != None and y_label != None: + fig.set_axis_labels(x_label, y_label) + fig.set_titles(row_template=row+": "+"{row_name}") + fig.tight_layout(w_pad=1) + return fig + +def plot_single_fraction(data, x, y, hue, x_label, y_label, y_twin_label, y_twin=None, palette='Set2', plot_sets=None): + if plot_sets: + rcParams.update(plot_sets) + fig, ax = plt.subplots() + ax = sns.scatterplot(data=data, x=x, y=y, hue=hue, palette=palette, marker='^') + if y_twin != None: + ax2 = ax.twinx() + sns.lineplot(data=data, x=x, y=y_twin, hue=hue, palette=palette, linestyle='-') + ax2.get_legend().remove() + ax2.invert_yaxis() + ax2.set_ylabel(y_twin) + ax2.set_ylabel(y_twin_label) + ax.set(xlabel=x_label, ylabel=y_label) + ax.legend_.set_title(hue) + return fig + +def plot_twin_grid(data, x, y, y_twin, hue, row, x_label=None, y_label=None, y_twin_label=None, palette='Set2', plot_sets=None): + if plot_sets: + rcParams.update(plot_sets) + [length, height] = plot_sets["figure.figsize"] + aspect = length/height + unique_palette = sns.color_palette(palette, n_colors=len(data[hue].unique())) + palette_dict = dict(zip(data[hue].unique(), unique_palette)) + + fig = sns.relplot(data=data, x=x, y=y, hue=hue, row=row, palette=palette_dict, style=hue, markers=['^'], height=height, aspect=aspect, alpha=0.85) + min_twin, max_twin = data[y_twin].min(), data[y_twin].max() + twin_diff = 0.05 * (max_twin - min_twin) + y_min, y_max = min_twin - twin_diff, max_twin + twin_diff + for index, [oar, ax] in enumerate(fig.axes_dict.items()): + ax2 = ax.twinx() + data_sub = data.loc[data[row]==oar] + sns.lineplot(data=data_sub, x=x, y=y_twin, hue=hue, palette=palette_dict, linestyle='-') + ax2.get_legend().remove() + ax2.invert_yaxis() + ax2.set_ylabel(y_twin_label) + ax2.set_ylim([y_min, y_max]) + fig.set_axis_labels(x_label, y_label) + fig.legend.set_title(hue) + fig.set_titles(row_template=row+": "+"{row_name}") + fig.tick_params(axis='x', which='minor', bottom=False) + fig.tight_layout(w_pad=1) + return fig \ No newline at end of file diff --git a/work/example_0.json b/work/example_0.json new file mode 100644 index 0000000..549dfb1 --- /dev/null +++ b/work/example_0.json @@ -0,0 +1,36 @@ +{ +"algorithm": "frac", +"level": 1, +"log": 0, +"keys": + { + "number_of_fractions": 6, + "fraction": 0, + "sparing_factors": [0.946, 0.906, 0.905, 0.796 , 0.868, + 0.916, 0.954], + "prob_update": 2, + "fixed_mean": 0.9, + "fixed_std": 0.04, + "shape": 10.61, + "scale": 0.003, + "shape_inv": 2.27, + "scale_inv": 0.002, + "accumulated_oar_dose": 0, + "accumulated_tumor_dose": 0, + "min_dose": 4, + "max_dose": 14, + "tumor_goal": 72, + "c": 1.5 + }, +"settings": + { + "dose_stepsize": 0.5, + "state_stepsize": 0.5, + "sf_stepsize": 0.005, + "plot_policy": 1, + "plot_remains": 0, + "plot_values": 0, + "plot_probability": 0, + "save_plot": 0 + } +} \ No newline at end of file diff --git a/work/testing.json b/work/example_1.json similarity index 69% rename from work/testing.json rename to work/example_1.json index a54198f..0a576fb 100644 --- a/work/testing.json +++ b/work/example_1.json @@ -6,8 +6,8 @@ { "number_of_fractions": 6, "fraction": 0, - "sparing_factors": [0.94657944, 0.90661544, 0.90529964, 0.7963188 , 0.8688744 , - 0.91669474, 0.95471553], + "sparing_factors": [0.946, 0.906, 0.905, 0.796 , 0.868, + 0.916, 0.954], "prob_update": 2, "fixed_mean": 0.9, "fixed_std": 0.04, @@ -20,16 +20,17 @@ "min_dose": 0, "max_dose": -1, "tumor_goal": 72, - "c": 0 + "c": 1.5 }, "settings": { "dose_stepsize": 0.5, "state_stepsize": 0.5, "sf_stepsize": 0.005, - "sf_low": 0.7, - "sf_high": 1.1, - "sf_prob_threshold": 0, - "plot_policy": 1 + "plot_policy": 1, + "plot_remains": 0, + "plot_values": 1, + "plot_probability": 1, + "save_plot": 0 } } \ No newline at end of file diff --git a/work/example_2.json b/work/example_2.json new file mode 100644 index 0000000..ace4c07 --- /dev/null +++ b/work/example_2.json @@ -0,0 +1,39 @@ +{ +"algorithm": "frac", +"level": 1, +"log": 0, +"keys": + { + "number_of_fractions": 6, + "fraction": 0, + "sparing_factors": [0.946, 0.906, 0.905, 0.796 , 0.868, + 0.916, 0.954], + "prob_update": 1, + "fixed_mean": 0.9, + "fixed_std": 0.04, + "shape": 10.61, + "scale": 0.003, + "shape_inv": 2.27, + "scale_inv": 0.002, + "accumulated_oar_dose": 0, + "accumulated_tumor_dose": 0, + "min_dose": 0, + "max_dose": -1, + "tumor_goal": 72, + "c": 1.5 + }, +"settings": + { + "dose_stepsize": 0.5, + "state_stepsize": 0.5, + "sf_stepsize": 0.005, + "sf_low": 0.3, + "sf_high": 1.4, + "sf_prob_threshold": 1e-3, + "plot_policy": 1, + "plot_remains": 1, + "plot_values": 0, + "plot_probability": 1, + "save_plot": 0 + } +} diff --git a/work/example_3.json b/work/example_3.json new file mode 100644 index 0000000..214f490 --- /dev/null +++ b/work/example_3.json @@ -0,0 +1,35 @@ +{ +"algorithm": "oar", +"level": 1, +"log": 0, +"keys": + { + "number_of_fractions": 6, + "fraction": 0, + "sparing_factors": [0.946, 0.906, 0.905, 0.796 , 0.868, + 0.916, 0.954], + "prob_update": 0, + "fixed_mean": 0.85, + "fixed_std": 0.1, + "shape": 10.61, + "scale": 0.003, + "shape_inv": 2.27, + "scale_inv": 0.002, + "accumulated_oar_dose": 0, + "accumulated_tumor_dose": 0, + "min_dose": 0, + "max_dose": -1, + "tumor_goal": 72 + }, +"settings": + { + "dose_stepsize": 0.5, + "state_stepsize": 0.5, + "sf_stepsize": 0.005, + "plot_policy": 1, + "plot_remains": 0, + "plot_values": 0, + "plot_probability": 1, + "save_plot": 0 + } +} \ No newline at end of file diff --git a/work/example_4.json b/work/example_4.json new file mode 100644 index 0000000..6dbcf11 --- /dev/null +++ b/work/example_4.json @@ -0,0 +1,35 @@ +{ +"algorithm": "oar", +"level": 1, +"log": 0, +"keys": + { + "number_of_fractions": 6, + "fraction": 0, + "sparing_factors": [0.946, 0.906, 0.905, 0.796 , 0.868, + 0.916, 0.954], + "prob_update": 0, + "fixed_mean": 0.85, + "fixed_std": 0.1, + "shape": 10.61, + "scale": 0.003, + "shape_inv": 2.27, + "scale_inv": 0.002, + "accumulated_oar_dose": 0, + "accumulated_tumor_dose": 0, + "min_dose": 0, + "max_dose": -1, + "tumor_goal": 72 + }, +"settings": + { + "dose_stepsize": 1, + "state_stepsize": 3, + "sf_stepsize": 0.01, + "plot_policy": 1, + "plot_remains": 0, + "plot_values": 0, + "plot_probability": 1, + "save_plot": 0 + } +} \ No newline at end of file diff --git a/work/example_5.json b/work/example_5.json new file mode 100644 index 0000000..9a5c6dc --- /dev/null +++ b/work/example_5.json @@ -0,0 +1,36 @@ +{ +"algorithm": "oar", +"level": 1, +"log": 0, +"keys": + { + "number_of_fractions": 6, + "fraction": 0, + "sparing_factors": [0.946, 0.906, 0.905, 0.796 , 0.868, + 0.916, 0.954], + "prob_update": 0, + "fixed_mean": 0.85, + "fixed_std": 0.1, + "shape": 10.61, + "scale": 0.003, + "shape_inv": 2.27, + "scale_inv": 0.002, + "accumulated_oar_dose": 0, + "accumulated_tumor_dose": 0, + "min_dose": 0, + "max_dose": -1, + "tumor_goal": 72 + }, +"settings": + { + "dose_stepsize": 3, + "state_stepsize": 3, + "sf_stepsize": 0.01, + "sf_prob_threshold": 1e-4, + "plot_policy": 1, + "plot_remains": 0, + "plot_values": 0, + "plot_probability": 1, + "save_plot": 0 + } +} \ No newline at end of file diff --git a/work/example_6.json b/work/example_6.json new file mode 100644 index 0000000..bae5237 --- /dev/null +++ b/work/example_6.json @@ -0,0 +1,36 @@ +{ +"algorithm": "oar_old", +"level": 1, +"log": 0, +"keys": + { + "number_of_fractions": 6, + "fraction": 0, + "sparing_factors": [0.946, 0.906, 0.905, 0.796 , 0.868, + 0.916, 0.954], + "prob_update": 0, + "fixed_mean": 0.85, + "fixed_std": 0.1, + "shape": 10.61, + "scale": 0.003, + "shape_inv": 2.27, + "scale_inv": 0.002, + "accumulated_oar_dose": 0, + "accumulated_tumor_dose": 0, + "min_dose": 0, + "max_dose": -1, + "tumor_goal": 72 + }, +"settings": + { + "dose_stepsize": 3, + "state_stepsize": 3, + "sf_stepsize": 0.01, + "sf_prob_threshold": 1e-4, + "plot_policy": 0, + "plot_remains": 0, + "plot_values": 0, + "plot_probability": 0, + "save_plot": 0 + } +} diff --git a/work/example_7.json b/work/example_7.json new file mode 100644 index 0000000..eb8dc69 --- /dev/null +++ b/work/example_7.json @@ -0,0 +1,38 @@ +{ +"algorithm": "tumor", +"level": 1, +"log": 0, +"keys": + { + "number_of_fractions": 6, + "fraction": 1, + "sparing_factors": [0.946, 0.906, 0.905, 0.796 , 0.868, + 0.916, 0.954], + "prob_update": 0, + "fixed_mean": 0.85, + "fixed_std": 0.1, + "shape": 10.61, + "scale": 0.003, + "shape_inv": 2.27, + "scale_inv": 0.002, + "accumulated_oar_dose": 0, + "accumulated_tumor_dose": 0, + "min_dose": 0, + "max_dose": -1, + "tumor_goal": 72, + "oar_limit": 90, + "c": 0 + }, +"settings": + { + "dose_stepsize": 15, + "state_stepsize": 15, + "sf_stepsize": 0.1, + "sf_prob_threshold": 1e-4, + "plot_policy": 1, + "plot_remains": 0, + "plot_values": 1, + "plot_probability": 0, + "save_plot": 0 + } +} diff --git a/work/frac_example.json b/work/frac_example.json deleted file mode 100644 index 949fd45..0000000 --- a/work/frac_example.json +++ /dev/null @@ -1,34 +0,0 @@ -{ -"algorithm": "frac", -"level": 1, -"log": 0, -"keys": - { - "number_of_fractions": 6, - "fraction": 0, - "sparing_factors": [0.98, 0.97, 0.8, 0.83, 0.8, 0.85, 0.94], - "prob_update": 0, - "fixed_mean": 0.9, - "fixed_std": 0.04, - "shape": "invalid", - "scale": "invalid", - "shape_inv": "invalid", - "scale_inv": "invalid", - "accumulated_oar_dose": 0, - "accumulated_tumor_dose": 0, - "min_dose": 0, - "max_dose": -1, - "tumor_goal": 72, - "c": 0 - }, -"settings": - { - "dose_stepsize": 0.5, - "state_stepsize": 0.5, - "sf_stepsize": 0.01, - "sf_low": 0.7, - "sf_high": 1.1, - "sf_prob_threshold": 0, - "plot_policy": 1 - } -} diff --git a/work/oar_example.json b/work/oar_example.json deleted file mode 100644 index 11790ed..0000000 --- a/work/oar_example.json +++ /dev/null @@ -1,34 +0,0 @@ -{ -"algorithm": "oar", -"level": 1, -"log": 0, -"keys": - { - "number_of_fractions": 6, - "fraction": 0, - "sparing_factors": [0.98, 0.97, 0.8, 0.83, 0.8, 0.85, 0.94], - "prob_update": 0, - "fixed_mean": 0.9, - "fixed_std": 0.04, - "shape": "invalid", - "scale": "invalid", - "shape_inv": "invalid", - "scale_inv": "invalid", - "accumulated_oar_dose": 0, - "accumulated_tumor_dose": 0, - "min_dose": 0, - "max_dose": -1, - "tumor_goal": 72, - "c": 0 - }, -"settings": - { - "dose_stepsize": 0.5, - "state_stepsize": 0.5, - "sf_stepsize": 0.01, - "sf_low": 0.7, - "sf_high": 1.1, - "sf_prob_threshold": 0, - "plot_policy": 1 - } -} \ No newline at end of file