# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Exercises to Bring it All Together

# %matplotlib inline
import matplotlib.pyplot as plt

# ## Exercise 0 -- Slicing and Dicing
#
# **This is a warm-up exercise; if you feel confident of your NumPy & SciPy skills, you can jump straight to Exercise 1**
#
# In the film industry, it is often necessary to impose actors on top of a rendered background.  To do that, the actors are filmed on a "green screen".  Here's an example shot (``images/greenscreen.jpg``):
#
# <img src="images/greenscreen.jpg" width="300px"/>
#
# Say we'd like to help these folks travel into a forest (``images/forest.jpg``):
#
# <img src="images/forest.jpg" width="300px"/>
#
# Can you transplant the foreground of the greenscreen onto the backdrop of the forest?

import sys
# !{sys.executable} -m pip install imageio

# +
import imageio

forest = imageio.imread('images/forest.jpg')
people = imageio.imread('images/greenscreen.jpg')
# -

f, (ax0, ax1) = plt.subplots(1, 2, figsize=(15, 10))
ax0.imshow(forest)
ax1.imshow(people);

# ## Solution for Exercise 0:

# the two images have different shapes:
forest.shape, people.shape

# resize the larger forest image to match the people image 
forest_crop = forest[:375, :500, :]

# Let's examine the red, green, and blue channels independently (note: each is a gray-level image).

# +
red = people[..., 0]
green = people[..., 1]
blue = people[..., 2]

f, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(15, 10))
ax0.imshow(red, cmap='gray')
ax1.imshow(green, cmap='gray')
ax2.imshow(blue, cmap='gray');
# -

# Let's examine the value distributions (histogram) to see where we can threshold the images to segment the foreground from the background.

plt.hist(red.ravel(), bins='auto');

plt.hist(green.ravel(), bins='auto');

plt.hist(blue.ravel(), bins='auto');

# +
import scipy
import scipy.ndimage

# set up matplotlib figure and plotting axis
f, (ax0, ax1) = plt.subplots(1, 2, figsize=(15, 10))

# generate a mask below a large value of green in the rgb data structure (second column)
mask = (green < 190) & (red > 5)

# fill in some of the selection holes in the mask using SciPy
mask = scipy.ndimage.morphology.binary_fill_holes(mask)

# assign the non-green region of people image to the same pixel regions in the forest image
forest_w_ppl = forest_crop.copy()
forest_w_ppl[mask] = people[mask]

ax0.imshow(mask, cmap='gray')
ax1.imshow(forest_w_ppl);

# there's still a bit of green leeking through, but that's not a bad first pass
# -

# ## Exercise 1 -- Demonstrate Diagram Duality in 2D
#
# For the random set of 2D input points selected below, demonstrate visually that the Delaunay triangulation and Voronoi diagram are duals. Specifically, you should be able to demonstrate that the circumcenters of the circumcircles of the Delaunay triangulation can be connected to form the edges of the Voronoi diagram.
#
# #### Delaunay triangulation / Voronoi tesselation (different data sets)
#
# <img src="images/delaunay.png" width="30%" style="float: left"/>
# <img src="images/Euclidean_Voronoi_diagram.svg.png" width="30%" style="float: left"/>
#
# <div style="clear: both"></div>
#
# (Voronoi tesselation by <NAME>, CC BY-SA 4.0, https://commons.wikimedia.org/wiki/File:Euclidean_Voronoi_diagram.svg)
#
# The initial setup is below.

# +
import numpy as np
import scipy
import scipy.spatial
import matplotlib
# %matplotlib inline
import matplotlib.pyplot as plt

np.random.seed(123)
generators = np.random.randn(15, 2)
generators.shape
# -

# Step 1: Plot the generators (input points)



# Step 2: Plot the Voronoi diagram along with its generators.



# Step 3: Plot the Delaunay triangulation alongside the Voronoi diagram and generators.



# Step 4: Plot the circumcircles and their centers over top of the Delaunay Triangulation & Voronoi diagram + generators.
#
# **Hint:** Plotting circles in Matplotlib:

# +
circle = plt.Circle((2, 5), 0.5, facecolor='None', edgecolor='red')

fig, ax = plt.subplots()

ax.add_artist(circle)
ax.scatter([1,2,3], [4,5,6])
ax.axis('equal');
# -

# Step 5: Connect the circumenters of circumcircles of neighboring triangles (hint: you may need to check the SciPy documentation for `Delaunay`--try looking for information about neighboring triangles)

# ## Solution for Exercise 1:

# +
# set up the matplotlib figure / axis objects as usual
fig_exercise_1 = plt.figure()
ax_exercise_1 = fig_exercise_1.add_subplot(111)

# scatter the generators (input points) on the plot
ax_exercise_1.scatter(generators[...,0],
                      generators[...,1],
                      c='black') # note that these generators actually get concealed by re-plotting in delaunay_plot_2d below

