<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Proof of concept
This is a proof of concept for the inclusion of positional uncertainty within the dwell time optimisation algorthim. The focus is on the optimisation method, and as such it is assumed that each seed radially deposits its energy purely based upon the inverse square law.
Two regions are defined, a central target region and a off centre avoid region.
Initialisation
Step1:  Create the calculation grid
Step2:  Define the target and avoid cubes
Step3:  Display the target and avoid cubes
Step4:  Create initial equidistant parrallel lines to represent catheters. Inbuild a slight random skew to emulate what occurs physically
Step5:  Display the lines overlayed
Step6:  Create a function to return x, y, z coords when a distance along a line is requested
Step7:  Pick dwell positons starting at a random position along the line
Step8:  Find all the dwell positions that are on the grid
Step9:  Plot the dwell positions
Step10:  Select the dwell positions that fall withing the target region. Only use these dwell positions.
Step11:  Calculate distance between each dwell and each grid point
Create an array containing the distance to each dwell position for each voxel and translate this to exposure at that voxel per unit dwell time at each dwell position. This is defined in this way so that a reasonable portion of the calculation of exposure can be done outside of the optimisation.
Step12:  Estimate uncertainty in distance
Estimate the uncertainty in distance given a shift along the catheter
Step13:  Create the PDF for exposure given a distance and uncertainty
Step14:  Create exposure calculation function
Step15:  Run a test of arbitrary dwell times
Step16:  Create function to display the results of the optimisation as it is being calculated.
Step17:  Create the optimisation cost function
The maximum exposure to a grid position not within the target is aimed to be less than 100. This cost function aims to achieve this.
Step18:  No grid point within the target is to be less than about 50. This cost function aims to achieve this.
Step19:  The avoid cost function aims to make no point within the avoid structure more than 45.
Step20:  Create the cost function to be used by the optimiser
Step21:  Create initial conditions
Step22:  Step noise
Step23:  Bounds
Step24:  Run the optimiser
Step25:  Presentation of results
Step26:  Display histogram of resulting dwell times
Step27:  Give overview of dwell times segmented by catheter
Step28:  Convert result into column vector for post analysis
Step29:  Create a custom resolution exposure calculation function for post analysis
Step30:  Analysis of results
Create analysis grid
Step31:  Define structures on new finer grid
Step32:  Make a function to plot DVHs given a reference
Step33:  Plotting of relevant DVHs
Step34:  Add a random shift to each catheter, see the effect it has on target DVH.
Step35:  With the shift the change in the DVH is the following. For the optimisation given here there is very little change in DVH.
Step36:  Overview of shift results
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import BoundaryNorm
from matplotlib.colorbar import ColorbarBase
%matplotlib inline
from utilities import BasinhoppingWrapper, create_green_cm
green_cm = create_green_cm()
Explanation: Proof of concept
This is a proof of concept for the inclusion of positional uncertainty within the dwell time optimisation algorthim. The focus is on the optimisation method, and as such it is assumed that each seed radially deposits its energy purely based upon the inverse square law.
Two regions are defined, a central target region and a off centre avoid region.
Initialisation
End of explanation
grid_spacing = 0.1 # Was 0.1, using 0.2 to speed up for testing
x_ = np.arange(-0.5, 0.5 + grid_spacing, grid_spacing)
y_ = np.arange(-0.5, 0.5 + grid_spacing, grid_spacing)
z_ = np.arange(-0.5, 0.5 + grid_spacing, grid_spacing)
x, y, z = np.meshgrid(x_, y_, z_)
x = np.ravel(x)
y = np.ravel(y)
z = np.ravel(z)
Explanation: Create the calculation grid
End of explanation
def target(x, y, z):
    return (
        (x < 0.45) & (x > -0.45) & 
        (y < 0.45) & (y > -0.45) & 
        (z < 0.45) & (z > -0.45))
def avoid(x, y, z):
    return (
        (x < 0.25) & (x > 0.05) & 
        (y < 0.15) & (y > -0.15) & 
        (z < 2) & (z > -2))
