<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Features for Trajectory Recommendation
Features
Load Data
Compute POI Info
Construct Travelling Sequences
Compute Some Sequence Statistics
Compute Transition Probabilities
Basic Definitions
Transition Probabilities between POI Categories
Transition Probabilities between POI Popularity Classes
Transition Probabilities between POI Pair Distance Classes
Compute Trajectory Likelihood (For tuning the discretization strategy)
Log Likelihood of Actual Trajectories
Log Likelihood of Enumerated Trajectories
Compare the Log Likelihood of Actual and Enumerated Trajectories
Compute the F1-score of Enumerated Trajectories
<a id='sec1'></a>
1. Features
POI Features ?
POI category
POI popularity
POI location
Transition Features
POI category (a transition matrix between different categories)
POI popularity (a transition matrix between different class of popularity)
POI pair distance (a transition matrix between different class of distance)
Computation Steps
First compute the above features using a set of travelling sequences,
Then compute the log likelihood of (actual/enumerated) sequences using POI category/popularity transition matrix and POI pair distance transition matrix and make comparison.
<a id='sec2'></a>
2. Load Data
Step1:  <a id='sec2.1'></a>
2.1 Compute POI Info
Compute POI (Longitude, Latitude) as the average coordinates of the assigned photos.
Step2:  Extract POI category and visiting frequency.
Step3:  <a id='sec2.2'></a>
2.2 Construct Travelling Sequences
Step4:  <a id='sec2.3'></a>
2.3 Compute Some Sequence Statistics
Step5:  Sequences with length {3, 4, 5}
Step6:  <a id='sec3'></a>
3. Compute Transition Probabilities
<a id='sec3.1'></a>
3.1 Basic Definitions
$\text{Pr}(\text{POI}_i \to \text{POI}_j)$
Step7:  <a id='sec3.3'></a>
3.3 Transition Probabilities between POI Popularity Classes
We model transition probabilities between POI popularities, i.e.
$\text{Pr}(\text{Pop}{\text{POI}_i} \to \text{Pop}{\text{POI}_j})$
after discretizing POI popularities.
3.3.1 Discretize POI Popularity
What the general criteria of data discretization? Discretization of continuous features from Wikipedia. 
~~TODO
Step8:  It could be seen from the above plot that discretization based on equal frequency (quantiles) performs better than that based on equal width, to balance the complexity and accuracy, we choose "quantile, nbins=9".
Step9:  Quantile based bins (equal frequency)
Step10:  Equal width bins
Step11:  Another Equal frequency bins
Step12:  3.3.2 Compute Transition Probabilities
Step14:  <a id='sec3.4'></a>
3.4 Transition Probabilities between POI Pair Distance Classes
We model transition probabilities between different POI pair distances, i.e.
$\text{Pr}(\text{Dist}{\text{POI}{i-1} \to \text{POI}i} \to \text{Dist}{\text{POI}_{i} \to \text{POI}_j})$
after discretize POI pair distances.
TODO
Step15:  3.4.2 Discretize POI Pair Distance
We use the same metrics as described above to choose a discretization strategy for POI pair distances.
<table>
<tr><td><b>descretization strategy</b></td><td><b>actSeqRank(Top%) Toronto (smaller is better)</b></td><td><b>rankedTop5(%) Toronto (larger is better)</b></td><td><b>actSeqRank(Top%) Glasgow</b></td><td><b>rankedTop5(%) Glasgow</b></td></tr>
<tr><td>quantile, nbins=2</td><td>mean
Step16:  Remove rows that contain NaN and plot the curve, we choose quantile, nbins=10 to balance the complexity and accuracy.
Step17:  Quantile based bins (equal frequency)
Step18:  Equal width bins
Step19:  Another Equal frequency bins
Step20:  3.4.3 Compute Transition Probabilities
Use POI pair that is observed in dataset to compute the transition matrix between different "class" of distances.
Step21:  <a id='sec4'></a>
4. Compute Trajectory Likelihood (For tuning the discretization strategy)
Log likelihood of trajectory $[\text{POI}1, \text{POI}_2, \dots, \text{POI}_i, ..., \text{POI}_N]$ is defined as
\begin{align}
\text{logl} =& 
\sum{i=1}^{N-1} \log(\text{Pr}(\text{Cat}{\text{POI}_i} \to \text{Cat}{\text{POI}{i+1}})) + 
\sum{i=1}^{N-1} \log(\text{Pr}(\text{Pop}{\text{POI}_i} \to \text{Pop}{\text{POI}{i+1}})) + 
\sum{i=2}^{N-1} \log(\text{Pr}(\text{Dist}{\text{POI}{i-1} \to \text{POI}i} \to 
\text{Dist}{\text{POI}{i} \to \text{POI}{i+1}})) \
& + \log(\text{Pr}(\text{POI}_1)) 
\end{align}
where $\text{Pr}(\text{POI}_1)$ is the prior of $\text{POI}_1$ and we assume $\text{Pr}(\text{POI}_1)=1.0$, 10-based logarithm is used here.
Step22:  Simple check.
Step23:  <a id='sec4.1'></a>
4.1 Log Likelihood of Actual Trajectories
To save computation power, we consider unique travelling sequences only (i.e. treat sequences with the same list of POI and visiting order of different users as one sequence) as no user-specific features are used here.
Step24:  <a id='sec4.2'></a>
4.2 Log Likelihood of Enumerated Trajectories
Compute log likelihood of enumerated trajectories for all unique actual sequences of length {3, 4, 5}.
Step25:  Enumerate trajectories of the same (start, end) and length (3, 4 or 5) with respect to an actual sequence.
Step26:  Compute the log likelihood of enumerated trajectories.
Step27:  <a id='sec4.3'></a>
4.3 Compare the Log Likelihood of Actual and Enumerated Trajectories
Compare the log likelihood between actual sequences $S_a$ and the one of highest log likelihood among enumerated sequences with respect to $S_a$ as well as the log likelihood rank of $S_a$.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import os
import re
import math
import random
import pickle
import pandas as pd
import numpy as np
import scipy.stats
#from numba import jit
from datetime import datetime
from joblib import Parallel, delayed
import matplotlib.pyplot as plt
nfeatures = 8 # number of features
EPS = 1e-12 # smooth, deal with 0 probability
random.seed(123456789) # control random choice when splitting training/testing set
data_dir = 'data/data-ijcai15'
#fvisit = os.path.join(data_dir, 'userVisits-Osak.csv')
#fcoord = os.path.join(data_dir, 'photoCoords-Osak.csv')
#fvisit = os.path.join(data_dir, 'userVisits-Glas.csv')
#fcoord = os.path.join(data_dir, 'photoCoords-Glas.csv')
#fvisit = os.path.join(data_dir, 'userVisits-Edin.csv')
#fcoord = os.path.join(data_dir, 'photoCoords-Edin.csv')
fvisit = os.path.join(data_dir, 'userVisits-Toro.csv')
fcoord = os.path.join(data_dir, 'photoCoords-Toro.csv')
suffix = fvisit.split('-')[-1].split('.')[0]
visits = pd.read_csv(fvisit, sep=';')
coords = pd.read_csv(fcoord, sep=';')
# merge data frames according to column 'photoID'
assert(visits.shape[0] == coords.shape[0])
traj = pd.merge(visits, coords, on='photoID')
traj.head()
num_photo = traj['photoID'].unique().shape[0]
num_user = traj['userID'].unique().shape[0]
num_poi = traj['poiID'].unique().shape[0]
num_seq = traj['seqID'].unique().shape[0]
pd.DataFrame({'#photo': num_photo, '#user': num_user, '#poi': num_poi, '#seq': num_seq, \
              '#photo/user': num_photo/num_user, '#seq/user': num_seq/num_user}, index=[str(suffix)])
