Homework 1

Homework 1¶

CS328 — Numerical Methods for Visual Computing and Machine Learning¶

Out on Tuesday 27/09/2022, due on Friday 14/10/2022.

This notebook contains literate code, i.e. brief fragments of Python surrounded by descriptive text. Please complete/extend this notebook for your homework submission:

Begin your response to each problem with a ## Solution markdown cell.
In addition to your code, please also provide a short description of what your solution is doing and how it works, either by adding comments or in an extra markdown cell.
Before handing in, please make sure that your notebook runs from top to bottom after selecting "Kernel->Restart & Run All" without causing any errors. To simplify the grading process, please do not clear the generated output.

Make sure to use the reference Python distribution so that project files can be opened by the TAs. In this course, we use Anaconda, specifically the version based on Python 3.8.

Homework assignments in CS328 count towards your final grade and must therefore be done individually.

Prelude¶

The following fragment imports NumPy and Matplotlib and configures the latter to produce nice graphics on modern high-resolution screens. The import statements at the end establish a shorthand notation for the most common integer and floating point formats.

In [ ]:

%matplotlib inline
%config InlineBackend.figure_format='retina'

import numpy as np
from matplotlib import pyplot as plt

from numpy import uint16  as u16
from numpy import uint32  as u32
from numpy import uint64  as u64
from numpy import float16 as f16
from numpy import float32 as f32
from numpy import float64 as f64

Two more definitions: the helper functions f2i and i2f below reinterpret floating point values as an integers and vice versa. We'll use these in Problem 3 to access the bit-level representation of an IEEE 754 floating point value.

In [ ]:

def f2i(value):
    ''' Reinterpret floating point value as integer '''
    return value.view('u%i' % value.itemsize)

def i2f(value):
    ''' Reinterpret integer as floating point value '''
    return value.view('f%i' % value.itemsize)

Problem 1: Images, slicing, and color balancing (15 pts)¶

The following two lines of code download a portrait photograph from the Wikipedia article on color balance and store it in a single precision NumPy array named image, whose pixels are in the range [0.0, 1.0].

In [ ]:

import PIL, requests
image = np.array(PIL.Image.open(requests.get('https://goo.gl/VcYNGN', stream=True, headers = {'User-agent': 'Mozilla/5.0'}).raw), dtype = np.float32) / 255

Try plotting this image using plt.imshow.

In [ ]:

You will notice that the colors are distorted by a strong blue tint. Such intense color casts can lead to unnatural appearance and are usually not desired in portrait photography. In such cases, color balancing can be used to post-process the colors of an existing image to e.g. ensure that white objects indeed appear white in the image. One of the simplest kinds of color balancing is known as white balancing. Here, we select a region of the image that is known to contain a neutral color, and we then determine the reciprocal of the associated red, green, and blue color values. Afterwards, the color channels of the entire image are scaled by these reciprocals, ensuring that the selected region becomes neutral (i.e. it has similar amounts of red, green, and blue).

Note that the woman holds a color checker in her hands. A color checker consists of an arrangement of colored patches with known color values, and the last row usually contains neutral colors. We will use the second neutral patch on the bottom left to white-balance this image.

TODO: Use NumPy slicing operations to crop out the second neutral patch in the last row (counting from the left side) and plot it using plt.imshow. The cropped region should be as large as possible without including the black frame or other patches.

In [ ]:

TODO: Compute the mean of the R, G, and B color channels over the patch. You should be able to accomplish this with a single function call to np.mean. Print the computed mean using print().

In [ ]:

TODO: Now, use broadcasting to multiply the image R, G, B channels by the reciprocals (inverses) of the values computed above and visualize the result using plt.imshow. Note: You may need to scale the image by a small amount to prevent pixels overflowing to a value greater than 1.0, which will lead to strange color artifacts. (Our reference implementation uses a value of 0.85.)

In [ ]:

Problem 2: Triple Birthday Paradox (25 pts)¶

Let us again consider the Birthday Paradox that was introduced during the first week's exercise session. In this task, we will study a slightly modified problem that can be summarized with the following question:

What is the probability that at least three people in a random set of $n$ people share the same birthday?

In the following, you will develop a simple numerical Python program to solve this problem and then make use of vectorization with NumPy arrays to speed it up. For this reason, you should use the "magic" Jupyter function %%time to profile the execution time of your implementation at every step.

Problem 2.1: Basic Python Implementation (10 pts)¶

TODO: Write a basic Python program^[1] that estimates of the "3-birthdays probability" based on random sampling.

Feel free to reuse code from the first weeks notebook to get some inspiration on how you would do this. The structure should stay exactly the same:

Choose $n$ random birthdays (i.e. integers)
Check if three of the birthdays are the same
Repeat this process $K$ times and keep track of how many times it was true
Return the ratio between "number successful events" and "total number events" ($K$)

When you chose $K$ large enough, a simulation like this should give you a good estimate of the true probability.

Remember that you can generate random integers by using randint from the random module:

import random as rnd
a = rnd.randint(0, 10) # Random integer between 0 and 10 (including)

