Rocket Fuel

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.