#plt.figure(figsize=[15, 5])
#plt.xlabel('Longitude')
#plt.ylabel('Latitude')
#plt.scatter(traj['photoLon'], traj['photoLat'], marker='+')
Explanation: Features for Trajectory Recommendation
Features
Load Data
Compute POI Info
Construct Travelling Sequences
Compute Some Sequence Statistics
Compute Transition Probabilities
Basic Definitions
Transition Probabilities between POI Categories
Transition Probabilities between POI Popularity Classes
Transition Probabilities between POI Pair Distance Classes
Compute Trajectory Likelihood (For tuning the discretization strategy)
Log Likelihood of Actual Trajectories
Log Likelihood of Enumerated Trajectories
Compare the Log Likelihood of Actual and Enumerated Trajectories
Compute the F1-score of Enumerated Trajectories
<a id='sec1'></a>
1. Features
POI Features ?
POI category
POI popularity
POI location
Transition Features
POI category (a transition matrix between different categories)
POI popularity (a transition matrix between different class of popularity)
POI pair distance (a transition matrix between different class of distance)
Computation Steps
First compute the above features using a set of travelling sequences,
Then compute the log likelihood of (actual/enumerated) sequences using POI category/popularity transition matrix and POI pair distance transition matrix and make comparison.
<a id='sec2'></a>
2. Load Data
End of explanation
poi_coords = traj[['poiID', 'photoLon', 'photoLat']].groupby('poiID').mean()
poi_coords.reset_index(inplace=True)
poi_coords.rename(columns={'photoLon':'poiLon', 'photoLat':'poiLat'}, inplace=True)
Explanation: <a id='sec2.1'></a>
2.1 Compute POI Info
Compute POI (Longitude, Latitude) as the average coordinates of the assigned photos.
End of explanation
poi_catfreq = traj[['poiID', 'poiTheme', 'poiFreq']].groupby('poiID').first()
poi_catfreq.reset_index(inplace=True)
poi_all = pd.merge(poi_catfreq, poi_coords, on='poiID')
poi_all.set_index('poiID', inplace=True)
#poi_all.to_csv(fpoi, index=True)
Explanation: Extract POI category and visiting frequency.
End of explanation
seq_all = traj[['userID', 'seqID', 'poiID', 'dateTaken']].copy().groupby(['userID', 'seqID', 'poiID'])\
          .agg([np.min, np.max, np.size])
seq_all.columns = seq_all.columns.droplevel()
seq_all.reset_index(inplace=True)
seq_all.rename(columns={'amin':'arrivalTime', 'amax':'departureTime', 'size':'#photo'}, inplace=True)
seq_all['poiDuration(sec)'] = seq_all['departureTime'] - seq_all['arrivalTime']
#seq_all.head()
seq_user = seq_all[['userID', 'seqID', 'poiID']].copy().groupby(['userID', 'seqID']).agg(np.size)
seq_user.reset_index(inplace=True)
seq_user.rename(columns={'size':'seqLen'}, inplace=True)
seq_user.set_index('seqID', inplace=True)
#seq_user.head()
Explanation: <a id='sec2.2'></a>
2.2 Construct Travelling Sequences
End of explanation
seq_len = seq_all[['seqID', 'poiID']].copy().groupby('seqID').agg(np.size)
seq_len.reset_index(inplace=True)
seq_len.rename(columns={'poiID':'seqLen'}, inplace=True)
#seq_len.head()
seq_stats = seq_all[['seqID', '#photo', 'poiDuration(sec)']].copy().groupby('seqID').agg(np.sum)
seq_stats.reset_index(inplace=True)
#seq_stats.rename(columns={'poiDuration(sec)':'totalPoiDuration(sec)'}, inplace=True)
seq_stats = pd.merge(seq_len, seq_stats, on='seqID')
seq_stats['poiDuration(sec)'] /= 60
seq_stats.rename(columns={'poiDuration(sec)':'totalPoiDuration(min)'}, inplace=True)
seq_stats.set_index('seqID', inplace=True)
#seq_stats.head()
#ax = seq_stats['seqLen'].hist(bins=50)
#ax.set_xlabel('sequence length')
#ax.set_ylim([0.1, 1e4])
#ax.set_yscale('log')
#ax = seq_stats['#photo'].hist(bins=50)
#ax.set_xlabel('#photo for sequence')
#ax.set_ylim([0.1, 1e4])
#ax.set_yscale('log')
#ax = seq_stats['totalPoiDuration(min)'].hist(bins=100)
#ax.set_xlabel('totalPoiDuration(min)')
#ax.set_ylim([0.1, 1e4])
#ax.set_yscale('log')
#ax.set_xscale('log')
Explanation: <a id='sec2.3'></a>
2.3 Compute Some Sequence Statistics
End of explanation
#seq_stats = seq_stats[seq_stats['seqLen'].isin({3, 4, 5})]
#ax = seq_stats['totalPoiDuration(min)'].hist(bins=50)
#ax.set_xlabel('totalPoiDuration(min)')
#ax.set_ylim([0.1, 1e4])
#ax.set_yscale('log')
def extract_seq(seqid, seq_all):
    seqi = seq_all[seq_all['seqID'] == seqid].copy()
    seqi.sort(columns=['arrivalTime'], ascending=True, inplace=True)
    return seqi['poiID'].tolist()
Explanation: Sequences with length {3, 4, 5}
End of explanation
def calc_poi_cat_transmat(seqid_set, poi_all, seq_all):
    poi_cats = poi_all['poiTheme'].unique().tolist()
    poi_cats.sort()
    poi_cat_transmat = pd.DataFrame(data=np.zeros((len(poi_cats), len(poi_cats)), dtype=np.float), \
                                    index=poi_cats, columns=poi_cats)
    for seqid in seqid_set:
        seq = extract_seq(seqid, seq_all)
        for j in range(len(seq)-1):
            poi1 = seq[j]
            poi2 = seq[j+1]
            cat1 = poi_all.loc[poi1, 'poiTheme']
            cat2 = poi_all.loc[poi2, 'poiTheme']
            poi_cat_transmat.loc[cat1, cat2] += 1
    return poi_cat_transmat
def normalise_transmat(transmat):
    assert(isinstance(transmat, pd.DataFrame))
    for row in range(transmat.index.shape[0]):
        nonzeroidx = np.nonzero(transmat.iloc[row])[0].tolist()
        if len(nonzeroidx) < transmat.columns.shape[0]:
            minv = np.min(transmat.iloc[row, nonzeroidx])
            EPS = 0.1 * minv  # row-specific smooth factor
            #zeroidx = list(set(range(len(transmat.columns))) - set(nonzeroidx))
            #transmat.iloc[row, zeroidx] = EPS
            transmat.iloc[row] += EPS
        rowsum = np.sum(transmat.iloc[row])
        assert(rowsum > 0)
        transmat.iloc[row] /= rowsum
    return transmat        
poi_cat_transmat = calc_poi_cat_transmat(seq_all['seqID'].unique(), poi_all, seq_all)
poi_cat_transmat
poi_cat_transmat = normalise_transmat(poi_cat_transmat)
poi_cat_transmat
poi_cat_transmat_log = np.log10(poi_cat_transmat)
poi_cat_transmat_log
Explanation: <a id='sec3'></a>
3. Compute Transition Probabilities
<a id='sec3.1'></a>
3.1 Basic Definitions
$\text{Pr}(\text{POI}_i \to \text{POI}_j)$: 
   the transition probability from $\text{POI}_i$ to $\text{POI }_j$