[1] Do not use any of the NumPy, SciPy, ... libraries in this part of the exercise—only pure Python is allowed.

In [ ]:

## TODO... Fill in the inner loop of the following function

# For n people, compute an approximation of the probability that at
# least three people share the same birthday, using K iterations
def birthday_paradox_basic(n, K):
    counter = 0 # Keep track of how often the statement is true
    
    # Perform K iterations of the same experiment.
    for it in range(K):
        success = True # Replace!
        
        if success:
            counter += 1  
    
    # Return average probability
    return counter / K

To test if your program is operating correctly, try running it for $n=88$ people using a very high value of $K$ (e.g. 100000). This should give you (in expectation) a probability slightly above $0.5$:

In [ ]:

TODO: Finally, profile the running time of the implementation and report how much time it took for $K=100000$. You can do this simply by putting the "magic command" %%time at the top of the cell that should be timed:

In [ ]:

%%time
# Your code below ..

Problem 2.2: Vectorized implementation using NumPy (15 pts)¶

As you likely noticed, an implementation in pure Python is rather slow! As a frame of reference: an implementation with list comprehensions (as was done during the exercise session) can easily take half a minute to complete for large values of $K$ like above. Creating lists of random integers is a particularly slow operation that can be considerably accelerated using vectorized NumPy functions that process entire NumPy arrays at a once.

TODO: Write another function that computes the same result, but this time, replace the inner part of the "for it in range(K)" loop with a combination of NumPy array functions.^[2]

Here are a number of NumPy functions that you may find helpful:

np.random.randint
np.bincount
np.unique
np.min, np.max
np.arange
np.any

You can read about these in the NumPy Documentation.

Hint: The output range of np.random.randint differs slightly from random.randint provided by the standard library that you used above. Make sure not to introduce any subtle off-by-one errors!

[2] In this part of the exercise, do not use Python lists and list comprehensions, and do not explicitly iterate over the days of the year. Do not use np.vectorize, np.apply_along_axis (or similar constructs) that merely emulate vectorization using a slow Python for loop internally.

In [ ]:

## TODO... Fill in the inner loop of the following function

def birthday_paradox_numpy(n, K):
    counter = 0 # Keep track of how often the statement is true
    
    # Perform K iterations of the same experiment.
    for it in range(K):
        success = True # Replace!
        
        if success:
            counter += 1  
    
    # Return average probability
    return counter / K

TODO: Profile the running time of the implementation vectorized implementation for $K=100000$ like above and report the (approximate) speedup in a markdown cell or comment.

In [ ]:

Problem 2.3: Optional Exercise (0 points, not graded)¶

Although the vectorized implementation is a tremendous improvement, it still contains a loop over a potentially large number of virtual experiments. We can do better!

Implement a version that fully relies on vectorized NumPy functions and contains no more loops. As above, do not use np.vectorize, np.apply_along_axis (or similar constructs) that merely emulate vectorization using a slow Python for loop internally.
Again, profile your code and report the speedup compared to both the "pure Python" and previous NumPy implementation. We expect to see at least some additional speedup here. (As a frame of reference, our solution to this problem achieves ~4x speedup compared to problem 2.2.)

In [ ]:

Problem 3: Measuring error using ULPs (25 pts)¶

A special property of the IEEE754 floating point specification is that contiguous ranges of floating point numbers are also contiguous when re-interpreted as unsiged integers, e.g. by after applying the function f2i defined above. See the next figure for an illustration of this.

Use this property along with the functions f2i and i2f and the aliases f16 etc. defined above to answer the following questions:

How many floating point values are located between the values 1 and 2, including the endpoints? How about 1001 and 1002 ? Give answers for 16, 32, and 64 bit precision variants.
How large is 1 ULP for the value $\pi$ expressed in 64 bit arithmetic? Remember that one ULP was defined as the jump that occurred when changing the last mantissa bit from a 1 to a 0 (or vice versa).
Suppose that we compute the surface area of the earth from its radius of $\approx 6.353\cdot 10^6m$ using the expression $A=4\pi r^2$. Assuming that the earth is perfectly spherical, and that the discretization of $\pi$ is the only source of error in this computation, bound the absolute error in square meters.

In [ ]:

Problem 4: Series approximations (25 pts)¶

In Lecture 1, we saw how various numerical issues could occur due to the floating point number representation.
Consider the power series representation of the sine function:

$$ \sin(x)=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\cdots $$

Hint: the first 20 coefficients of the sine power series representation look like this:

Evaluate the power series at (x=20) using the first 100 terms and compute the absolute and relative error (you can assume that np.sin is accurate). Hint: the factorial function is provided in scipy.special.
Express the error as an integer number of ULPs (see also Problem 3) -- in other words, how many floating point numbers are between the true and the approximate answer?
Now, evaluate the exponential function at x=30 using the same approach: using its corresponding power series, and assuming that np.exp is accurate. Specify the absolute and relative error. What do you observe compared to the sine power series? Why is this the case?
List three different kinds of numerical problems that can arise when evaluating the above two power series with very many terms (e.g. thousands), and what parts of the expression specifically cause them.

In [ ]: