text:
# Neural Sequence Distance Embeddings

[](https://colab.research.google.com/github/gcorso/NeuroSEED/blob/master/tutorial/NeuroSEED.ipynb)

The improvement of data-dependent heuristics and representation for biological sequences is a critical requirement to fully exploit the recent technological and scientific advancements for human microbiome analysis. This notebook presents Neural Sequence Distance Embeddings (NeuroSEED), a novel framework to embed biological sequences in geometric vector spaces that unifies recently proposed approaches. We demonstrate its capacity by presenting different ways it can be applied to the tasks of edit distance approximation, closest string retrieval, hierarchical clustering and multiple sequence alignment. In particular, the hyperbolic space is shown to be a key component to embed biological sequences and obtain competitive heuristics. Benchmarked with common bioinformatics and machine learning baselines, the proposed approaches display significant accuracy and/or runtime improvements on real-world datasets formed by sequences from samples of the human microbiome.

Figure 1: On the left, a diagram of the NeuroSEED underlying idea: embed sequences in vector spaces preserving the edit distance between them and then extract information from the vector space. On the right, an example of the hierarchical clustering produced on the Poincarè disk from the P53 tumour protein from 20 different organisms.

## Introduction and Motivation

### Motivation

Dysfunctions of the human microbiome (Morgan & Huttenhower, 2012) have been linked to many serious diseases ranging from diabetes and antibiotic resistance to inflammatory bowel disease. Its usage as a biomarker for the diagnosis and as a target for interventions is a very active area of research. Thanks to the advances in sequencing technologies, modern analysis relies on sequence reads that can be generated relatively quickly. However, to fully exploit the potential of these advances for personalised medicine, the computational methods used in the analysis have to significantly improve in terms of speed and accuracy.

Figure 2: Traditional approach to the analysis of the 16S rRNA sequences from the microbiome.

### Problem

While the number of available biological sequences has been growing exponentially over the past decades, most of the problems related to string matching have not been addressed by the recent advances in machine learning. Classical algorithms are data-independent and, therefore, cannot exploit the low-dimensional manifold assumption that characterises real-world data. Exploiting the available data to produce data-dependent heuristics and representations would greatly accelerate large-scale analyses that are critical to microbiome analysis and other biological research.

Unlike most tasks in computer vision and NLP, string matching problems are typically formulated as combinatorial optimisation problems. These discrete formulations do not fit well with the current deep learning approaches causing these problems to be left mostly unexplored by the community. Current supervised learning methods also suffer from the lack of labels that characterises many downstream applications with biological sequences. On the other hand, common self-supervised learning approaches, very successful in NLP, are less effective in the biological context where relations tend to be per-sequence rather than per-token (McDermott et al. 2021).

### Neural Sequence Distance Embedding

In this notebook, we present Neural Sequence Distance Embeddings (NeuroSEED), a general framework to produce representations for biological sequences where the distance in the embedding space is correlated with the evolutionary distance between sequences. This control over the geometric interpretation of the representation space enables the use of geometrical data processing tools for the analysis of the spectrum of sequences.

Figure 3: The key idea of NeuroSEED is to learn an encoder function that preserves distances between the sequence and vector space.

Examining the task of embedding sequences to preserve the edit distance reveals the importance of data-dependent approaches and of using a geometry that matches well the underlying distribution in the data analysed. For biological datasets, that have an implicit hierarchical structure given by evolution, the hyperbolic space provides significant improvement.

We show the potential of the framework by analysing three fundamental tasks in bioinformatics: closest string retrieval, hierarchical clustering and multiple sequence alignment. For all tasks, relatively simple unsupervised approaches using NeuroSEED encoders significantly outperform data-independent heuristics in terms of accuracy and/or runtime. In the paper (preprint will be available soon) and the [complete repository](https://github.com/gcorso/NeuroSEED) we also present more complex geometrical approaches to hierarchical clustering and multiple sequence alignment.

## 2. Analysis

To improve readability and limit the size of the notebook we make use of some subroutines in the [official repository](https://github.com/gcorso/NeuroSEED) for the research project. The code in the notebook is our best effort to convey the promising application of hyperbolic geometry to this novel research direction.

Install and import the required packages.

```python

!pip install geoopt==0.3.1
!pip install torch==1.7.1+cu101 torchvision==0.8.2+cu101 torchaudio==0.7.2 -f https://download.pytorch.org/whl/torch_stable.html
!pip3 install geomstats==2.2
#!apt install clustalw
#pip install biopython
#!pip install python-Levenshtein
#!pip install Cython
#!pip install networkx
#!pip install tqdm
#!pip install gdown
#!cd hierarchical_clustering/relaxed/mst; python setup.py build_ext --inplace; cd ../unionfind; python setup.py build_ext --inplace; cd ..; cd ..; cd ..;
os.environ['GEOMSTATS_BACKEND'] = 'pytorch'
```

    Requirement already satisfied: geoopt==0.3.1 in ./neuroseed/lib/python3.8/site-packages (0.3.1)
    Requirement already satisfied: torch>=1.4.0 in ./neuroseed/lib/python3.8/site-packages (from geoopt==0.3.1) (1.7.1+cu101)
    Requirement already satisfied: numpy in ./neuroseed/lib/python3.8/site-packages (from geoopt==0.3.1) (1.22.3)
    Requirement already satisfied: typing-extensions in ./neuroseed/lib/python3.8/site-packages (from torch>=1.4.0->geoopt==0.3.1) (4.1.1)
    WARNING: You are using pip version 21.1.1; however, version 22.0.4 is available.
    You should consider upgrading via the '/home/noga/NeuroSEED/neuroseed/bin/python3 -m pip install --upgrade pip' command.
    Looking in links: https://download.pytorch.org/whl/torch_stable.html
    Requirement already satisfied: torch==1.7.1+cu101 in ./neuroseed/lib/python3.8/site-packages (1.7.1+cu101)
    Requirement already satisfied: torchvision==0.8.2+cu101 in ./neuroseed/lib/python3.8/site-packages (0.8.2+cu101)
    Requirement already satisfied: torchaudio==0.7.2 in ./neuroseed/lib/python3.8/site-packages (0.7.2)
    Requirement already satisfied: typing-extensions in ./neuroseed/lib/python3.8/site-packages (from torch==1.7.1+cu101) (4.1.1)
    Requirement already satisfied: numpy in ./neuroseed/lib/python3.8/site-packages (from torch==1.7.1+cu101) (1.22.3)
    Requirement already satisfied: pillow>=4.1.1 in ./neuroseed/lib/python3.8/site-packages (from torchvision==0.8.2+cu101) (9.0.1)
    WARNING: You are using pip version 21.1.1; however, version 22.0.4 is available.
    You should consider upgrading via the '/home/noga/NeuroSEED/neuroseed/bin/python3 -m pip install --upgrade pip' command.
    Requirement already satisfied: geomstats==2.2 in ./neuroseed/lib/python3.8/site-packages (2.2.0)
    Requirement already satisfied: scikit-learn>=0.22.1 in ./neuroseed/lib/python3.8/site-packages (from geomstats==2.2) (1.0.2)
    Requirement already satisfied: scipy>=1.4.1 in ./neuroseed/lib/python3.8/site-packages (from geomstats==2.2) (1.8.0)
    Requirement already satisfied: matplotlib in ./neuroseed/lib/python3.8/site-packages (from geomstats==2.2) (3.5.1)
    Requirement already satisfied: joblib==0.14.1 in ./neuroseed/lib/python3.8/site-packages (from geomstats==2.2) (0.14.1)
    Requirement already satisfied: autograd==1.3 in ./neuroseed/lib/python3.8/site-packages (from geomstats==2.2) (1.3)
    Requirement already satisfied: jupyter in ./neuroseed/lib/python3.8/site-packages (from geomstats==2.2) (1.0.0)
    Requirement already satisfied: numpy>=1.18.1 in ./neuroseed/lib/python3.8/site-packages (from geomstats==2.2) (1.22.3)
    Requirement already satisfied: future>=0.15.2 in ./neuroseed/lib/python3.8/site-packages (from autograd==1.3->geomstats==2.2) (0.18.2)
    Requirement already satisfied: threadpoolctl>=2.0.0 in ./neuroseed/lib/python3.8/site-packages (from scikit-learn>=0.22.1->geomstats==2.2) (3.1.0)
    Requirement already satisfied: notebook in ./neuroseed/lib/python3.8/site-packages (from jupyter->geomstats==2.2) (6.4.8)
    Requirement already satisfied: ipykernel in ./neuroseed/lib/python3.8/site-packages (from jupyter->geomstats==2.2) (6.9.1)
    Requirement already satisfied: nbconvert in ./neuroseed/lib/python3.8/site-packages (from jupyter->geomstats==2.2) (6.4.2)
    Requirement already satisfied: qtconsole in ./neuroseed/lib/python3.8/site-packages (from jupyter->geomstats==2.2) (5.2.2)
    Requirement already satisfied: ipywidgets in ./neuroseed/lib/python3.8/site-packages (from jupyter->geomstats==2.2) (7.6.5)
    Requirement already satisfied: jupyter-console in ./neuroseed/lib/python3.8/site-packages (from jupyter->geomstats==2.2) (6.4.3)
    Requirement already satisfied: nest-asyncio in ./neuroseed/lib/python3.8/site-packages (from ipykernel->jupyter->geomstats==2.2) (1.5.4)
    Requirement already satisfied: jupyter-client<8.0 in ./neuroseed/lib/python3.8/site-packages (from ipykernel->jupyter->geomstats==2.2) (7.1.2)
    Requirement already satisfied: ipython>=7.23.1 in ./neuroseed/lib/python3.8/site-packages (from ipykernel->jupyter->geomstats==2.2) (8.1.1)
    Requirement already satisfied: matplotlib-inline<0.2.0,>=0.1.0 in ./neuroseed/lib/python3.8/site-packages (from ipykernel->jupyter->geomstats==2.2) (0.1.3)
    Requirement already satisfied: traitlets<6.0,>=5.1.0 in ./neuroseed/lib/python3.8/site-packages (from ipykernel->jupyter->geomstats==2.2) (5.1.1)
    Requirement already satisfied: tornado<7.0,>=4.2 in ./neuroseed/lib/python3.8/site-packages (from ipykernel->jupyter->geomstats==2.2) (6.1)
    Requirement already satisfied: debugpy<2.0,>=1.0.0 in ./neuroseed/lib/python3.8/site-packages (from ipykernel->jupyter->geomstats==2.2) (1.5.1)
    Requirement already satisfied: backcall in ./neuroseed/lib/python3.8/site-packages (from ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (0.2.0)
    Requirement already satisfied: prompt-toolkit!=3.0.0,!=3.0.1,<3.1.0,>=2.0.0 in ./neuroseed/lib/python3.8/site-packages (from ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (3.0.28)
    Requirement already satisfied: jedi>=0.16 in ./neuroseed/lib/python3.8/site-packages (from ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (0.18.1)
    Requirement already satisfied: pickleshare in ./neuroseed/lib/python3.8/site-packages (from ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (0.7.5)
    Requirement already satisfied: stack-data in ./neuroseed/lib/python3.8/site-packages (from ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (0.2.0)
    Requirement already satisfied: pexpect>4.3 in ./neuroseed/lib/python3.8/site-packages (from ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (4.8.0)
    Requirement already satisfied: decorator in ./neuroseed/lib/python3.8/site-packages (from ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (5.1.1)
    Requirement already satisfied: pygments>=2.4.0 in ./neuroseed/lib/python3.8/site-packages (from ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (2.11.2)
    Requirement already satisfied: setuptools>=18.5 in ./neuroseed/lib/python3.8/site-packages (from ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (56.0.0)
    Requirement already satisfied: parso<0.9.0,>=0.8.0 in ./neuroseed/lib/python3.8/site-packages (from jedi>=0.16->ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (0.8.3)
    Requirement already satisfied: entrypoints in ./neuroseed/lib/python3.8/site-packages (from jupyter-client<8.0->ipykernel->jupyter->geomstats==2.2) (0.4)
    Requirement already satisfied: pyzmq>=13 in ./neuroseed/lib/python3.8/site-packages (from jupyter-client<8.0->ipykernel->jupyter->geomstats==2.2) (22.3.0)
    Requirement already satisfied: python-dateutil>=2.1 in ./neuroseed/lib/python3.8/site-packages (from jupyter-client<8.0->ipykernel->jupyter->geomstats==2.2) (2.8.2)
    Requirement already satisfied: jupyter-core>=4.6.0 in ./neuroseed/lib/python3.8/site-packages (from jupyter-client<8.0->ipykernel->jupyter->geomstats==2.2) (4.9.2)
    Requirement already satisfied: ptyprocess>=0.5 in ./neuroseed/lib/python3.8/site-packages (from pexpect>4.3->ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (0.7.0)
    Requirement already satisfied: wcwidth in ./neuroseed/lib/python3.8/site-packages (from prompt-toolkit!=3.0.0,!=3.0.1,<3.1.0,>=2.0.0->ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (0.2.5)
    Requirement already satisfied: six>=1.5 in ./neuroseed/lib/python3.8/site-packages (from python-dateutil>=2.1->jupyter-client<8.0->ipykernel->jupyter->geomstats==2.2) (1.16.0)
    Requirement already satisfied: nbformat>=4.2.0 in ./neuroseed/lib/python3.8/site-packages (from ipywidgets->jupyter->geomstats==2.2) (5.2.0)
    Requirement already satisfied: widgetsnbextension~=3.5.0 in ./neuroseed/lib/python3.8/site-packages (from ipywidgets->jupyter->geomstats==2.2) (3.5.2)
    Requirement already satisfied: ipython-genutils~=0.2.0 in ./neuroseed/lib/python3.8/site-packages (from ipywidgets->jupyter->geomstats==2.2) (0.2.0)
    Requirement already satisfied: jupyterlab-widgets>=1.0.0 in ./neuroseed/lib/python3.8/site-packages (from ipywidgets->jupyter->geomstats==2.2) (1.0.2)
    Requirement already satisfied: jsonschema!=2.5.0,>=2.4 in ./neuroseed/lib/python3.8/site-packages (from nbformat>=4.2.0->ipywidgets->jupyter->geomstats==2.2) (4.4.0)
    Requirement already satisfied: importlib-resources>=1.4.0 in ./neuroseed/lib/python3.8/site-packages (from jsonschema!=2.5.0,>=2.4->nbformat>=4.2.0->ipywidgets->jupyter->geomstats==2.2) (5.4.0)
    Requirement already satisfied: pyrsistent!=0.17.0,!=0.17.1,!=0.17.2,>=0.14.0 in ./neuroseed/lib/python3.8/site-packages (from jsonschema!=2.5.0,>=2.4->nbformat>=4.2.0->ipywidgets->jupyter->geomstats==2.2) (0.18.1)
    Requirement already satisfied: attrs>=17.4.0 in ./neuroseed/lib/python3.8/site-packages (from jsonschema!=2.5.0,>=2.4->nbformat>=4.2.0->ipywidgets->jupyter->geomstats==2.2) (21.4.0)
    Requirement already satisfied: zipp>=3.1.0 in ./neuroseed/lib/python3.8/site-packages (from importlib-resources>=1.4.0->jsonschema!=2.5.0,>=2.4->nbformat>=4.2.0->ipywidgets->jupyter->geomstats==2.2) (3.7.0)
    Requirement already satisfied: Send2Trash>=1.8.0 in ./neuroseed/lib/python3.8/site-packages (from notebook->jupyter->geomstats==2.2) (1.8.0)
    Requirement already satisfied: terminado>=0.8.3 in ./neuroseed/lib/python3.8/site-packages (from notebook->jupyter->geomstats==2.2) (0.13.3)
    Requirement already satisfied: prometheus-client in ./neuroseed/lib/python3.8/site-packages (from notebook->jupyter->geomstats==2.2) (0.13.1)
    Requirement already satisfied: argon2-cffi in ./neuroseed/lib/python3.8/site-packages (from notebook->jupyter->geomstats==2.2) (21.3.0)
    Requirement already satisfied: jinja2 in ./neuroseed/lib/python3.8/site-packages (from notebook->jupyter->geomstats==2.2) (3.0.3)
    Requirement already satisfied: argon2-cffi-bindings in ./neuroseed/lib/python3.8/site-packages (from argon2-cffi->notebook->jupyter->geomstats==2.2) (21.2.0)
    Requirement already satisfied: cffi>=1.0.1 in ./neuroseed/lib/python3.8/site-packages (from argon2-cffi-bindings->argon2-cffi->notebook->jupyter->geomstats==2.2) (1.15.0)
    Requirement already satisfied: pycparser in ./neuroseed/lib/python3.8/site-packages (from cffi>=1.0.1->argon2-cffi-bindings->argon2-cffi->notebook->jupyter->geomstats==2.2) (2.21)
    Requirement already satisfied: MarkupSafe>=2.0 in ./neuroseed/lib/python3.8/site-packages (from jinja2->notebook->jupyter->geomstats==2.2) (2.1.0)
    Requirement already satisfied: cycler>=0.10 in ./neuroseed/lib/python3.8/site-packages (from matplotlib->geomstats==2.2) (0.11.0)
    Requirement already satisfied: kiwisolver>=1.0.1 in ./neuroseed/lib/python3.8/site-packages (from matplotlib->geomstats==2.2) (1.3.2)
    Requirement already satisfied: pillow>=6.2.0 in ./neuroseed/lib/python3.8/site-packages (from matplotlib->geomstats==2.2) (9.0.1)
    Requirement already satisfied: fonttools>=4.22.0 in ./neuroseed/lib/python3.8/site-packages (from matplotlib->geomstats==2.2) (4.30.0)
    Requirement already satisfied: pyparsing>=2.2.1 in ./neuroseed/lib/python3.8/site-packages (from matplotlib->geomstats==2.2) (3.0.7)
    Requirement already satisfied: packaging>=20.0 in ./neuroseed/lib/python3.8/site-packages (from matplotlib->geomstats==2.2) (21.3)
    Requirement already satisfied: bleach in ./neuroseed/lib/python3.8/site-packages (from nbconvert->jupyter->geomstats==2.2) (4.1.0)
    Requirement already satisfied: mistune<2,>=0.8.1 in ./neuroseed/lib/python3.8/site-packages (from nbconvert->jupyter->geomstats==2.2) (0.8.4)
    Requirement already satisfied: defusedxml in ./neuroseed/lib/python3.8/site-packages (from nbconvert->jupyter->geomstats==2.2) (0.7.1)
    Requirement already satisfied: nbclient<0.6.0,>=0.5.0 in ./neuroseed/lib/python3.8/site-packages (from nbconvert->jupyter->geomstats==2.2) (0.5.12)
    Requirement already satisfied: jupyterlab-pygments in ./neuroseed/lib/python3.8/site-packages (from nbconvert->jupyter->geomstats==2.2) (0.1.2)
    Requirement already satisfied: testpath in ./neuroseed/lib/python3.8/site-packages (from nbconvert->jupyter->geomstats==2.2) (0.6.0)
    Requirement already satisfied: pandocfilters>=1.4.1 in ./neuroseed/lib/python3.8/site-packages (from nbconvert->jupyter->geomstats==2.2) (1.5.0)
    Requirement already satisfied: webencodings in ./neuroseed/lib/python3.8/site-packages (from bleach->nbconvert->jupyter->geomstats==2.2) (0.5.1)
    Requirement already satisfied: qtpy in ./neuroseed/lib/python3.8/site-packages (from qtconsole->jupyter->geomstats==2.2) (2.0.1)
    Requirement already satisfied: executing in ./neuroseed/lib/python3.8/site-packages (from stack-data->ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (0.8.3)
    Requirement already satisfied: pure-eval in ./neuroseed/lib/python3.8/site-packages (from stack-data->ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (0.2.2)
    Requirement already satisfied: asttokens in ./neuroseed/lib/python3.8/site-packages (from stack-data->ipython>=7.23.1->ipykernel->jupyter->geomstats==2.2) (2.0.5)
    WARNING: You are using pip version 21.1.1; however, version 22.0.4 is available.
    You should consider upgrading via the '/home/noga/NeuroSEED/neuroseed/bin/python3 -m pip install --upgrade pip' command.

```python
import torch
import os
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
import time
from geomstats.geometry.poincare_ball import PoincareBall

from edit_distance.train import load_edit_distance_dataset
from util.data_handling.data_loader import get_dataloaders
from util.ml_and_math.loss_functions import AverageMeter

```

### Dataset description

As microbiome analysis is one of the most critical applications where the methods presented could be applied, we chose to use a dataset containing a portion of the 16S rRNA gene widely used in the biological literature to analyse microbiome diversity. Qiita (Clemente et al. 2015) contains more than 6M sequences of up to 152 bp that cover the V4 hyper-variable region collected from skin, saliva and faeces samples of uncontacted Amerindians. The full dataset can be found on the [European Nucleotide Archive](https://www.ebi.ac.uk/ena/browser/text-search?query=ERP008799), but, in this notebook, we will only use a subset of a few tens of thousands that have been preprocessed and labelled with pairwise distances. We also provide results on the RT988 dataset (Zheng et al. 2019), another dataset of 16S rRNA that contains slightly longer sequences (up to 465 bp).

```python
!gdown --id 1yZTOYrnYdW9qRrwHSO5eRc8rYIPEVtY2 # for edit distance approximation
#!gdown --id 1hQSHR-oeuS9bDVE6ABHS0SoI4xk3zPnB # for closest string retrieval
#!gdown --id 1ukvUI6gUTbcBZEzTVDpskrX8e6EHqVQg # for hierarchical clustering
```

    /home/noga/NeuroSEED/neuroseed/lib/python3.8/site-packages/gdown/cli.py:127: FutureWarning: Option `--id` was deprecated in version 4.3.1 and will be removed in 5.0. You don't need to pass it anymore to use a file ID.
      warnings.warn(
    Downloading...
    From: https://drive.google.com/uc?id=1yZTOYrnYdW9qRrwHSO5eRc8rYIPEVtY2
    To: /home/noga/NeuroSEED/edit_qiita_large.pkl
    100%|████████████████████████████████████████| 218M/218M [00:03<00:00, 62.5MB/s]

### Edit distance approximation

**Edit distance**  The task of finding the distance or similarity between two strings and the related task of global alignment lies at the foundation of bioinformatics. Due to the resemblance with the biological mutation process, the edit distance and its variants are typically used to measure similarity between sequences. Given two string $s_1$ and $s_2$, their edit distance $ED(s_1, s_2)$ is defined as the minimum number of insertions, deletions or substitutions needed to transform $s_1$ in $s_2$. We always deal with the classical edit distance where the same weight is given to every operation, however, all the approaches developed can be applied to any distance function of choice.

**Task and loss function** As represented in Figure 3, the task is to learn an encoding function $f$ such that given any pair of sequences from the domain of interest $s_1$ and $s_2$:
\begin{equation}ED(s_1, s_2) \approx n \; d(f(s_1), f(s_2)) \end{equation}

where $n$ is the maximum sequence length and $d$ is a distance function over the vector space. In practice this is enforced in the model by minimising the mean squared error between the actual and the predicted edit distance. To make the results more interpretable and comparable across different datasets, we report results using \% RMSE defined as:
\begin{equation}
\text{% RMSE}(f, S) = \frac{100}{n} \, \sqrt{L(f, S)} = \frac{100}{n} \, \sqrt{\sum_{s_1, s_2 \in S} (ED(s_1, s_2) - n \; d(f(s_1), f(s_2)))^2}
\end{equation}

which can be interpreted as an approximate average error in the distance prediction as a percentage of the size of the sequences.

In this notebook, we only show the code to run a simple linear layer on the sequence which, in the hyperbolic space, already gives particularly good results. Later we will also report results for more complex models whose implementation can be found in the [NeuroSEED repository](https://github.com/gcorso/NeuroSEED).

```python
class LinearEncoder(nn.Module):
    """  Linear model which simply flattens the sequence and applies a linear transformation. """

    def __init__(self, len_sequence, embedding_size, alphabet_size=4):
        super(LinearEncoder, self).__init__()
        self.encoder = nn.Linear(in_features=alphabet_size * len_sequence,
                                 out_features=embedding_size)

    def forward(self, sequence):
        # flatten sequence and apply layer
        B = sequence.shape[0]
        sequence = sequence.reshape(B, -1)
        emb = self.encoder(sequence)
        return emb

class PairEmbeddingDistance(nn.Module):
    """ Wrapper model for a general encoder, computes pairwise distances and applies projections """

    def __init__(self, embedding_model, embedding_size, scaling=False):
        super(PairEmbeddingDistance, self).__init__()
        self.hyperbolic_metric = PoincareBall(embedding_size).metric.dist
        self.embedding_model = embedding_model
        self.radius = nn.Parameter(torch.Tensor([1e-2]), requires_grad=True)
        self.scaling = nn.Parameter(torch.Tensor([1.]), requires_grad=True)

    def normalize_embeddings(self, embeddings):
        """ Project embeddings to an hypersphere of a certain radius """
        min_scale = 1e-7
        max_scale = 1 - 1e-3
        return F.normalize(embeddings, p=2, dim=1) * self.radius.clamp_min(min_scale).clamp_max(max_scale)

    def encode(self, sequence):
        """ Use embedding model and normalization to encode some sequences. """
        enc_sequence = self.embedding_model(sequence)
        enc_sequence = self.normalize_embeddings(enc_sequence)
        return enc_sequence

    def forward(self, sequence):
        # flatten couples
        (B, _, N, _) = sequence.shape
        sequence = sequence.reshape(2 * B, N, -1)

        # encode sequences
        enc_sequence = self.encode(sequence)

        # compute distances
        enc_sequence = enc_sequence.reshape(B, 2, -1)
        distance = self.hyperbolic_metric(enc_sequence[:, 0], enc_sequence[:, 1])
        distance = distance * self.scaling

        return distance
```

General training and evaluation routines used to train the models:

```python

def train(model, loader, optimizer, loss, device):
    avg_loss = AverageMeter()
    model.train()

    for sequences, labels in loader:
        # move examples to right device
        sequences, labels = sequences.to(device), labels.to(device)

        # forward propagation
        optimizer.zero_grad()
        output = model(sequences)

        # loss and backpropagation
        loss_train = loss(output, labels)
        loss_train.backward()
        optimizer.step()

        # keep track of average loss
        avg_loss.update(loss_train.data.item(), sequences.shape[0])

    return avg_loss.avg

def test(model, loader, loss, device):
    avg_loss = AverageMeter()
    model.eval()

    for sequences, labels in loader:
        # move examples to right device
        sequences, labels = sequences.to(device), labels.to(device)

        # forward propagation and loss computation
        output = model(sequences)
        loss_val = loss(output, labels).data.item()
        avg_loss.update(loss_val, sequences.shape[0])

    return avg_loss.avg
```

The linear model is trained on 7000 sequences (+700 of validation) and tested on 1500 different sequences:

```python
EMBEDDING_SIZE = 128

device = 'cuda' if torch.cuda.is_available() else 'cpu'
torch.manual_seed(2021)
if device == 'cuda':
    torch.cuda.manual_seed(2021)

# load data
datasets = load_edit_distance_dataset('./edit_qiita_large.pkl')
loaders = get_dataloaders(datasets, batch_size=128, workers=1)

# model, optimizer and loss
encoder = LinearEncoder(152, EMBEDDING_SIZE)
model = PairEmbeddingDistance(embedding_model=encoder, embedding_size=EMBEDDING_SIZE)
model.to(device)

optimizer = optim.Adam(model.parameters(), lr=0.001)
loss = nn.MSELoss()

# training
for epoch in range(0, 21):
    t = time.time()
    loss_train = train(model, loaders['train'], optimizer, loss, device)
    loss_val = test(model, loaders['val'], loss, device)

    # print progress
    if epoch % 5 == 0:
        print('Epoch: {:02d}'.format(epoch),
              'loss_train: {:.6f}'.format(loss_train),
              'loss_val: {:.6f}'.format(loss_val),
              'time: {:.4f}s'.format(time.time() - t))

# testing
for dset in loaders.keys():
    avg_loss = test(model, loaders[dset], loss, device)
    print('Final results {}: loss = {:.6f}'.format(dset, avg_loss))

```

Therefore, our linear model after only 50 epochs has a $\% RMSE \approx 2.6$ that, as we will see, is significantly better than any data-independent baseline.

### Closest string retrieval

This task consists of finding the sequence that is closest to a given query among a large number of reference sequences and is very commonly used to classify sequences. Given a set of reference strings $R$ and a set of queries $Q$, the task is to identify the string $r_q \in R$ that minimises $ED(r_q, q)$ for each $q \in Q$. This task is performed in an unsupervised setting using models trained for edit distance approximation. Therefore, given a pretrained encoder $f$, its prediction is taken to be the string $r_q \in R$ that minimises $d(f(r_q), f(q))$ for each $q \in Q$. This allows for sublinear retrieval (via locality-sensitive hashing or other data structures) which is critical in real-world applications where databases can have billions of reference sequences. As performance measures, we report the top-1, top-5 and top-10 scores, where top-$k$ indicates the percentage of times the model ranks the closest string within its top-$k$ predictions.

```python
from closest_string.test import closest_string_testing

closest_string_testing(encoder_model=model, data_path='./closest_qiita_large.pkl',
                       batch_size=128, device=device, distance='hyperbolic')
```

    Results: accuracy 0.441 0.524 0.573 0.605 0.633 0.650 0.674 0.686 0.698 0.707
    Top1: 0.441  Top5: 0.633  Top10: 0.707
    Total time elapsed: 108.2532s

Evaluated on a dataset composed of 1000 reference and 1000 query sequences (disjoint from the edit distance training set) the simple model we trained is capable of detecting the closest sequence correctly 44\% of the time and in approximately 3/4 of the cases it places the real closest sequence in its top-10 choices.

### Hierarchical clustering

Hierarchical clustering (HC) consists of constructing a hierarchy over clusters of data by defining a tree with internal points corresponding to clusters and leaves to datapoints. The goodness of the tree can be measured using Dasgupta's cost (Dasgupta 2016).

One simple approach to use NeuroSEED to speed up hierarchical clustering is similar to the one adopted in the previous section: estimate the pairwise distance matrix with a model pretrained for *edit distance approximation* and then use the matrix as the basis for classical agglomerative clustering algorithms (e.g. Single, Average and Complete Linkage). The computational cost to generate the matrix goes from $O(N^2M^2)$ to $O(N(M+N))$ and by using optimisations like locality-sensitive hashing the clustering itself can be accelerated.

The following code computes the pairwise distance matrix and then runs a series of agglomerative clustering heuristics (Single, Average, Complete and Ward Linkage) on it.

```python
from hierarchical_clustering.unsupervised.unsupervised import hierarchical_clustering_testing

hierarchical_clustering_testing(encoder_model=model, data_path='./hc_qiita_large_extr.pkl',
                                batch_size=128, device=device, distance='hyperbolic')
```

    Hierarchical torch.Size([10000, 152])
    {'single': {'DC': 334868934852.27704}, 'complete': {'DC': 334732900213.30273}, 'average': {'DC': 333694281572.618}, 'ward': {'DC': 334340801844.171}}

An alternative approach to performing hierarchical clustering we propose uses the continuous relaxation of Dasgupta's cost (Chami et al. 2020) to embed sequences in the hyperbolic space. In comparison to Chami et al. (2020), we show that it is possible to significantly decrease the number of pairwise distances required by directly mapping the sequences.
This allows to considerably speed up the construction especially when dealing with a large number of sequences without requiring any pretrained model. Figure 1 shows an example of this approach when applied to a small dataset of proteins and the code for it is in the NeuroSEED repository.

### Multiple Sequence Alignment

Multiple Sequence Alignment is another very common task in bioinformatics and there are several ways of using NeuroSEED to accelerate heuristics. The most commonly used programs such as the Clustal series and MUSCLE are based on a phylogenetic tree estimation phase from the pairwise distances which produces a guide tree, which is then used to guide a progressive alignment phase.

In Clustal algorithm for MSA on a subset of RT988 of 1200 sequences, the construction of the distance matrix and the tree takes 99\% of the total running time (the rest takes 24s out of 35 minutes). Therefore, one obvious improvement that NeuroSEED can bring is to speed up this phase using the hierarchical clustering techniques seen in the previous section.

The following code uses the model pretrained for edit distance to approximate the neighbour joining tree construction and the runs clustalw using that guide tree:

```python
from multiple_alignment.guide_tree.guide_tree import approximate_guide_trees

# performs neighbour joining algorithm on the estimate of the pairwise distance matrix
approximate_guide_trees(encoder_model=model, dataset=datasets['test'],
                        batch_size=128, device=device, distance='hyperbolic')

# Command line clustalw using the tree generated with the previous command.
# The substitution matrix and gap penalties are set to simulate the classical edit distance used to train the model
!clustalw -infile="sequences.fasta" -dnamatrix=multiple_alignment/guide_tree/matrix.txt -transweight=0 -type='DNA' -gapopen=1 -gapext=1 -gapdist=10000 -usetree='njtree.dnd'  | grep 'Alignment Score'
```

    Constructed approximate distance matrix in 151.0619s
    Constructed matrix
    Constructed nj tree in 3213.0527s
    /bin/bash: clustalw: command not found

An alternative method we propose for the MSA uses an autoencoder to convert the Steiner string approximation problem in a continuous optimisation task. More details on this in our paper and repository.

## 3. Benchmark

In this section, we compare the NeuroSEED approach to classical baseline alignment-free approaches such as k-mer and contrast the performance of neural models with different architectures and on different geometric spaces.

### Edit distance approximation

Figure 4: \% RMSE test set results on the Qiita and RT988 datasets. The first five models are the k-mer baselines and, in parentheses, we indicate the dimension of the embedding space. The remaining are encoder models trained with the NeuroSEED framework and they all have an embedding space dimension equal to 128. - indicates that the model did not converge.

Figure 4 highlights the advantage provided by data-dependent methods when compared to the data-independent baseline approaches. Moreover, the results show that it is critical for the geometry of the embedding space to reflect the structure of the low dimensional manifold on which the data lies. In these biological datasets, there is an implicit hierarchical structure given by the evolution process which is well reflected by the *hyperbolic* plane. Thanks to this close correspondence, even relatively simple models like the linear regression and MLP perform very well with this distance function.

Figure 5: \% RMSE on Qiita dataset for a Transformer with different distance functions.

The clear benefit of using the hyperbolic space is evident when analysing the dimension required for the embedding space (Figure 5). In these experiments, we run the Transformer model tuned on the Qiita dataset with an embedding size of 128 on a range of dimensions. The hyperbolic space provides significantly more efficient embeddings, with the model reaching the 'elbow' at dimension 32 and matching the performance of the other spaces with dimension 128 with only 4 to 16. Given that the space to store the embedding and the time to compute distances between them scale linearly with the size of the space, this would provide a significant improvement in downstream tasks over other NeuroSEED approaches.

**Running time** A critical step behind most of the algorithms analysed in the rest of the paper is the computation of the pairwise distance matrix of a set of sequences. Taking as an example the RT988 dataset (6700 sequences of length up to 465 bases), optimised C code computes on a CPU approximately 2700 pairwise distances per second and takes 2.5 hours for the whole matrix. In comparison, using a trained NeuroSEED model, the same matrix can be approximated in 0.3-3s (similar value for the k-mer baseline) on the same CPU. The computational complexity for $N$ sequences of length $M$ is reduced from $O(N^2\; M^2)$ to $O(N(M + N))$ (assuming the model is linear w.r.t. the length and constant embedding size). The training process takes typically 0.5-3 hours on a GPU. However, in applications such as microbiome analysis, biologists typically analyse data coming from the same distribution (e.g. the 16S rRNA gene) for multiple individuals, therefore the initial cost would be significantly amortised.

### Closest string retrieval

Figure 6 shows that also in this task the data-dependent models outperform the baselines even when these operate on larger spaces. In terms of distance function, the *cosine* distance achieves performances on par with the *hyperbolic*. This can be explained by the fact that for a set of points on the same hypersphere, the ones with the smallest *cosine* or *hyperbolic* distance are the same. So the *cosine* distance is capable of providing good orderings of sequence similarity but inferior approximations of their distance.

Figure 6: Accuracy of different models in the *closest string retrieval* task on the Qiita dataset.

### Hierarchical clustering

The results (Figure 7) show that the difference in performance between the most expressive models and the round truth distances is not statistically significant. The *hyperbolic* space achieves the best performance and, although the relative difference between the methods is not large in terms of percentage Dasgupta's cost (but still statistically significant), it results in a large performance gap when these trees are used for tasks such as MSA. The total CPU time taken to construct the tree is reduced from more than 30 minutes to less than one in this dataset and the difference is significantly larger when scaling to datasets of more and longer sequences.

Figure 7: Average Linkage \% increase in Dasgupta's cost of NeuroSEED models compared to the performance of clustering on the ground truth distances. Average Linkage was the best performing clustering heuristic across all models.

### Multiple Sequence Alignment

The results reported in Figure 8 show that the alignment scores obtained when using the NeuroSEED heuristics with models such as GAT are not statistically different from those obtained with the ground truth distances. Most of the models show a relatively large variance in performance across different runs. This has positive and negative consequences: the alignment obtained using a single run may not be very accurate, but, by training an ensemble of models and applying each of them, we are likely to obtain a significantly better alignment than the one from the ground truth matrix while still only taking a fraction of the time.

Figure 8: Percentage change in the alignment cost (- alignment score) returned by Clustal when using the heuristics to generate the tree as opposed to using NJ on real distances. The alignment was done on 1.2k unseen sequences from the RT988 dataset.

## 4. Limitations

As mentioned in the introduction, we believe that the NeuroSEED framework has the potential to be applied to numerous problems and, therefore, this project constitutes only an initial analysis of its geometrical properties and applications. Below we list some of the limitations of the current analysis and potential directions of research to cover them.

**Type of sequences** Both the datasets analysed consist of sequence reads of the same part of the genome. This is a very common set-up for sequence analysis (for example for microbiome analysis) and it is enabled by biotechnologies that can amplify and sequence certain parts of the genome selectively, but it is not ubiquitous. Shotgun metagenomics consists of sequencing random parts of the genome. This would, we believe, generate sequences lying on a low-dimensional manifold where the hierarchical relationship of evolution is combined with the relationship based on the specific position in the whole genome. Therefore, more complex geometries such as product spaces might be best suited.

**Type of labels** In this project, we work with edit distances between strings, these are very expensive when large scale analysis is required, but it is feasible to produce several thousand exact pairwise distance values from which the model can learn. For different definitions of distance, however, this might not be the case. If it is only feasible to determine which sequences are closest, then encoders can be trained using triplet loss and then most of the approaches presented would still apply. Future work could explore the robustness of this framework to inexact estimates of the distances as labels and whether NeuroSEED models, once trained, could provide more accurate predictions than its labels.

**Architectures** Throughout the project we used models that have been shown to work well for other types of sequences and tasks. However, the correct inductive biases that models should have to perform SEED are likely to be different to the ones used for other tasks and even dependent on the type of distance it tries to preserve. Moreover, the capacity of the hyperbolic space could be further exploited using models that directly operate in the hyperbolic space (Peng et al. 2021).

**Self-supervised embeddings** One potential application of NeuroSEED that was not explored in this project is the direct use of the embedding produced by NeuroSEED for downstream tasks. This would enable the use of a wide range of geometric data processing tools for the analysis of biological sequences.

## References

(Morgan & Huttenhower, 2012) Xochitl C Morgan and Curtis Huttenhower. Human microbiome analysis. PLoS Comput Biol, 2012.

(McDermott et al. 2021) Matthew McDermott, Brendan Yap, Peter Szolovits, and Marinka Zitnik. Rethinking relational encoding in language model: Pre-training for general sequences. arXiv preprint, 2021.

(Clemente et al. 2015) Jose C Clemente, Erica C Pehrsson, Martin J Blaser, Kuldip Sandhu, Zhan Gao, Bin Wang, Magda Magris, Glida Hidalgo, Monica Contreras, Oscar Noya-Alarcon, et al. ´The microbiome of uncontacted amerindians. Science advances, 2015.

(Zheng et al. 2019)Wei Zheng, Le Yang, Robert J Genco, Jean Wactawski-Wende, Michael Buck, and Yijun Sun. Sense: Siamese neural network for sequence embedding and alignment free comparison. Bioinformatics, 2019.

(Dasgupta 2016) Sanjoy Dasgupta. A cost function for similarity-based hierarchical clustering. In Proceedings of the forty-eighth annual ACM symposium on Theory of Computing, 2016.

(Chami et al. 2020) Ines Chami, Albert Gu, Vaggos Chatziafratis, and Christopher Re. From trees to continuous embeddings and back: Hyperbolic hierarchical clustering. Advances in Neural Information Processing Systems 33, 2020.

max_stars_repo_name: nongiga/NeuroSEED

max_issues_repo_name: nongiga/NeuroSEED

max_forks_repo_name: nongiga/NeuroSEED


text:
# Travelling Salesman Problem

#### Importación de libraría necesarias para la creación del problema

```python
#Import libraries
import pandas
import itertools
import numpy
import random
import networkx

#Import graphing libraries
%matplotlib inline

import matplotlib as mpl
import matplotlib.pyplot as plt
```

#### Creación del modelo con DOcplex
Parámetros
\begin{equation}
c_ij = \text{distance from i to j}
\end{equation}
Variables de decisión
\begin{equation}
x_{ij} = \begin{cases} 0, & \text{if the path from i to j is selected}  \\1, & \text{otherwise}\end{cases} \text{   }\forall i = 1,..,n \text{ , } j = 1,..,n\end{equation}

Función Objetivo
\begin{equation}
\text{Min} \sum\limits_{i=1}^n \sum\limits_{i\neq j,j=1}^n c_{ij}x_{ij}

\end{equation}

Restricciones
\begin{equation}
\sum\limits_{i=1,i\neq j}^n x_{ij} = 1,\forall j = 1..n
\end{equation}
\begin{equation}
\sum\limits_{j=1,j\neq i}^n x_{ij} = 1,\forall i = 1..n
\end{equation}
\begin{equation}
\sum\limits_{i\in S,j\neq i,j\in S} x_{ij} \leq \lvert S \rvert - 1,\text{   } \forall S\subset {1,..,n},2\leq \lvert S \rvert \leq n-2
\end{equation}

```python
from docplex.mp.model import Model

#create the base model
def tsp(data, max_length):
    m = Model(name='tsp',log_output=False)
    m.parameters.threads.set(1)

    all_locations = set(data.point_from.unique()).union(set(data.point_to.unique()))

    ####VARIABLES####
    data['travel_on'] = m.binary_var_list([(row.point_from, row.point_to) for row in data.itertuples()], name = 'travel_on')

    ####OBJECTIVE####
    distance_traveled = m.sum([row.travel_on * row.distance for row in data.itertuples()])
    m.add_kpi(distance_traveled, "distance")

    m.minimize(distance_traveled)
    ####CONSTRAINTS####
    #force each location to only have one outgoing edge
    for start_point in all_locations:
        m.add_constraint(m.sum([row.travel_on for row in data[data.point_from == start_point].itertuples()]) == 1, \
            'outgoing_edge_%s' % start_point)

    #force each location to only have one incoming edge
    for end_point in all_locations:
        m.add_constraint(m.sum([row.travel_on for row in data[data.point_to == end_point].itertuples()]) == 1, \
            'incoming_edge_%s' % end_point)

    return m

#update the parameters, solve the model, and return the results
def update_and_solve(m, data, mipgap = 0.0001):
    m.parameters.mip.tolerances.mipgap.set(mipgap)
    m.parameters.timelimit = 600
    m.solve()
    data['travel_on_val'] = [var.solution_value for var in data.travel_on.values]

def add_cycle_breaking_constraints(m,subtour_id, data):
    g = networkx.DiGraph()
    g.add_edges_from([(row.point_from,row.point_to) for row in data[data.travel_on_val > 0.5].itertuples()])
    num_cycles = 0
    longest_cycle = None
    subtours = {}
    for cycle in networkx.simple_cycles(g):
        cycleLength = len(cycle)
        num_cycles += 1
        subtours[subtour_id] = cycle
        subtour_id += 1
    if num_cycles == 1: #finished
        return m.kpis_as_dict()['distance'], cycle, subtour_id
    else:
        for new_subtour in subtours.keys():
            cycle = subtours[new_subtour]
            idxs = data[(data.point_from.isin(cycle)) & (data.point_to.isin(cycle)) & (data.travel_on_val > 0.5)].index.values
            m.add_constraint(m.sum(data.loc[idxs,'travel_on'].values) <= len(cycle) - 1,\
                'subtour_elimintaion_%d' % subtour_id)
        return m.kpis_as_dict()['distance'], [], subtour_id

def distance_tsp(data, max_length):
    m = tsp(data, max_length)
    mipgap = 0.1
    subtour_id = 0
    iteracion = 0
    max_iter = 500
    estado_iter = ''
    update_and_solve(m, data, mipgap = mipgap)
    distance, tour, subtour_id = add_cycle_breaking_constraints(m,subtour_id, data)
    if len(tour) == max_length:
        mipgap = mipgap - 0.01
    while (len(tour) < max_length or mipgap > 0.009) and iteracion < max_iter:
        update_and_solve(m, data, mipgap = mipgap)
        distance, tour, subtour_id = add_cycle_breaking_constraints(m,subtour_id, data)
        if len(tour) == max_length:
            mipgap = 0
        iteracion = iteracion + 1
    if iteracion == max_iter:
        estado_iter = 'SIN SOLUCION'
    else:
        estado_iter='OPTIMO'
    return distance, tour, estado_iter
```

#### Creación de coordenadas aleatorias
Se definen p puntos con coordenadas aleatorias. Estos puntos serán los lugares a visitar

```python
p = 20
df = pandas.DataFrame({'X':[random.randrange(100) for i in range(p)],'Y':[random.randrange(100) for i in range(p)]})
```

```python
df.plot(kind='scatter', x='X', y='Y', figsize=(6, 6), color='darkblue')
```

Se definen todas las posibles combinaciones entre dos puntos

```python
perm = [x for x in itertools.permutations(range(p), 2)]
```

Se definen las distancias para cada par de puntos

```python
data = pandas.DataFrame(columns=['point_from','point_to','distance'])

for i in range(len(perm)):
    p1 = perm[i][0]
    p2 = perm[i][1]
    data.loc[i] = p1 , p2 , \
    round(numpy.lib.scimath.sqrt((df.loc[p1,'X'] - df.loc[p2,'X'])**2+(df.loc[p1,'Y'] - df.loc[p2,'Y'])**2),2)

data.head(5)
```

      point_from
      point_to
      distance

      0
      0.0
      1.0
      49.20

      1
      0.0
      2.0
      14.42

      2
      0.0
      3.0
      35.34

      3
      0.0
      4.0
      26.68

      4
      0.0
      5.0
      3.16

#### Ejecución de modelo - DOcplex
Por medio del modelo definido en tsp_model, se obtiene la ruta óptima para los puntos definidos

```python
max_length = len(set(data.point_from.unique()).union(set(data.point_to.unique())))
distance, tour, estado = distance_tsp(data, max_length)
```

#### Solución de modelo - DOcplex
La solución del modelo nos retornará la distancia recorrida, el orden de recorrido de los puntos y el estado del modelo (óptimo o sin solución óptima)

```python
distance
```

    362.94999999999993

```python
tour
```

    [0.0,
     4.0,
     12.0,
     3.0,
     1.0,
     7.0,
     16.0,
     6.0,
     10.0,
     9.0,
     15.0,
     13.0,
     17.0,
     8.0,
     14.0,
     19.0,
     11.0,
     18.0,
     2.0,
     5.0]

```python
estado
```

    'OPTIMO'

#### Solución gráfica
Para evaluar la solución, graficamos las rutas definidas por la soluci´´on del modelo

```python
def plot_tsp(df,points, style='bo-'):
    "Plot lines to connect a series of points."
    plt.plot([df.loc[p,'X'] for p in points], [df.loc[p,'Y'] for p in points], style)
    plt.axis('scaled'); plt.axis('off')
```

```python
first_point= tour[0]
tour.append(first_point)
plot_tsp(df,tour)
```

ext: ipynb


text:
# Exercise 2
* Due: 20.1. at 14:00
* Max points: 8

## General rules
* **Submit your answers by using Moodle.** The link is https://moodle.jyu.fi/mod/assign/view.php?id=673612&forceview=1.
* Name your file as *lastname_ex2.xxx*. If you have multiple files, use *lastname_ex2_1.xxx, lastname_ex2_2.xxx* etc.
* **Remember to add comments to your answers and include also some testing part!**
* You will get feedback about your answers about a week from the due date
* The exercises will be given on the previous Wednesday's lecture

## Exercises

1. **(2 points)** Assume that you are buying and selling electricity. The more electricity you buy, the better price you get **per kwH** and the **buying prize** of electricity follows function $f(x) = 1-0.01x$, where $0\leq x\leq50$ is the amount of electricity you buy. On the other hand, the **price** that you get from selling electricity **follows** function $g(x) = 2-0.01x^2$ with $x$ again the amount of electricity that you sell. Formulate an optimization problem that maximizes the profit.
2. **(2 points)** Use bisection search to optimize problem
$$
\begin{align}
\min & (1-x)^2+x\\
\text{s.t. }&x\in [0,2].
\end{align}
$$
Code the search algorithm by yourself, not use the one shown in the lecture material.
3. **(2 points)** Use golden section search to optimize the problem from above exercise. Code the search algorithm by yourself, not use the one shown in the lecture material.
4. **(2 points)** Use differentiation/derivation to optimize the optimization problem of task 1 (above) and to verify the answer of optimization problem in task 2 (above).

max_stars_repo_path: Exercise 2.ipynb

max_issues_repo_path: Exercise 2.ipynb

max_forks_repo_path: Exercise 2.ipynb


text:
## Computer Vision Interpret

[`vision.interpret`](/vision.interpret.html#vision.interpret) is the module that implements custom [`Interpretation`](/train.html#Interpretation) classes for different vision tasks by inheriting from it.

```python
from fastai.gen_doc.nbdoc import *
from fastai.vision import *
from fastai.vision.interpret import *
```

```python
show_doc(SegmentationInterpretation)
```

class SegmentationInterpretation[source][test]

> SegmentationInterpretation(**`learn`**:[`Learner`](/basic_train.html#Learner), **`preds`**:`Tensor`, **`y_true`**:`Tensor`, **`losses`**:`Tensor`, **`ds_type`**:[`DatasetType`](/basic_data.html#DatasetType)=***`<DatasetType.Valid: 2>`***) :: [`Interpretation`](/train.html#Interpretation)

×
No tests found for SegmentationInterpretation. To contribute a test please refer to this guide and this discussion.

Interpretation methods for segmenatation models.

```python
show_doc(SegmentationInterpretation.top_losses)
```

top_losses[source][test]

> top_losses(**`sizes`**:`Tuple`, **`k`**:`int`=***`None`***, **`largest`**=***`True`***)

×
No tests found for top_losses. To contribute a test please refer to this guide and this discussion.

Reduce flatten loss to give a single loss value for each image

```python
show_doc(SegmentationInterpretation._interp_show)
```

_interp_show[source][test]

> _interp_show(**`ims`**:[`ImageSegment`](/vision.image.html#ImageSegment), **`classes`**:`Collection`\[`T_co`\]=***`None`***, **`sz`**:`int`=***`20`***, **`cmap`**=***`'tab20'`***, **`title_suffix`**:`str`=***`None`***)

×
No tests found for _interp_show. To contribute a test please refer to this guide and this discussion.

Show ImageSegment with color mapping labels

```python
show_doc(SegmentationInterpretation.show_xyz)
```

show_xyz[source][test]

> show_xyz(**`i`**, **`classes`**:`list`=***`None`***, **`sz`**=***`10`***)

×
No tests found for show_xyz. To contribute a test please refer to this guide and this discussion.

show (image, true and pred) from self.ds with color mappings, optionally only plot

```python
show_doc(SegmentationInterpretation._generate_confusion)
```

_generate_confusion[source][test]

> _generate_confusion()

×
No tests found for _generate_confusion. To contribute a test please refer to this guide and this discussion.

Average and Per Image Confusion: intersection of pixels given a true label, true label sums to 1

```python
show_doc(SegmentationInterpretation._plot_intersect_cm)
```

_plot_intersect_cm[source][test]

> _plot_intersect_cm(**`cm`**, **`title`**=***`'Intersection with Predict given True'`***)

×
No tests found for _plot_intersect_cm. To contribute a test please refer to this guide and this discussion.

Plot confusion matrices: self.mean_cm or self.single_img_cm generated by `_generate_confusion`

Let's show how [`SegmentationInterpretation`](/vision.interpret.html#SegmentationInterpretation) can be used once we train a segmentation model.

### train

```python
camvid = untar_data(URLs.CAMVID_TINY)
path_lbl = camvid/'labels'
path_img = camvid/'images'
```

```python
codes = np.loadtxt(camvid/'codes.txt', dtype=str)
get_y_fn = lambda x: path_lbl/f'{x.stem}_P{x.suffix}'
```

```python
data = (SegmentationItemList.from_folder(path_img)
        .split_by_rand_pct()
        .label_from_func(get_y_fn, classes=codes)
        .transform(get_transforms(), tfm_y=True, size=128)
        .databunch(bs=16, path=camvid)
        .normalize(imagenet_stats))
```

```python
data.show_batch(rows=2, figsize=(7,5))
```

```python
learn = unet_learner(data, models.resnet18)
learn.fit_one_cycle(3,1e-2)
learn.save('mini_train')
```

      epoch
      train_loss
      valid_loss
      time

      0
      10.024513
      3.442348
      00:15

      1
      6.325253
      2.343699
      00:03

      2
      4.759998
      2.108100
      00:02

```python
jekyll_warn("Following results will not make much sense with this underperforming model but functionality will be explained with ease")
```

 Warning: Following results will not make much sense with this underperforming model but functionality will be explained with ease

### interpret

```python
interp = SegmentationInterpretation.from_learner(learn)
```

Since `FlattenedLoss of CrossEntropyLoss()` is used we reshape and then take the mean of pixel losses per image. In order to do so we need to pass `sizes:tuple` to `top_losses()`

```python
top_losses, top_idxs = interp.top_losses(sizes=(128,128))
```

```python
top_losses, top_idxs
```

    (tensor([3.3195, 3.1692, 2.6574, 2.5976, 2.4910, 2.3759, 2.3710, 2.2064, 2.0871,
             2.0834, 2.0479, 1.8645, 1.8412, 1.7956, 1.7013, 1.6126, 1.6015, 1.5470,
             1.4495, 1.3423]),
     tensor([12,  4, 17, 13, 19, 18,  7,  8, 10,  1, 15,  0,  2,  9, 16, 11, 14,  5,
              6,  3]))

```python
plt.hist(to_np(top_losses), bins=20);plt.title("Loss Distribution");
```

Next, we can generate a confusion matrix similar to what we usually have for classification. Two confusion matrices are generated: `mean_cm` which represents the global label performance and `single_img_cm` which represents the same thing but for each individual image in dataset.

Values in the matrix are calculated as:

\begin{align}
\ CM_{ij} & = IOU(Predicted , True | True) \\
\end{align}

Or in plain english: ratio of pixels of predicted label given the true pixels

```python
learn.data.classes
```

    array(['Animal', 'Archway', 'Bicyclist', 'Bridge', 'Building', 'Car', 'CartLuggagePram', 'Child', 'Column_Pole',
           'Fence', 'LaneMkgsDriv', 'LaneMkgsNonDriv', 'Misc_Text', 'MotorcycleScooter', 'OtherMoving', 'ParkingBlock',
           'Pedestrian', 'Road', 'RoadShoulder', 'Sidewalk', 'SignSymbol', 'Sky', 'SUVPickupTruck', 'TrafficCone',
           'TrafficLight', 'Train', 'Tree', 'Truck_Bus', 'Tunnel', 'VegetationMisc', 'Void', 'Wall'], dtype='<U17')

```python
mean_cm, single_img_cm = interp._generate_confusion()
```

```python
mean_cm.shape, single_img_cm.shape
```

    ((32, 32), (20, 32, 32))

`_plot_intersect_cm` first displays a dataframe showing per class score using the IOU definition we made earlier. These are the diagonal values from the confusion matrix which is displayed after.

`NaN` indicate that these labels were not present in our dataset, in this case validation set. As you can imagine it also helps you to maybe construct a better representing validation set.

```python
df = interp._plot_intersect_cm(mean_cm, "Mean of Ratio of Intersection given True Label")
```

Next let's look at the single worst prediction in our dataset. It looks like this dummy model just predicts everything as `Road` :)

```python
i = top_idxs[0]
df = interp._plot_intersect_cm(single_img_cm[i], f"Ratio of Intersection given True Label, Image:{i}")
```

Finally we will visually inspect this single prediction

```python
interp.show_xyz(i, sz=15)
```

```python
jekyll_warn("""With matplotlib colormaps the max number of unique qualitative colors is 20.
So if len(classes) > 20 then close class indexes may be plotted with the same color.
Let's fix this together :)""")
```

 Warning: With matplotlib colormaps the max number of unique qualitative colors is 20.
So if len(classes) > 20 then close class indexes may be plotted with the same color.
Let's fix this together :)

```python
interp.c2i
```

    {'Animal': 0,
     'Archway': 1,
     'Bicyclist': 2,
     'Bridge': 3,
     'Building': 4,
     'Car': 5,
     'CartLuggagePram': 6,
     'Child': 7,
     'Column_Pole': 8,
     'Fence': 9,
     'LaneMkgsDriv': 10,
     'LaneMkgsNonDriv': 11,
     'Misc_Text': 12,
     'MotorcycleScooter': 13,
     'OtherMoving': 14,
     'ParkingBlock': 15,
     'Pedestrian': 16,
     'Road': 17,
     'RoadShoulder': 18,
     'Sidewalk': 19,
     'SignSymbol': 20,
     'Sky': 21,
     'SUVPickupTruck': 22,
     'TrafficCone': 23,
     'TrafficLight': 24,
     'Train': 25,
     'Tree': 26,
     'Truck_Bus': 27,
     'Tunnel': 28,
     'VegetationMisc': 29,
     'Void': 30,
     'Wall': 31}

```python
show_doc(ObjectDetectionInterpretation)
```

class ObjectDetectionInterpretation[source][test]

> ObjectDetectionInterpretation(**`learn`**:[`Learner`](/basic_train.html#Learner), **`preds`**:`Tensor`, **`y_true`**:`Tensor`, **`losses`**:`Tensor`, **`ds_type`**:[`DatasetType`](/basic_data.html#DatasetType)=***`<DatasetType.Valid: 2>`***) :: [`Interpretation`](/train.html#Interpretation)

×
No tests found for ObjectDetectionInterpretation. To contribute a test please refer to this guide and this discussion.

Interpretation methods for classification models.

```python
jekyll_warn("ObjectDetectionInterpretation is not implemented yet. Feel free to implement it :)")
```

 Warning: ObjectDetectionInterpretation is not implemented yet. Feel free to implement it :)

## Undocumented Methods - Methods moved below this line will intentionally be hidden

## New Methods - Please document or move to the undocumented section


text:
```python
1+1
```

    2

```python
import re
from functools import partial
```

## Preparing the List of Prime Gaps

The file `raw_prime_gaps_list.txt` contains the list of known lead prime gaps for all gaps up to $9998$. The data is uniformly stored. A typical row looks like the following, and the meaning of each element of the row is described below.

    6* CFC Glaisher 1877   1.91     2  23
    || |   |        |      |        |  |
    || |   |        |      |        |  ---> first prime in the prime gap
    || |   |        |      |        ------> number of digits of the first prime
    || |   |        |      ---------------> "merit", which is bigger if the gap is bigger wrt prime size
    || |   |        ----------------------> year of discovery
    || |   -------------------------------> discoverer
    || -----------------------------------> C if "conventional gap" (all here are)
    ||                                      F if known first occurrence of gap size, N if not, ? if unknown
    ||                                      C if certified primes, P if probabilistic primes
    |-------------------------------------> * if a "maximal gap", larger than all previous gaps
    --------------------------------------> size of gap

Sometimes the first prime is stored not as a number, but as an expression in the form $a \cdot b\# / c - d$, where $b\#$ denotes a *primorial*, or as an expression in the form $2^a + b$. We first read in the list. We extract those rows which are maximal (i.e. has * following the gapsize). Then we keep the gap size and the first prime.

```python
with open("raw_prime_gaps_list.txt", "r") as f:
    lines = f.readlines()
print(lines[3])
```

         6* CFC Glaisher 1877   1.91     2  23

```python
def is_maximal_gap(splitlist):
    if "*" in splitlist[0]:
        return True
    return False

def extract_gap_and_prime(line):
    splitlist = line.split(None, 6) # Split on whitespace (None), and make at most 6 splits
                                    # Splitting at most 6 times keeps the expression for the first prime together
    gap_str = splitlist[0]
    first_prime_str = splitlist[6]
    first_prime_str = first_prime_str.rstrip()

    #assert gap_str[-1] == "*"
    if gap_str[-1] == "*":
        gap = gap_str[:-1]
    else:
        gap = gap_str
    gap = int(gap)

    if "*" in first_prime_str:
        first_prime = extract_primorial(first_prime_str)
    elif "^" in first_prime_str:
        first_prime = extract_two_power_prime(first_prime_str)
    else:
        first_prime = int(first_prime_str)

    return gap, first_prime

def memoize(f):
    memo = {}
    def helper(x):
        if x not in memo:
            memo[x] = f(x)
        return memo[x]
    return helper

primes_to_541 = primes_first_n(100)
@memoize
def primorial(n):
    return prod([p for p in primes_to_541 if p<= n])

def extract_primorial(expression):
    pattern = re.compile(r"(\w+)\W*\*\W*(\w+)\W*\#\W*/\W*(\w+)\W*-\W*(\w+)")
    m = pattern.match(expression)
    a,b,c,d = map(int, m.groups())
    return a*primorial(b)/c - d

def extract_two_power_prime(expression):
    pattern = re.compile(r"2\^(\w+)\W*\+\W*(\w+)")
    m = pattern.match(expression)
    exp, addend = map(int, m.groups())
    return 2**exp + addend
```

```python
# First we test the method on just the first 10 prime gaps
gaps = []
gaps_dict = dict()
for line in lines[:10]:
    gap, first_prime = extract_gap_and_prime(line)
    gaps.append(gap)
    gaps_dict[gap] = [first_prime, first_prime + gap]
```

```python
# These are the gaps
gaps
```

    [1, 2, 4, 6, 8, 10, 12, 14, 16, 18]

```python
# Clear prime pair storage
gaps_dict
```

    {1: [2, 3],
     2: [3, 5],
     4: [7, 11],
     6: [23, 29],
     8: [89, 97],
     10: [139, 149],
     12: [199, 211],
     14: [113, 127],
     16: [1831, 1847],
     18: [523, 541]}

```python
# First we test the method on just the first 10 prime gaps
gaps = []
gaps_dict = dict()
for line in lines:
    gap, first_prime = extract_gap_and_prime(line)
    gaps.append(gap)
    gaps_dict[gap] = [first_prime, first_prime + gap]
```

```python
# Choose random representative to check
gaps_dict[116]
```

    [5845193, 5845309]

## Setting up the model functions

We will model the distribution of primes in a few different ways. The nature of the distribution of primes is uncertain. By the prime number theorem, the average prime gap looks like
$$ X, X + \log X.$$
We expect instead that there are many gaps that are much larger. Cramer thought that prime gaps were never larger than
$$ X, X + \log^2 X.$$
It turns out that he wasn't correct, but it may be true that there is a constant $B$ such that gaps are never larger than
$$ X, X + B \log^2 X.$$
It is known that this constant, if it exists, is at least $2e^{-\gamma} \approx 1.1$.
It is also conjectured that no such constant exists, but that for any $A > 2$ there exists a corresponding constant (which I also denote by $B$) such that gaps are never larger than
$$ X, X + B \log^A X.$$

Recall in the background that we conjecture that for each $\alpha$ there exists a constant $C(\alpha)$ such that
$$ p_{n+1}^\alpha - p_n^\alpha \leq C(\alpha).$$

We will model how gaps behave according to these models and compare them to the actual behavior.
To do this, we define the function
$$ f(A, B, \alpha, x) := (X + B \log^A X)^\alpha - X^\alpha.$$
For each fixed $A, B$, and $\alpha$, maximizing $f_{A, B, \alpha}(x)$ gives a largest-possible estimate of the value of $C(\alpha)$. We say largest-possible because this presupposes that there are primes at exactly the maximizing value of $x$, which corresponds to $x$ and $x + B \log^A x$ being primes themselves.

Note that if Cramer's Conjecture were true (for instance), then it may be that the $C(\alpha)$ predictions we get from the Cramer model might be slightly larger than reality, and in general each model errs on the side of overestimating rather than underestimating.

We will set up model functions with values of $(A, B)$ being $(1,1)$ (PNT), $(2,1)$ (Naive Cramer), $(2,2e^{-\gamma})$ (Current Cramer), $(2, 10)$ (Cramer), and $(2.5, 1)$ (Adleman-McCurley).

### Maximizing the model functions

Maximizing $f_{A,B,\alpha}(x)$ is done naively. Compute the derivative and find when it's zero. There is a unique zero greater than $1$, which is our maximum point. A first-order estimation shows that $f_{A,B,\alpha}$ is maximized near $m(A,B,\alpha) = e^{A/(1-\alpha)}$. Note that $B$ does not affect this first-order estimation, so I do not expect the three variants of Cramer to behave significantly differently.

We have that
$$ \begin{align}
f(A, B, \alpha, x) &= (x + B \log^A x)^\alpha - x^\alpha, \\
f'(A, B, \alpha, x) &= \alpha(x + B \log^A x)^{\alpha - 1} \left( 1 + AB \frac{\log^{A-1} x}{x} \right) - \alpha x^{\alpha - 1}.
\end{align}$$

Note that $f' = 0$ (for $x > 1$) is the same as $f'/(\alpha x^{\alpha - 1}) = 0$, which can be rearranged into
$$ g_{A, B, \alpha}(x) := \left( 1 + AB \frac{\log^{A-1} x}{x} \right) - \left( 1 + B \frac{\log^A x}{x} \right)^{1 - \alpha} = 0.$$

This is a pretty well-behaved function, and numerical root solvers should be able to handle it. This is what we implement now.

```python
def F(A, B, alpha, x):
    return (x + B*log(x)**A)**alpha - x**alpha
```

```python
def g(A, B, alpha, x):
    return (1 + A * B * log(x)**(A-1)/x) - (1 + B*log(x)**A / x)**(1-alpha)
```

```python
find_root(partial(g,2,1,.931), 2, 10**17)
```

    3875382310049.5977

Actually, this caps out at 17 digit precision and won't suffice for $\alpha > 0.93$. But for really large $x$, the binomial expansion should become more and more accurate, leading to $e^{A/(1-\alpha)}$ being a better and better estimate of the maximizing $x$. For $\alpha > 0.93$, I need to use a different technique. This runs into 64-bit precision arithmetic problems. To solve this, for $\alpha > 0.9$, I will use a homegrown bisection rootfinder in a high precision arithmetic sage environment.

```python
RF = RealField(1000)
def RF_F(A, B, alpha, x):
    return (RF(x) + B*RF(log(x))**A)**alpha - RF(x)**alpha
```

```python
start = exp(2/.005)
```

```python
RF_F(2, 1, .995, start)
```

    21545.3770912687021748108986403427084075846916931792212525384307640395688240237337048227691696804791326950035399102608956974569724349727991073964722581240384122762528927904117476919639364452675336235413841573393011675588938510358312010132614449997916430297316353906271329030935986082731144694523430359

```python
def RF_g(A, B, alpha, x):
    return (RF(1) + A * B * RF(log(x))**(A-1)/RF(x)) - (RF(1) + B*RF(log(x))**A / RF(x))**(1-alpha)
```

```python
def RF_gsim(x):
    return RF_g(2,1,0.8,x)
```

```python
start = exp(2/.2)
```

```python
RF_gsim(start)
```

    1.64444546432196045785420540155893726063200274572267336535350345517820986134145980935118162918424805451405780491393962735635632068351784816308968439342270391699018760425327496611550914706761553155070977740717893810316453496034348416368377098158251454957552350910050505616367388235938291760692342296640e-6

```python
find_root(RF_gsim, 2, 10**17)
```

    22423.973931885273

```python
# Note that this requires that func(2*start) < 0 to work properly
# In our case, this means that the `start` guess can't be too bad
def bisection_rootfinding(func, start, iterations=1, goal=10**(-7)):
    if start*2**(-iterations) < goal:
        return start
    if func(start) > 0:
        return bisection_rootfinding(func, start + start * 2**(-iterations), iterations=iterations+1, goal=goal)
    return bisection_rootfinding(func, start - start * 2**(-iterations), iterations=iterations+1, goal=goal)
```

```python
bisection_rootfinding(RF_gsim, start)
```

    22423.9739319633

Notice in particular that this is essentially the same as scipy's find_root result, so my method seems to work pretty well.

Now let's try it on something larger.

```python
def RF_gsim2(x):
    return RF_g(2,1,0.95,x)
```

```python
start = RF(exp(2/.05))
```

```python
find_root(RF_gsim2, 2, 10**16)
```

    9999999710581220.0

```python
bisection_rootfinding(RF_gsim2, start)
```

    2.35385266837042022843327073581445209550189717806195566505913760542398087081961498860906821406821821384500066808118082097402013219336066587283757576744855862678903728952660301215698355426390611932122470916392998635276106218830903304888702478640727370013260828376486514621410432957225533987162136976562e17

```python
fr = find_root(RF_gsim2, 2, 10**16)
br = bisection_rootfinding(RF_gsim2, start)
```

```python
br/fr
```

    23.538527364953389

```python
RF_F(2,1,0.95,br)
```

    205.709630519650881691924521792765193253506704662129143200009919811578164796694847391279127980866813238127336707532494092877547455355110417039290138055484169795377093169784536561532277570158710114624712809121313730522096004967342309066003863302396648990201591070015928728914879882408923972678419871216

```python
RF_F(2,1,0.95,fr)
```

    204.359554823890368403363980346311818555686040610577329335810851723660578683578255065736205720488148424074134614592141667900060298707050938806333994911863866233402912888249520944989763989818618341662786563588774804722056915716107185579765431262102230328894524070859044620196971152617523013227880281789

This is interesting --- there is a region in which find_root fails due to precision reasons, but not in an obviously detectable way. But as this clearly shows, my bisection method suffices.

## Data Creation

The plan is to loop through $\alpha$ from $0.005$ to $0.995$, compute everything, and compare everything. But for many initial $\alpha$, everything is boring. Let's first find exactly when $\alpha$ is boring.

```python
def estimate_c_alpha(alpha):
    MAX = 0
    MAXIMIZING_PRIME = 0
    MAXIMIZING_GAP = 0
    for gap in gaps:  # Note: It is possible to optimize this selection based on vals for other alpha
        a,b = gaps_dict[gap]
        diff = RF(b)**alpha - RF(a)**alpha
        if diff > MAX:
            MAX = diff
            MAXIMIZING_PRIME = a
            MAXIMIZING_GAP = gap
    return MAX, MAXIMIZING_GAP, MAXIMIZING_PRIME
```

```python
for a in range(20):
    alpha = 0.05 + 0.04 * a
    print alpha, estimate_c_alpha(alpha)
```

    0.0500000000000000 (0.0273310781848302800658590880398522815347590177971289917771832998742533538405129347101901809408899390154774514006487716615848197146224848388750750143388527967917247425289525849885984701315952214994461518869467334113477417254612430110708981668255496105519249618121865605658797754277991125986546424926508, 2, 3)
    0.0900000000000000 (0.0519370809347252520685872649136154113889575896783660909311486774384967763294410100201853785718930335741505701554636256358169356949220163741754562328432777159041140235202911046252875832381116056987398716939257761619006205913299310880723598291675100978615964192750372577866405938486631637439508364320858, 2, 3)
    0.130000000000000 (0.0792030225749460145408332846210646211349709553727122695589775991379642815865105358693876664136993146934495792581811221651860043234387972830839584009007088431244195480057541946715422880580136690530288636745487275691029720040890241179152279971140525502557180885918952288021915731111539843731655078009650, 2, 3)
    0.170000000000000 (0.111181263499968631320136517388785040808714392234047453641809201804708639290279254938652323470002276346079710485543325649046090275258826662034012458243570622205524629138788775066951458339047395559511537749810240846304395803420508311453056782745024686529577255723133918595809482942348237543967846160003, 4, 7)
    0.210000000000000 (0.149826273160475994962243083701533895900842934406639563544922783813024602459916508959343188251692967568113386862740598596871970918514700800403775752958239520903512136772589541047142530312391410824945291551508017637758651939264125623653517866914330922645399426878943623423072451500026730493789321506636, 4, 7)
    0.250000000000000 (0.194583725140086155235410867344262555198847337413419396460571687675616808191366806740479594064082037831924757984080937834304785500452520906853755534860850295188536908645060113128466605462778750838482301927791913256899345519152556683960954609815854565747156297795965378511888898199097339590595739963184, 4, 7)
    0.290000000000000 (0.246248813819103170253557250449355251524182358470425963859211259303798615289095890248478723442933301746121823090879497185382980418432530720592316569672647132718753876297567208875018396019594127419469968573820605380516809172280817095923124075568117875266835223189694717618200541375491919546679627098467, 4, 7)
    0.330000000000000 (0.305711343206456535065347320087438603464924475308724063528038459208099421322271293811067461182623734565140057698384504321768168205679408867081213453787162538160335263784242371458734357856213201079531522786702433012164371829170292653473361351573652392532724606923212418835000408796683348784457264111052, 4, 7)
    0.370000000000000 (0.373966429543365074817695190057943133748849296316293016605715723578219057277796879857315638543939944511494607661966046966136173460500005119867241883996572229561550012814088131390300062064484716155148271213476199228912371725988542658561793549101843469396474674653770477560665905219304663796828624873265, 4, 7)
    0.410000000000000 (0.452126376701729070189230306846564166572718811750208095294640063601049113790103772310658620779822483542466716249347768517447112764806799712472932886358915375650933223964845094374376513000612507694391310907365957562029322741530071576739195666770017192298200796425941866160240244584347982547800294338410, 4, 7)
    0.450000000000000 (0.541433850160674856517448019977648670845387468609120824704039687062082365058862967523552961325586337836348173605227514038998558396631077325859615480296792512772133696891677065764305911703916867104804927875994085360984285663433909652187556181498138595496488010425855281207978804879931061852800799756917, 4, 7)
    0.490000000000000 (0.643276488631938665905352690029550515258841518006625633245321825184230438817001007398166180674071026761014493558595088627772662382841185013505094743174300363527368853334206213217354029111657745924317354222327336246191131256221501182273919522940639882386609873185598418596322399256320888450121570394212, 4, 7)
    0.530000000000000 (0.782274375222570659961116752409324914028545518530568244320104431746104885735461769616383026384078842399462241244581615655826225425861047007355692369799948454253659291398201594313709988207531393383294538319815566762354200807118161016244669292398154414432780957805186418730929458670499355583639366807496, 14, 113)
    0.570000000000000 (1.01884142574031551509021338668081500059651084095766586110657227080467560613268941623508491511136916158752831640137800592349753586674611323376787825530400982051690086491024022896463721852472889392683539307104501345623645811741468826667843990073735076008410738714094795624533929203128446699606735615606, 14, 113)
    0.610000000000000 (1.32041627446091249684010282433791853388208830739656554504138875808096002722343829956112453368302449982309552879586646808802152034816717620334622823023810628771492706007214793761849531553602264771273615622992631932925771633549693623574091063218430747286446606981048690064876572566265298842449766093359, 14, 113)
    0.650000000000000 (1.77606150694455767148184879493035887223003417395280279720283805453747271273934846538180951981092679591381731044243677198786743664042472880933235210451979387660865985682252547374928040543948333428791944282726895450499622924451542431454744301895332923950015913854758146588145372911146101685185764793123, 34, 1327)
    0.690000000000000 (2.51494758180812497524004824914538864215457135214415764156856021955949292849607775801439209422484553476893581354288954398402859547490062743970891503842286272825447055396110162612007004241446273913524824849127987779554312507146909153549001997795834074088281582457832963082371713882142574398014399748687, 34, 1327)
    0.730000000000000 (3.54926107221611583152257524364752677913475564845984515579046097207485619605424351240800845901580940996317144955978790853632298886397526616816851581340281228865107002973080947704010636374478688987465249159849949548386146640449395044703904284644139242341066041035459498716689367467428185462776529563374, 34, 1327)
    0.770000000000000 (5.12182230348182730576342401351876551208525168860138306791801121718534370070661507654496625833636423338279349496907637710728402502004783185640546392968476514453708497086637272112232436492441948446521882798746778382139248643241899662976156695772601498121414397776266212838948354134427604190418086056576, 72, 31397)
    0.810000000000000 (8.15293966895682191988905194405319109615825254285457443436102887983737009629372030499598020896742552803941325831592633996341564615822777249900088377876688336896985842163801399806797855932027735326506916167614518727288169134555219040681812889496872163795527757013603506051256509476820385589025707472966, 72, 31397)

It appears that starting at $\alpha = 0.005$ will be fine, though not a lot happens until $\alpha \approx 0.5$. The semiimportant data there is really just the value of $C(\alpha)$.

Let us now set up the data structures to store the data. These will be a set of dictionaries, contained in a "model" dictionary.

```python
experimental_model = dict()
PNT_model = dict()
naive_cramer_model = dict()
current_cramer_model = dict()
weak_cramer_model = dict()
AM_model = dict()

models = dict()
models["experimental"] = experimental_model
models["PNT"] = PNT_model
models["naive"] = naive_cramer_model
models["current"] = current_cramer_model
models["weak"] = weak_cramer_model
models["AM"] = AM_model
```

In each model dictionary, the keys will be values of $\alpha$ from $0.005$ to $0.995$, and the values will be `(C(alpha), gap, lower_prime)` tuples. That's that.

I'll store the data as pickled dictionaries, and as json dumps (for external portability).

```python
RF = RealField(1000)
def RF_F(A, B, alpha, x):
    return (RF(x) + B*RF(log(x))**A)**alpha - RF(x)**alpha

def RF_g(A, B, alpha, x):
    return (RF(1) + A * B * RF(log(x))**(A-1)/RF(x)) - (RF(1) + B*RF(log(x))**A / RF(x))**(1-alpha)

def maximize_model(A, B, alpha):
    fun = lambda x: RF_F(A, B, alpha, x)
    deriv = lambda x: RF_g(A, B, alpha, x)
    estimate_lower_prime = 2*RF(exp(A/(1 - alpha)))
    lower_prime = bisection_rootfinding(deriv, start)
    Calpha = fun(lower_prime)
    gap = B*log(lower_prime)**A
    return (Calpha, gap, lower_prime)
```

```python
def make_experimental_data(alpha):
    Calpha, gap, lower_prime = estimate_c_alpha(alpha)
    return (Calpha, gap, lower_prime)

def make_model_data(A, B, alpha):
    return maximize_model(A,B,alpha)

def make_PNT_model_data(alpha):
    return make_model_data(1,1,alpha)

def make_naive_model_data(alpha):
    return make_model_data(2,1,alpha)

def make_current_model_data(alpha):
    return make_model_data(2,1.1229189671337703396, alpha)

def make_weak_model_data(alpha):
    return make_model_data(2, 10, alpha)

def make_AM_model_data(alpha):
    return make_model_data(2.1, 1, alpha)
```

```python
for a in range(1000):
    alpha = 0.000 + a/1000
    models["experimental"][alpha] = make_experimental_data(alpha)
    models["PNT"][alpha] = make_PNT_model_data(alpha)
    models["naive"][alpha] = make_naive_model_data(alpha)
    models["current"][alpha] = make_current_model_data(alpha)
    models["weak"][alpha] = make_weak_model_data(alpha)
    models["AM"][alpha] = make_AM_model_data(alpha)
    if a%50 == 0:
        print a
```

    0
    50
    100
    150
    200
    250
    300
    350
    400
    450
    500
    550
    600
    650
    700
    750
    800
    850
    900
    950

```python
#t = models.keys()[0]
#for k in models[t].keys():
#    print k
#    for key in models.keys():
#        print key, "{}".format(float(models[key][k][0]))
#    print
```

```python
alist = models['experimental'].keys()
```

```python
alist = sorted(alist)
alist[:10]
```

    [0.000000000000000,
     0.00100000000000000,
     0.00200000000000000,
     0.00300000000000000,
     0.00400000000000000,
     0.00500000000000000,
     0.00600000000000000,
     0.00700000000000000,
     0.00800000000000000,
     0.00900000000000000]

```python
models['experimental'][alist[-50]]
```

    (186.266539364655933626490915855618094964571041849497438999997826389079814853291919149871340027107619054174485385512758967203039410422848728780559014112645466156493093420011269068962647210601068265065442552646499480648773915312419366944230198770879264447877675770996717230251871662312084943372625852009,
     1132,
     1693182318746371)

```python
exp_plotys = [models['experimental'][a][0] for a in alist]
cur_plotys = [models['current'][a][0] for a in alist]
weak_plotys = [models['weak'][a][0] for a in alist]
naive_plotys = [models['naive'][a][0] for a in alist]
pnt_plotys = [models['PNT'][a][0] for a in alist]
am_plotys = [models['AM'][a][0] for a in alist]
```

```python
exp_plotys[:10]
```

    [0,
     0.000511517768525288490431299564019530566521268528293360224090648225671358406177645208427960618154537414232201941218446425306818783533501147565303119802041928147556730502807612639983069318331337966366989636147947059894660635945989823381565534871988347392491415512712789530278446363938930989574701960971573,
     0.00102442172449282884332508157534822491563469058472756719379918434497150349164072269796738895556658217479601021664620351202402928247792955884906721383928526283258088661094103784976857592343906755172954042356731975418539998880830592392230883140406784021001485475079790186985430084990582641020616554607448,
     0.00153871471873605261609594567254478753444623540399582289121122633476589317565823494066369723602547767528506869164326737509390842483499917786372877323718544377009408577779985644905224201012604546334079646010266208402550801858445699818597717680152286230410860527428679105037060870377855630342605362083545,
     0.00205439960735970041704645998947656177812210389930506653845904193767730286099795243878248319312669702013552227025385042847607239098634986006808941199232325363397106263902487262040597462454010379778333332955140235088061999175180377117872404786751318950384094489121442812986482906623907565586558415282809,
     0.00257147925174905936109663014134097336834442393609311620079119771972767194500766471612512274531782092337622866174440543746300567206334808347961436348689967356555975464533599100596026029910740394817019762901927392208024221858943093784678802392647510480682890612867455763120157088180379124648103287551767,
     0.00308995651857921622546796559200596882746424951676754671589696002353232888475187180207023465138402128824739465016627718117073708463068539941891599637570574693047123871730515453760687505233732008693440140525073857761664831298381288803400485841363604104199122966657955365153196391365677979590315607030346,
     0.00360983427982432633151290676273828177065418823632371050144958102026278975990287331278230039244296151938098161366750880745556329436894367278846820410099015999510215259570924904252430315259755884537232301410584943365318499431178635391607582265002963417075775796112480314734163131768771943838626823279178,
     0.00413111541276689817892355442974284642043215842675225341315076577505076312445243624385504765217391148512516335857815162425101960944675790086942557310550369600595084936465873338520976161756438902253630154479581694290320640417906763915305024844757572628498939719652771807118350659671983418068743491981217,
     0.00465380280000709340462715100328122372167500192805071772269677766578089916306974130382072879733511866037963297784029049269126452696333254096221820907436597556468574696079353308055141116896900016154860397150953509265025577577698312210423534852671473702537036156003845806675868441170706672772677146830196]

```python
P = scatter_plot(zip(alist, exp_plotys), xmin=0, xmax=0.995, ymin=0, ymax=1000, markersize=2, facecolor='black', edgecolor='')
P += scatter_plot(zip(alist, cur_plotys), markersize=2, facecolor='red', edgecolor='')
P += scatter_plot(zip(alist, weak_plotys), markersize=2, facecolor='green', edgecolor='')
P += scatter_plot(zip(alist, pnt_plotys), markersize=2, facecolor='blue', edgecolor='')
P += scatter_plot(zip(alist, am_plotys), markersize=2, facecolor='orange', edgecolor='')
P += scatter_plot(zip(alist, naive_plotys), markersize=2, facecolor='purple', edgecolor='')
P
```

```python
alpha_in_8_to_10 = []
for a in range(200):
    alpha = 0.8 + 0.001 * a
    print "\r" + str(alpha),
    alpha_in_8_to_10.append((alpha, estimate_c_alpha(alpha)[0]))
```

    0.999000000000000

```python
alpha_in_8_to_10[:10]
```

    [(0.800000000000000,
      7.26014575499710583007752218476427663828926179611157953282726566695049107265377270133560692997036040530500477768867195065477012290917719308708782018245669078582044688607873551418902354119005826714341041787179613542691696812731368111577302778025140783938715803646904732042037141424376103751080011998710),
     (0.801000000000000,
      7.34488929791419664639472074261502192281957152607915337406435419916315200585249998380204802464158172422854707558003676849509272095842610708260861187394075385302864085424482425636215460071693157771352179388596036196527196233733066745050675507537225353771965938807629237803255305635045814031962805504380),
     (0.802000000000000,
      7.43061042252953118360499889585345215018025636058791613024779180241959349552833645769717749958099242955401006554977231966370929595157074587453403662539834509389164431402822851597631998341664334362076793184715469597355362017466344732729542069893132901288343947966821105638284381796933618371826846838619),
     (0.803000000000000,
      7.51732029848634200607433174510508482594122453398141176966268473529751351184646813784800405777531754550192381100544366450785598014002895673006395936770365346627588379366606543940945237679157583436394422090354325265492733676752398517946425200692808124042667052044810697692986390060537707486364599143525),
     (0.804000000000000,
      7.60503022204000734801611994175910058752179636206651996147193043335584596511723929814175288234159427739125200776567201269508854261086012056659928873030403103232887807707163451358405790873269189649627354113518436152672779483523449395753734423912537899862066995252477954706024219583312729407926179911563),
     (0.805000000000000,
      7.69375161748367094826638210122414507148859404247720099871848487689468214713202082463994969124267991718825786287112188912678100536191574433336629908447710876306600467603662495511710583014759510914202253333621545659830430946106714917813648256271009955928561611769164808920368003266509206523186762233538),
     (0.806000000000000,
      7.78349603858982244721779300896340370719651892262359702987062149173543437496559797298960185871544210722387598361933245542587426523378352635215594880061012627819890365697867503048798428438187426312900797894815467141758355394778799680454552418211658510893648886609371746014810249794042746393852202824007),
     (0.807000000000000,
      7.87427517006801615180146452499849061927596528455360826721739657562644794302766771087638979639602449116574992077357141394603976062924856101812872569383559596854294631178413183904768796450844864789738475061388541959034575255714643195898000598548584414383123038816412410548231650309595727849497961986431),
     (0.808000000000000,
      7.96610082903890794669462114539672129478305173791675241547842421561455113734040619186986942287577193152223869834403630008542757211046233366303597315991416245196747104402192182372699339898663234679669488422469870090098190800507409854987061382181736808634586083222536568831129544352980099058547400185061),
     (0.809000000000000,
      8.05898496652479212401601208637250723827102374629869062269780036497368626661197110233759851953267391453351914055670274814265889392222687741155299567683025316275204849480434304438911646122416249013640276960552932240889258471885514419724654940837698763355241134565033813625655292442352013103404866524715)]

```python
alphaprime = []
yprime = []
for tup in alpha_in_8_to_10:
    alpha, y = tup
    alphaprime.append(log(1-alpha))
    yprime.append(log(y))
```

```python
alphaprime[:5]
```

    [-1.60943791243410,
     -1.61445045425764,
     -1.61948824828760,
     -1.62455155024415,
     -1.62964061975162]

```python
yprime[:5]
```

    [1.98239990508105792432914206584034425361291820164801420571371447279489000805005017210551940329748997949158276953758402139559706347895794635989738106647433102557324635006924478567083315201989122813678874733027136858749688442475437260290841954181034549727611229670749280412218485884946091336845146816553,
     1.99400473775108917242607592672379876282962602758486599176148576772489434639025636624056896545690094280224195809733448464124202508093422016580455567072652267572463293430462118712221270594276058468607769040072042105653811588064721872113903111728451820492918543610398905495968245410433690439953389342567,
     2.00560801181928100231139813531008896336733947147889255572229770276025059061492979832016543643104245292920992912477217573183288125389017314482489589303875686728543697919670636242183802728761179540256596711814791281490102327649928484622757922109228990643321882068823404558914950808859897896899546791065,
     2.01720973117000196837249208594601178602418005180730883439142747157620936053011217274044740680743596917763962963449407257327502987553886810995997876898715267399577354310049293809033854648554959973418779247307890041118469773630291760102404488143229295135303350509732825944937643716389433950777297807287,
     2.02880989967311768574223509850438224081778861713460381994857172328879726399734062796210317666818748089330102932762383113852809191013114520866613652369138698179850038119705032754770196695094060071695316927452740911833479460919282968072815444250088745055725905599645333532839537328901190948844738968695]

```python
zip(alphaprime[:5], yprime[:5])
```

    [(-1.60943791243410,
      1.98239990508105792432914206584034425361291820164801420571371447279489000805005017210551940329748997949158276953758402139559706347895794635989738106647433102557324635006924478567083315201989122813678874733027136858749688442475437260290841954181034549727611229670749280412218485884946091336845146816553),
     (-1.61445045425764,
      1.99400473775108917242607592672379876282962602758486599176148576772489434639025636624056896545690094280224195809733448464124202508093422016580455567072652267572463293430462118712221270594276058468607769040072042105653811588064721872113903111728451820492918543610398905495968245410433690439953389342567),
     (-1.61948824828760,
      2.00560801181928100231139813531008896336733947147889255572229770276025059061492979832016543643104245292920992912477217573183288125389017314482489589303875686728543697919670636242183802728761179540256596711814791281490102327649928484622757922109228990643321882068823404558914950808859897896899546791065),
     (-1.62455155024415,
      2.01720973117000196837249208594601178602418005180730883439142747157620936053011217274044740680743596917763962963449407257327502987553886810995997876898715267399577354310049293809033854648554959973418779247307890041118469773630291760102404488143229295135303350509732825944937643716389433950777297807287),
     (-1.62964061975162,
      2.02880989967311768574223509850438224081778861713460381994857172328879726399734062796210317666818748089330102932762383113852809191013114520866613652369138698179850038119705032754770196695094060071695316927452740911833479460919282968072815444250088745055725905599645333532839537328901190948844738968695)]

```python
var('A')
var('B')
model(x) = B - A*x
```

```python
find_fit(zip(alphaprime, yprime), model)
```

    [A == 1.6670411794480962, B == -0.28990653254233467]

This indicates that the exponent from the linear fit is 1.66, which is not so good. I suspect what's happening is the lack of sufficiently large prime gaps means that the experimental data for alpha very near $1$ is simply incorrect. So let's cut out the highest portion of the range.

```python
# From 0.8 to 0.95
LOWER = 0
UPPER = 150
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
find_fit(zip(alphaprime_9, yprime_9), model)
```

    [A == 2.38952780953096, B == -1.8741114946012887]

```python
# From 0.85 to 0.95
LOWER = 50
UPPER = 150
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
find_fit(zip(alphaprime_9, yprime_9), model)
```

    [A == 2.3838016290728183, B == -1.8595471688135894]

Wow, look at the stability of that result.

```python
# From 0.8 to 0.9
LOWER = 0
UPPER = 100
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
find_fit(zip(alphaprime_9, yprime_9), model)
```

    [A == 2.381841660510989, B == -1.8594294921692536]

```python
# From 0.8 to 0.85
LOWER = 0
UPPER = 50
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
find_fit(zip(alphaprime_9, yprime_9), model)
```

    [A == 2.3446674868637274, B == -1.7982001515169852]

The sheer stability of this in these ranges indicates to me that maybe there is something here. What if we eek just a little bit more out of the upper range?

```python
# From 0.8 to 0.96
LOWER = 0
UPPER = 160
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
find_fit(zip(alphaprime_9, yprime_9), model)
```

    [A == 2.3517182353898307, B == -1.7986473277851818]

```python
# From 0.8 to 0.97
LOWER = 0
UPPER = 170
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
find_fit(zip(alphaprime_9, yprime_9), model)
```

    [A == 2.2790797480846114, B == -1.6504264401594178]

The guess really plummets when we increase from 0.96 to 0.97, so I would guess that this is around where the prime gaps peter off.

Let's use the model data on the range $\alpha \leq 0.8 < 0.95$. On this range, we have

    A' = 2.3895
    B' = -1.874

I use $A'$ and $B'$ to denote that this is what the model estimates. The estimate for $A$ is exactly $A'$, so we guess that $A \approx 2.3895$. The estimate for $B$ is a bit more complicated, and work in my notebook RN5.28 shows that $B \approx \exp(B' - A' \log(A'/e))$.

```python
# From 0.8 to 0.95
var('A')
var('B')
model(x) = B - A*x
LOWER = 0
UPPER = 150
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
res = find_fit(zip(alphaprime_9, yprime_9), model, solution_dict=True)

```

```python
Aprime = res[A]
Bprime = res[B]
print Aprime, Bprime
```

    2.38952780953 -1.8741114946

```python
exper_A = Aprime
exper_B = exp(Bprime - Aprime * log(Aprime / e))
exper_B = exper_B.n()
print exper_A, exper_B
```

    2.38952780953 0.208859966269026

Hmm, interesting. The growth estimate appears to be overcompensating and tossing more into the exponent and taking away from the coefficient.

As an additional experiment, shat if we force the leading coefficient to be $-1$? I'm not quite sure what this would mean, but sometimes one must play a trombone to petunias.

```python
# From 0.8 to 0.95
var('A')
var('B')
model(x) = -1 - A*x
LOWER = 0
UPPER = 150
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
res = find_fit(zip(alphaprime_9, yprime_9), model, solution_dict=True)
Aprime = res[A]
print Aprime
```

    1.99402560961

```python
exper_A = Aprime
exper_B = exp(-1 - Aprime * log(Aprime / e))
exper_B = exper_B.n()
print exper_A, exper_B
```

    1.99402560961 0.682384387756181

It's amusingly close to causing the exponent to be $2$, but otherwise it is not particularly meaningful.

```python

```

max_stars_repo_path: PrimeGaps/Prime Gaps.ipynb

max_stars_repo_name: davidlowryduda/notebooks

max_issues_repo_path: PrimeGaps/Prime Gaps.ipynb

max_issues_repo_name: davidlowryduda/notebooks

max_forks_repo_path: PrimeGaps/Prime Gaps.ipynb

max_forks_repo_name: davidlowryduda/notebooks