$\text{Pr}(\text{Cat}_i \to \text{Cat}_j)$:
   the transition probability from a POI of category $\text{Cat}_i$ to a POI of category $\text{Cat}_j$
$\text{Pr}(\text{Pop}_i \to \text{Pop}_j)$:
   the transition probability from a POI of Popularity class $\text{Pop}_i$ to a POI of Popularity class $\text{Pop}_j$
$\text{Pr}(\text{Dist}_i \to \text{Dist}_j)$:
   the transition probability from a POI-POI pair with distance (between the two) class $\text{Dist}_i$ to a POI-POI pair with distance (between the two) class $\text{Dist}_j$
By design, $\text{Pr}(\text{POI}_i \to \text{POI}_j)$ should be bigger 
if any of $\text{Pr}(\text{Cat}_i \to \text{Cat}_j)$, $\text{Pr}(\text{Pop}_i \to \text{Pop}_j)$ and
$\text{Pr}(\text{Dist}_i \to \text{Dist}_j)$ becomes bigger (if other factors stay the same).
So how to combine these probabilities?
Both addition and multiplication seems to be able to serve this purpose, what is the difference?
The Addition Case
For the addtion case, 
\begin{equation}
\text{Pr}(\text{POI}i \to \text{POI}_j) = \frac{
\text{Pr}(\text{Cat}{\text{POI}i} \to \text{Cat}{\text{POI}j}) + 
\text{Pr}(\text{Pop}{\text{POI}i} \to \text{Pop}{\text{POI}j}) + 
\text{Pr}(\text{Dist}{\text{POI}{i-1} \to \text{POI}_i} \to \text{Dist}{\text{POI}_{i} \to \text{POI}_j})
}
{Z_i}
\end{equation}
where $\text{POI}_{i-1}$ is the direct predecessor of $\text{POI}_i$ in a trajectory and $Z_i$ is a normalizing constant.
The Multiplication Case
For the multiplication case, 
\begin{equation}
\text{Pr}(\text{POI}i \to \text{POI}_j) = \frac{
\text{Pr}(\text{Cat}{\text{POI}i} \to \text{Cat}{\text{POI}j}) \times
\text{Pr}(\text{Pop}{\text{POI}i} \to \text{Pop}{\text{POI}j}) \times
\text{Pr}(\text{Dist}{\text{POI}{i-1} \to \text{POI}_i} \to \text{Dist}{\text{POI}_{i} \to \text{POI}_j})
}
{Z_i}
\end{equation}
similarly, $\text{POI}_{i-1}$ is the direct predecessor of $\text{POI}_i$ in a trajectory and $Z_i$ is again a normalizing constant.
The Difference between Addition and Multiplication
It is important to note the fact that, by design, $\text{Pr}(\text{POI}_i \to \text{POI}_j)$ should be very small
if any of $\text{Pr}(\text{Cat}_i \to \text{Cat}_j)$, $\text{Pr}(\text{Pop}_i \to \text{Pop}_j)$ and
$\text{Pr}(\text{Dist}_i \to \text{Dist}_j)$ is very small, in the extreme case, if any of the three probabilities is $0$, then $\text{Pr}(\text{POI}_i \to \text{POI}_j)$ should be $0$
because the event "Transition from POI$_i$ to POI$_j$" is impossible.
From the equation of the addition case, it is clear that the addition rule contradicts the above fact while the multiplication rule is consistent with it.
Intuitively, the addition rule could make an unlikely event become much more likely, specifically, 
make an impossible event become possible.
<a id='sec3.2'></a>
3.2 Transition Probabilities between POI Categories
We model transition probabilities between POI categories, i.e.
$\text{Pr}(\text{Cat}{\text{POI}_i} \to \text{Cat}{\text{POI}_j})$.
We count the number of transition first, then normalise each row while taking care of zero by adding each cell a small number (i.e. $0.1$ times the minimum value of that row) if there exists a zero cell.
End of explanation
rank_mean_Toro = [18.496, 16.049, 16.478, 16.811, 15.049, 15.831, 15.567, 14.556, 14.398, 14.874, 13.491, 13.528, \
                  13.279, 12.784, 24.761, 21.841, 23.636, 20.154, 19.129, 16.922, 18.240, 17.507, 18.196, 17.711, \
                  17.389, 15.820, 15.681, 15.712, 15.977]
rank_std_Toro = [24.746, 19.873, 21.002, 22.159, 18.722, 20.140, 19.811, 17.937, 17.694, 18.915, 16.760, 16.964, \
                 15.960, 15.796, 29.356, 28.366, 29.252, 26.713, 25.530, 20.619, 22.739, 22.659, 23.321, 23.398, \
                 22.718, 20.856, 19.559, 19.794, 20.373]
rank_mean_Glas = [16.328, 16.188, 15.130, 14.316, 14.581, 14.777, 15.019, 14.255, 13.771, 13.568, 11.464, 12.416, \
                  12.596, 12.209, 23.705, 23.225, 19.416, 19.201, 19.907, 17.030, 19.977, 18.183, 18.158, 15.555, \
                  15.174, 15.184, 12.922, 14.274, 20.427]
rank_std_Glas = [19.763, 18.117, 18.643, 17.294, 17.522, 17.175, 17.101, 16.184, 16.043, 15.522, 13.896, 13.081, \
                 14.259, 13.527, 24.304, 25.065, 22.059, 23.250, 23.132, 19.898, 23.118, 22.388, 21.773, 19.722, \
                 17.188, 18.837, 15.119, 16.828, 21.596]
rank_top5_Toro = [49.25, 49.25, 49.62, 47.74, 48.12, 48.87, 50.75, 51.88, 50.75, 51.88, 54.51, 54.14, \
                  53.38, 56.77, 43.23, 47.74, 45.11, 46.99, 46.99, 47.37, 45.49, 48.50, 46.99, 49.25, \
                  49.62, 53.01, 51.13, 51.13, 49.62]
rank_top5_Glas = [55, 52, 58, 61, 58, 59, 56, 64, 62, 63, 66, 63, 65, 66, 47, 47, 53, 54, 53, 57, 54, 55, 52, 60, 60, \
                  61, 63, 61, 50]
#xlabels = ['qbins='+str(x) for x in range(2, 16)]
#xlabels.extend(['ebins='+str(x) for x in range(2, 16)])
#xlabels.append('another')
xlabels = [x for x in range(2, len(rank_mean_Toro)+2)]
plt.figure(figsize=[15, 10])
plt.xlim([0, len(rank_mean_Toro)+2])
plt.ylim([-20, 100])
plt.errorbar(xlabels, rank_mean_Toro, rank_std_Toro, linestyle='--', marker='s', label='errorbar_Toronto')
plt.errorbar(xlabels, rank_mean_Glas, rank_std_Glas, linestyle='--', marker='s', label='errorbar_Glasgow')
plt.plot(xlabels, rank_top5_Toro, linestyle='--', marker='s', label='top5_Toronto')
plt.plot(xlabels, rank_top5_Glas, linestyle='--', marker='s', label='top5_Glasgow')
plt.legend()
#idx = 10
idx = 7
plt.annotate('choose', xy=(xlabels[idx], rank_top5_Glas[idx]), xytext=(xlabels[idx], rank_top5_Glas[idx]+15), \
             arrowprops=dict(facecolor='green', shrink=0.1))
