<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Ants in Space! An introduction to the code in beam_paco__gtoc5
Luís F. Simões
2017-04
<h1 id="tocheading">Table of Contents</h1>
<div id="toc"></div>
Step1:  Taking a look at our Python environment.
Step2:  Solving a TSPLIB problem with P-ACO
In this section we show the steps to solve a Travelling Salesman Problem (TSP) instance with the Population-based Ant Colony Optimization (P-ACO) algorithm.
We'll use a TSP instance downloaded from the TSPLIB
Step3:  Load each city's (x, y) coordinates.
Step4:  Calculate distances matrix.
Step5:  Instantiate the TSP "path handler" with this distances matrix, and P-ACO, with its default parameters.
Step6:  Solve it.
Step7:  Continue refining the solution for a few more generations.
Step8:  Let's see what we found.
Step9:  Basic steps for assembling GTOC5 trajectories
The two primary functions for assembling a GTOC5 trajectory are here mission_to_1st_asteroid() and add_asteroid(). The first initializes the mission's data structure with the details of the Earth launch leg, that takes the spacecraft towards the mission's first asteroid. Subsequently, via multiple calls to add_asteroid(), the mission is extended with additional exploration targets. Each call to add_asteroid() creates a rendezvous leg towards the specified asteroid, immediately followed by a flyby of that same asteroid, and so increases the mission's overall score by 1.
Here's an example of a mission that launches towards asteroid 1712, and moves next to asteroid 4893. The True value returned by add_asteroid() indicates that a feasible transfer leg was found, and asteroid 4893 was therefore added to the mission.
Step10:  We can evaluate this trajectory with respect to its score (number of asteroids fully explored), final mass (in kg), and time of flight (here converted from days to years).
Step11:  An aggregation of the mission's mass and time costs can be obtained with resource_rating(). It measures the extent to which the mass and time budgets available for the mission have been depleted by the trajectory. It produces a value of 1.0 at the start of the mission, and a value of 0.0 when the mission has exhausted its 3500 kg of available mass, or its maximum duration of 15 years.
Step12:  As the score increments discretely by 1.0 with each added asteroid, and the resource rating evaluates mass and time available in a range of [0, 1], both can be combined to give a single-objective evaluation of the trajectory, that should be maximized
Step13:  Calling seq(), we can see either the full sequence of asteroids visited in each leg, or just the distinct asteroids visited in the mission. In this example, we see that the mission starts on Earth (id 0), performs a rendezvous with asteroid 1712, followed by a flyby of the same asteroid, and then repeats the pattern at asteroid 4893.
Step14:  The trajectory data structure built by mission_to_1st_asteroid() and add_asteroid() is a list of tuples summarizing the evolution of the spacecraft's state. It provides the minimal sufficient information from which a more detailed view can be reproduced, if so desired. Each tuple contains
Step15:  Epochs are given here as Modified Julian Dates (MJD), and can be converted as
Step16:  Greedy search
In this section we perform a Greedy search for a GTOC5 trajectory. We'll start by going to asteroid 1712. Then, and at every following step, we attempt to create legs towards all still available asteroids. Among equally-scored alternatives, we greedily pick the one with highest resource rating to adopt into the trajectory, and continue from there. Search stops when no feasible legs are found that can take us to another asteroid. This will happen either because no solutions were found that would allow for a leg to be created, or because adding a found solution would require the spacecraft to exceed the mission's mass or time budgets.
Step17:  Greedy search gave us a trajectory that is able to visit 14 distinct asteroids. However, by the 14th, the spacecraft finds itself unable to find a viable target to fly to next, even though it has 84.8 kg of mass still available (the spacecraft itself weighs 500 kg, so the mission cannot go below that value), and 2 years remain in its 15 year mission.
Phasing indicators
A big disadvantage of the approach followed above is the high computational cost of deciding which asteroid to go to next. It entails the optimization of up to 7075 legs, only to then pick a single one and discard all the other results.
An alternative is to use one of the indicators available in gtoc5/phasing.py. They can provide an indication of how likely a specific asteroid is to be an easily reachable target.
Step18:  We use here the (improved) orbital phasing indicator to rate destinations with respect to the estimated ΔV of hypothetical legs that would depart from dep_ast, at epoch dep_t, towards each possible asteroid, arriving there within leg_dT days. We don't know exactly how long the transfer time chosen by add_asteroid() would be, but we take leg_dT=125 days as reference transfer time.
Step19:  Below are the 5 asteroids the indicator estimates would be most easily reachable. As we've seen above in the results from the greedy search, asteroid 4893, here the 2nd best rated alternative, would indeed be the target reachable with lowest ΔV.
Step20:  The indicator is however not infallible. If we attempt to go from asteroid 1712 towards each of these asteroids, we find that none of them are actually reachable, except for 4893! Still, the indicator allows us to narrow our focus considerably.
Step21:  Armed with the indicator, we can reimplement the greedy search, so it will only optimize legs towards a number of top rated alternatives, and then proceed with the best out of those.
Step22:  We were able to find another score 14 trajectory, but this time it took us ~1 second, whereas before it was taking us 2 and a half minutes.
Finding a GTOC5 trajectory of score 17 with Beam Search
Step23:  Generate the Table of Contents
<ASSISTANT_TASK:>
Python Code:
# https://esa.github.io/pykep/
# https://github.com/esa/pykep
# https://pypi.python.org/pypi/pykep/
import PyKEP as pk
import numpy as np
from tqdm import tqdm, trange
import matplotlib.pylab as plt
%matplotlib inline
import seaborn as sns
plt.rcParams['figure.figsize'] = 10, 8
from gtoc5 import *
from gtoc5.multiobjective import *
from gtoc5.phasing import *
from paco import *
from paco_traj import *
from experiments import *
from experiments__paco import *
Explanation: Ants in Space! An introduction to the code in beam_paco__gtoc5
Luís F. Simões
2017-04
<h1 id="tocheading">Table of Contents</h1>
<div id="toc"></div>
End of explanation
%load_ext watermark
%watermark -v -m -p PyKEP,numpy,scipy,tqdm,pandas,matplotlib,seaborn
# https://github.com/rasbt/watermark
Explanation: Taking a look at our Python environment.
End of explanation
from urllib.request import urlretrieve
import gzip
urlretrieve('http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/tsp/eil101.tsp.gz', filename='eil101.tsp.gz');
Explanation: Solving a TSPLIB problem with P-ACO
In this section we show the steps to solve a Travelling Salesman Problem (TSP) instance with the Population-based Ant Colony Optimization (P-ACO) algorithm.
We'll use a TSP instance downloaded from the TSPLIB: eil101 (symmetric, 101 cities, total distance in best known solution: 629).
End of explanation
with gzip.open('eil101.tsp.gz') as f:
    xy_locs = np.loadtxt(f, skiprows=6, usecols=(1,2), comments='EOF', dtype=np.int)