target_cube = target(x, y, z)
avoid_cube = avoid(x, y, z)
target_cube = target_cube & ~avoid_cube
Explanation: Define the target and avoid cubes
End of explanation
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(
    x[target_cube], y[target_cube], z[target_cube], 
    alpha=0.2, s=2, color='blue')
ax.scatter(
    x[avoid_cube], y[avoid_cube], z[avoid_cube], 
    alpha=0.2, s=2, color='red')
ax.set_xlim([np.min(x_), np.max(x_)])
ax.set_ylim([np.min(y_), np.max(y_)])
ax.set_zlim([np.min(z_), np.max(z_)])
ax.view_init(azim=-90, elev=90)
fig
Explanation: Display the target and avoid cubes
End of explanation
number_of_lines = 16 # Was 16, reduced to 9 for testing
line_start = np.meshgrid(
    [-0.39, -0.13, 0.13, 0.39],
    [-0.39, -0.13, 0.13, 0.39],
    [1])
line_finish = np.array([
    line_start[0] + np.random.normal(scale=0.02, size=[4, 4, 1]),
    line_start[1] + np.random.normal(scale=0.02, size=[4, 4, 1]),
    -line_start[2]])
line_start = np.array([np.ravel(mesh) for mesh in line_start])
line_finish = np.array([np.ravel(mesh) for mesh in line_finish])
Explanation: Create initial equidistant parrallel lines to represent catheters. Inbuild a slight random skew to emulate what occurs physically
End of explanation
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(
    x[target_cube], y[target_cube], z[target_cube], 
    alpha=0.2, s=2, color='blue')
ax.scatter(
    x[avoid_cube], y[avoid_cube], z[avoid_cube], 
    alpha=0.2, s=2, color='red')
ax.scatter(*line_start)
ax.scatter(*line_finish)
for i in range(len(line_start[0])):
    plt_coords = [
        [line_start[j][i], line_finish[j][i]]
        for j in range(len(line_start))]
    ax.plot(*plt_coords, color='black', alpha=0.5)
ax.set_xlim([-1, 1])
ax.set_ylim([-1, 1])
ax.set_zlim([-1, 1])
Explanation: Display the lines overlayed
End of explanation
diff = (line_finish - line_start)
line_length = np.sqrt(diff[0]**2 + diff[1]**2 + diff[2]**2)
def find_distance_coords(line_num=None, distance=None):
    relative_dist = distance / line_length[line_num]
    
    if (relative_dist > 1) | (relative_dist < 0):
        return np.array([np.nan]*3)
    
    x = (
        line_start[0][line_num] * (1 - relative_dist) + 
        line_finish[0][line_num] * relative_dist)
    
    y = (
        line_start[1][line_num] * (1 - relative_dist) + 
        line_finish[1][line_num] * relative_dist)
        
    z = (
        line_start[2][line_num] * (1 - relative_dist) + 
        line_finish[2][line_num] * relative_dist)
    
    coords = np.array([x, y, z])
    
    return coords
Explanation: Create a function to return x, y, z coords when a distance along a line is requested
End of explanation
dwell_spacing = 0.1 # Was 0.1, increased to 0.3 for testing
dwell_distances_from_initial = np.arange(0, 2, dwell_spacing)
number_of_dwells = len(dwell_distances_from_initial)
inital_dwell_position = np.random.uniform(
    low=0, high=dwell_spacing, size=number_of_lines)
inital_dwell_position
dwell_distances = np.reshape(inital_dwell_position, (-1, 1)) + np.reshape(dwell_distances_from_initial, (1, -1))
def find_dwell_coords(line_num=None, dwell_num=None):
    distance = dwell_distances[line_num, dwell_num]
    
    coords = find_distance_coords(
        line_num=line_num, distance=distance)
    
    return coords
Explanation: Pick dwell positons starting at a random position along the line
End of explanation
dwell_positions = np.array([
    [
        find_dwell_coords(
            line_num=line_num, dwell_num=dwell_num)
        for dwell_num in range(number_of_dwells)]
 for line_num in range(number_of_lines)])