Explanation: <a id='sec3.3'></a>
3.3 Transition Probabilities between POI Popularity Classes
We model transition probabilities between POI popularities, i.e.
$\text{Pr}(\text{Pop}{\text{POI}_i} \to \text{Pop}{\text{POI}_j})$
after discretizing POI popularities.
3.3.1 Discretize POI Popularity
What the general criteria of data discretization? Discretization of continuous features from Wikipedia. 
~~TODO: Improve the discritization using Fayyad's minimum description length principle (MDLP)  described by this paper~~.
NOTE: MDLP is a supervised discretization method which requires some sort of class labels that are not available here. We need an unsupervised discretization method, well-known ones including equal width, equal frequency, clustering based, etc.
Try different discretization strategies and choose the best one using the following metrics:
actSeqRank(Top%) = $100 \times \frac{\text{rank of actual sequence}~S_a}{\text{number of enumerated sequence w.r.p}~S_a}$
rankedTop5(%) = $100 \times \frac{\sum{\delta(\text{rank of actual sequence} \le 5)}}
{\text{number of actual sequence}}$, where $\delta(\text{True}) = 1$ and $\delta(\text{False}) = 0$
i.e. the percentage of actual sequence $S_a$ ranked top 5 among all enumerated sequences with respect to $S_a$.
<table>
<tr><td><b>descretization strategy</b></td><td><b>actSeqRank(Top%) Toronto (smaller is better)</b></td><td><b>rankedTop5(%) Toronto (larger is better)</b></td><td><b>actSeqRank(Top%) Glasgow</b></td><td><b>rankedTop5(%) Glasgow</b></td></tr>
<tr><td>quantile, nbins=2</td><td>mean: 18.496, std: 24.746</td><td>131/266 = 49.25%</td><td>mean: 16.328, std: 19.763</td><td>55/100 = 55.00%</td></tr>
<tr><td>quantile, nbins=3</td><td>mean: 16.049, std: 19.873</td><td>131/266 = 49.25%</td><td>mean: 16.188, std: 18.117</td><td>52/100 = 52.00%</td></tr>
<tr><td>quantile, nbins=4</td><td>mean: 16.478, std: 21.002</td><td>132/266 = 49.62%</td><td>mean: 15.130, std: 18.643</td><td>58/100 = 58.00%</td></tr>
<tr><td>quantile, nbins=5</td><td>mean: 16.811, std: 22.159</td><td>127/266 = 47.74%</td><td>mean: 14.316, std: 17.294</td><td>61/100 = 61.00%</td></tr>
<tr><td>quantile, nbins=6</td><td>mean: 15.049, std: 18.722</td><td>128/266 = 48.12%</td><td>mean: 14.581, std: 17.522</td><td>58/100 = 58.00%</td></tr>
<tr><td>quantile, nbins=7</td><td>mean: 15.831, std: 20.140</td><td>130/266 = 48.87%</td><td>mean: 14.777, std: 17.175</td><td>59/100 = 59.00%</td></tr>
<tr><td>quantile, nbins=8</td><td>mean: 15.567, std: 19.811</td><td>135/266 = 50.75%</td><td>mean: 15.019, std: 17.101</td><td>56/100 = 56.00%</td></tr>
<tr><td>quantile, nbins=9</td><td>mean: 14.556, std: 17.937</td><td>138/266 = 51.88%</td><td>mean: 14.255, std: 16.184</td><td>64/100 = 64.00%</td></tr>
<tr><td>quantile, nbins=10</td><td>mean: 14.398, std: 17.694</td><td>135/266 = 50.75%</td><td>mean: 13.771, std: 16.043</td><td>62/100 = 62.00%</td></tr>
<tr><td>quantile, nbins=11</td><td>mean: 14.874, std: 18.915</td><td>138/266 = 51.88%</td><td>mean: 13.568, std: 15.522</td><td>63/100 = 63.00%</td></tr>
<tr><td>quantile, nbins=12</td><td>mean: 13.491, std: 16.760</td><td>145/266 = 54.51%</td><td>mean: 11.464, std: 13.896</td><td>66/100 = 66.00%</td></tr>
<tr><td>quantile, nbins=13</td><td>mean: 13.528, std: 16.964</td><td>144/266 = 54.14%</td><td>mean: 12.416, std: 13.081</td><td>63/100 = 63.00%</td></tr>
<tr><td>quantile, nbins=14</td><td>mean: 13.279, std: 15.960</td><td>142/266 = 53.38%</td><td>mean: 12.596, std: 14.259</td><td>65/100 = 65.00%</td></tr>
<tr><td>quantile, nbins=15</td><td>mean: 12.784, std: 15.796</td><td>151/266 = 56.77%</td><td>mean: 12.209, std: 13.527</td><td>66/100 = 66.00%</td></tr>
<tr><td>equalWidth, nbins=2</td><td>mean: 24.761, std: 29.356</td><td>115/266 = 43.23%</td><td>mean: 23.705, std: 24.304</td><td>47/100 = 47.00%</td></tr>
<tr><td>equalWidth, nbins=3</td><td>mean: 21.841, std: 28.366</td><td>127/266 = 47.74%</td><td>mean: 23.225, std: 25.065</td><td>47/100 = 47.00%</td></tr>
<tr><td>equalWidth, nbins=4</td><td>mean: 23.636, std: 29.252</td><td>120/266 = 45.11%</td><td>mean: 19.416, std: 22.059</td><td>53/100 = 53.00%</td></tr>
<tr><td>equalWidth, nbins=5</td><td>mean: 20.154, std: 26.713</td><td>125/266 = 46.99%</td><td>mean: 19.201, std: 23.250</td><td>54/100 = 54.00%</td></tr>
<tr><td>equalWidth, nbins=6</td><td>mean: 19.129, std: 25.530</td><td>125/266 = 46.99%</td><td>mean: 19.907, std: 23.132</td><td>53/100 = 53.00%</td></tr>
<tr><td>equalWidth, nbins=7</td><td>mean: 16.922, std: 20.619</td><td>126/266 = 47.37%</td><td>mean: 17.030, std: 19.898</td><td>57/100 = 57.00%</td></tr>
<tr><td>equalWidth, nbins=8</td><td>mean: 18.240, std: 22.739</td><td>121/266 = 45.49%</td><td>mean: 19.977, std: 23.118</td><td>54/100 = 54.00%</td></tr>
<tr><td>equalWidth, nbins=9</td><td>mean: 17.507, std: 22.659</td><td>129/266 = 48.50%</td><td>mean: 18.183, std: 22.388</td><td>55/100 = 55.00%</td></tr>
<tr><td>equalWidth, nbins=10</td><td>mean: 18.196, std: 23.321</td><td>125/266 = 46.99%</td><td>mean: 18.158, std: 21.773</td><td>52/100 = 52.00%</td></tr>
<tr><td>equalWidth, nbins=11</td><td>mean: 17.711, std: 23.398</td><td>131/266 = 49.25%</td><td>mean: 15.555, std: 19.722</td><td>60/100 = 60.00%</td></tr>
<tr><td>equalWidth, nbins=12</td><td>mean: 17.389, std: 22.718</td><td>132/266 = 49.62%</td><td>mean: 15.174, std: 17.188</td><td>60/100 = 60.00%</td></tr>
<tr><td>equalWidth, nbins=13</td><td>mean: 15.820, std: 20.856</td><td>141/266 = 53.01%</td><td>mean: 15.184, std: 18.837</td><td>61/100 = 61.00%</td></tr>
<tr><td>equalWidth, nbins=14</td><td>mean: 15.681, std: 19.559</td><td>136/266 = 51.13%</td><td>mean: 12.922, std: 15.119</td><td>63/100 = 63.00%</td></tr>
<tr><td>equalWidth, nbins=15</td><td>mean: 15.712, std: 19.794</td><td>136/266 = 51.13%</td><td>mean: 14.274, std: 16.828</td><td>61/100 = 61.00%</td></tr>
<tr><td>another, bins=[0, 500, 1500, 10000]</td><td>mean: 15.977, std: 20.373</td><td>132/266 = 49.62%</td><td>mean: 20.427, std: 21.596</td><td>50/100 = 50.00%</td></tr>
</table>
End of explanation
poi_all['poiFreq'].get_values()
poi_all['poiFreq'].describe()
poi_all['poiFreq'].quantile([.25, .5, .75]).tolist()
ax = poi_all['poiFreq'].hist(bins=10)
ax.set_xlabel('POI Popularity')
ax.set_ylabel('#POI')
#plt.plot(np.ones(poi_all.index.shape[0]), np.sqrt(poi_all['poiFreq']), marker='+')
Explanation: It could be seen from the above plot that discretization based on equal frequency (quantiles) performs better than that based on equal width, to balance the complexity and accuracy, we choose "quantile, nbins=9".
End of explanation
nbins = 9
quantiles = np.round(np.linspace(0, 1, nbins+1), 2)[1:-1]
quantiles
bins_qt = [0]
bins_qt.extend(poi_all['poiFreq'].quantile(quantiles))
bins_qt.append(poi_all['poiFreq'].max() + 1)
bins_qt
Explanation: Quantile based bins (equal frequency)
End of explanation
#nbins = 15
#inter = round((poi_all['poiFreq'].max() + 1) / nbins)
#bins_ew = [x*inter for x in range(nbins)]
#bins_ew.append(poi_all['poiFreq'].max() + 1)
#bins_ew
Explanation: Equal width bins
End of explanation
#bins = np.linspace(0, 10000, 11)
#bins = np.logspace(0, 4, 5)
#bins = [1, 100, 500, 1000, 2000, 5000]
#bins_ef = [0, 500, 1500, 10000]
bins_pop = bins_qt
#bins_pop = bins_ew
#bins_pop = bins_ef
ax = poi_all['poiFreq'].hist(bins=bins_pop)
ax.set_xlabel('POI Popularity')
ax.set_ylabel('#POI')
ax.set_xscale('log')
poi_all['popClass'] = np.digitize(poi_all['poiFreq'].get_values(), bins_pop)
#poi_all
Explanation: Another Equal frequency bins
End of explanation
def calc_poi_pop_transmat(seqid_set, poi_all, seq_all):
    pop_class = poi_all['popClass'].unique().tolist()
    pop_class.sort()
    poi_pop_transmat = pd.DataFrame(data=np.zeros((len(pop_class), len(pop_class)), dtype=np.float), \
                                    index=pop_class, columns=pop_class)
    for seqid in seqid_set:
        seq = extract_seq(seqid, seq_all)
        for j in range(len(seq)-1):
            poi1 = seq[j]
            poi2 = seq[j+1]
            pc1 = poi_all.loc[poi1, 'popClass']
            pc2 = poi_all.loc[poi2, 'popClass']
            poi_pop_transmat.loc[pc1, pc2] += 1
    return poi_pop_transmat