nr_cities = len(xy_locs)
xy_locs[:5]
Explanation: Load each city's (x, y) coordinates.
End of explanation
distances = np.zeros((nr_cities, nr_cities))
for a in range(nr_cities):
    for b in range(a, nr_cities):
        distances[a,b] = distances[b,a] = np.linalg.norm(xy_locs[a] - xy_locs[b])
distances[:4, :4]
Explanation: Calculate distances matrix.
End of explanation
rng, seed = initialize_rng(seed=None)
print('Seed:', seed)
path_handler = tsp_path(distances, random_state=rng)
aco = paco(path_handler.nr_nodes, path_handler, random_state=rng)
Explanation: Instantiate the TSP "path handler" with this distances matrix, and P-ACO, with its default parameters.
End of explanation
%time (quality, best) = aco.solve(nr_generations=100)
quality
Explanation: Solve it.
End of explanation
%time (quality, best) = aco.solve(nr_generations=400, reinitialize=False)
quality
Explanation: Continue refining the solution for a few more generations.
End of explanation
xy = np.vstack([xy_locs[best], xy_locs[best][0]]) # to connect back to the start
line, = plt.plot(xy[:,0], xy[:,1], 'go-')
Explanation: Let's see what we found.
End of explanation
t = mission_to_1st_asteroid(1712)
add_asteroid(t, 4893)
Explanation: Basic steps for assembling GTOC5 trajectories
The two primary functions for assembling a GTOC5 trajectory are here mission_to_1st_asteroid() and add_asteroid(). The first initializes the mission's data structure with the details of the Earth launch leg, that takes the spacecraft towards the mission's first asteroid. Subsequently, via multiple calls to add_asteroid(), the mission is extended with additional exploration targets. Each call to add_asteroid() creates a rendezvous leg towards the specified asteroid, immediately followed by a flyby of that same asteroid, and so increases the mission's overall score by 1.
Here's an example of a mission that launches towards asteroid 1712, and moves next to asteroid 4893. The True value returned by add_asteroid() indicates that a feasible transfer leg was found, and asteroid 4893 was therefore added to the mission.
End of explanation
score(t), final_mass(t), tof(t) * DAY2YEAR
Explanation: We can evaluate this trajectory with respect to its score (number of asteroids fully explored), final mass (in kg), and time of flight (here converted from days to years).
End of explanation
resource_rating(t)
Explanation: An aggregation of the mission's mass and time costs can be obtained with resource_rating(). It measures the extent to which the mass and time budgets available for the mission have been depleted by the trajectory. It produces a value of 1.0 at the start of the mission, and a value of 0.0 when the mission has exhausted its 3500 kg of available mass, or its maximum duration of 15 years.
End of explanation
score(t) + resource_rating(t)
Explanation: As the score increments discretely by 1.0 with each added asteroid, and the resource rating evaluates mass and time available in a range of [0, 1], both can be combined to give a single-objective evaluation of the trajectory, that should be maximized:
End of explanation
print(seq(t))
print(seq(t, incl_flyby=False))
Explanation: Calling seq(), we can see either the full sequence of asteroids visited in each leg, or just the distinct asteroids visited in the mission. In this example, we see that the mission starts on Earth (id 0), performs a rendezvous with asteroid 1712, followed by a flyby of the same asteroid, and then repeats the pattern at asteroid 4893.
End of explanation
t[-1]
Explanation: The trajectory data structure built by mission_to_1st_asteroid() and add_asteroid() is a list of tuples summarizing the evolution of the spacecraft's state. It provides the minimal sufficient information from which a more detailed view can be reproduced, if so desired. Each tuple contains:
asteroid ID
spacecraft mass
epoch
the leg's $\Delta T$
the leg's $\Delta V$
The mass and epoch values correspond to the state at the given asteroid, at the end of a rendezvous or self-fly-by leg, after deploying the corresponding payload. The $\Delta T$ and $\Delta V$ values refer to that leg that just ended.
End of explanation
pk.epoch(t[-1][2], 'mjd')
Explanation: Epochs are given here as Modified Julian Dates (MJD), and can be converted as:
End of explanation
import os
from copy import copy
def greedy_step(traj):
    traj_asts = set(seq(traj, incl_flyby=False))
    progress_bar_args = dict(leave=False, file=os.sys.stdout, desc='attempting score %d' % (score(traj)+1))
    
    extended = []
    for a in trange(len(asteroids), **progress_bar_args):
        if a in traj_asts:
            continue
        tt = copy(traj)
        if add_asteroid(tt, next_ast=a, use_cache=False):
            extended.append(tt)
    
    return max(extended, key=resource_rating, default=[])