Explanation: Find all the dwell positions that are on the grid
End of explanation
line_colours = np.random.uniform(size=(number_of_lines,3))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(
    x[target_cube], y[target_cube], z[target_cube], 
    alpha=0.2, s=2, color='blue')
ax.scatter(
    x[avoid_cube], y[avoid_cube], z[avoid_cube], 
    alpha=0.2, s=2, color='red')
for line_num in range(number_of_lines):
    ax.scatter(*np.transpose(dwell_positions[line_num]), 
               c=line_colours[line_num], alpha=0.7)
ax.set_xlim([-1, 1])
ax.set_ylim([-1, 1])
ax.set_zlim([-1, 1])
Explanation: Plot the dwell positions
End of explanation
dwell_positions_to_be_filtered = np.reshape(dwell_positions, (-1, 3))
distance_to_dwell_pos_initial = np.array([
    np.sqrt(
        (x[i] - dwell_positions_to_be_filtered[:,0])**2 + 
        (y[i] - dwell_positions_to_be_filtered[:,1])**2 + 
        (z[i] - dwell_positions_to_be_filtered[:,2])**2
    )
    for i in range(len(x))
])
closest_voxel_to_dwell = np.reshape(np.argmin(distance_to_dwell_pos_initial, axis=0), [1, -1])
target_cube_voxels = np.reshape(np.where(target_cube)[0], [-1, 1])
is_dwell_closest_to_target = np.any(closest_voxel_to_dwell == target_cube_voxels, axis=0)
relevant_dwell_positions = dwell_positions_to_be_filtered[is_dwell_closest_to_target]
line_number_index = np.reshape(np.arange(number_of_lines), (number_of_lines, 1))
line_number_index = np.ones([number_of_lines, number_of_dwells]) * line_number_index
line_number_index = np.ravel(line_number_index)[is_dwell_closest_to_target]
line_number_index = line_number_index.astype(int)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(
    x[target_cube], y[target_cube], z[target_cube], 
    alpha=0.2, s=2, color='blue')
ax.scatter(
    x[avoid_cube], y[avoid_cube], z[avoid_cube], 
    alpha=0.2, s=2, color='red')
for line_num in range(number_of_lines):
    ref = line_number_index == line_num
    ax.scatter(*np.transpose(relevant_dwell_positions[ref]), 
               c=line_colours[line_num], alpha=0.7)
ax.set_xlim([np.min(x_), np.max(x_)])
ax.set_ylim([np.min(y_), np.max(y_)])
ax.set_zlim([np.min(z_), np.max(z_)])
Explanation: Select the dwell positions that fall withing the target region. Only use these dwell positions.
End of explanation
distance_to_dwell_pos = np.array([
    np.sqrt(
        (x[i] - relevant_dwell_positions[:,0])**2 + 
        (y[i] - relevant_dwell_positions[:,1])**2 + 
        (z[i] - relevant_dwell_positions[:,2])**2
    )
    for i in range(len(x))
])
exposure_per_unit_time = 1 / distance_to_dwell_pos**2
Explanation: Calculate distance between each dwell and each grid point
Create an array containing the distance to each dwell position for each voxel and translate this to exposure at that voxel per unit dwell time at each dwell position. This is defined in this way so that a reasonable portion of the calculation of exposure can be done outside of the optimisation.
End of explanation
relevant_dwell_distances = np.ravel(dwell_distances)[is_dwell_closest_to_target]
shift_uncertainty = 0.05
shifted_up = np.array([
    find_distance_coords(
        distance=relevant_dwell_distances[i] + shift_uncertainty, 
        line_num=line_number_index[i])
    for i in range(len(relevant_dwell_distances))
])
shifted_down = np.array([
    find_distance_coords(
        distance=relevant_dwell_distances[i] - shift_uncertainty, 
        line_num=line_number_index[i])
    for i in range(len(relevant_dwell_distances))
])
distance_to_dwell_shift_up = np.array([
    np.sqrt(
        (x[i] - shifted_up[:,0])**2 + 
        (y[i] - shifted_up[:,1])**2 + 
        (z[i] - shifted_up[:,2])**2
    )
    for i in range(len(x))
])
distance_to_dwell_shift_down = np.array([
    np.sqrt(
        (x[i] - shifted_down[:,0])**2 + 
        (y[i] - shifted_down[:,1])**2 + 
        (z[i] - shifted_down[:,2])**2
    )
    for i in range(len(x))
])
dwell_distance_uncertainty = (
    np.abs(distance_to_dwell_shift_up - distance_to_dwell_pos) + 
    np.abs(distance_to_dwell_shift_down - distance_to_dwell_pos)
) / 2
Explanation: Estimate uncertainty in distance
Estimate the uncertainty in distance given a shift along the catheter
End of explanation
num_relevant_dwells = np.sum(is_dwell_closest_to_target)
def create_pdf_given_time(exposure, distance, uncertainty):   
    exposure = np.reshape(exposure, [1, 1, -1])
    distance = np.expand_dims(distance, axis=2)
    uncertainty = np.expand_dims(uncertainty, axis=2)    
    
    bot = 2*np.sqrt(2 * np.pi) * uncertainty * exposure**1.5
    
    def pdf(time):
        time = np.reshape(time, [1, -1, 1])
        top = (
            np.sqrt(time) * 
            np.exp(
                -(np.sqrt(time/exposure) - distance)**2 / 
                (2*uncertainty**2)))
    
        return top/bot
    
    return pdf