poi_pop_transmat = calc_poi_pop_transmat(seq_all['seqID'].unique(), poi_all, seq_all)
poi_pop_transmat
poi_pop_transmat = normalise_transmat(poi_pop_transmat)
poi_pop_transmat
poi_pop_transmat_log = np.log10(poi_pop_transmat)
poi_pop_transmat_log
Explanation: 3.3.2 Compute Transition Probabilities
End of explanation
def calc_dist(longitude1, latitude1, longitude2, latitude2):
    Calculate the distance (unit: km) between two places on earth
    # convert degrees to radians
    lon1 = math.radians(longitude1)
    lat1 = math.radians(latitude1)
    lon2 = math.radians(longitude2)
    lat2 = math.radians(latitude2)
    radius = 6371.009 # mean earth radius is 6371.009km, en.wikipedia.org/wiki/Earth_radius#Mean_radius
    # The haversine formula, en.wikipedia.org/wiki/Great-circle_distance
    dlon = math.fabs(lon1 - lon2)
    dlat = math.fabs(lat1 - lat2)
    return 2 * radius * math.asin( math.sqrt( \
               (math.sin(0.5*dlat))**2 + math.cos(lat1) * math.cos(lat2) * (math.sin(0.5*dlon))**2 ))
def calc_obs_poipair_distmat(seqid_set, poi_all, seq_all):
    poi_distmat = pd.DataFrame(data=np.full((poi_all.shape[0], poi_all.shape[0]), np.nan, dtype=np.float), \
                               index=poi_all.index, columns=poi_all.index)
    for seqid in seqid_set:
        seq = extract_seq(seqid, seq_all)
        if len(seq) < 2: continue
        for j in range(len(seq)-1):
            poi1 = seq[j]
            poi2 = seq[j+1]
            if np.isnan(poi_distmat.loc[poi1, poi2]):
                dist = calc_dist(poi_all.loc[poi1, 'poiLon'], poi_all.loc[poi1, 'poiLat'], \
                                 poi_all.loc[poi2, 'poiLon'], poi_all.loc[poi2, 'poiLat'])
                poi_distmat.loc[poi1, poi2] = dist
                poi_distmat.loc[poi2, poi1] = dist
    return poi_distmat
poi_distmat = calc_obs_poipair_distmat(seq_all['seqID'].unique(), poi_all, seq_all)
#poi_distmat
Explanation: <a id='sec3.4'></a>
3.4 Transition Probabilities between POI Pair Distance Classes
We model transition probabilities between different POI pair distances, i.e.
$\text{Pr}(\text{Dist}{\text{POI}{i-1} \to \text{POI}i} \to \text{Dist}{\text{POI}_{i} \to \text{POI}_j})$
after discretize POI pair distances.
TODO: Improve the distance calculation using Google maps distance API with different travel modes demonstrated here.
3.4.1 Compute POI Pair Distance
Compute POI-pair distance if the pair is observed in dataset.
End of explanation
rank_mean_Toro = [15.156, 14.917, 14.389, 13.645, 14.299, 12.689, 12.996, 12.510, 12.467, 12.548, 11.980, \
                  12.170, 11.384, 11.444, 10.932, 10.991, 10.836, 15.110]
rank_std_Toro = [16.911, 17.484, 17.527, 16.550, 17.550, 14.674, 15.606, 14.946, 14.549, 14.758, 13.883, \
                 13.983, 12.787, 12.888, 12.621, 12.950, 12.383, 16.204]
rank_mean_Glas = [14.354, 14.450, 14.255, 14.085, 13.156, 12.755, 11.716, 11.355, 11.181, 10.214, 10.041, \
                  9.345, 9.008, 8.613, 8.553, 8.025, 7.922, 14.937]
rank_std_Glas = [16.541, 16.173, 16.184, 16.159, 14.436, 13.681, 13.003, 12.893, 12.107, 10.872, 10.895, \
                 9.759, 9.444, 9.052, 9.336, 8.502, 8.702, 16.630]