# measure time taken at one level to attempt legs towards all asteroids (that aren't already in the traj.)
%time _ = greedy_step(mission_to_1st_asteroid(1712))
def greedy_search(first_ast):
    t = mission_to_1st_asteroid(first_ast)
    while True:
        tt = greedy_step(t)
        if tt == []:
            # no more asteroids could be added
            return t
        t = tt
%time T = greedy_search(first_ast=1712)
score(T), resource_rating(T), final_mass(T), tof(T) * DAY2YEAR
print(seq(T, incl_flyby=False))
Explanation: Greedy search
In this section we perform a Greedy search for a GTOC5 trajectory. We'll start by going to asteroid 1712. Then, and at every following step, we attempt to create legs towards all still available asteroids. Among equally-scored alternatives, we greedily pick the one with highest resource rating to adopt into the trajectory, and continue from there. Search stops when no feasible legs are found that can take us to another asteroid. This will happen either because no solutions were found that would allow for a leg to be created, or because adding a found solution would require the spacecraft to exceed the mission's mass or time budgets.
End of explanation
t = mission_to_1st_asteroid(1712)
Explanation: Greedy search gave us a trajectory that is able to visit 14 distinct asteroids. However, by the 14th, the spacecraft finds itself unable to find a viable target to fly to next, even though it has 84.8 kg of mass still available (the spacecraft itself weighs 500 kg, so the mission cannot go below that value), and 2 years remain in its 15 year mission.
Phasing indicators
A big disadvantage of the approach followed above is the high computational cost of deciding which asteroid to go to next. It entails the optimization of up to 7075 legs, only to then pick a single one and discard all the other results.
An alternative is to use one of the indicators available in gtoc5/phasing.py. They can provide an indication of how likely a specific asteroid is to be an easily reachable target.
End of explanation
r = rate__orbital_2(dep_ast=t[-1][0], dep_t=t[-1][2], leg_dT=125)
r[seq(t)] = np.inf  # (exclude bodies already visited)
Explanation: We use here the (improved) orbital phasing indicator to rate destinations with respect to the estimated ΔV of hypothetical legs that would depart from dep_ast, at epoch dep_t, towards each possible asteroid, arriving there within leg_dT days. We don't know exactly how long the transfer time chosen by add_asteroid() would be, but we take leg_dT=125 days as reference transfer time.
End of explanation
r.argsort()[:5]
Explanation: Below are the 5 asteroids the indicator estimates would be most easily reachable. As we've seen above in the results from the greedy search, asteroid 4893, here the 2nd best rated alternative, would indeed be the target reachable with lowest ΔV.
End of explanation
[add_asteroid(copy(t), a) for a in r.argsort()[:5]]
Explanation: The indicator is however not infallible. If we attempt to go from asteroid 1712 towards each of these asteroids, we find that none of them are actually reachable, except for 4893! Still, the indicator allows us to narrow our focus considerably.
End of explanation
def narrowed_greedy_step(traj, top=10):
    traj_asts = set(seq(traj, incl_flyby=False))
    
    extended = []
    ratings = rate__orbital_2(dep_ast=traj[-1][0], dep_t=traj[-1][2], leg_dT=125)
    for a in ratings.argsort()[:top]:
        if a in traj_asts:
            continue
        tt = copy(traj)
        if add_asteroid(tt, next_ast=a, use_cache=False):
            extended.append(tt)
    
    return max(extended, key=resource_rating, default=[])