n = 2**9
exposure_spacing = np.linspace(0.01, 100, n)
pdf = create_pdf_given_time(exposure_spacing, distance_to_dwell_pos, dwell_distance_uncertainty)
time = np.ones(num_relevant_dwells)
check = pdf(time)
for i in range(num_relevant_dwells):
    plt.plot(exposure_spacing, check[6,i,:])
    plt.ylim([0, 0.4])
np.shape(check[0,:,:])
fft_of_rows = np.fft.rfft(check[2,:,:])
fft_of_convolution = fft_of_rows.prod(axis=0)
convolution_not_normalised = np.fft.irfft(fft_of_convolution)
# convolution = convolution_not_normalised[0:n] / np.trapz(convolution_not_normalised[0:n], x=exposure_spacing)
# plt.plot(exposure_spacing, convolution)
plt.plot(exposure_spacing, convolution_not_normalised)
Explanation: Create the PDF for exposure given a distance and uncertainty
End of explanation
def calculate_exposure(dwell_times):
    exposure = np.sum(dwell_times * exposure_per_unit_time, axis=1)
    
    return exposure
Explanation: Create exposure calculation function
End of explanation
num_relevant_dwells = len(relevant_dwell_positions)
random_pick = np.random.uniform(
    size=2, high=num_relevant_dwells, low=0).astype(int)
dwell_times = np.zeros([1, num_relevant_dwells])
dwell_times[0, random_pick] = 10
exposure = calculate_exposure(dwell_times)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
reference = exposure > 80
colour = exposure[reference]
colour[colour > 200] = 200
ax.scatter(
    x[reference], y[reference], z[reference], 
    alpha=0.2, s=20, c=colour, cmap=green_cm)
ax.set_xlim([np.min(x_), np.max(x_)])
ax.set_ylim([np.min(y_), np.max(y_)])
ax.set_zlim([np.min(z_), np.max(z_)])
Explanation: Run a test of arbitrary dwell times
End of explanation
def display_results(dwell_times):
    dwell_times = np.reshape(dwell_times, (1, num_relevant_dwells))
    exposure = calculate_exposure(dwell_times)
    
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    reference = exposure > 25
    colour = exposure[reference]
    colour[colour > 100] = 100
    
    small = exposure[reference] < 50
    large = ~small
    ax.scatter(
        x[reference][small], y[reference][small], z[reference][small], 
        alpha=0.2, s=3, c=colour[small], cmap=green_cm)
    
    ax.scatter(
        x[reference][large], y[reference][large], z[reference][large], 
        alpha=0.4, s=20, c=colour[large], cmap=green_cm)
    ax.set_xlim([np.min(x_), np.max(x_)])
    ax.set_ylim([np.min(y_), np.max(y_)])
    ax.set_zlim([np.min(z_), np.max(z_)])
    
    cost_function(dwell_times, debug=True)
    
    plt.show()