rank_top5_Toro = [45.11, 48.87, 52.26, 55.64, 51.88, 52.63, 53.76, 52.63, 54.89, 53.76, 55.26, 52.63, 57.52, \
                  56.77, 58.65, 58.65, 57.14, 43.98]
rank_top5_Glas = [60, 60, 64, 63, 63, 62, 64, 66, 68, 68, 73, 73, 73, 75, 74, 79, 79, 60]
xlabels = [x for x in range(2, len(rank_mean_Toro)+2)]
plt.figure(figsize=[15, 10])
plt.xlim([0, len(rank_mean_Toro)+2])
plt.ylim([-20, 100])
plt.errorbar(xlabels, rank_mean_Toro, rank_std_Toro, linestyle='--', marker='s', label='errorbar_Toronto')
plt.errorbar(xlabels, rank_mean_Glas, rank_std_Glas, linestyle='--', marker='s', label='errorbar_Glasgow')
plt.plot(xlabels, rank_top5_Toro, linestyle='--', marker='s', label='top5_Toronto')
plt.plot(xlabels, rank_top5_Glas, linestyle='--', marker='s', label='top5_Glasgow')
plt.legend()
#idx = 10
idx = 8
plt.annotate('choose', xy=(xlabels[idx], rank_top5_Glas[idx]), xytext=(xlabels[idx], rank_top5_Glas[idx]+15), \
             arrowprops=dict(facecolor='green', shrink=0.1))
Explanation: 3.4.2 Discretize POI Pair Distance
We use the same metrics as described above to choose a discretization strategy for POI pair distances.
<table>
<tr><td><b>descretization strategy</b></td><td><b>actSeqRank(Top%) Toronto (smaller is better)</b></td><td><b>rankedTop5(%) Toronto (larger is better)</b></td><td><b>actSeqRank(Top%) Glasgow</b></td><td><b>rankedTop5(%) Glasgow</b></td></tr>
<tr><td>quantile, nbins=2</td><td>mean: 15.156, std: 16.911</td><td>120/266 = 45.11%</td><td>mean: 14.354, std: 16.541</td><td>60/100 = 60.00%</td></tr>
<tr><td>quantile, nbins=3</td><td>mean: 14.917, std: 17.484</td><td>130/266 = 48.87%</td><td>mean: 14.450, std: 16.173</td><td>60/100 = 60.00%</td></tr>
<tr><td>quantile, nbins=4</td><td>mean: 14.389, std: 17.527</td><td>139/266 = 52.26%</td><td>mean: 14.255, std: 16.184</td><td>64/100 = 64.00%</td></tr>
<tr><td>quantile, nbins=5</td><td>mean: 13.645, std: 16.550</td><td>148/266 = 55.64%</td><td>mean: 14.085, std: 16.159</td><td>63/100 = 63.00%</td></tr>
<tr><td>quantile, nbins=6</td><td>mean: 14.299, std: 17.550</td><td>138/266 = 51.88%</td><td>mean: 13.156, std: 14.436</td><td>63/100 = 63.00%</td></tr>
<tr><td>quantile, nbins=7</td><td>mean: 12.689, std: 14.674</td><td>140/266 = 52.63%</td><td>mean: 12.755, std: 13.681</td><td>62/100 = 62.00%</td></tr>
<tr><td>quantile, nbins=8</td><td>mean: 12.996, std: 15.606</td><td>143/266 = 53.76%</td><td>mean: 11.716, std: 13.003</td><td>64/100 = 64.00%</td></tr>
<tr><td>quantile, nbins=9</td><td>mean: 12.510, std: 14.946</td><td>140/266 = 52.63%</td><td>mean: 11.355, std: 12.893</td><td>66/100 = 66.00%</td></tr>
<tr><td>quantile, nbins=10</td><td>mean: 12.467, std: 14.549</td><td>146/266 = 54.89%</td><td>mean: 11.181, std: 12.107</td><td>68/100 = 68.00%</td></tr>
<tr><td>quantile, nbins=11</td><td>mean: 12.548, std: 14.758</td><td>143/266 = 53.76%</td><td>mean: 10.214, std: 10.872</td><td>68/100 = 68.00%</td></tr>
<tr><td>quantile, nbins=12</td><td>mean: 11.980, std: 13.883</td><td>147/266 = 55.26%</td><td>mean: 10.041, std: 10.895</td><td>73/100 = 73.00%</td></tr>
<tr><td>quantile, nbins=13</td><td>mean: 12.170, std: 13.983</td><td>140/266 = 52.63%</td><td>mean: 9.345, std: 9.759</td><td>73/100 = 73.00%</td></tr>
<tr><td>quantile, nbins=14</td><td>mean: 11.384, std: 12.787</td><td>153/266 = 57.52%</td><td>mean: 9.008, std: 9.444</td><td>73/100 = 73.00%</td></tr>
<tr><td>quantile, nbins=15</td><td>mean: 11.444, std: 12.888</td><td>151/266 = 56.77%</td><td>mean: 8.613, std: 9.052</td><td>75/100 = 75.00%</td></tr>
<tr><td>quantile, nbins=16</td><td>mean: 10.932, std: 12.621</td><td>156/266 = 58.65%</td><td>mean: 8.553, std: 9.336</td><td>74/100 = 74.00%</td></tr>
<tr><td>quantile, nbins=17</td><td>mean: 10.991, std: 12.950</td><td>156/266 = 58.65%</td><td>mean: 8.025, std: 8.502</td><td>79/100 = 79.00%</td></tr>
<tr><td>quantile, nbins=18</td><td>mean: 10.728, std: 12.285</td><td>155/266 = 58.27%</td><td>EMPTY ROW 13</td><td>NaN</td></tr>
<tr><td>quantile, nbins=19</td><td>mean: 10.754, std: 12.368</td><td>150/266 = 56.39%</td><td>EMPTY ROW 11</td><td>NaN</td></tr>
<tr><td>quantile, nbins=20</td><td>mean: 10.836, std: 12.383</td><td>152/266 = 57.14%</td><td>mean: 7.922, std: 8.702</td><td>79/100 = 79.00%</td></tr>
<tr><td>quantile, nbins=21</td><td>mean: 10.579, std: 12.301</td><td>162/266 = 60.90%</td><td>EMPTY ROW 15</td><td>NaN</td></tr>
<tr><td>equalWidth, nbins=2</td><td>EMPTY LAST BIN</td><td>NaN</td><td>EMPTY LAST BIN</td><td>NaN</td></tr>
<tr><td>equalWidth, nbins=3</td><td>EMPTY LAST BIN</td><td>NaN</td><td>EMPTY LAST TWO BIN</td><td>NaN</td></tr>
<tr><td>equalWidth, nbins=4</td><td>EMPTY LAST TWO BIN</td><td>NaN</td><td>EMPTY LAST TWO BIN</td><td>NaN</td></tr>
<tr><td>another, bins=[0, 2, 5, 100]</td><td>mean: 15.110, std: 16.204</td><td>117/266 = 43.98%</td><td>mean: 14.937, std: 16.630</td><td>60/100 = 60.00%</td></tr>
</table>
End of explanation
#distdata = pd.Series([x for x in np.unique(poi_distmat.get_values().flatten()) if not np.isnan(x)])
distdata = pd.Series([poi_distmat.iloc[x, y] \
                      for x in range(poi_distmat.index.shape[0]) \
                      for y in range(x+1, poi_distmat.index.shape[0]) \
                      if not np.isnan(poi_distmat.iloc[x, y])])
