# Rocket Fuel

June 17, 2020Each winter on those long and chilly evenings we test our wits on another edition of the Advent of Code. Competitive programming is not my cup of tea, but I find it quite satisfying to test my skills each year. However, the exhilarating effect of a flawless couple of days fades, as certainly, life will mess up our plans and the calendar unexpectedly fills up. There are also those days where the solution takes longer to appear than you had hoped for – when you suddenly realize you never got right that tiny detail in the problem description.

You miss a day or two and begin to backlog exercises. You promise yourself that on the weekend you will try to keep up. However, the odds are not in your favour, as each day we have a probability of around 15% of giving up the chase. Every year we expect it to become our victory this time around, but that success is unlikely to come.

The longest I have lasted was until day 20 back in 2017, after realising it was quite silly of me to spend Christmas Eve trying to catch up and further commit, instead of spending time with my loved ones.

Last year I only got to do a handful of problems. I skipped some problems and only did those that seem more fun. The last I did was Day 14 – Space Stoichiometry. I found it quite interesting and it entered directly to the top of my favorite Advent of Code puzzles.

## We are out of fuel

As we travel throughout the Solar System our spaceship begins to run low on `FUEL`

. We stop at the rings of Saturn and are tasked to convert `ORE`

into `FUEL`

using a nanofactory. As input we are given a list of reactions that the factory can do.
For example:

```
5 ORE => 7 A
3 A => 10 B
1 A, 3 B => 1 C
5 C, 2 B => 1 FUEL
```

In this problem:

- A reaction has a list of the required ingredients and which product it produces;
- Each reaction creates only one type of product;
- There is only one way of producing a given product;
- Reactions cannot be partially run;
- The output of a reaction can be used as the input of another reaction;
`ORE`

is the only raw-material available to us at the beginning;- All chemicals, except
`ORE`

, must be produced using a reaction; `FUEL`

is the end-product the we need to create.

The key to the problem is to establish the relationship between `ORE`

and `FUEL`

.

In Part I, we need to calculate the quantity of `ORE`

required to produce a certain quantity of `FUEL`

.
$$ \mathtt{ORE} = f(\mathtt{FUEL}) $$
In Part II, we need to calculate its inverse and find how many `FUEL`

can we generate with a given quantity of `ORE`

.
$$\mathtt{FUEL} = f^-1(\mathtt{ORE}). $$

## PART I - Making rocket fuel

### Stockpiling chemicals

We can, starting with the 1 `FUEL`

, replace each chemical by its ingredients until we turn all ingredients into their `ORE`

equivalents – like seeing a culinary video on reverse.

To do this, we need to keep a list of `needs`

which we seed with 1 `FUEL`

.
For each iteration of our algorithm, we process one chemical in the `needs`

list. The chemical is replaced by its ingredients. We must calculate how many times we need to run a reaction to create enough of the required chemical. We also need to keep a count of the excess chemicals produced – a stock.

If the ingredient being processed is `ORE`

we remove it from our list and add it to the total `ORE`

required.

In this implementation, we choose to process the last item of our `needs`

list at each iteration. Running this example on our input would tell us that we need 10 `ORE`

to generate 1 `FUEL`

. It took 6 iterations deplete our `needs`

list.

```
def ore_required_stock(reactions, qty):
ore_required = 0
needs = defaultdict(int, {'FUEL': qty})
stock = defaultdict(int)
def get_from_stock(qty, chemical):
in_stock = stock[chemical]
qty -= in_stock
stock[chemical] = abs(qty) if qty < 0 else 0
return in_stock - stock[chemical]
iterations = 0
while needs:
iterations += 1
chemical, qty_required = needs.popitem()
from_store = get_from_stock(qty_required, chemical)
qty_required -= from_store
qty_produced, ingredients = reactions[chemical]
n = math.ceil(qty_required/qty_produced)
stock[chemical] = qty_produced*n - qty_required
for qty_ingredient, ingredient in ingredients:
if ingredient == 'ORE':
ore_required += qty_ingredient*n
else:
needs[ingredient] += qty_ingredient*n
return (ore_required, iterations)
```

### Ordering chemicals