Explanation: Create function to display the results of the optimisation as it is being calculated.
End of explanation
def hot_exterior_cost_function(max_target_exterior):
    return ((max_target_exterior)/50)**4
testx = np.linspace(0, 120)
testy = hot_exterior_cost_function(testx)
plt.plot(testx, testy)
plt.ylim([0, 20])
Explanation: Create the optimisation cost function
The maximum exposure to a grid position not within the target is aimed to be less than 100. This cost function aims to achieve this.
End of explanation
def cold_target_cost_function(min_target):
    return ((min_target-80)/15)**4
testx = np.linspace(0, 120)
testy = cold_target_cost_function(testx)
plt.plot(testx, testy)
plt.ylim([0, 20])
Explanation: No grid point within the target is to be less than about 50. This cost function aims to achieve this.
End of explanation
def hot_avoid_cost_function(max_avoid):
    return ((max_avoid)/25)**4
testx = np.linspace(0, 120)
testy = hot_avoid_cost_function(testx)
plt.plot(testx, testy)
plt.ylim([0, 20])
Explanation: The avoid cost function aims to make no point within the avoid structure more than 45.
End of explanation
def cost_function(dwell_times, debug=False):
    dwell_times = np.reshape(dwell_times, (1, num_relevant_dwells))
    exposure = calculate_exposure(dwell_times)
    
    min_target = np.min(exposure[target_cube])
    max_avoid = np.max(exposure[avoid_cube])
    max_target_exterior = np.max(exposure[~target_cube])
    
    cold_target_cost = cold_target_cost_function(min_target)
    hot_exterior_cost = hot_exterior_cost_function(max_target_exterior)
    hot_avoid_cost = hot_avoid_cost_function(max_avoid)
    
    total_cost = hot_exterior_cost + cold_target_cost + hot_avoid_cost
    
    if debug:
        print("Minimum target = %0.4f, resulting cost = %0.4f" %
              (min_target, cold_target_cost))
        print("Maximum exterior = %0.4f, resulting cost = %0.4f" %
              (max_target_exterior, hot_exterior_cost))
        print("Maximum avoid = %0.4f, resulting cost = %0.4f" %
              (max_avoid, hot_avoid_cost))
        print("Total cost = %0.4f" % (total_cost))
        
    
    return total_cost
Explanation: Create the cost function to be used by the optimiser
End of explanation
num_relevant_dwells
initial_conditions = np.ones(num_relevant_dwells)*0.1
Explanation: Create initial conditions
End of explanation
step_noise = np.ones(num_relevant_dwells) * 0.3
Explanation: Step noise
End of explanation
bounds = ((0, None),)*num_relevant_dwells
Explanation: Bounds
End of explanation
optimisation = BasinhoppingWrapper(
    to_minimise=cost_function,
    initial=initial_conditions,
    step_noise=step_noise,
    basinhopping_confidence=5,
    optimiser_confidence=0.0001,
    n=5,
    debug=display_results,
    bounds=bounds
)
Explanation: Run the optimiser
End of explanation
display_results(optimisation.result)
Explanation: Presentation of results
End of explanation
plt.hist(optimisation.result)
Explanation: Display histogram of resulting dwell times
End of explanation
for line_num in range(number_of_lines):
    
    print("Line number %d" % (line_num))
    
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(
        x[target_cube], y[target_cube], z[target_cube], 
        alpha=0.2, s=2, color='blue')
    ax.scatter(
        x[avoid_cube], y[avoid_cube], z[avoid_cube], 
        alpha=0.2, s=2, color='red')
    
    ref = line_number_index == line_num
    ax.scatter(*np.transpose(relevant_dwell_positions[ref]), 
               c=line_colours[line_num], alpha=0.7)
    ax.set_xlim([np.min(x_), np.max(x_)])
    ax.set_ylim([np.min(y_), np.max(y_)])
    ax.set_zlim([np.min(z_), np.max(z_)])
    
    plt.show() 
    
    
    plt.scatter(np.arange(np.sum(ref)), optimisation.result[ref],
               c=line_colours[line_num])
    plt.xlabel("Dwell number")
    plt.ylabel("Units of dwell time")
    
    plt.ylim(bottom=0)
    plt.show()