distdata.describe()
ax = distdata.hist(bins=20)
ax.set_xlabel('POI-Pair Distance (km)')
ax.set_ylabel('#POI-Pair')
Explanation: Remove rows that contain NaN and plot the curve, we choose quantile, nbins=10 to balance the complexity and accuracy.
End of explanation
nbins = 10
quantiles = np.round(np.linspace(0, 1, nbins+1), 2)[1:-1]
quantiles
bins_qt = [0]
bins_qt.extend(distdata.quantile(quantiles))
bins_qt.append(10*round(distdata.max()))
bins_qt
Explanation: Quantile based bins (equal frequency)
End of explanation
#nbins = 4
#inter = round((round(distdata.max()) + 1) / nbins)
#maxdist = 30 # Toronto, maximum distance among all POI pairs
#maxdist = 46 # Glasgow
#inter = round(maxdist / nbins) 
#bins_ew = [x*inter for x in range(nbins)]
#bins_ew.append(maxdist)
#bins_ew
Explanation: Equal width bins
End of explanation
#bins = np.linspace(0, 10, 7)
#bins = np.logspace(0, 2, 4)
#bins = [0, 1, 2, 3, 10]
#bins_a = [0, 2, 5, 100] # walk, ride, drive
#bins_ef = [0, 1.15, 2.25, 100]
bins_dist = bins_qt
#bins_dist = bins_ew
#bins_dist = bins_ef
#bins_dist = bins_a
ax = distdata.ix[np.nonzero(distdata)].hist(bins=bins_dist)
ax.set_xlabel('POI-Pair Distance (km)')
ax.set_ylabel('#POI-Pair')
ax.set_xscale('log')
poi_distclass_mat = pd.DataFrame(data=np.zeros((poi_all.shape[0], poi_all.shape[0]), dtype=np.int), \
                                 index=poi_all.index, columns=poi_all.index)
for i in range(poi_all.index.shape[0]):
    poi1 = poi_all.index[i]
    for j in range(i+1, poi_all.index.shape[0]):
        poi2 = poi_all.index[j]
        dc = None
        if np.isnan(poi_distmat.loc[poi1, poi2]):
            dist = calc_dist(poi_all.loc[poi1, 'poiLon'], poi_all.loc[poi1, 'poiLat'], \
                             poi_all.loc[poi2, 'poiLon'], poi_all.loc[poi2, 'poiLat'])
            dc = np.digitize([dist], bins_dist)[0]
        else:
            dc = np.digitize([poi_distmat.loc[poi1, poi2]], bins_dist)[0]
        assert(dc is not None)
        poi_distclass_mat.loc[poi1, poi2] = dc
        poi_distclass_mat.loc[poi2, poi1] = dc
poi_distclass_mat
Explanation: Another Equal frequency bins
End of explanation
def calc_poipair_dist_transmat(seqid_set, poi_all, seq_all, poi_distclass_mat, bins_dist):
    dist_class = list(range(1, len(bins_dist)))
    poipair_dist_transmat = pd.DataFrame(data=np.zeros((len(dist_class), len(dist_class)), dtype=np.float), \
                                         index=dist_class, columns=dist_class)
    for seqid in seqid_set:
        seq = extract_seq(seqid, seq_all)
        if len(seq) < 3: continue
        for j in range(1, len(seq)-1):
            poi1 = seq[j-1]
            poi2 = seq[j]
            poi3 = seq[j+1]
            dc1 = poi_distclass_mat.loc[poi1, poi2]
            dc2 = poi_distclass_mat.loc[poi2, poi3]
            poipair_dist_transmat.loc[dc1, dc2] += 1
    return poipair_dist_transmat
poipair_dist_transmat = calc_poipair_dist_transmat(seq_all['seqID'].unique(), poi_all, seq_all, \
                                                   poi_distclass_mat, bins_dist)
poipair_dist_transmat
poipair_dist_transmat = normalise_transmat(poipair_dist_transmat)
poipair_dist_transmat
poipair_dist_transmat_log = np.log10(poipair_dist_transmat)
poipair_dist_transmat_log
Explanation: 3.4.3 Compute Transition Probabilities
Use POI pair that is observed in dataset to compute the transition matrix between different "class" of distances.
End of explanation
def calc_seq_loglikelihood(seq, poi_all, poi_cat_transmat_log, poi_pop_transmat_log, \
                           poi_distclass_mat, poipair_dist_transmat_log):
    assert(len(seq) > 1)
    cat1 = poi_all.loc[seq[0], 'poiTheme']
    cat2 = poi_all.loc[seq[1], 'poiTheme']
    pc1  = poi_all.loc[seq[0], 'popClass']
    pc2  = poi_all.loc[seq[1], 'popClass']
    logL = poi_cat_transmat_log.loc[cat1, cat2] + poi_pop_transmat_log.loc[pc1, pc2]
    for j in range(1, len(seq)-1):
        poi1 = seq[j-1]
        poi2 = seq[j]
        poi3 = seq[j+1]
        cat2 = poi_all.loc[poi2, 'poiTheme']
        cat3 = poi_all.loc[poi3, 'poiTheme']
        pc2  = poi_all.loc[poi2, 'popClass']
        pc3  = poi_all.loc[poi3, 'popClass']
        dc12 = poi_distclass_mat.loc[poi1, poi2]
        dc23 = poi_distclass_mat.loc[poi2, poi3]
        logL += poi_cat_transmat_log.loc[cat2, cat3] + poi_pop_transmat_log.loc[pc2, pc3]
        #print(seq, dc12, dc23)
        logL += poipair_dist_transmat_log.loc[dc12, dc23]
    return logL
Explanation: <a id='sec4'></a>
4. Compute Trajectory Likelihood (For tuning the discretization strategy)
Log likelihood of trajectory $[\text{POI}1, \text{POI}_2, \dots, \text{POI}_i, ..., \text{POI}_N]$ is defined as
\begin{align}
\text{logl} =& 
\sum{i=1}^{N-1} \log(\text{Pr}(\text{Cat}{\text{POI}_i} \to \text{Cat}{\text{POI}{i+1}})) + 
\sum{i=1}^{N-1} \log(\text{Pr}(\text{Pop}{\text{POI}_i} \to \text{Pop}{\text{POI}{i+1}})) + 
\sum{i=2}^{N-1} \log(\text{Pr}(\text{Dist}{\text{POI}{i-1} \to \text{POI}i} \to 
\text{Dist}{\text{POI}{i} \to \text{POI}{i+1}})) \
& + \log(\text{Pr}(\text{POI}_1)) 
\end{align}
where $\text{Pr}(\text{POI}_1)$ is the prior of $\text{POI}_1$ and we assume $\text{Pr}(\text{POI}_1)=1.0$, 10-based logarithm is used here.
End of explanation
seq1 = [10, 21, 28, 22]
d12 = calc_dist(poi_all.loc[10,'poiLon'], poi_all.loc[10,'poiLat'], poi_all.loc[21,'poiLon'], poi_all.loc[21, 'poiLat'])
d23 = calc_dist(poi_all.loc[21,'poiLon'], poi_all.loc[21,'poiLat'], poi_all.loc[28,'poiLon'], poi_all.loc[28, 'poiLat'])
d34 = calc_dist(poi_all.loc[28,'poiLon'], poi_all.loc[28,'poiLat'], poi_all.loc[22,'poiLon'], poi_all.loc[22, 'poiLat'])
print(d12, d23, d34)
print(bins_dist)
s1 = poi_cat_transmat_log.loc[poi_all.loc[10, 'poiTheme'], poi_all.loc[21, 'poiTheme']] + \
     poi_cat_transmat_log.loc[poi_all.loc[21, 'poiTheme'], poi_all.loc[28, 'poiTheme']] + \
     poi_cat_transmat_log.loc[poi_all.loc[28, 'poiTheme'], poi_all.loc[22, 'poiTheme']] + \
     poi_pop_transmat_log.loc[poi_all.loc[10, 'popClass'], poi_all.loc[21, 'popClass']] + \
     poi_pop_transmat_log.loc[poi_all.loc[21, 'popClass'], poi_all.loc[28, 'popClass']] + \
     poi_pop_transmat_log.loc[poi_all.loc[28, 'popClass'], poi_all.loc[22, 'popClass']]