# overlay the Voronoi diagram using SciPy
vor = scipy.spatial.Voronoi(generators)
scipy.spatial.voronoi_plot_2d(vor=vor,
                              ax=ax_exercise_1, # let SciPy use the axis object that we've already intialized with generators present
                              show_points=False) # we've already plotted the generators above 

# likewise for adding in the Delaunay triangles
tri = scipy.spatial.Delaunay(generators)
scipy.spatial.delaunay_plot_2d(tri=tri,
                               ax=ax_exercise_1)

# write a function for calculating the circumcircles of the triangles (this is admittedly rather difficult)
import math

def calc_circumcircle(array_triangle_vertices):
    '''The input array of triangle vertices should have shape (N_triangles, 3, 2) and this function should use the vertex information to calculate
       and return numpy arrays containing the circumcenters and circumradii for the circumcircles of each triangle.'''
    
    def sum_squares_vertex(coordinates):
        square_coordinates = [coordinate**2 for coordinate in coordinates]
        sum_squares = sum(square_coordinates)
        return sum_squares
    
    def calc_sum_squares_each_vertex(vertex_1_coords,vertex_2_coords,vertex_3_coords):
        vertex_1_sum_squares = sum_squares_vertex(vertex_1_coords)
        vertex_2_sum_squares = sum_squares_vertex(vertex_2_coords)
        vertex_3_sum_squares = sum_squares_vertex(vertex_3_coords)
        return (vertex_1_sum_squares, vertex_2_sum_squares, vertex_3_sum_squares)
    
    list_circumcenters = [] # an ordered list of circumcenter coords
    list_circumradii = [] # a correspondingly ordered list of circumradii
    
    for triangle in array_triangle_vertices:
        vertex_1_sum_squares, vertex_2_sum_squares, vertex_3_sum_squares = calc_sum_squares_each_vertex(triangle[0,...],triangle[1,...],triangle[2,...])
        x_determinant = np.array([[vertex_1_sum_squares, triangle[0,1], 1],
                                  [vertex_2_sum_squares, triangle[1,1], 1],
                                  [vertex_3_sum_squares, triangle[2,1], 1]])
        y_determinant = np.array([[vertex_1_sum_squares, triangle[0,0], 1],
                                  [vertex_2_sum_squares, triangle[1,0], 1],
                                  [vertex_3_sum_squares, triangle[2,0], 1]])
        a_determinant = np.array([[triangle[0,0], triangle[0,1], 1],
                                  [triangle[1,0], triangle[1,1], 1],
                                  [triangle[2,0], triangle[2,1], 1]])
        c_determinant = np.array([[vertex_1_sum_squares,triangle[0,0],triangle[0,1]],
                                  [vertex_2_sum_squares,triangle[1,0],triangle[1,1]],
                                  [vertex_3_sum_squares,triangle[2,0],triangle[2,1]]])
        
        denominator = 2.0 * np.linalg.det(a_determinant)
        circumcenter_x_coordinate = np.linalg.det(x_determinant) / denominator
        circumcenter_y_coordinate = -1.0 * (np.linalg.det(y_determinant)/ denominator)
        
        #adjuting for the negative coefficient of c:
        circumradius = math.sqrt(np.linalg.det(x_determinant)**2 + np.linalg.det(y_determinant)**2 + (4 * np.linalg.det(a_determinant) * np.linalg.det(c_determinant))) / (2 * abs(np.linalg.det(a_determinant)))
        #append the values of interest to their respective lists:
        list_circumcenters.append([circumcenter_x_coordinate,circumcenter_y_coordinate])
        list_circumradii.append([circumradius])
    #return a tuple of numpy arrays:
    return (np.array(list_circumcenters), np.array(list_circumradii))

#deal with calculating and plotting the circumcircles
from matplotlib.collections import PatchCollection
circumcenters, circumradii = calc_circumcircle(tri.points[tri.simplices])
patches = []

for circumcenter_coordinates, circumradius in zip(circumcenters,circumradii):
    patches.append(matplotlib.patches.Circle((circumcenter_coordinates[0], circumcenter_coordinates[1]), circumradius[0], fill=False, alpha=1.0, fc='none', ec='none',))
p = PatchCollection(patches, alpha=0.1, match_original = True)
ax_exercise_1.add_collection(p)

# use an equal aspect ratio (circles may appear stretched otherwise)
ax_exercise_1.set_aspect('equal')

# adjust the size of the figure object
fig_exercise_1.set_size_inches(10, 10)
# -

# It should now be clear without any additional effort that the Voronoi vertices are actually
# the circumcenters of the circumcircles of the Delaunay triangles, and that the connection
# of these circumcenters forms the Voronoi edges.