def narrowed_greedy_search(first_ast, **kwargs):
    t = mission_to_1st_asteroid(first_ast)
    while True:
        tt = narrowed_greedy_step(t, **kwargs)
        if tt == []:
            # no more asteroids could be added
            return t
        t = tt
# measure time taken at one level to attempt legs towards the best `top` asteroids
%time _ = narrowed_greedy_step(mission_to_1st_asteroid(1712), top=10)
%time T = narrowed_greedy_search(first_ast=1712, top=10)
score(T), resource_rating(T), final_mass(T), tof(T) * DAY2YEAR
print(seq(T, incl_flyby=False))
Explanation: Armed with the indicator, we can reimplement the greedy search, so it will only optimize legs towards a number of top rated alternatives, and then proceed with the best out of those.
End of explanation
gtoc_ph = init__path_handler(multiobj_evals=True)
# configuring Beam P-ACO to behave as a deterministic multi-objective Beam Search
_args = {
    'beam_width': 20,
    'branch_factor': 250,
    'alpha': 0.0,        # 0.0: no pheromones used
    'beta': 1.0,
    'prob_greedy': 1.0,  # 1.0: deterministic, greedy branching decisions
    }
bpaco = beam_paco_pareto(nr_nodes=len(asteroids), path_handler=gtoc_ph, random_state=None, **_args)
# start the search
# given we're running the algoritm in deterministic mode, we execute it for a single generation
%time best_pf = bpaco.solve(nr_generations=1)
# being this a `_pareto` class, .best returns a Pareto front
# pick the first solution from the Pareto front
best_eval, best = best_pf[0]
# Evaluation of the best found solution
# (score, mass consumed, time of flight)
best_eval
# sequence of asteroids visited (0 is the Earth)
print(seq(best, incl_flyby=False))
# mission data structure, up to the full scoring of the first two asteroids
best[:5]
Explanation: We were able to find another score 14 trajectory, but this time it took us ~1 second, whereas before it was taking us 2 and a half minutes.
Finding a GTOC5 trajectory of score 17 with Beam Search
End of explanation
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')
// https://github.com/kmahelona/ipython_notebook_goodies
Explanation: Generate the Table of Contents
End of explanation
<END_TASK>