We can improve our approach by choosing which chemical we process at each iteration. Instead of removing the last item of our needs list, we should try to extract an item that minimizes the number of iterations.

But how could we decide on criteria for this?

In our previous implementation, we needed to process `B`

twice. However, if instead we processed `C`

we would generate more `B`

. Then, we can convert all the `B`

chemicals in a single iteration. This example encourages us to process first the chemicals that are not needed later as ingredients.

As our reactions list is a directed acyclic graph (DAG), we can do this is by finding its topological ordering. This will return the chemical process order that guarantees we process a product before its sub-products. In our example, this would guarantee that we process `C`

before its sub-product `B`

. Depending on the reaction list the topological order can be unique or not, but this as no practical impact on our solution.

We can find the topological ordering using depth-first search (DFS). After finding the deepest node and while returning backwards we list the nodes we visit. The topological ordering will then be the reverse of this list.

```
def topological_order(reactions):
visited_nodes = []
order = []
def dfs(reactions, node):
visited_nodes.append(node)
if node == 'ORE':
return
_, ingredients = reactions[node]
for _, ingredient in ingredients:
if ingredient not in visited_nodes:
dfs(reactions, ingredient)
order.append(node)
dfs(reactions, 'FUEL')
order.reverse()
order.append('ORE')
ordering = {element: index for index, element in enumerate(order)}
return ordering
```

We now need to make some changes to our `ore_required`

function to process the chemicals in our newly found ordering.

```
def ore_required_topological(reactions, ordering, qty):
ore_required = 0
needs = defaultdict(int, {'FUEL': qty})
iterations = 0
while needs:
iterations += 1
chemical = sorted(needs.keys(), key=lambda k: ordering[k])[0]
qty_required = needs.pop(chemical)
qty_produced, ingredients = reactions[chemical]
n = math.ceil(qty_required/qty_produced)
for qty_ingredient, ingredient in ingredients:
if ingredient == 'ORE':
ore_required += qty_ingredient*n
else:
needs[ingredient] += qty_ingredient*n
return (ore_required, iterations)
```

Using the topological ordering we do not need to consider any extra chemical produced in a reaction and add them to the stock, as each chemical will only be processed once.

Thus, the number of iterations in this solution will always be the same as the number of reactions given in the input. For our example reaction list we reduced the number of iterations from 6 to 4. However, this comes at the cost of calculating the topological order before processing any chemicals.

## PART 2 - We found a trillion ore nuggets

Part 2 task us to find how much `FUEL`

can we produce with one trillion `ORE`

.

As we know how many `ORE`

we require to produce `FUEL`

, but do not know the inverse, we can use successive approximation to find the solution.

We start by calculating the `ORE`

required for 1 `FUEL`

. If the result is less that 1 trillion we try again, but this time calculating the `ORE`

required to make 2 `FUEL`

. A so on until we reach our target. However, this approach could take a long time to compute.

We could speed-up the search by providing an initial estimation. As an initial estimate, we can divide a trillion by the amount of `ORE`

required to make 1 `FUEL`

.

$$\mathtt{estimation_{\nobreakspace under}} = \frac{\mathtt{1\nobreakspace trillion}}{f(1 \nobreakspace \mathtt{FUEL})}$$

This underestimates the solution, as extra byproducts generated by creating one `FUEL`

could be used as ingredients when making another. We can see this by calculating the `ORE/FUEL`

ratio when producing different amounts of `FUEL`

.

We follow the same procedure as before, evaluating the function and incrementing `FUEL`

on each iteration. Although this shortens the search, it can still be unacceptably long.

### Bisecting the search

For our search to finish quickly, we can use the bisect method, which has a time complexity of \( \log_2(\mathtt{N}) \), where \(\mathtt{N}\) is the size of our search interval. We can do this because `ore_required`

is a monotonic function.

The method works like a binary search and in each iteration, we split the search interval in two. We do this by calculating a midpoint between a low and high bounds and computing \( f( \mathtt{mid}) \). If it is higher than our target, we change our upper bound to the newly found midpoint. If it is lower, we change the lower bound to the midpoint.

As we are always using integers we need to increment or decrement the midpoint when assigning it to be the lower and upper bound respectively. This prevents an infinite loop.