# ## Exercise 2 -- Algorithm Time Complexity
#
# A common concern in computational geometry is the time complexity of an algorithm (the growth of code execution time as a function of algorithm input size).
#
# Compare the execution times of `numpy.sum()` and `scipy.spatial.Voronoi()` for the same input data sets, collecting the results in NumPy arrays and reporting average and standard deviation values of the execution times over a set of trials with increasingly large data sets.

# ## Solution for Exercise 2

# Step 1: Initialize some NumPy arrays for storing the results of the timings for sum and Voronoi

# +
timing_string = 'time_values'
dimensions = 2 # use 2D data
list_sizes = [10, 100, 500, 10**3, 10**4, 10**5] # use a range of data input sizes
num_trials = 5

# using a dictionary to organize / categorize the NumPy arrays is helpful here
dict_timing_results = {'sum':{}, 'vor':{}}

# build up the dictionary data structure
for benchmark_name in dict_timing_results.keys():
    for size in list_sizes:
        benchmark_string = timing_string + str(size)
        dict_timing_results[benchmark_name][benchmark_string] = np.zeros(num_trials)
    
# visually confirm a sensible dict & NumPy array data structure    
dict_timing_results
# -

# Step 2: Write a for loop that runs the timing trials and stores the timing results in the NumPy arrays that are appropriately named for the data inputs being stored within them. This will also involve generating matching data sets that work with both sum() and Voronoi()

# +
import time # useful for timing

for trial in range(num_trials):
    for function, bench_name in zip([np.sum, scipy.spatial.Voronoi], ['sum', 'vor']):
        for size in list_sizes:
            # pin down the random number seed for consistency
            # each algorithm will get the same data structure & values
            np.random.seed(123) # any integer number should work here
            input_data = np.random.randn(size, dimensions)
            
            # only time the operation of interest, in seconds
            start = time.time()
            # the actual benchmarked function call for the current data input size:
            function(input_data)
            end = time.time()
            
            # store the timing results in the dictionary object within the appropriate array index
            dict_timing_results[bench_name]['time_values' + str(size)][trial] = end - start
            
# -

# Step 3: Calculate the average and standard deviation of the results and compare the increases in execution time for sum() and Voronoi() as data input size increases in a plot. Include the standard deviation "error bars" in your plot. Hint: matplotlib has an `errorbar` plot function, for example.

# +
# iterate through the dictionary entries and calculate / plot / label accordingly

# first calculate average / std. dev values and store in the dict
for bench_name in dict_timing_results.keys():
    for size in list_sizes:
        key = 'time_values' + str(size)
        timings = dict_timing_results[bench_name][key]
        dict_timing_results[bench_name][str(size) + 'avg'] = np.average(timings)
        dict_timing_results[bench_name][str(size) + 'std'] = np.std(timings)
    
# then plot the avg / std. dev. results

fig_exercise_2 = plt.figure()
ax_exercise_2 = fig_exercise_2.add_subplot(111)

for bench_name, color in zip(dict_timing_results.keys(), ['orange', 'green']):
    for size in list_sizes:
        avg_val = dict_timing_results[bench_name][str(size) + 'avg']
        std_val = dict_timing_results[bench_name][str(size) + 'std']
        if size != 10: # avoid duplicate legend entries
            ax_exercise_2.errorbar(size, avg_val, yerr=std_val, c=color, alpha=0.5, fmt='o')
        else: # add the legend initially
            ax_exercise_2.errorbar(size, avg_val, yerr=std_val, c=color, alpha=0.5, label=bench_name, fmt='o')

# use log scale on x, y axes for clarity
ax_exercise_2.set_xscale("log")
ax_exercise_2.set_yscale("log")

# adjust figure size
fig_exercise_2.set_size_inches(10, 10)

# set various labels / axis limits
ax_exercise_2.set_title("Timing 5 trials of np.sum and scipy.spatial.Voronoi\n using the same input data", fontsize=25)
ax_exercise_2.set_xlabel("Log input size", fontsize=24)
ax_exercise_2.set_ylabel("Log execution time (sec)", fontsize=25)
ax_exercise_2.set_ylim(1e-6, 1)

# add a legend
ax_exercise_2.legend(fontsize=24, loc=2)
    

# -

# If we wanted to formally assign a time complexity to the algorithms we would typically assess this directly from the source code, asking questions such as "how many times do we need to pass through all the points in this algorithm?". 
#
# However, we could also perform an empirical estimate
# by using i.e., `scipy.optimize.curve_fit` with linear, loglinear, and quadratic functions to see which one produces the most suitable fit to the observed timings data.
#
# Informally, we might expect a linear time complexity for sum and a loglinear time complexity for Voronoi. Certainly it seems to be the case that Voronoi has a larger time complexity
# than sum here since the time consumed relative to sum increases as input size increases.

# ## Course feedback
#
# On a blank piece of paper, please write **(a)** one thing you enjoyed about the course and **(b)** one thing that we can improve.
#
# Thank you!

StarcoderdataJupyterScriptsDedupFiltered