<jupyter_start><jupyter_text>Ants in Space! An introduction to the code in [beam_paco__gtoc5](https://github.com/lfsimoes/beam_paco__gtoc5)<NAME>2017-04 Table of Contents ---<jupyter_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 *<jupyter_output><empty_output><jupyter_text>--- Taking a look at our Python environment.<jupyter_code>%load_ext watermark
%watermark -v -m -p PyKEP,numpy,scipy,tqdm,pandas,matplotlib,seaborn
# https://github.com/rasbt/watermark<jupyter_output>CPython 3.5.3
IPython 5.3.0

PyKEP {'major': 1, 'minor': 3, 'bugfix': 9}
numpy 1.12.1
scipy 0.19.0
tqdm 4.11.2
pandas 0.19.2
matplotlib 2.0.0
seaborn 0.7.1

compiler   : MSC v.1900 64 bit (AMD64)
system     : Windows
release    : 10
machine    : AMD64
processor  : Intel64 Family 6 Model 94 Stepping 3, GenuineIntel
CPU cores  : 8
interpreter: 64bit<jupyter_text>---  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](http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/): `eil101` (symmetric, 101 cities, total distance in best known solution: [629](http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/STSP.html)).<jupyter_code>from urllib.request import urlretrieve
import gzip
urlretrieve('http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/tsp/eil101.tsp.gz', filename='eil101.tsp.gz');<jupyter_output><empty_output><jupyter_text>Load each city's (x, y) coordinates.<jupyter_code>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]<jupyter_output><empty_output><jupyter_text>Calculate distances matrix.<jupyter_code>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]<jupyter_output><empty_output><jupyter_text>Instantiate the TSP "path handler" with this distances matrix, and P-ACO, with its default parameters.<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Solve it.<jupyter_code>%time (quality, best) = aco.solve(nr_generations=100)
quality<jupyter_output>Wall time: 6.88 s<jupyter_text>Continue refining the solution for a few more generations.<jupyter_code>%time (quality, best) = aco.solve(nr_generations=400, reinitialize=False)
quality<jupyter_output>Wall time: 27.2 s<jupyter_text>Let's see what we found.<jupyter_code>xy = np.vstack([xy_locs[best], xy_locs[best][0]]) # to connect back to the start
line, = plt.plot(xy[:,0], xy[:,1], 'go-')<jupyter_output><empty_output><jupyter_text>---  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.<jupyter_code>t = mission_to_1st_asteroid(1712)
add_asteroid(t, 4893)<jupyter_output><empty_output><jupyter_text>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).<jupyter_code>score(t), final_mass(t), tof(t) * DAY2YEAR<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>resource_rating(t)<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>score(t) + resource_rating(t)<jupyter_output><empty_output><jupyter_text>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`.<jupyter_code>print(seq(t))
print(seq(t, incl_flyby=False))<jupyter_output>[0, 1712, 1712, 4893, 4893]
[0, 1712, 4893]<jupyter_text>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:1. asteroid ID2. spacecraft mass3. [epoch][epoch]4. the leg's $\Delta T$5. 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.[epoch]: https://en.wikipedia.org/wiki/Epoch_(astronomy)<jupyter_code>t[-1]<jupyter_output><empty_output><jupyter_text>Epochs are given here as [Modified Julian Dates (MJD)](https://en.wikipedia.org/wiki/Julian_dayVariants), and can be converted as:<jupyter_code>pk.epoch(t[-1][2], 'mjd')<jupyter_output><empty_output><jupyter_text>Greedy search In this section we perform a [Greedy search](https://en.wikipedia.org/wiki/Greedy_algorithm) 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.<jupyter_code>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))<jupyter_output>[0, 1712, 4893, 4028, 6939, 1059, 3295, 4006, 4981, 6437, 5344, 2090, 3993, 4186, 1130]<jupyter_text>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.<jupyter_code>t = mission_to_1st_asteroid(1712)<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>r.argsort()[:5]<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>[add_asteroid(copy(t), a) for a in r.argsort()[:5]]<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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))<jupyter_output>[0, 1712, 4893, 4028, 6939, 1059, 505, 6060, 3907, 5051, 2413, 6819, 913, 6638, 6723]<jupyter_text>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<jupyter_code>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]<jupyter_output><empty_output><jupyter_text>--- Generate the Table of Contents<jupyter_code>%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')
// https://github.com/kmahelona/ipython_notebook_goodies<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09491773308957953, 'fraction_numerical': 0.031078610603290677, 'mean_word_length': 5.2276867030965395, 'pattern_counts': {'":': 0, '<': 92, '<?xml version=': 0, '>': 92, 'https://': 10, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '346750', 'license': 'MIT', 'n_tokens_mistral': 4761, 'n_tokens_neox': 4351, 'n_words': 1520, 'path': 'usage_demos.ipynb', 'repo_name': 'AdityaChaudhary/beam_paco__gtoc5', 'size': '13683'}