```
def bisect(f, target, interval, args=()):
iterations = 0
low, high = interval
while low <= high:
iterations += 1
mid = low + (high-low)//2
f_mid = f(mid, *args)
if f_mid == target:
return (mid, iterations)
if f_mid < target:
low = mid + 1
else:
high = mid - 1
return low-1, iterations
```

To use this method, we need to estimate a lower and upper bound. For the lower bound we can use the underestimation we have calculated previously.

For the upper bound we can calculate the `ORE/FUEL`

ratio for several orders of magnitude of `FUEL`

, starting at log10(estimation_over) and stop when we produced more than a trillion `ORE`

.

Using this ratio we can then make another estimation. This time overestimating the solution.

$$\mathtt{estimation_{\nobreakspace over}} = \frac{\mathtt{1\nobreakspace trillion}}{\mathtt{lowest \nobreakspace ORE/FUEL \nobreakspace ratio}} $$

```
def estimate_fuel_produced(reactions, ore_qty):
ordering = topological_order(reactions)
ratio, _ = ore_required_topological(reactions, ordering, 1)
under_estimation = ore_qty//ratio
order = math.log10(under_estimation)
qty = 10**order
ore, _ = ore_required_topological(reactions, ordering, qty)
while ore < ore_qty:
print(qty)
order += 1
qty = 10**order
ore, _ = ore_required_topological(reactions, ordering, qty)
ratio = math.floor(ore/qty)
over_estimation = ore_qty//ratio
return (under_estimation,over_estimation), order
```

With our input example, for 1 trillion `ORE`

, we estimate to produce between 100000000000 and 142857142857 `FUEL`

.

We can then use our estimates as the interval of the bisect method.

```
TRILLION = int(1E12)
def f(x, reactions, ordering):
ore, _ = ore_required_topological(reactions, ordering, x)
return ore
ordering = topological_order(reactions)
args = (reactions, ordering)
estimations, _ = estimate_fuel_produced(reactions, TRILLION)
result, i = bisect(f, TRILLION, estimations, args)
```

For our example it will return 138613861385 `FUEL`

produced with 1 trillion `ORE`

. It took 35 iterations of the bisect method to find the solution.

### Intersecting the search

Another way to approximate our target value, its to use the secant method. This method is reminiscent of the Newton method but uses the slope between two values instead of the derivative.

Instead of evaluating one point and seeing if it is higher or lower than our target value like in the bisect method, the secant method evaluates the two extremes of the search interval and calculates the slope \(m\) between them – the secant.

$$m = \frac{f(x_1)-f(x_0)}{x_1-x_0}$$

Then it intersects the secant with our target.

$$x_2 = \frac{\mathtt{target} -f(x_1) + m x_1}{m}$$

If the value is below or above the target it changes the search interval accordingly.

This method converges faster to a solution than the bisect method. Also, in this problem the slope of our function is almost constant for high values of `FUEL`

produced. This makes the secant line between our two initial guesses almost coincident with the function. Thus, the search algorithm only needs to pinpoint the value of `FUEL`

which required over 1 trillion `ORE`

.

```
def secant(f, target, interval, args=()):
x0, x1 = interval
iterations = 0
while(x0 != x1):
iterations += 1
fx0 = f(x0, *args)
fx1 = f(x1, *args)
m = (fx1-fx0) / (x1-x0)
x2 = math.floor((target - fx1 + m*x1)/m)
x0, x1 = x1, x2
return x0, iterations
```

Using the secant method took 2 iteration instead of the 35 took by the bisect method. However, each iteration of the secant method runs the evaluation function twice.

## Liftoff

In summary we:

- Applied two unique approaches to solve each part of the problem;
- Answered part I by keeping a stock of the excess byproducts created on each reaction;
- Solved part I using topological ordering;
- Solved part II using two successive approximation algorithms: the bisect and the secant methods.

The complete code is available on Github.

### Further reading

- Scipy - Optimization and Root Finding documentation
- WilliamFiset - Topological Sort Algorithm - video

If you liked this article you can follow me on Twitter: @alexdsmartins or subscribe to the feed.