# Longest Common Subsequence

## The Longest Common Subsequence Problem

A subsequence of a sequence is the original sequence with zero or more elements removed from it. A common subsequence of two sequences is a subsequence that belongs to both sequences. With the Longest Common Subsequence (LCS) problem you are given two sequences - $X = \langle x_1, x_2,\ldots,x_m\rangle$ and $Y = \langle y_1, y_2,\ldots,y_n\rangle$ and you want to find a common subsequence of $X$ and $Y$ that has the maximum length (there might be more than one common subsequence that have the maximum length).

Note: The elements of the subsequence have to be in both sequences in the same relative order, but don't need to be consecutive. "AK" is a common subsequence of "ARK" and "BACK" because "A" and "K" appears in both sequences and the "A" comes before the "K," even though there are letters between them in both of the original sequences.

## Finding The LCS

We basically have to look at three cases.

$c[i,j] = \begin{cases} 0 & \text{if } i = 0 \text{ or } j = 0\\ c[i - 1, j - 1] + 1 & \text{if } i,j > 0 \text{ and } x_i = y_j\\ \text{Max}(c[i, j-1], c[i - 1, j) & \text{if } i, j > 0 \text{ and } x_i \neq y_j \end{cases}$

• $c$ is the lookup table that we will build
• $i$ is the index of an element in $X$ and a row in $c$
• $j$ is the index of an element in $Y$ and a column in $c$
• $X$ and $Y$ use 1-based indices, not 0-based

### The Length-Finder

As with the other Dynamic Programming algorithms, our function is going to find the length of the longest sequence, not the actual sequence (this is left to another function).

\begin{algorithm}
\caption{Longest Common Subsequence Length Finder}
\begin{algorithmic}
\INPUT $X$: Sequence of elements
\INPUT $Y$: Sequence of elements
\OUTPUT Table of common subsequence lengths, table to reconstruct the best subsequence.
\PROCEDURE{LCS-Length}{$X, Y$}
\STATE $m \gets X.length$
\STATE $n \gets Y.length$
\STATE $b[1\ldots m, 1 \ldots n]$ and $c[0 \ldots m, 0 \ldots n]$ are new tables.
\FOR {$i \in \{1 \ldots m\}$}
\STATE $c[i, 0] \gets 0$
\ENDFOR

\FOR {$j \in \{0 \ldots n\}$}
\STATE $c[0, j] \gets 0$
\ENDFOR

\FOR {$i \in \{1 \ldots m\}$}
\FOR {$j \in \{1 \ldots n\}$}
\IF {$X_i = Y_j$}
\STATE $c[i, j] \gets c[i - 1, j - 1] + 1$
\STATE $b[i, j] \gets "\nwarrow"$
\ELSEIF {$c[i - 1, j] \geq c[i, j-1]$}
\STATE $c[i, j] \gets c[i - 1, j]$
\STATE $b[i, j] \gets "\uparrow"$
\ELSE
\STATE $c[i, j] \gets c[i, j - 1]$
\STATE $b[i,j] \gets "\leftarrow"$
\ENDIF
\ENDFOR
\ENDFOR
\RETURN $c, b$
\ENDPROCEDURE
\end{algorithmic}
\end{algorithm}


### Implementation

# python
from collections import namedtuple

# pypi
from attrs import define
from expects import contain_only, equal, expect

Direction = namedtuple("Direction", ["up_left", "up", "left"],
defaults=["\\", "|", "-"])()

@define
class LongestCommonSubsequence:
"""Finds the longest common subsequence of two sequences

Args:
x: first sequence
y: second sequence
"""
x: list
y: list
_table: list=None
_trace_table: list=None
_length: int=None

@property
def table(self) -> list:
"""The memoization table"""
if self._table is None:
self._table = [[0] + [None] * len(self.y)
for row in range(len(self.x) + 1)]
for column in range(len(self.y) + 1):
self._table[0][column] = 0
return self._table

@property
def trace_table(self) -> list:
"""Table to reconstruct longest common subsequence"""
if self._trace_table is None:
self._trace_table = [[None] * len(self.y)
for row in range(len(self.x))]
return self._trace_table

@property
def length(self) -> int:
"""Number of characters in longest common subsequence"""
if self._length is None:
for row in range(1, len(self.x) + 1):
for column in range(1, len(self.y) + 1):
previous_row, previous_column = row - 1, column - 1
unpadded_row, unpadded_column = row - 1, column - 1
if self.x[previous_row] == self.y[previous_column]:
self.table[row][column] = self.table[previous_row][previous_column] + 1
self.trace_table[unpadded_row][unpadded_column] = Direction.up_left
elif self.table[previous_row][column] >= self.table[row][previous_column]:
self.table[row][column] = self.table[previous_row][column]
self.trace_table[unpadded_row][unpadded_column] = Direction.up
else:
self.table[row][column] = self.table[row][previous_column]
self.trace_table[unpadded_row][unpadded_column] = Direction.left
self._length = self.table[-1][-1]
return self._length

def print_lcs(self, row: int, column: int) -> None:
"""Prints the elements of the longest common subsequence

Note:
to start row and column should match the last cell in the trace table.

Args:
row: row in the trace_table to start with
column: column in the trace_table to start with
"""
if row < 0 or column < 0:
return

if self.trace_table[row][column] == Direction.up_left:
self.print_lcs(row - 1, column - 1)
print(self.x[row])
elif self.trace_table[row][column] == Direction.up:
self.print_lcs(row - 1, column)
else:
self.print_lcs(row, column - 1)
return

def print_longest(self):
self.print_lcs(len(self.trace_table) - 1, len(self.trace_table[0]) - 1)
return

x = "A B C B D A B".split()
y = "B D C A B A".split()
lcs = LongestCommonSubsequence(x, y)

expect(len(lcs.table)).to(equal(len(x) + 1))
expect(len(lcs.table[0])).to(equal(len(y) + 1))
for row in lcs.table:
expect(row[0]).to(equal(0))

expect(lcs.table[0]).to(contain_only(*([0] * (len(y) + 1))))
expect(len(lcs.trace_table)).to(equal(len(x)))
expect(len(lcs.trace_table[0])).to(equal(len(y)))
expect(lcs.length).to(equal(4))


## Printing The Sequence

lcs.print_longest()

B
C
B
A


# Comparing Energy Prices

# from python
from fractions import Fraction

def print_energy(megajoules_per_dollar: float, source: str) -> None:
"""Print the energy provided by the source.

Args:
megajoules_per_dollar: how many megajoules provided by the source for each dollar
source: where the energy is coming from
"""
print(f"{source} provides {megajoules_per_dollar:0.2f} megajoules per dollar.")




## How Many Megajoules For a Dollar Of Natural Gas?

• A standard cubic foot of natural gas has about $1.1 \times 10^6$ joules of energy.
• You can get about $5\times 10^5$ BTUs of gas for a dollar.
• There are about 1,030 BTUs in a cubic foot.
\begin{align} \left(1.1 \times 10^6 \frac{joules}{\cancel{cubic foot}}\right)\left(\frac{1\ \cancel{cubic foot}}{1.03 \times 10^3\ BTU}\right) &= \left(\frac{1.1 \times 10^3\ joules}{1.03\ \cancel{BTU}}\right)\left(\frac{5 \times 10^5\ \cancel{BTU}}{dollar}\right)\\ &= \left(\frac{5.5 \times 10^8\ \cancel{joules}}{1.03\ dollars}\right)\left(\frac{1\ Megajoule}{10^6 \cancel{joules}}\right)\\ &= \left(\frac{5.5 \times 10^2\ Megajoules}{1.03 dollars}\right)\\ &\approx 533.98\ \frac{Megajoules}{dollar}\\ \end{align}

### Check The Math

joules_per_cubic_foot = 1.1E6
cubic_feet_per_btu = 1/1030
btus_per_dollar = 5e5

joules_per_btu = joules_per_cubic_foot * cubic_feet_per_btu

joules_per_dollar = joules_per_btu * btus_per_dollar
natural_gas_megajoules_per_dollar = joules_per_dollar * MEGAJOULES_PER_JOULE
natural_gas = 1/natural_gas_megajoules_per_dollar
print_energy(natural_gas_megajoules_per_dollar, "Natural Gas")
COSTS.append((natural_gas, "Natural Gas"))

Natural Gas provides 533.98 megajoules per dollar.
Natural Gas costs $0.0019 per megajoule.  ## How Many Megajoules For a Dollar of Coal? • A Ton of coal holds about $3.2 × 10^{10}$ joules of energy. • A Ton of coal costs$40.

How many Megajoules of energy in the form of coal can you get for a dollar?

\begin{align} \left(\frac{3.2 \times 10^{10}\ joules}{40\ dollars}\right) &= \left(\frac{8 \times 10^{8} joules}{dollar}\right)\\ &= \left(\frac{8 \times 10^{8}\ \cancel{joules}}{dollar}\right)\left(\frac{1\ Megajoule}{1 \times 10^6\ \cancel{joules}}\right)\\ &=8 \times 10^2 \frac{Megajoules}{dollar} \end{align}

### Checking the Math

coal_megajoules_per_dollar = MEGAJOULES_PER_JOULE * 3.2e10/40
coal = 1/coal_megajoules_per_dollar
print_energy(coal_megajoules_per_dollar, "Coal")
COSTS.append((coal, "Coal"))

Coal provides 800.00 megajoules per dollar.
Coal costs $0.0013 per megajoule.  ## How Many Megajoules Per Dollar In Corn (Biodiesel)? • Corn oil costs about$0.10 per fluid ounce (wholesale).
• A fluid ounce carries about 240 dietary calories (kilo-calories).
• A calorie is about 4.2 joules.

How many Megajoules of energy in the form of corn oil can you get for a dollar?

\begin{align} \left(\frac{1\cancel{fluid-ounce}}{0.10\ dollar}\right)\left(\frac{240\ kilocalories}{\cancel{fluid- ounce}}\right) &= \left(\frac{240\ \cancel{kilocalories}}{0.10\ dollar}\right)\left(\frac{10^3\ calories}{1\ \cancel{kilocalorie}}\right)\\ &=\left(\frac{24\times 10^4\ \cancel{calories}}{0.10\ dollar}\right)\left(\frac{4.2\ joules}{1 \cancel{calorie}}\right)\\ &=\left(\frac{100.8\times 10^4\ joules}{1 \times 10^{-1} dollar}\right)\\ &=\left(\frac{1.008 \times 10^7 joules}{dollar}\right)\left(\frac{1\ megajoule}{10^6\ joules}\right)\\ &= 1.008 \times 10^{1}\ \frac{megajoules}{dollar}\\ &= 10.008 \end{align}

### Checking the Math

ounces_per_dollar = 1/0.10
kilocalories_per_ounce = 240
joules_per_calorie = 4.2
calories_per_kilocalorie = 10**3

corn_megajoules_per_dollar = (ounces_per_dollar
* kilocalories_per_ounce
* calories_per_kilocalorie
* joules_per_calorie
* MEGAJOULES_PER_JOULE)
corn = 1/corn_megajoules_per_dollar
print_energy(corn_megajoules_per_dollar, "Corn")
COSTS.append((corn, "Corn"))

Corn provides 10.08 megajoules per dollar.
Corn costs $0.0992 per megajoule.  ## What Is the Ratio of the Cost of the Most Expensive to the Cheapest? most_expensive, cheapest, cost, name = max(COSTS), min(COSTS), 0, 1 print(f"The most expensive form of energy is {most_expensive[name]} " f"at${most_expensive[cost]:0.2f} per megajoule.")
print(f"The cheapest form of energy is {cheapest[name]} at "
f"${cheapest[cost]:0.4f} per megajoule.")  The most expensive form of energy is Corn at$0.10 per megajoule.
The cheapest form of energy is Coal at $0.0013 per megajoule.  print(f"{most_expensive[name]} is {most_expensive[cost]/cheapest[cost]:0.2f} " f"times more expensive than {cheapest[name]}.")  Corn is 79.37 times more expensive than Coal.  ## Source # Number of Subsets In a Set ## Short Answer The number of subsets for a set of n elements is $2^n$. ## One Way To Remember It Let's say we have a set of n things and we want to know how many different ways we can pick subsets of those things (including all of them or none of them). If we have three things, for instance, we might pick only the first item, or the first and the last item, or only the second item, etc. We can represent the items as boxes in a row and if an item is in the subset we put a one in its box and if it isn't then we put a zero there. So for the three item example we can represent the subset that has the first and third item as 101 (the boxes are invisible and elastic). Writing it this way we can think of each element of the set as a binary digit and the number of subsets you can create is the same as the number of binary numbers with the same number of digits as the elements in our set. To see all the possible subsets you can write out the binary numbers. For the case of three elements we represent them as three digits so finding the subsets will be the equivalent of counting from 0 to 7 in binary. Item 1 Item 2 Item 3 0 0 0 0 1 0 0 1 2 0 1 0 3 0 1 1 4 1 0 0 5 1 0 1 6 1 1 0 7 1 1 1 In this case we get 8 possibilities. More generally, for n digits we can count from $0$ to $2^n - 1$ for a total of $2^n$ numbers so the total number of subsets of a set with n elements is $2^n$. # Algorithms Illuminated - Part III: Greedy Algorithms and Dynamic Programming ## Citation • [AI3] Roughgarden T. Algorithms illuminated. Part 3: Greedy algorithms and dynamic programming. First edition. San Francisco, CA: Soundlikeyourself Publishing, LLC; 2019. 217 p. ## On The Web ## Other Books in the Series # The Knapsack Problem ## The Knapsack Problem(s) The basic premise of the Knapsack Problem is that we have a knapsack with a maximum capacity (C) and a selection of n items, each item having a weight ($w_i$) and a value ($v_i$), and we want to pick the items that give us the most total value without exceeding the capacity of our knapsack. There are two forms of the problem. The first form is what CLRS calls the $\textit{0-1 Knapsack Problem}$ because we can either take the item ($1$) or not take the item ($0$) as opposed to the $\textit{Fractional Knapsack Problem}$ where the "items" are things that we can take fractional portions of (if we're loading up on spices the $\textit{0-1}$ problem might be pre-packaged spices while the $\textit{fractional}$ problem involves scooping spices out of the bulk bins). Although the problems look more or less the same, the fractional problem is one that can be solved using a greedy algorithm while the 0-1 problem can fail if you use a greedy approach so the all-or-nothing version is the more interesting case and the one we'll look at here. ### Imports # python from collections import namedtuple from functools import partial import random # pypi from attrs import define from expects import (be_below_or_equal, contain_exactly, equal, expect, raise_error) from joblib import Parallel, delayed import altair import numpy import pandas # other monkey stuff from graeae.visualization.altair_helpers import output_path, save_chart from graeae import Timer  ### Setup Solution = namedtuple("Solution", "value inventory count") TableSolution = namedtuple("TableSolution", "value inventory count table") NEGATIVE_INFINITY = float("-inf") TIMER = Timer() SLUG = "the-knapsack-problem" OUTPUT_PATH = output_path(SLUG) save_it = partial(save_chart, output_path=OUTPUT_PATH)  ## The Brute Force Solution One way to find the optimal load for our knapsack is to find all the possible loads and taking the best load from those that will fit in the knapsack. This requires us to calculate the values for every subset of items which means we'll have to do $2^n$ calculations. ### An Iterative Brute-Force Version To have a reference solution I'll make a brute-force iterative version of the knapsack solver. This will have a greater runtime ($n 2^n$) but sometimes starting with something easy is useful, even if it's not the solution we ultimately want. def brute_knapsack(capacity: int, values: list, weights: list) -> Solution: """Finds the best selection of items to maximize value Args: capacity: maximum weight allowed values: value for each item available weights: how much each item weighs Returns: best-value, count of each item, count of loops """ assert len(values) == len(weights) number_of_items = len(values) items = [0] * number_of_items best_value = 0 count = 0 for combination in range(2**(number_of_items)): value, weight, carry = 0, 0, 1 for item in range(number_of_items): increment = items[item] + carry keep = increment % 2 carry = increment//2 if keep: value += values[item] weight += weights[item] items[item] = keep count += 1 if weight <= capacity and value > best_value: best_value = value solution = items[:] return Solution(value=best_value, inventory=solution, count=count)  ### Checking the Brute def totals(inventory: list, values: list, weights: list) -> tuple: """Reduce the inventory values and weights to a total value and a total weight Args: inventory: list of item counts values: list of values for the items list of weights for the items Returns: value, weight for the inventory """ value = sum((value for index, value in enumerate(values) if inventory[index])) weight = sum((weight for index, weight in enumerate(weights) if inventory[index])) return value, weight  def check_solution(solution: Solution, expected_inventory: list, values: list, weights: list, capacity: int): """Check that the solution matches the expected Args: solution: namedtuple with the knapsack solution expected_inventory: list of 0's and 1's representing which items to keep expected_value: total value expected for the solution values: values for the items weights: weights for the items capacity: maximum weight for the knapsack Raises: AssertionError if something isn't the expected """ expect(solution.inventory).to(contain_exactly(*expected_inventory)) value, weight = totals(solution.inventory, values, weights) expected_value, expected_weight = totals(expected_inventory, values, weights) expect(weight).to(be_below_or_equal(capacity)) expect(weight).to(equal(expected_weight)) expect(value).to(equal(expected_value)) return  def check_examples(solver: object) -> None: """Check the toy examples Args: solver: function to find the optimal knapsack load """ # values and weights don't match # broken = lambda : solver(5, [0, 1], [2, 1, 3]) # expect(broken).to(raise_error(AssertionError)) capacity = 10 values = [42, 12, 40, 25] weights = [7, 3, 4, 5] expected = [0, 0, 1, 1] solution = solver(capacity, values, weights) check_solution(solution, expected, values, weights, capacity) capacity = 6 values = [3, 2, 4, 4] weights = [4, 3, 2, 3] expected = [0, 0, 1, 1] solution = solver(capacity, values, weights) check_solution(solution, expected, values, weights, capacity) capacity = 18 values = [0, 3, 7, 7, 2, 5, 3, 0] weights = [4, 4, 6, 6, 1, 5, 2, 5] expected = [0, 0, 1, 1, 1, 1, 0, 0] solution = solver(capacity, values, weights) check_solution(solution, expected, values, weights, capacity) # this won't work for greedy algorithms capacity = 10 values = [42, 20, 25, 6] weights = [7, 4, 5, 6] expected = [0, 1, 1, 0] return check_examples(brute_knapsack)  Let's look at a particular solution. values = [3, 4, 2, 4] weights = [4, 2, 3, 3] capacity = 6 solution = brute_knapsack(capacity=capacity, values=values, weights=weights) print(f"Call Count: {solution.count}") print(f"Chosen knapsack value {solution.value}") print(f"Item inventory: {solution.inventory}") expect(solution.count).to(equal(len(values) * 2**len(values))) expect(solution.value).to(equal(8)) expect(solution.inventory).to(contain_exactly(0, 1, 0, 1))  Call Count: 64 Chosen knapsack value 8 Item inventory: [0, 1, 0, 1]  We have a solution that works, but the runtime is $n2^n$ so let's make a version that does a little better. ### A Recursive Exhaustive Search def exhausted(capacity: int, values: list, weights: list, this_item: int=0) -> Solution: """Find the optimal knapsack using an exhaustive search Args: capacity: how much weight the knapsack can hold values: how much the items are worth weights: hom much the items weigh this_item: index of the current item in the values and weights count: number of times this function is called """ assert len(values) == len(weights) next_item = this_item + 1 # quit this branch if the knapsack is already out of space if capacity == 0: return Solution(0, [0] * (len(weights) - this_item), 1) # to save on an extra base-case call handle the last item separately here if next_item == len(weights): skip_this_item = Solution(0, [0], 1) if weights[this_item] > capacity: return skip_this_item use_this_item = Solution(value=values[this_item], inventory=[1], count=1) return max((skip_this_item, use_this_item), key=lambda x: x.value) # now on to the recursive cases descendant_solution = exhausted(this_item=next_item, capacity=capacity, values=values, weights=weights) skip_count = descendant_solution.count + 1 skip_this_item = Solution(value=descendant_solution.value, inventory=[0] + descendant_solution.inventory, count=skip_count) if capacity < weights[this_item]: solution = skip_this_item count = skip_count else: capacity_after_this_item_is_added = capacity - weights[this_item] descendant_solution = exhausted( this_item=next_item, capacity=capacity_after_this_item_is_added, values=values, weights=weights) check_count = skip_count + descendant_solution.count include_this_item = Solution(value=values[this_item] + descendant_solution.value, inventory=[1] + descendant_solution.inventory, count=check_count) skip_this_item = Solution(value=skip_this_item.value, inventory=skip_this_item.inventory, count=check_count) solution = max((skip_this_item, include_this_item), key=lambda x: x.value) count = check_count return solution check_examples(exhausted)  ### Checking The Exhaustive Let's look at that example that we looked at for the iterative brute-force version. values = [3, 4, 2, 4] weights = [4, 2, 3, 3] capacity = 6 solution = exhausted(capacity=capacity, values=values, weights=weights) brute_solution = brute_knapsack(capacity=capacity, values=values, weights=weights) print(f"Call Count: {solution.count}") print(f"Chosen knapsack value {solution.value}") print(f"Item inventory: {solution.inventory}") expect(solution.value).to(equal(brute_solution.value)) expect(solution.inventory).to(contain_exactly(*brute_solution.inventory))  Call Count: 12 Chosen knapsack value 8 Item inventory: [0, 1, 0, 1]  So now the number calls has gone down to $\approx 2^n$, which is better, but not what we want just yet. ## Levitin's Memory Function This is a memoized function that is in Levitin's book. It looks slightly different from the other memoized functions in the other books (but they all look slightly different from each other anyway) but it's only cosmetic. I've been creating the final solution list of items to use in the functions themselves but I'm going to try doing it the way the books do and separate out the solution using a re-creation function afterwards. ### Some Pseudocode #### The Memoizer Note: Levitin keeps the weights, values, and solution table in the global space so it doesn't appear in the pseudocode. I'm going to copy that here but change it when I get to implementing it. I'm also going to change the variables a little to get them a little closer to the names I use. I'll call the eternal collections $\textit{Table, Weights}$, and $\textit{Values}$. The $Table$ is an $items \times capacity$ table, with from 0 to number of items rows and 0 to the capacity columns. The 0 row and 0 column get initialized with 0 and the other cells with -1. If we have 4 items and a knapsack capacity of 5 we'd have an initial table like this. 0 1 2 3 4 5 0 0 0 0 0 0 0 1 0 -1 -1 -1 -1 -1 2 0 -1 -1 -1 -1 -1 3 0 -1 -1 -1 -1 -1 4 0 -1 -1 -1 -1 -1 Where the rows are the items and the columns are the used-capacities for the knapsack. \begin{algorithm} \caption{Memory Function Knapsack Solver} \begin{algorithmic} \INPUT$i$: the number of the first items to consider. \INPUT$c$: the knapsack's capacity. \OUTPUT Value of the optimal subset of the first$i$items that fit in the knapsack. \PROCEDURE{MFKnapsack}{$i, c$} \IF {\textit{Table}$[i, c] < 0$} \IF {$c < \textit{Weights}[i]$} \STATE$v \gets $\textsc{MFKnapsack}($i - 1, c$) \ELSE \STATE$v \gets $\textsc{Max}(\textsc{MFKnapsack}($i - 1, c$),$\textit{Values}[i] + $\textsc{MFKnapsack}($i - 1, c - \textit{Weights}[i]$)) \ENDIF \STATE$\textit{Table}[i, c] \gets v$\ENDIF \RETURN$\textit{Table}[i, c]$\ENDPROCEDURE \end{algorithmic} \end{algorithm}  To start the function you would pass in the total number of items as the argument for $i$. Since we initialized the cells (other than the zero row and column) with -1 the initial if is a check to see if the item and capacity passed to the function is already in the table and if it isn't we run the body but if it is we can just return the value from the table. In the body if the weight of the current item is beyond the remaining capacity of the knapsack we pick the value for the previous item using the current capacity. If the current item will fit in the knapsack then we pick the larger of the previous item's entry with the current capacity and the value of the current item plus the previous item's entry for the current capacity minus the weight of the current item - meaning we pick the bigger of the values we get if we skip this item or keep it. 1. If the item and capacity aren't in the table: • If the item's weight is greater than the remaining capacity use the previous item's value for the current capacity. • Otherwise use the greater of the previous item's value and this item's value plus the previous item's value for the current capacity minus the current item's weight (the capacity if you use the current item) • Whichever value you use, set it to the table's entry for this item and the current capacity 2. Return the table entry for this item and the current capacity #### A Reconstructor The main algorithm builds a memo-table and returns the value of the optimal solution but doesn't tell you which items are actually taken. For that we'll need a separate function. The pseudocode assumes that the weights ($w_i$) and values ($v_i$) are global variables. \begin{algorithm} \caption{Knapsack Inventory Reconstructor} \begin{algorithmic} \INPUT$A$: The table created to solve the knapsack problem \INPUT$c$: the knapsack's capacity. \OUTPUT An optimal knapsack solution \PROCEDURE{KnapsackReconstruction}{$A, C$} \STATE $S \gets \emptyset$ \STATE $c \gets C$ \FOR {$i = n \ldots 1$} \IF {$w_i \leq c$ and $A[i - 1][c - w_i] + v_i \geq A[i - 1][c]$} \STATE $S \gets S \cup \{i\}$ \STATE $c \gets c - w_i$ \ENDIF \ENDFOR \RETURN $S$ \ENDPROCEDURE \end{algorithmic} \end{algorithm}  ### Memory-Function Knapsack The counts and such are cluttering up the function so I'm going to make this class-based. @define class Memorizer: """Dynamic Programming solution to the knapsack problem Args: capacity: total capacity (by weight) of knapsack values: values for items to put in knapsack weights: weights for items to put in knapsack """ capacity: int values: list weights: list _items: int=None _table: list=None count: int=0 _value: int=None _inventory: list=None @property def items(self) -> int: """The number of items available for the knapsack Raises: AssertionError: the number of values and weights don't match Returns: number of items """ if self._items is None: assert len(self.values) == len(self.weights) self._items = len(self.values) return self._items @property def value(self) -> int: """The total value of the optimal knapsack""" if self._value is None: self._value = self.find_value(self.items, self.capacity) return self._value @property def table(self) -> list: """The memo table Returns: items + 1 x capacity + 1 list of lists: 0's in 0 column/row, -1 elsewhere """ if self._table is None: first_row = [0] * (self.capacity + 1) row = [0] + [-1] * self.capacity table = [row[:] for item in range(self.items)] self._table = [first_row] + table return self._table def find_value(self, item: int, capacity: int) -> int: """Find the best total value for the knapsack Args: item: the number of the item to use (0...item) capacity: maximum weight allowed Returns: best-value """ self.count += 1 # the table is padded # so we need to adjust the item index for weights, values this_item = item - 1 if self.table[item][capacity] < 0: previous_item = item - 1 previous_value = self.find_value(previous_item, capacity) if capacity < self.weights[this_item]: value = previous_value else: value = max(previous_value, self.values[this_item] + self.find_value( previous_item, capacity - self.weights[this_item])) self.table[item][capacity] = value return self.table[item][capacity] @property def inventory(self) -> list: """Reconstructs the optimal knapsack load using the table Returns: inventory of items in the optimal knapsack """ if self._inventory is None: # make sure that the problem has already been solved self() self._inventory = [0] * self.items remaining_capacity = self.capacity for table_item in range(self.items, 0, -1): # table_item is padded by one this_item = previous_table_item = table_item - 1 if (self.weights[this_item] <= remaining_capacity and self.table[previous_table_item][ remaining_capacity - self.weights[this_item]] + self.values[this_item] >= self.table[previous_table_item][remaining_capacity]): self._inventory[this_item] = 1 remaining_capacity -= self.weights[this_item] return self._inventory def __call__(self) -> int: """Finds the best solution: As a side effect this also sets self.value Returns: value for optimal knapsack """ return self.value  #### Check the table maker capacity, items = 5, 4 values = weights = [0] * items table = Memorizer(capacity=capacity, weights=weights, values = values).table # one row per item plus a zero row expect(len(table)).to(equal(items + 1)) # columns from 0...capacity expect(len(table[0])).to(equal(capacity + 1)) # first row should be 0's expect(sum(table[0])).to(equal(0)) # first column should be 0's expect(sum(row[0] for row in table)).to(equal(0)) # everything else should be -1 (items x capacity sub-array) expect(sum(sum(row) for row in table)).to(equal(-1 * (items * capacity)))  #### Check the Final Table weights = [2, 1, 3, 2] values = [12, 10, 20, 15] capacity = 5 memoizer = Memorizer(weights=weights, values=values, capacity=capacity) memoizer() expect(memoizer.value).to(equal(37)) expected_table = [[0, 0, 0, 0, 0, 0], [0, 0, 12, 12, 12, 12], [0, -1, 12, 22, -1, 22], [0, -1, -1, 22, -1, 32], [0, -1, -1, -1, -1, 37]] for row_index, row in enumerate(memoizer.table): expect(row).to(contain_exactly(*expected_table[row_index]))  #### Check the Recovered Solution Although knowing what the optimal value is for the knapsack is somewhat informative in that it tells us what we can expect to achieve, it isn't really the solution since we don't know what items actually give us this value, so we're going to need to reconstruct it from the table. weights = [2, 1, 3, 2] values = [12, 10, 20, 15] capacity = 5 solution = Memorizer(capacity=capacity, values=values, weights=weights) expect(solution.inventory).to(contain_exactly(1, 1, 0, 1))  #### Check It Against The Examples values = [3, 4, 2, 4] weights = [4, 2, 3, 3] capacity = 6 solution = Memorizer(capacity, values, weights) print(f"Chosen knapsack value {solution.value}") print(f"Item inventory: {solution.inventory}") print(f"Call Count: {solution.count}") check_examples(Memorizer)  Chosen knapsack value 8 Item inventory: [0, 1, 0, 1] Call Count: 17  Our solution is correct, but if you count all the function calls, not just the calls where the solution isn't in the table yet, it takes more calls than our exhaustive function. ### Compared to the Exhaustive Search sizes = list(range(2, 61)) # 0 weights break the Memorizer so make sure everything weighs at least 1 values = [random.choices(list(range(1, size)), k=size) for size in sizes] weights = [random.choices(list(range(1, random.randint(1, size) * size)), k=size) for size in sizes] capacities = [sum(random.choices(weight, k=4)) for weight in weights] capacities_values_weights = list(zip(capacities, values, weights))  with TIMER: exhaustive_output = Parallel(n_jobs=-1)( delayed(exhausted)(capacity, values, weights) for capacity,values,weights in capacities_values_weights)  Started: 2022-07-11 02:35:10.337731 Ended: 2022-07-11 02:50:47.862952 Elapsed: 0:15:37.525221  def memorizer_knapsack(capacity, values, weights): memorizer = Memorizer(capacity, values, weights) memorizer() return memorizer with TIMER: memorized_output = Parallel(n_jobs=-1)( delayed(memorizer_knapsack)(capacity, values, weights) for capacity,values,weights in capacities_values_weights) for index, output in enumerate(exhaustive_output): try: expect(output.value).to(equal(memorized_output[index].value)) capacity, value, weight = capacities_values_weights[index] b_value, b_weight = totals(output.inventory, value, weight) m_value, m_weight = totals(memorized_output[index].inventory, value, weight) expect(m_value).to(equal(b_value)) expect(m_weight).to(be_below_or_equal(capacity)) except AssertionError as error: c, v, w = capacities_values_weights[index] print(f"Index: {index}") print(error) print(f"Brute: {brute_knapsack(c, v, w)}") raise  Started: 2022-07-11 05:00:50.158584 Ended: 2022-07-11 05:00:52.849549 Elapsed: 0:00:02.690965  frame = pandas.DataFrame({"Items": sizes, "Exhaustive": [ solution.count for solution in exhaustive_output], "Memoized": [ solution.count for solution in memorized_output]}) melted = frame.melt(id_vars=["Items"], value_vars=["Exhaustive", "Memoized"], var_name="Algorithm", value_name="Calls") chart = altair.Chart(melted).mark_line(point=True).encode( x="Items", y="Calls", color="Algorithm", tooltip=[altair.Tooltip("Items", format=","), altair.Tooltip("Calls", format=","), "Algorithm"], ).properties( title="Exhaustive vs Memoized Knapsack Solution", width=800, height=525 ) save_it(chart, "exhaustive-vs-memoized")  The height of those last points squashes the previous points down to make it look like the two algorithms do about the same until you hit 47 items, but if you trim off those end points you'll see that the exhaustive algorithm generally requires much more calls than the memoized version. The points aren't on a smooth line as a function of the number of items because whenever an item won't fit in the remaining capacity of the knapsack we skip the second recursive call. UPPER_BOUND = 47 trimmed = melted[melted.Items < UPPER_BOUND] chart = altair.Chart(trimmed).mark_line(point=True).encode( x="Items", y="Calls", color="Algorithm", tooltip=[altair.Tooltip("Items", format=","), altair.Tooltip("Calls", format=","), "Algorithm"], ).properties( title=f"Exhaustive vs Memoized Knapsack Solution (< {UPPER_BOUND})", width=800, height=525 ) save_it(chart, "exhaustive-vs-memoized-trimmed")  ## Dynamic Programming This is taken from Algorithms Illuminated Part 3. ### Some Pseudocode \begin{algorithm} \caption{Dynamic Programming Knapsack Solver} \begin{algorithmic} \INPUT Item Values: $v_1, v_2, \ldots, v_n$ \INPUT Item Weights: $w_1, w_2, \ldots, w_n$ \INPUT Knapsack Capacity $C$ \OUTPUT Subset $S$ of items with maximum possible sum of values and size at most $C$ \PROCEDURE{DynamicKnapsack}{$v, w, C$} \STATE $A \gets (n + 1) \times (c + 1)$ two dimensional array. \FOR {$c \in \{0 \ldots C\}$} \STATE $A[0][c] \gets 0$ \ENDFOR \FOR {$i \in \{1 \ldots n\}$} \FOR {$c \in \{0 \ldots C \}$} \IF {$w_i > C$} \STATE $A[i][c] \gets A[i - 1][c]$ \ELSE \STATE $A[i][c] \gets$\textsc{Max}($A[i - 1][c], A[i - 1][c - w_i] + v_i$) \ENDIF \ENDFOR \ENDFOR \RETURN $A[n][C]$ \ENDPROCEDURE \end{algorithmic} \end{algorithm}  This looks pretty much like the memoized-recursive version so it shouldn't be too hard to understand. ### In Python @define class CaptainDynamic(Memorizer): """Dynamic Programming solution to the knapsack problem Args: capacity: total capacity (by weight) of knapsack values: values for items to put in knapsack weights: weights for items to put in knapsack """ @property def value(self) -> int: """The total value of the optimal knapsack""" if self._value is None: self._value = self.find_value() return self._value def find_value(self) -> int: """Finds the optimal value""" for item_row in range(1, self.items + 1): previous_item = this_item = item_row - 1 for capacity in range(self.capacity + 1): skip_this_item = self.table[previous_item][capacity] if self.weights[this_item] > capacity: self.table[item_row][capacity] = skip_this_item else: use_this_item = ( self.table[previous_item][ capacity - self.weights[this_item]] + self.values[this_item]) self.table[item_row][capacity] = max( (skip_this_item, use_this_item) ) self.count += 1 return self.table[self.items][self.capacity]  #### Check the table maker capacity, items = 5, 4 values = weights = [0] * items table = CaptainDynamic(capacity=capacity, weights=weights, values = values).table # one row per item plus a zero row expect(len(table)).to(equal(items + 1)) # columns from 0...capacity expect(len(table[0])).to(equal(capacity + 1)) # first row should be 0's expect(sum(table[0])).to(equal(0)) # first column should be 0's expect(sum(row[0] for row in table)).to(equal(0)) # everything else should be -1 (items x capacity sub-array) expect(sum(sum(row) for row in table)).to(equal(-1 * (items * capacity)))  #### Check the Recovered Solution Although knowing what the optimal value is for the knapsack is somewhat informative in that it tells us what we can expect to achieve, it isn't really the solution since we don't know what items actually give us this value, so we're going to need to reconstruct it from the table. weights = [2, 1, 3, 2] values = [12, 10, 20, 15] capacity = 5 solution = CaptainDynamic(capacity=capacity, values=values, weights=weights) expect(solution.inventory).to(contain_exactly(1, 1, 0, 1)) m_solution = Memorizer(capacity=capacity, values=values, weights=weights) expect(m_solution.inventory).to(contain_exactly(1, 1, 0, 1))  #### Check It Against The Examples values = [3, 4, 2, 4] weights = [4, 2, 3, 3] capacity = 6 d_solution = CaptainDynamic(capacity, values, weights) print("Captain Dynamic") print(f"Chosen knapsack value {d_solution.value}") print(f"Item inventory: {d_solution.inventory}") print(f"Call Count: {d_solution.count}") check_examples(CaptainDynamic) print("\nMemorizer") m_solution = Memorizer(capacity, values, weights) print(f"Chosen knapsack value: {m_solution.value}") print(f"Item inventory: {m_solution.inventory}") print(f"Call Count: {m_solution.count}") check_examples(Memorizer)  Captain Dynamic Chosen knapsack value 8 Item inventory: [0, 1, 0, 1] Call Count: 28 Memorizer Chosen knapsack value: 8 Item inventory: [0, 1, 0, 1] Call Count: 17  Our solution is correct, but if you count all the function calls, not just the calls where the solution isn't in the table yet, it takes more calls than our exhaustive function. ### Compared def dynamic_knapsack(capacity, values, weights): captain = CaptainDynamic(capacity, values, weights) captain() return captain with TIMER: dynamic_output = Parallel(n_jobs=-1)( delayed(dynamic_knapsack)(capacity, values, weights) for capacity,values,weights in capacities_values_weights) for index, output in enumerate(exhaustive_output): try: expect(output.value).to(equal(dynamic_output[index].value)) capacity, value, weight = capacities_values_weights[index] b_value, b_weight = totals(output.inventory, value, weight) m_value, m_weight = totals(dynamic_output[index].inventory, value, weight) expect(m_value).to(equal(b_value)) expect(m_weight).to(be_below_or_equal(capacity)) except AssertionError as error: c, v, w = capacities_values_weights[index] print(f"Index: {index}") print(error) print(f"Brute: {brute_knapsack(c, v, w)}") raise  Started: 2022-07-11 05:08:19.593276 Ended: 2022-07-11 05:08:21.585963 Elapsed: 0:00:01.992687  frame = pandas.DataFrame({"Items": sizes, "Dynamic": [ solution.count for solution in dynamic_output], "cN": [ solution.capacity * solution.items for solution in dynamic_output], "Memoized": [ solution.count for solution in memorized_output]}) melted = frame.melt(id_vars=["Items"], value_vars=["Dynamic", "Memoized", "cN"], var_name="Algorithm", value_name="Calls") chart = altair.Chart(melted).mark_line(point=True).encode( x="Items", y="Calls", color="Algorithm", tooltip=[altair.Tooltip("Items", format=","), altair.Tooltip("Calls", format=","), "Algorithm"], ).properties( title="Dynamic vs Memoized Knapsack Solution", width=800, height=525 ) save_it(chart, "dynamic-vs-memoized")  ## Sources # Learning Algorithms Through Programming and Problem Solving ## Table Of Contents (Partial) ### Algorithms and Complexity • What is an Algorithm? • Pseudocode • Problem Versus Problem Instance • Correct Versus Incorrect Algorithms • Fast Versus Slow Algorithms • Big-O Notation ### Algorithm Design Techniques • Exhaustive Search Algorithms • Branch-and-Bound Algorithms • Greedy Algorithms • Dynamic Programming Algorithms • Recursive Algorithms • Divide-and-Conquer Algorithms • Randomized Algorithms ### Programming Challenges • Sum of Two Digits • Maximum Pairwise Product • Naive Algorithm • Fast Algorithm • Testing and Debugging • Can You Tell Me What Error I Have Made? • Stress Testing • Even Faster Algorithm • A More Compact Algorithm • Solving a Programming Challenge in Five Easy Steps • Reading the Problem Statement • Designing an Algorithm • Implementing an Algorithm • Testing and Debugging • Submitting the Grading System • Good Programming Practices ### Algorithmic Warm Up • Fibonacci Numbers • Last Digit of a Fibonacci Number • Greatest Common Divisor • Least Common Multiple • Fibonacci Numbers Again • Last Digit of the Sum of Fibonacci Numbers • Last Digit of the Sum of Fibonacci Numbers Again ### Greedy Algorithms • Money Change • Maximum Value of the Loot • Maximum Advertisement Revenue • Collecting Signatures • Maximum Number of Prizes • Maximum Salary ### Divide and Conquer • Binary Search • Majority Element • Improving QuickSort • Number of Inversions • Organizing a Lottery • Closest Points ### Dynamic Programming • Money Change Again • Primitive Calculator • Edit Distance • Longest Common Subsequence of Two Sequences • Longest Common Subsequence of Three Sequences • Maximum Amount of Gold • Partitioning Souvenirs • Maximum Value of an Arithmetic Expression ## Citation • Kulikov A. Learning algorithms through programming and puzzle solving. La Jolla, CA: Active Learning Technologies; 2018. # The Coin Changing Problem ## The Change Problem The Change Problem (Coin Changing Problem) asks us to convert some amount of money into denominations using the fewest coins. Input: An integer (money) and an array of d denominations ($c = c_1, c_2, \ldots, c_d$) in decreasing order of value ($c_1 > c_2> \cdots > c_d$). Output: A list of d integers $i_1, i_2, \ldots, i_d$ such that $c_1 \cdot i_1 + c_2 \cdot i_2 + \cdots + c_d \cdot i_d = money$ and $i_1 + i_2 + \cdots + i_d$ is as small as possible. ## SetUp ### Imports # python from collections import namedtuple from functools import partial # pypi from expects import contain_exactly, expect, equal import altair import pandas # my stuff from graeae.visualization.altair_helpers import output_path, save_chart  ### Setup ChangeCounts = namedtuple("ChangeCounts", ["coin_counts", "run_count"]) ChangeCountsTable = namedtuple("ChangeCounts", ["coin_counts", "run_count", "table"]) SLUG = "the-coin-changing-problem" OUTPUT_PATH = output_path(SLUG) save_it = partial(save_chart, output_path=OUTPUT_PATH) PlotParameters = namedtuple("PlotParameters", ["height", "width"], defaults=[525, 750])() Denominations = namedtuple("Denominations", "US double_dime obsolete", defaults=[[25, 10, 5, 1], [25, 20, 10, 5, 1], [25, 20, 10, 5, 3, 2, 1]])()  ## The Cashier's Algorithm The way many cashier's make change for customers is using a greedy approach. ### Pseudocode \begin{algorithm} \caption{CashiersChange} \begin{algorithmic} \INPUT Amount owed ($money$) and a supply of coins ($c$). \PROCEDURE{CashiersChange}{$money, c$} \WHILE {$money > 0$} \STATE $coin \gets$ Coin with the largest denomination that doesn't exceed $money$. \STATE Give$coin$to the customer. \STATE $money \gets money - coin$ \ENDWHILE \ENDPROCEDURE \end{algorithmic} \end{algorithm}  ### Python Implementation def cashier_changer(amount: int, denominations: list) -> ChangeCounts: """Implements the cashier's method for making change Args: amount: the value of the change needed denominations: values for coins (in descending order) Returns: list of coin-counts, number of loops ran """ coins = [0] * len(denominations) remainder = amount loop_count = 0 for index, denomination in enumerate(denominations): coin_count = 0 while denomination <= remainder: coin_count += 1 loop_count += 1 remainder = remainder - denomination coins[index] = coin_count return ChangeCounts(coin_counts=coins, run_count=loop_count)  #### A First Test Examples from LATPAPS. def check_coins(answer: list, denominations: list, expected_coins: int, expected_total: int) -> None: """Check that our answer matches the expected Args: answer: list of counts for each denomination denominations: list of coin denominations available expected_coins: the expected list of counts for each denomination expected_total: what the value our coins should add up to Raises: AssertionError: something about our answer doesn't match the expected """ expect(sum(answer)).to(equal(sum(expected_coins))) expect(answer).to(contain_exactly(*expected_coins)) expect(sum(count * coin for count, coin in zip(answer, denominations))).to( equal(expected_total) ) return  # case 1 money = 2 coins = [10, 5, 1] actual = cashier_changer(money, coins) check_coins(actual.coin_counts, coins, [0, 0, 2], money) print(f"case 1: {actual.run_count}") # case 2 money = 28 actual = cashier_changer(money, coins) check_coins(actual.coin_counts, coins, [2, 1, 3], money) print(f"case 2: {actual.run_count}") # case 3 money = 99 actual = cashier_changer(money, coins) check_coins(actual.coin_counts, coins, [9, 1, 4], money) print(f"case 2: {actual.run_count}")  case 1: 2 case 2: 6 case 2: 14  The number of loops is dependent on the change owed, with an upper limit based on the denominations. In this case if we assume you wouldn't give out more than 99 cents in change then we would cap out at $9 \times 10 + 1 \times 5 + 4 \times 1$ = 14 coins/loops. #### U.S. Denominations This is a quick check using the U.S. coins most commonly used by cashiers. By adding a 25 cent piece we reduce the upper limit on the amount of coins needed to $3 \times 25 + 2 \times 10 + 4 \times 1$ = 9 coins. print(Denominations.US)  [25, 10, 5, 1]  money = 28 actual = cashier_changer(money, Denominations.US) check_coins(actual.coin_counts, Denominations.US, [1, 0, 0, 3], money) money = 14 actual = cashier_changer(money, Denominations.US) check_coins(actual.coin_counts, Denominations.US, [0, 1, 0, 4], money)  ### Plotting the Coin Counts/Loops PLOT_AMOUNTS = list(range(1, 100)) CASHIER_COUNTS = [cashier_changer(amount, Denominations.US).run_count for amount in PLOT_AMOUNTS] PLOT_FRAME = pandas.DataFrame({"Amount to Change": PLOT_AMOUNTS, "Coins": CASHIER_COUNTS}) chart = altair.Chart(PLOT_FRAME).mark_bar().encode( x=altair.X("Amount to Change", axis=altair.Axis(tickMinStep=1)), y=altair.Y("Coins", axis=altair.Axis(tickMinStep=1)), tooltip=["Amount to Change", "Coins"]).properties( title="Cashier's Change Coin Counts (Common U.S. Coins)", width=PlotParameters.width, height=PlotParameters.height, ) save_it(chart, "cashiers-change-us-coins")  #### The Double-Dime Our cashier algorithm turns out to work for the specific set of coins that cashiers commonly use, but will it work for other coins as well? According to Wikipedia, there was at one time a proposal in the United States for a twenty-cent piece (and at one time there were half-cent, two-cent, and three-cent coins).If we include the twenty cent piece amongst our denominations, we find that there are cases where the cashier algorithm will miss the optimal solution. print(Denominations.double_dime)  [25, 20, 10, 5, 1]  money = 40 actual = cashier_changer(money, Denominations.double_dime) try: check_coins(actual.coin_counts, Denominations.double_dime, [0, 2, 0, 0, 0], money) except AssertionError as error: print(f"AssertionError: {error}")  AssertionError: expected: 3 to equal 2  Because the cashier algorithm always takes the largest possible coins first, it ends up using 25¢ + 10¢ + 5¢ as the solution instead of the optimal 20¢ + 20¢. While it might seem artificial, given the characterization of this as a solution to making change, it's important to note that generalizing the cashier algorithm beyond the curated denominations or even beyond coins specifically leaves it vulnerable to cases where it will fail. ## A Greedy Version As noted, the Cashier's Algorithm is a greedy algorithm, but let's translate it into a silghtly smarter version which uses a little arithmetic to get to the same point. ### In Pseudocode \begin{algorithm} \caption{GreedyChange} \begin{algorithmic} \INPUT An integer ($money$) and an array of$d$denominations ($c = c_1, c_2, \ldots, c_d$) in decreasing order of value ($c_1 > c_2> \cdots > c_d$). \OUTPUT A list of$d$integers $i_1, i_2, \ldots, i_d$ such that $c_1 \cdot i_1 + c_2 \cdot i_2 + \cdots + c_d \cdot i_d = money$ and $i_1 + i_2 + \cdots + i_d$ is as small as possible. \PROCEDURE{CashiersChange}{$money, c$} \STATE$remainder \gets money$\STATE$d \gets $\textsc{Length}($c$) \FOR {$k \in \{1 \ldots d\}$} \STATE$i_k \gets \lfloor \frac{remainder}{c_k} \rfloor$\STATE$remainder \gets remainder - c_k \cdot i_k$\ENDFOR \RETURN$(i_1, i_2, \ldots, i_d)\ENDPROCEDURE \end{algorithmic} \end{algorithm}  ### In Python def greedy_changer(money: int, denominations: list) -> ChangeCounts: """Make change using the fewest coins Args: money: the amount to change denominations: list of coin denominations in decreasing order Returns: Number of each denomination needed to make the change, loop-count """ coins = [0] * len(denominations) remainder = money count = 0 for location, denomination in enumerate(denominations): coins[location] = remainder // denomination remainder = remainder - denomination * coins[location] count += 1 return ChangeCounts(coins, count)  ### Testing It Out #### A First Test # case 1 money = 2 coins = [10, 5, 1] actual = greedy_changer(money, coins) check_coins(actual.coin_counts, coins, [0, 0, 2], money) # case 2 money = 28 actual = greedy_changer(money, coins) check_coins(actual.coin_counts, coins, [2, 1, 3], money)  #### U.S. Denominations money = 28 actual = greedy_changer(money, Denominations.US) check_coins(actual.coin_counts, Denominations.US, [1, 0, 0, 3], money) money = 14 actual = greedy_changer(money, Denominations.US) check_coins(actual.coin_counts, Denominations.US, [0, 1, 0, 4], money)  #### The Double-Dime Looking at the greedy-algorithm you can see that it only has one loop that traverses the denominations of coins - so it is a very quick algorithm, but while our greedy algorithm turns out to work for the specific set of coins that cashiers use, it also falls prey to sometimes missing the optimal solution. money = 40 actual = greedy_changer(money, Denominations.double_dime) try: check_coins(actual.coin_counts, Denominations.double_dime, [0, 2, 0, 0, 0], money) except AssertionError as error: print(f"AssertionError: {error}")  AssertionError: expected: 3 to equal 2  ### Plotting greedy_counts = [greedy_changer(amount, Denominations.US).run_count for amount in PLOT_AMOUNTS] GREEDY_FRAME = PLOT_FRAME.rename(columns={"Coins": "Cashier"}) GREEDY_FRAME["Greedy"] = greedy_counts melted = GREEDY_FRAME.melt(id_vars=["Amount to Change"], value_vars=["Greedy", "Cashier"], var_name="Change Method", value_name="Loop Count") chart = altair.Chart(melted).mark_line(point=True).encode( x="Amount to Change", y=altair.Y("Loop Count", axis=altair.Axis(tickMinStep=1)), tooltip=["Amount to Change", "Loop Count", "Change Method"], color="Change Method").properties( title="Cashier's and Greedy Change Coin Counts (Common U.S. Coins)", width=PlotParameters.width, height=PlotParameters.height, ) save_it(chart, "cashiers-greedy-change-us-coins")  Since the Cashier works by repeated subtraction instead of division it goes up with the amount of change owed, and so generally does worse than the Greedy changer, except in the cases where it takes fewer coins to make the change than the number of denominations there are. If I had included the fifty cent piece it would have done a little better in the second half of the plot. ## A Brute Force Changer To get around the limitations of the greedy-style of change making we'll figure out all the ways we can make change for a particular amount and pick the solution with the fewest coins. def brute_changer(money: int, denominations: list) -> ChangeCounts: """Make change using the fewest coins Args: money: the amount to change denominations: list of coin denominations in decreasing order Returns: Number of each denomination needed to make the change, loop_count """ best = float("inf") number_of_denominations = len(denominations) counter = 0 for first in range(number_of_denominations): remainder = money counts = [0] * number_of_denominations total = 0 for next_location in range(first, number_of_denominations): denomination = denominations[next_location] count = remainder//denomination remainder = remainder - count * denomination counts[next_location] = count total += count counter += 1 if total < best: best_counts = counts best = total return ChangeCounts(best_counts, counter)  # case 1 money = 2 coins = [10, 5, 1] actual = brute_changer(money, coins) check_coins(actual.coin_counts, coins, [0, 0, 2], money) # case 2 money = 28 actual = brute_changer(money, coins) check_coins(actual.coin_counts, coins, [2, 1, 3], money) coins = [25, 10, 5, 1] money = 28 actual = brute_changer(money, coins) check_coins(actual.coin_counts, coins, [1, 0, 0, 3], money) money = 14 actual = brute_changer(money, coins) check_coins(actual.coin_counts, coins, [0, 1, 0, 4], money)  And now our double-dime case. money = 40 actual = brute_changer(money, Denominations.double_dime) check_coins(actual.coin_counts, Denominations.double_dime, [0, 2, 0, 0, 0], money) print(actual)  ChangeCounts(coin_counts=[0, 2, 0, 0, 0], run_count=15)  The Brute-Force Changer fixes our double-dimes case, and has a constant runtime of: $T(n) = \sum_{i = 1}^{n} = \frac{n(n + 1)}{2}$ Where $n$ is the number of denominations, not the input amount. Let's make sure our functions all agree. for amount in range(100): cashiers_change = cashier_changer(amount, Denominations.US).coin_counts greedy_change = greedy_changer(amount, Denominations.US).coin_counts brute_change = brute_changer(amount, Denominations.US).coin_counts check_coins(answer=cashiers_change, denominations=Denominations.US, expected_coins=greedy_change, expected_total=amount) check_coins(answer=brute_change, denominations=Denominations.US, expected_coins=greedy_change, expected_total=amount)  ### Plotting Let's compare the three methods we've created so far using the typical U.S. Coins. brute_counts = [brute_changer(amount, Denominations.US).run_count for amount in PLOT_AMOUNTS] GREEDY_FRAME["Brute-Force"] = brute_counts melted = GREEDY_FRAME.melt(id_vars=["Amount to Change"], value_vars=["Greedy", "Cashier", "Brute-Force"], var_name="Change Method", value_name="Loop Count") chart = altair.Chart(melted).mark_line(point=True).encode( x="Amount to Change", y=altair.Y("Loop Count", axis=altair.Axis(tickMinStep=1)), tooltip=["Amount to Change", "Loop Count", "Change Method"], color="Change Method").properties( title="Greedy Change & Brute-Force Coin Counts (Common U.S. Coins)", width=PlotParameters.width, height=PlotParameters.height, ) save_it(chart, "brute-force-change-us-coins")  The Greedy and Brute-Force methods both have a loop-count based on the number of coin denominations we have. Since we have four denominations (25, 10, 5, and 1) the Greedy method has four loops and the Brute-Force method has ten loops $\left(^{4(4 + 1)}/_{2} = 10\right)$, so it does a little more calculating but not much, although in this case the benefit it has is wasted since we shouldn't find sub-optimal solutions with the greedy approach and using the coins that we did. ## Memoized Changer def memoized_changer(money: int, denominations: list, table: dict=None, counter: int=0) -> ChangeCounts: """Make change using the fewest coins Args: money: the amount to change denominations: list of coin denominations in decreasing order table: memoization table (largest denomination, amount to change): best coin counts counter: number of times the function is called Returns: Number of each denomination needed to make the change, call_count """ table = dict() if table is None else table counter += 1 if money == 0: return ChangeCounts([0] * len(denominations), counter) if (denominations[0], money) in table: return ChangeCounts(table[(denominations[0], money)], counter) last_denomination = len(denominations) - 1 best = float("inf") best_counts = None for current, denomination in enumerate(denominations): count = money//denomination remaining = money - count * denomination if current == last_denomination: if remaining > 0: count = float("inf") counts = [count] counter += 1 else: counts, counter = memoized_changer( remaining, denominations[current + 1:], table, counter) counts = [count] + counts total_counts = sum(counts) if total_counts < best: best = total_counts best_counts = [0] * current + counts table[(denomination, money)] = counts return ChangeCounts(best_counts, counter)  ### Checking the Memoized Changer def check_changers(changer: object, amount: int, denominations: list): """Check that the changer matches the brute-force solution Arguments: changer: change-function to compare to brute-force amount: cents to change denominations: descending list of coin-denominations Raises: AssertionError: there's a discrepancy somewhere """ brute_counts = brute_changer(amount, denominations) check_counts = changer(amount, denominations) check_coins(check_counts.coin_counts, denominations, brute_counts.coin_counts, amount) return  def check_cases(changer: object) -> None: """Check the pre-picked cases are okay Arguments: changer: the dynamic-programming changer to check Raises: AssertionError if something doesn't check out """ for amount in range(100): check_changers(changer, amount, Denominations.US) check_changers(changer, amount, Denominations.double_dime) check_changers(changer, amount, Denominations.obsolete) money = 2 coins = [10, 5, 1] check_changers(memoized_changer, money, coins) money = 28 check_changers(memoized_changer, money, coins) coins = [25, 10, 5, 1] check_changers(memoized_changer, money, coins) money = 14 check_changers(memoized_changer, money, coins) # the double-dimes coins = [25, 20, 10, 5, 1] money = 40 check_changers(memoized_changer, money, coins) return  check_cases(memoized_changer)  ### Plotting #### All the Methods So Far memoized_counts = [memoized_changer(amount, Denominations.US, {}).run_count for amount in PLOT_AMOUNTS] GREEDY_FRAME["Memoized"] = memoized_counts melted = GREEDY_FRAME.melt(id_vars=["Amount to Change"], value_vars=["Greedy", "Cashier", "Brute-Force", "Memoized"], var_name="Change Method", value_name="Loop Count") chart = altair.Chart(melted).mark_line(point=True).encode( x="Amount to Change", y=altair.Y("Loop Count", axis=altair.Axis(tickMinStep=1)), tooltip=["Amount to Change", "Loop Count", "Change Method"], color="Change Method").properties( title="Greedy, Brute-Force & Memoized Change Counts (Common U.S. Coins)", width=PlotParameters.width, height=PlotParameters.height, ) save_it(chart, "memoized-change-us-coins")  With a Brute-Force ceiling of 10 loops it's hard to see a marked improvement, given how much easier it is to write the Brute-Force method rather than the memoized one, but let's get rid of the two greedy methods and compare the brute force and memoized with more obsolete coins added in. #### Obsolete Coins First, we'll check to make sure that the brute-force and memoized methods agree on the best coins for making change with our new denominations, which include some that were once used in the United States but aren't any more. print(Denominations.obsolete)  [25, 20, 10, 5, 3, 2, 1]  Now let's see how the run-times compare. BRUTE_OBSOLETE = [brute_changer(amount, Denominations.obsolete).run_count for amount in PLOT_AMOUNTS] MEMOIZED_OBSOLETE = [memoized_changer(amount, Denominations.obsolete, {}).run_count for amount in PLOT_AMOUNTS] OBSOLETE_FRAME = pandas.DataFrame({"Amount to Change": PLOT_AMOUNTS, "Brute-Force": BRUTE_OBSOLETE, "Memoized": MEMOIZED_OBSOLETE}) melted = OBSOLETE_FRAME.melt(id_vars=["Amount to Change"], value_vars=["Brute-Force", "Memoized"], var_name="Change Method", value_name="Loop Count") chart = altair.Chart(melted).mark_line(point=True).encode( x="Amount to Change", y=altair.Y("Loop Count", axis=altair.Axis(tickMinStep=1)), tooltip=["Amount to Change", "Loop Count", "Change Method"], color="Change Method").properties( title="Brute-Force & Memoized Change Counts (Including Obsolete Coins)", width=PlotParameters.width, height=PlotParameters.height, ) save_it(chart, "memoized-change-obsolete")  This time we have seven coin denominations so the brute-force method takes $^7(8)/_2 = 28$ loops. Now that we have more coins the memoized does slightly (very slightly) worse than the brute-force for larger amounts, but in the examples so far I've been creating a new memo-table for every amount - what happens if we build one table and re-use it as we find the change for new amounts? OBSOLETE_TABLE = {} MEMOIZED_OBSOLETE_2 = [memoized_changer(amount, Denominations.obsolete, OBSOLETE_TABLE).run_count for amount in PLOT_AMOUNTS] OBSOLETE_FRAME_2 = pandas.DataFrame({"Amount to Change": PLOT_AMOUNTS, "Brute-Force": BRUTE_OBSOLETE, "Memoized": MEMOIZED_OBSOLETE_2}) melted = OBSOLETE_FRAME_2.melt(id_vars=["Amount to Change"], value_vars=["Brute-Force", "Memoized"], var_name="Change Method", value_name="Loop Count") chart = altair.Chart(melted).mark_line(point=True).encode( x="Amount to Change", y=altair.Y("Loop Count", axis=altair.Axis(tickMinStep=1)), tooltip=["Amount to Change", "Loop Count", "Change Method"], color="Change Method").properties( title="Brute-Force & Memoized Change Counts (Obsolete Coins, Re-Used Table)", width=PlotParameters.width, height=PlotParameters.height, ) save_it(chart, "memoized-change-obsolete-keep-table")  Re-using the table allows us to avoid some re-calculating and the number of calls eventually drops down to the number of coins we have, the same as the number of loops for the greedy method, although in this case we can handle cases where the greedy method fails. Even better, now that we have the table built with all our allowed values, we can just look things up and don't need to run the function (although I will here just to make sure nothing funky is going on), so it has a cost of one. MEMOIZED_OBSOLETE_3 = [memoized_changer(amount, Denominations.obsolete, OBSOLETE_TABLE).run_count for amount in PLOT_AMOUNTS] OBSOLETE_FRAME_3 = pandas.DataFrame({"Amount to Change": PLOT_AMOUNTS, "Brute-Force": BRUTE_OBSOLETE, "Memoized": MEMOIZED_OBSOLETE_3}) melted = OBSOLETE_FRAME_3.melt(id_vars=["Amount to Change"], value_vars=["Brute-Force", "Memoized"], var_name="Change Method", value_name="Loop Count") chart = altair.Chart(melted).mark_line(point=True).encode( x="Amount to Change", y=altair.Y("Loop Count", axis=altair.Axis(tickMinStep=1)), tooltip=["Amount to Change", "Loop Count", "Change Method"], color="Change Method").properties( title="Brute-Force & Memoized Change Counts (Obsolete Coins, Pre-Filled Table)", width=PlotParameters.width, height=PlotParameters.height, ) save_it(chart, "memoized-change-obsolete-pre-filled-table")  ## Iterative Dynamic Programming The idea behind this problem (making change) is to show how Dynamic Programming can eleminate redundant calculations, so this is the point where I'll switch to an iterative approach that uses Dynamic Programming. The iterative, top-down approach works in the opposite direction of the memoized approach, starting with the smallest values and then moving to ever increasing values, storing the outcomes as we go so that we can re-use the lower solutions as the values increase. This is where I was stumped for a while - since the memoized function I used is based on the remainders left over after subtracting a coin's value from it (e.g. 26 cents minus a quarter has a remainder of 1 cent so we use the stored 1-cent solution from the memo-table rather than re-calculating it) in order to create an iterative version we somehow have to store the values in the table before they are needed, but the scheme I used where we iterate over the coin-denominations doesn't work (or it makes it really convoluted to force it to work). From what I can tell this is because I didn't really make a dynamic programming solution with the memoized version. With Dynamic Programming our memo-table should be built around the change we need to give out, so you would have a table with all the possible amounts, not just the ones needed as remainders. This works with our depth-first version because we're starting at the leaf nodes and working back up the tree, but it doesn't work with the iterative version because we don't know ahead of time which remainders we need later on in the table. For example, if we can only use 1-cent and 10-cent coins and we need to make change for 11 cents, then we have to have an entry for 1-cent in the table before we get to figuring out the solution that includes the 10-cent piece. If we need to break 12 cents, though, we need to have an entry for 2 cents before we get to calculating the entry for 10 cents. But how do we know which values we need to enter into the table for the 1-cent piece before we get to the 10 cent piece? The solution is to fill out the table completely from 1 to the amount to make change for each coin before moving on to the next larger denomination. So, for our 12-cent example we need entries for $1 \ldots 12$ cents for the 1 cent piece filled out before we move on to calculating the entries for the 10 cent piece so that no matter what remainder we have when we reach the 10-cent piece we can look it up in the table (we actually only need up to 9 cents, but since the denominations aren't fixed doing the trimming once again adds to the complexity more than I'd like). So, I'll just fill out the table and we can see how it goes. def centaur(amount: int, denominations: list, table: dict=None) -> ChangeCountsTable: """Make change using the fewest coins Warning: To keep the answers matching the other functions it's assumed that the denominations are in decreasing order and the solution will also be match a decreasing order Args: money: the amount to change denominations: list of coin denominations in decreasing order table: memoization table (largest denomination, amount to change): best coin counts Returns: Number of each denomination needed to make the change, call_count, lookup table """ # the denominations need to be in ascending order, the opposite of the other methods denominations = list(sorted(denominations)) counter = 0 if not table: table = {denomination: [None] * (amount + 1) for denomination in denominations} for index, denomination in enumerate(denominations): table[denomination][0] = [0] * (index + 1) best = float("inf") best_counts = None for owed in range(amount + 1): for index, denomination in enumerate(denominations): counter += 1 count = owed//denomination remainder = owed - count * denomination counts = [count] if index > 0: counts = table[denominations[index - 1]][remainder] + counts table[denomination][owed] = counts if amount == owed and sum(counts) < best: best = sum(counts) best_counts = counts + [0] * (len(denominations) - index - 1) # put solution in descending order to match other functions best_counts.reverse() return ChangeCountsTable(best_counts, counter, table)  check_cases(centaur)  ### Plotting Since the iterative version works differently from the others the loop-counts aren't really easy to compare. centaur_counts = [centaur(amount, Denominations.US).run_count for amount in PLOT_AMOUNTS] GREEDY_FRAME["Iterative"] = centaur_counts melted = GREEDY_FRAME.melt(id_vars=["Amount to Change"], value_vars=["Greedy", "Cashier", "Brute-Force", "Memoized", "Iterative"], var_name="Change Method", value_name="Loop Count") chart = altair.Chart(melted).mark_line(point=True).encode( x="Amount to Change", y=altair.Y("Loop Count", axis=altair.Axis(tickMinStep=1)), tooltip=["Amount to Change", "Loop Count", "Change Method"], color="Change Method").properties( title="Coin Change Counts (Common U.S. Coins)", width=PlotParameters.width, height=PlotParameters.height, ) save_it(chart, "all-change-us-coins")  Since I'm filling out the whole table for the iterative function the loop is dependent on the amount of change owed, not the number of coin-denominations as with the other functions so it increases linearly. There's something I'm missing here as to why either the memoized-version is incorrect or there is an advantage to using the iterative version, but I haven't figured it out yet. The solutions seem to be the same in either case, though, and the linear growth for the number of loops is still tractable so I'll move on for now. ## Sources # Global Warming I: The Science and Modeling of Climate Change ## Syllabus ### Overview ### Heat, Light, and Energy ### First Climate Model ### Greenhouse Gasses and the Atmosphere ## Source # How much coal does it take to power a light bulb? ## Table of Contents ## How many Joules of energy does it take to power a 100 Watt light bulb for one day? Given: $1 Watt = 1 \frac{Joule}{Second}$ For a 100 Watt Bulb: $\require{cancel} (100 \cancel{Watt})\left(1 \frac{\frac{Joule}{Second}}{\cancel{Watt}}\right) = 100 \frac{Joule}{Second}$ For One Minute: $\left(100 \frac{Joule}{\cancel{Second}}\right) \left(60 \frac{\cancel{Second}}{Minute}\right) = 6{,}000 \frac{Joule}{Minute}$ For One Hour: $\left(6 \times 10^3 \frac{Joule}{\cancel{Minute}}\right)\left(60 \frac{\cancel{Minute}}{Hour}\right) = 3.6 \times 10^5 \frac{Joule}{Hour}$ For One Day: $\left(3.6 \times 10^5 \frac{Joule}{\cancel{Hour}}\right)\left(24 \frac{\cancel{Hour}}{Day}\right) = 8.64 \times 10^6 \frac{Joule}{Day}$ ## Double-Check The Math bulb = 100 # Watt (Joule/Second) per_minute = bulb * 60 # Second/Minute per_hour = per_minute * 60 # Minute/Hour per_day = per_hour * 24 # Hour/Day print(f"A 100 Watt Light bulb uses {per_day:0.2e} Joules per day.")  A 100 Watt Light bulb uses 8.64e+06 Joules per day.  ## Source # The Rod Cutting Problem ## The Rod Cutting Problem The Rod Cutting Problem begins with the premise that we have a metal rod (or some other kind of rod that's sellable) of a given length (n) and the price of a rod depends on the length so cutting the rod into shorter pieces will yield different total values - sometimes cutting a rod up make it more valuable than keeping it whole. Goal: Cut the rod in such a way that you maximize the total value of the pieces. To make this simpler we'll assume that lengths are always whole numbers (non-negative integers). ### An Example Say we have the following prices for different rod lengths.  Length 1 2 3 4 5 6 7 8 Value 1 5 8 9 10 17 17 20 Now suppose we have a rod that's 10 units in length. what combination of cuts will give us the best outcome? cuts values by piece total value 10 pieces of length 1 10 x 1 10 2 pieces of length 5 2 x 10 20 1 piece 8, 1 piece 2 20 + 5 25 1 piece 6, 1 piece 4 17 + 9 26 Out of these four possible combinations, 26 is the most we could get. ### Imports # from pypi from expects import contain_exactly, equal, expect  ## Brute Force In the previous example we could choose to make one cut, splitting our rod into 6 units and 4 units, and it would give us the highest value of the four combinations that we looked at. but what if we want to know the highest value out of all the combinations? Let's start by looking at a Brute-force approach to solving the problem. A Brute Force approach to finding the maximal cuts might do the following: • Try every possible combination of cuts • Evaluate the total values for each of the combinations • Pick the combination with the highest total value What will the runtime for this approach be? Let's start with two ideas: • Each cut location can have two values (cut, don't cut) so it is a binary value • Given a rod of length n, there are n - 1 possible cuts you can make (e.g. if your rod has length three, you can make two cuts (at length 1 and length 2) to give you three pieces) Although we are concerned with cuts, to get the total number of combinations you can think of each cut as a binary digit - (e.g. cut at length 1, don't cut at length 2 could be written as two digits: 10) and since the number of combinations of binary digits is $2^{\textit{number of digits}}$, this means that if we do an exhaustive exploration of the combinations we need to check \begin{align} 2^{n - 1} &= \left(2^{-1}\right)\left(2^n\right) \\ &= \frac{1}{2}2^n\\ T(n) &= \Theta\left(2^n\right) \end{align} So we have a problem with exponential growth. ### The CLRS Brute-Force Solution First, let's define some variables. Variable Representing $i$ Number of units of length of the rod $p_i$ Price for rod of length $i$ $r_n$ Revenue for rod segments with total length $n$ $q$ Maximum possible revenue for a particular length And now a little pseudocode. \begin{algorithm} \caption{Brute-Force Rod-Cutting} \begin{algorithmic} \INPUT Arrayp[0 \ldots n]$of prices mapped to length by index \INPUT Integer$n$the length of the rod \OUTPUT$q$the highest value for the combinations of rods totaling$n$in length \PROCEDURE{CutRodBruteForce}{$p, n$} \IF {$n = 0$} \RETURN 0 \ENDIF \STATE$q \gets -\inf$\FOR {$i \in \{1 \ldots n\}$} \STATE$q \gets $\textsc{Max}($q, p[i]$+ \textsc{CutRodBruteForce}($p, n - i$)) \ENDFOR \RETURN$q$\ENDPROCEDURE \end{algorithmic} \end{algorithm}  #### In Python def cut_rod_brute_force(prices: list, length: int, count: int=0) -> int: """Finds the maximum value you can extract from the rod Args: prices: map of length to price length: total length of the rod before cutting count: number of calls made Returns: the best total price you would get from cutting and selling the rod """ count += 1 if length == 0: return 0, count best_total = float("-inf") for next_cut in range(1, length + 1): next_total, count = cut_rod_brute_force( prices, length - next_cut, count=count) best_total = max(best_total, prices[next_cut] + next_total) return best_total, count  This naively assumes that there's an entry in prices for every length from 1 to the total length so price-lists need to be padded if there's missing lengths, as in the next example. #### The First Example The first thing we're going to do is to check the example given earlier, padding the price-list to make it have 10 entries. I originally had it just short-circuit if the list was shorter but then the counts were off by a little bit so I decided to get rid of that. It might make it slightly more efficient in certain cases, but brute-force isn't really what we're going for anyway. PRICES = [0, 1, 5, 8, 9, 10, 17, 17, 20, 0, 0] best_total, count = cut_rod_brute_force(PRICES, 10) expect(best_total).to(equal(27)) expect(count).to(equal(2**10)) print(f"Best Total Value: {best_total}") print(f"Count: {count:,}") print(f"Combinations: {2**10:,}")  Best Total Value: 27 Count: 1,024 Combinations: 1,024  So our actual best value is 27, not the 26 from the sub-set of combinations I used in the earlier example. ## Memoized Cut Rod The main reason why our brute-force version is so expensive is that it does all the calculations for every length over and over again when we test the different combinations. One way to get around this is by storing the values as they're calculated so that we can just look them up instead of repeating the calculations. This first version is very similar to the brute-force version except that the brute-force version makes a recursive call for every length we check, while for this memoized version we maintain an array to store previously calculated values and if the next one we want is in it we pull it from the array instead of making another recursive call. ### Cut Rod Memoized This first function is sort of a mask to make it look like the brute-force version. It sets up an empty memo table (as an array) and then passes it to the CutRodMemoizedAuxiliary function to do the actual calculations. \begin{algorithm} \caption{Memoized Rod-Cutting} \begin{algorithmic} \INPUT Array$p[0 \ldots n]$of prices mapped to length by index \INPUT Integer$n$the length of the rod to cut \OUTPUT$q$the highest value for the combinations of rods totaling$n$in length \PROCEDURE{CutRodMemoized}{$p, n$} \STATE Let$r[0 \ldots n]$be a new array. \FOR {$i \in \{0\ldots n\}$} \STATE$r[i]\gets -\infty$\ENDFOR \RETURN \textsc{CutRodMemoizedAuxiliary}($p, n, r$) \ENDPROCEDURE \end{algorithmic} \end{algorithm}  ### Cut Rod Memoized Auxiliary \begin{algorithm} \caption{Memoized Rod-Cutting Auxiliary} \begin{algorithmic} \INPUT Array$p[0 \ldots n]$of prices mapped to length by index \INPUT Integer$n$the length of the rod to cut \INPUT Array$r$of previously calculated values \OUTPUT$q$the highest value for the combinations of rods totaling$n$in length \PROCEDURE{CutRodMemoizedAuxiliary}{$p, n, r$} \IF {$r[n] \geq 0$} \RETURN$r[n]$\ENDIF \IF {$n=0$} \STATE$q \gets 0$\ELSE \STATE$q \gets -\infty$\FOR {$i \in \{1 \ldots n\}$} \STATE$q \gets$\textsc{Max}($q, p[i] + $\textsc{CutRodMemoizedAuxiliary}($p, n-i, r$)) \ENDFOR \ENDIF \STATE$r[n] \gets q\RETURN q \ENDPROCEDURE \end{algorithmic} \end{algorithm}  If you squint at CutRodMemoizedAuxiliary you might notice that it looks similar to the brute-force version except that there's an initial check to see if the value we want is already in our lookup-table and only makes the recursive call if it isn't. ### Python Version def cut_rod_memoized(prices: list, length: int) -> int: """Finds the maximum value for a rod after it has been cut up Args: prices: map of length to price length: the length of the rod to be cut up Returns: the maximum value that can be gained by cutting up and selling the rod """ table = [float("-inf")] * (length + 1) return cut_rod_memoized_auxiliary(prices, length, table)  def cut_rod_memoized_auxiliary(prices: list, length: int, best_values: list, count: int=0) -> int: """Find the maximum value from cutting up and selling rod Args: prices: map of length to price length: the length of the rod to be cut up best_values: lookup-table for previously calculated values (index is starting length) Returns: the maximum value that can be gained by cutting up and selling the rod """ count += 1 if best_values[length] >= 0: return best_values[length], count if length == 0: best_total = 0 else: best_total = float("-inf") for next_cut in range(1, length + 1): leftover = length - next_cut next_total, count = cut_rod_memoized_auxiliary(prices, leftover, best_values, count) best_total = max(best_total, prices[next_cut] + next_total) best_values[length] = best_total return best_total, count  best_total, count = cut_rod_memoized(PRICES, 4) print(f"Best Total Value: {best_total}") print(f"Count: {count}")  Best Total Value: 10 Count: 11  best_total, count = cut_rod_memoized(PRICES, 5) print(f"Best Total Value: {best_total}") print(f"Count: {count}")  Best Total Value: 13 Count: 16  ### Example best_total, count = cut_rod_memoized(PRICES, 10) expect(best_total).to(equal(27)) expect(count).to(equal(1 + (10 * 11)/2)) print(f"Best Total Value: {best_total}") print(f"Count: {count}")  Best Total Value: 27 Count: 56  We've gone from 1,024 calls to 56 calls, a pretty good improvement. The number of calls comes from the for loop plus one for the initial call. The for loop goes from 1 through the length of the rod, but passes in the difference between the starting length and the loop value. So if we start with a length of 4, the for-loop makes recursive calls using lengths of 4-1=3, 4-2=2, 4-3=1, 4-4=0. But then each of the calls goes through the for-loop as well (except for the base-case of 0). Since the first call of the for-loop is always one less than the starting length, we end up memoizing the values for all the starting lengths as we go so the subsequent calls don't need to go into the for-loop. So the number of calls we make equals $1 + 2 + \cdots + n$ plus one for the first call before the recursion starts. This means the runtime is $1 + \sum_{i=1}^n i = 1 + \frac{n(n+1)}{2} \Rightarrow O(n^2)$ So for our case with length 10, we have \begin{align} T(10) &= \frac{10(10 + 1)}{2} + 1\\ &= 56 \end{align} ## Non-Recursive Solution The memoized cut-rod solution is a top-down, depth-first search solution that uses recursion. We can eliminate the recursion altogether using a for-loop along with our look-up array. The trick is to make it a bottoms-up approach - that is to say that we start with the solutions for the smaller lengths and work up to the longer ones so that the look-up array always has the sub-problem values that we need to look up. \begin{algorithm} \caption{Bottoms-Up Rod-Cutting} \begin{algorithmic} \INPUT Arrayp[0 \ldots n]$of prices mapped to length by index \INPUT Integer$n$the length of the rod to cut \OUTPUT$q$the highest value for the combinations of rods totaling$n$in length \PROCEDURE{CutRodBottomsUp}{$p, n$} \STATE Let$r[0\ldots n]$be a new array. \STATE$r[0] \gets 0$\FOR {$j \in \{1 \ldots n\}$} \STATE$q \gets -\infty$\FOR {$i \in \{1 \ldots j\}$} \STATE$q \gets$\textsc{Max}($q, p[i] + r[j - i]$) \ENDFOR \STATE$r[j] \gets q$\ENDFOR \RETURN$r[n]$\ENDPROCEDURE \end{algorithmic} \end{algorithm}  ### Python Version def cut_rod_bottom_up(prices: list, length: int) -> tuple: """Find the maximum value for a rod after cutting Args: prices: map of length to price length: total length of rod to cut up Returns: best-value, count """ count = 1 best_values = [0] * (length + 1) for rod_length in range(1, length + 1): best_value_this_length = float("-inf") for cut_length in range(1, rod_length + 1): count += 1 leftover = rod_length - cut_length best_value_this_length = max( best_value_this_length, prices[cut_length] + best_values[leftover]) best_values[rod_length] = best_value_this_length return best_values[length], count  ### The Example Again best_total, count = cut_rod_bottom_up(PRICES, 10) expect(best_total).to(equal(27)) expect(count).to(equal(1 + 110/2)) print(f"Best Total Value: {best_total}") print(f"Count: {count}")  Best Total Value: 27 Count: 56  The runtime for this version is the same as the memoized version. It's sort of the backwards case - the inner for-loop runs 1 then 1, 2, then 1, 2, 3 up to the length of the rod, so the number of times it runs is $1 + 2 + 3 + \cdots + n$ while the memoized count goes $n + \cdots + 3 + 2 + 1$. In any case, the count ends up the same. $1 + \sum_{i=1}^n i = 1 + \frac{n(n+1)}{2} \Rightarrow O(n^2)$ ## Recovering the Cuts The previous functions all return the best value for a length but don't tell you the actual cuts that you need in order to get it. We can alter the function slightly to get both the best-revenue table and the cuts you need to use to get the best value. ### Pseudocode #### Extended Bottoms-Up \begin{algorithm} \caption{Extended Bottoms-Up Rod-Cutting} \begin{algorithmic} \INPUT Array$p[0 \ldots n]$of prices mapped to length by index \INPUT Integer$n$the length of the rod to cut \OUTPUT$r$the best revenue for each length \OUTPUT$s$list of next cut-lengths to use to get best revenue \PROCEDURE{CutRodBottomsUpExtended}{$p, n$} \STATE Let$r[0\ldots n]$and$s[0 \ldots n]$be new arrays. \STATE$r[0] \gets 0$\FOR {$j \in \{1 \ldots n\}$} \STATE$q \gets -\infty$\FOR {$i \in \{1 \ldots j\}$} \IF {$q < p[i] + r[j - i]$} \STATE$q \gets p[i] + r[j - i]$) \STATE$s[j] \gets i$\ENDIF \ENDFOR \STATE$r[j] \gets q$\ENDFOR \RETURN$(r, s)$\ENDPROCEDURE \end{algorithmic} \end{algorithm}  #### Printing the Solution To actually see the cuts we just need to retrieve the next cuts from s and calculate the remaining length after each cut. \begin{algorithm} \caption{Print Rod-Cutting Solution} \begin{algorithmic} \INPUT Array$p[0 \ldots n]$of prices mapped to length by index \INPUT Integer$n$the length of the rod to cut \PROCEDURE{Print-Cut-Rod-Solution}{$p, n$} \STATE$(r, s) \gets $\textsc{CutRodBottomsUpExtended}($p, n$) \WHILE {$n > 0$} \STATE \textsc{Print}($s[n]$) \STATE$n \gets n - s[n]\$
\ENDWHILE
\ENDPROCEDURE
\end{algorithmic}
\end{algorithm}


### In Python

def cut_rod_extended(prices: list, length: int) -> tuple:
"""Find the maximum values for a cut rod

Args:
prices: map of length to price
length: total length of rod to cut up

Returns:
best-revenues, best-lengths
"""
best_revenues = [0] * (length + 1)
best_cuts = [0] * (length + 1)

for next_length in range(1, length + 1):
best_revenue = float("-inf")
for next_cut in range(1, next_length + 1):
remaining_length = next_length - next_cut
next_revenue = prices[next_cut] + best_revenues[next_length - next_cut]
if best_revenue < next_revenue:
best_revenue = next_revenue
best_cuts[next_length] = next_cut
best_revenues[next_length] = best_revenue
return best_revenues, best_cuts

def print_solution(prices: list, length: int) -> tuple:
"""Solve and print the best cuts to maximize revenue

Args:
prices: list mapping length to price
length: the pre-cut length of the rod

Returns:
best_revenues, best_cuts
"""
revenues, cuts = cut_rod_extended(prices, length)
check, remaining_length = 0, length
output = []
while remaining_length > 0:
next_cut = cuts[remaining_length]
check += next_cut
output.append(str(next_cut))
remaining_length -= next_cut
print(f"Best Revenue for rod of length {length}: {revenues[length]}")
print(f"Cuts: ({', '.join(output)})")
expect(check).to(equal(length))
return revenues, cuts

print_solution(PRICES, 10)

Best Revenue for rod of length 10: 27
Cuts: (2, 2, 6)


### CLRS Example

CLRS_PRICES = [0, 1, 5, 8, 9, 10, 17, 17, 20, 24, 30]
revenues, cuts = print_solution(CLRS_PRICES, 10)
expect(revenues).to(contain_exactly(0, 1, 5, 8, 10, 13, 17, 18, 22, 25, 30))
expect(cuts).to(contain_exactly(0, 1, 2, 3, 2, 2, 6, 1, 2, 3, 10))

Best Revenue for rod of length 10: 30
Cuts: (10)

for length in range(1, 11):
print_solution(CLRS_PRICES, length)
print()

Best Revenue for rod of length 1: 1
Cuts: (1)

Best Revenue for rod of length 2: 5
Cuts: (2)

Best Revenue for rod of length 3: 8
Cuts: (3)

Best Revenue for rod of length 4: 10
Cuts: (2, 2)

Best Revenue for rod of length 5: 13
Cuts: (2, 3)

Best Revenue for rod of length 6: 17
Cuts: (6)

Best Revenue for rod of length 7: 18
Cuts: (1, 6)

Best Revenue for rod of length 8: 22
Cuts: (2, 6)

Best Revenue for rod of length 9: 25
Cuts: (3, 6)

Best Revenue for rod of length 10: 30
Cuts: (10)