Explanation: Give overview of dwell times segmented by catheter
End of explanation
dwell_times = np.reshape(optimisation.result, (1, num_relevant_dwells))
Explanation: Convert result into column vector for post analysis
End of explanation
def create_calculate_exposure(dwell_positions, 
                              x_new=None, y_new=None, z_new=None):
    x, y, z = np.meshgrid(
        x_new, y_new, z_new)
    x = np.ravel(x)
    y = np.ravel(y)
    z = np.ravel(z)
    distance_to_dwell_pos = np.array([
        np.sqrt(
            (x[i] - dwell_positions[:,0])**2 + 
            (y[i] - dwell_positions[:,1])**2 + 
            (z[i] - dwell_positions[:,2])**2
        )
        for i in range(len(x))
    ])
    exposure_per_unit_time = 1 / distance_to_dwell_pos**2
    
    def calculate_exposure(dwell_times, reshape=False):
        exposure = np.sum(dwell_times * exposure_per_unit_time, axis=1)
        if reshape:
            exposure = np.reshape(exposure, (
                len(x_new), len(y_new), len(z_new)))
        return exposure
        
    return calculate_exposure
dx = 0.01; dy = 0.01; dz = 0.1
contour_x_ = np.arange(-1, 1 + dx, dx)
contour_y_ = np.arange(-1, 1 + dy, dy)
contour_z_ = z_ #np.arange(-0.7, 0.7 + dz, dz)
contour_calculate_exposure = create_calculate_exposure(
    relevant_dwell_positions, x_new=contour_x_, y_new=contour_y_, z_new=contour_z_)
contour_exposure_init = contour_calculate_exposure(dwell_times, reshape=True)
vmin = 0
vmax = 200
num_contours = 41
fig = plt.figure(figsize=(10,0.3))
cbar_ax = fig.add_subplot(1,1,1)
bounds = np.linspace(vmin, vmax, num_contours)
norm = BoundaryNorm(bounds, green_cm.N)
cb = ColorbarBase(
    cbar_ax, norm=norm, cmap=green_cm,
    orientation='horizontal')
cb.ax.xaxis.set_ticks_position('top')
cb.ax.xaxis.set_label_position('top')
cb.set_clim([vmin, vmax])
contour_exposure = contour_exposure_init.copy()
contour_exposure[contour_exposure > vmax] = vmax
for i, z_value in enumerate(contour_z_):
    plt.figure(figsize=(6,6))
    
    c = plt.contourf(contour_x_, contour_y_,
        contour_exposure[:, :, i], num_contours,
        vmin=vmin, vmax=vmax, cmap=green_cm)
    
    reference = z[target_cube] == z_value
    x_target = x[target_cube][reference]
    y_target = y[target_cube][reference]
    plt.scatter(x_target, y_target, alpha=0.3, color='blue')
    
    reference = z[avoid_cube] == z_value
    x_avoid = x[avoid_cube][reference]
    y_avoid = y[avoid_cube][reference]
    plt.scatter(x_avoid, y_avoid, alpha=0.3, color='red')
    
    plt.title("Slice z = %0.1f" % (z_value))
    plt.xlim([np.min(contour_x_), np.max(contour_x_)])
    plt.ylim([np.min(contour_y_), np.max(contour_y_)])
    plt.show()
Explanation: Create a custom resolution exposure calculation function for post analysis
End of explanation
dx = 0.02
dy = 0.02
dz = 0.02
post_x_ = np.arange(-1, 1 + dx, dx)
post_y_ = np.arange(-1, 1 + dy, dy)
post_z_ = np.arange(-1, 1 + dz, dz)
post_x, post_y, post_z = np.meshgrid(
    post_x_, post_y_, post_z_)