s2 = poipair_dist_transmat_log.loc[np.digitize([d12], bins_dist)[0], np.digitize([d23], bins_dist)[0]] + \
     poipair_dist_transmat_log.loc[np.digitize([d23], bins_dist)[0], np.digitize([d34], bins_dist)[0]]
print(s1+s2)
calc_seq_loglikelihood([10, 21, 28, 22], poi_all, poi_cat_transmat_log, poi_pop_transmat_log, \
                       poi_distclass_mat, poipair_dist_transmat_log)
def parse_seqstr(seqstr):
    term = re.sub('[ \[\]]', '', seqstr).split(',')
    return [int(x) for x in term]
Explanation: Simple check.
End of explanation
unique_seq = dict() # seq -> [(seqid, userid)]
for seqid in sorted(seq_all['seqID'].unique().tolist()):
    seq = extract_seq(seqid, seq_all)
    if str(seq) not in unique_seq:
        unique_seq[str(seq)] = [(seqid, seq_user.loc[seqid])]
    else:
        unique_seq[str(seq)].append((seqid, seq_user.loc[seqid]))
unique_seq345 = [parse_seqstr(x) for x in sorted(unique_seq.keys()) if len(x.split(',')) in {3,4,5}]
unique_seq345_logL = pd.DataFrame(data=np.zeros((len(unique_seq345), 2), dtype=np.float), \
                                  index=[str(x) for x in unique_seq345], columns=['logLikelihood', 'seqLen'])
unique_seq345_logL.index.name = 'actSeq'
for seq in unique_seq345:
    assert(len(seq) in {3,4,5})
    logL = calc_seq_loglikelihood(seq, poi_all, poi_cat_transmat_log, poi_pop_transmat_log, \
                                  poi_distclass_mat, poipair_dist_transmat_log)
    unique_seq345_logL.loc[str(seq), 'logLikelihood'] = logL
    unique_seq345_logL.loc[str(seq), 'seqLen'] = len(seq)
    #print('Sequence %-20s Log likelihood: %.3f' % (str(seq), logL))
print(unique_seq345_logL.index.shape[0])
unique_seq345_logL.head()
unique_seq345_logL['seqLen'].hist(bins=10)
Explanation: <a id='sec4.1'></a>
4.1 Log Likelihood of Actual Trajectories
To save computation power, we consider unique travelling sequences only (i.e. treat sequences with the same list of POI and visiting order of different users as one sequence) as no user-specific features are used here.
End of explanation
poi_list = poi_all.index.tolist()
#poi_list
Explanation: <a id='sec4.2'></a>
4.2 Log Likelihood of Enumerated Trajectories
Compute log likelihood of enumerated trajectories for all unique actual sequences of length {3, 4, 5}.
End of explanation
def enum_seq345(seq, poi_list): 
    assert(len(seq) in {3, 4, 5})
    p0 = seq[0]
    pN = seq[-1]
    
    # enumerate sequences with length 3
    if len(seq) == 3:
        return [[p0, p, pN] \
                for p in poi_list if p not in {p0, pN}]
    
    # enumerate sequences with length 4
    if len(seq) == 4:
        return [[p0, p1, p2, pN] \
                for p1 in poi_list if p1 not in {p0, pN} \
                for p2 in poi_list if p2 not in {p0, p1, pN}]
    
    # enumerate sequences with length 5
    if len(seq) == 5:
        return [[p0, p1, p2, p3, pN] \
                for p1 in poi_list if p1 not in {p0, pN} \
                for p2 in poi_list if p2 not in {p0, p1, pN} \
                for p3 in poi_list if p3 not in {p0, p1, p2, pN}]
Explanation: Enumerate trajectories of the same (start, end) and length (3, 4 or 5) with respect to an actual sequence.
End of explanation
enum_logL_df = pd.DataFrame()
for seq in unique_seq345:
    enum_seqs = enum_seq345(seq, poi_list)
    df = pd.DataFrame(data=sorted([str(x) for x in enum_seqs]), columns=['enumSeq'])
    df.set_index('enumSeq', inplace=True)
    df['actSeq'] = str(seq)
    enum_logL_df = enum_logL_df.append(df)
print(enum_logL_df.shape)
enum_logL_df.head()
t1 = datetime.now()
logL = Parallel(n_jobs=-2)(delayed(calc_seq_loglikelihood)\
       (seq, poi_all, poi_cat_transmat_log, poi_pop_transmat_log, poi_distclass_mat, poipair_dist_transmat_log)\
       for seq in [parse_seqstr(x) for x in enum_logL_df.index])
print('%d seconds used' % (datetime.now()-t1).total_seconds()) # 930 seconds
enum_logL_df['enumSeqLogLikelihood'] = logL
#enum_logL_df.head(23)
Explanation: Compute the log likelihood of enumerated trajectories.
End of explanation
df = pd.DataFrame(data=sorted([str(x) for x in unique_seq345]), columns=['actSeq'])
df.set_index('actSeq', inplace=True)
df['actSeqLogLikelihood'] = unique_seq345_logL.loc[df.index, 'logLikelihood']
df['enumSeq'] = ''
df['enumSeqLogLikelihood'] = 0
df['actSeqRank'] = 0
df['#enumSeq'] = 0
for seqstr in df.index:
    sub_df = enum_logL_df[enum_logL_df['actSeq'] == seqstr].copy()
    sub_df.reset_index(inplace=True)
    sub_df.sort(columns=['enumSeqLogLikelihood'], ascending=False, inplace=True)
    df.loc[seqstr, 'enumSeq'] = sub_df.iloc[0]['enumSeq']
    df.loc[seqstr, 'enumSeqLogLikelihood'] = sub_df.iloc[0]['enumSeqLogLikelihood']
    df.loc[seqstr, 'actSeqRank'] = 1 + np.nonzero(sub_df['enumSeq'] == seqstr)[0][0] # rank of actual sequence
    df.loc[seqstr, '#enumSeq'] = sub_df.index.shape[0]
df['actSeqRank(Top%)'] = 100*df['actSeqRank']/df['#enumSeq']
#df
print('mean: %.3f, std: %.3f' % (round(df['actSeqRank(Top%)'].mean(),3), round(df['actSeqRank(Top%)'].std(),3)))
df['actSeqRank(Top%)'].describe()
ntop = np.nonzero(df['actSeqRank'] <= 5)[0].shape[0]
print('%d/%d = %.2f%%' % (ntop, df.index.shape[0], 100*ntop/df.index.shape[0]))
Explanation: <a id='sec4.3'></a>
4.3 Compare the Log Likelihood of Actual and Enumerated Trajectories
Compare the log likelihood between actual sequences $S_a$ and the one of highest log likelihood among enumerated sequences with respect to $S_a$ as well as the log likelihood rank of $S_a$.
End of explanation
<END_TASK>