post_x = np.ravel(post_x)
post_y = np.ravel(post_y)
post_z = np.ravel(post_z)
post_calculate_exposure = create_calculate_exposure(
    relevant_dwell_positions, x_new=post_x_, y_new=post_y_, z_new=post_z_)
post_exposure = post_calculate_exposure(dwell_times)
Explanation: Analysis of results
Create analysis grid
End of explanation
post_target_cube = target(post_x, post_y, post_z)
post_avoid_cube = avoid(post_x, post_y, post_z)
Explanation: Define structures on new finer grid
End of explanation
def create_dvh(reference, exposure, ax=None):
    if ax is None:
        fig = plt.figure()
        ax = fig.add_subplot(1,1,1)
    
    results = exposure[reference]
    results[results>150] = 150
    hist = np.histogram(results, 100)
    
    freq = hist[0]
    bin_edge = hist[1]
    bin_mid = (bin_edge[1::] + bin_edge[:-1:])/2
    
    cumulative = np.cumsum(freq[::-1])
    cumulative = cumulative[::-1]
    bin_mid = np.append([0], bin_mid)
    
    cumulative = np.append(cumulative[0], cumulative)
    percent_cumulative = cumulative / cumulative[0] * 100
    ax.plot(bin_mid, percent_cumulative)
Explanation: Make a function to plot DVHs given a reference
End of explanation
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
create_dvh(post_target_cube, post_exposure, ax=ax)
create_dvh(post_avoid_cube, post_exposure, ax=ax)
create_dvh(~post_target_cube, post_exposure, ax=ax)
plt.legend(["Target", "Avoid", "Patient"])
plt.xlim([0, 150])
Explanation: Plotting of relevant DVHs
End of explanation
shift_uncertainty = 0.05
line_displacement = np.random.normal(scale=shift_uncertainty, size=number_of_lines)
shifted_dwell_positions = np.array([
    find_distance_coords(
        distance=relevant_dwell_distances[i] + line_displacement[line_number_index[i]], 
        line_num=line_number_index[i])
    for i in range(len(relevant_dwell_distances))
])
shifted_calculate_exposure = create_calculate_exposure(
    shifted_dwell_positions, x_new=post_x_, y_new=post_y_, z_new=post_z_)
shifted_exposure = shifted_calculate_exposure(dwell_times)
Explanation: Add a random shift to each catheter, see the effect it has on target DVH.
End of explanation
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
create_dvh(post_target_cube, post_exposure, ax=ax)
create_dvh(post_target_cube, shifted_exposure, ax=ax)
plt.legend(["Original", "With position error"], loc='lower left')
plt.xlim([0, 110])
# plt.ylim([95, 100])
np.min(shifted_exposure[post_target_cube])
np.min(post_exposure[post_target_cube])
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
create_dvh(post_avoid_cube, post_exposure, ax=ax)
create_dvh(post_avoid_cube, shifted_exposure, ax=ax)
plt.legend(["Original", "With position error"], loc='upper right')
plt.xlim([0, 110])
Explanation: With the shift the change in the DVH is the following. For the optimisation given here there is very little change in DVH.
End of explanation
exposure = calculate_exposure(dwell_times)
np.min(exposure[target_cube])
def return_cold_spot_in_random_displace():
    line_displacement = np.random.normal(scale=shift_uncertainty, size=number_of_lines)
    shifted_dwell_positions = np.array([
        find_distance_coords(
            distance=relevant_dwell_distances[i] + line_displacement[line_number_index[i]], 
            line_num=line_number_index[i])
        for i in range(len(relevant_dwell_distances))
    ])
    shifted_calculate_exposure = create_calculate_exposure(
        shifted_dwell_positions, x_new=x_, y_new=y_, z_new=z_)
    shifted_exposure = shifted_calculate_exposure(dwell_times)
    return np.min(shifted_exposure[target_cube])
n = 1000
cold_spots = np.zeros(n)
for i in range(n):
    cold_spots[i] = return_cold_spot_in_random_displace()
plt.hist(cold_spots,20);
Explanation: Overview of shift results
End of explanation
<END_TASK>