Skip to content

Commit

Permalink
Add Hawkes earthquake aftershocks study example
Browse files Browse the repository at this point in the history
  • Loading branch information
Mbompr committed Mar 5, 2019
1 parent 6ad84c5 commit 4b8cbef
Show file tree
Hide file tree
Showing 7 changed files with 6,403 additions and 4 deletions.
6,253 changes: 6,253 additions & 0 deletions doc/_static/example_data/earthquakes.csv

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions examples/earthquakes.csv
111 changes: 111 additions & 0 deletions examples/plot_hawkes_earthquake_example.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
"""
==============================================================
Study earthquake aftershocks propagation with Hawkes processes
==============================================================
This experiments aims to detect how aftershocks propagate near Japan. A Hawkes
process has been fit on data from
Ogata, Y., 1988.
Statistical models for earthquake occurrences and residual analysis
for point processes.
Journal of the American Statistical association, 83(401), pp.9-27.
On the left we can see where earthquakes have occurred and on the right
the propagation matrix, i.e. how likely a earthquake in a given zone will
trigger an aftershock in another zone.
We can observe than zone 1, 2 and 3 are tightly linked while zone 4 and
5 are more self-excited.
Note that this Hawkes process does not take the location of an earthquake
to recover this result.
Dataset `'earthquakes.csv'` is available
`here <../_static/example_data/earthquakes.csv>`_.
"""

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from mpl_toolkits.basemap import Basemap

from tick.hawkes import HawkesEM
from tick.plot import plot_hawkes_kernel_norms
from tick.plot.plot_utilities import get_plot_color

df = pd.read_csv('earthquakes.csv')

lats = df.Latitude.values
lons = df.Longitude.values

# put timestamps in minutes and remove minimum
timestamps = pd.to_datetime(df.DateTime).values.astype(np.float64)
timestamps *= 1e-9 / 60
timestamps -= min(timestamps)

fig, ax_list = plt.subplots(1, 2, figsize=(8, 3.4))
ax_list[0].set_title('Earthquakes near Japan')

# Draw map
lon0, lat0 = (139, 34)
lon1, lat1 = (147, 44)
m = Basemap(projection='merc', llcrnrlon=lon0, llcrnrlat=lat0, urcrnrlon=lon1,
urcrnrlat=lat1, resolution='l', ax=ax_list[0])
m.drawmapboundary(fill_color='#99ffff')
m.fillcontinents(color='#cc9966', lake_color='#99ffff')

# Number of splits in the map
n_groups_lats = 5

group_lats = []
group_lons = []
group_timestamps = []
group_names = []

delta_lats = (lats.max() - lats.min()) / n_groups_lats
delta_lons = lons.max() - lons.min()

for group in range(n_groups_lats):
# Split events into groups
min_lat_group = lats.min() + group * delta_lats
max_lat_group = lats.min() + (group + 1) * delta_lats

mask = (min_lat_group <= lats) & (lats < max_lat_group)

group_lats += [lats[mask]]
group_lons += [lons[mask]]
group_timestamps += [timestamps[mask]]

# Draw earthquakes on map
x, y = m(group_lons[group], group_lats[group])
m.scatter(x, y, 0.3, marker='o', color=get_plot_color(group))

# Draw zone labels on map
group_name = 'Z{}'.format(group)
group_names += [group_name]

center_lat = 0.5 * (min_lat_group + max_lat_group)
x_zone, y_zone = m(lons.max() + 0.1 * delta_lons, center_lat)

zone_bbox_style = dict(boxstyle="round", ec=(0, 0, 0, 0.9), fc=(1., 1, 1,
0.6))
ax_list[0].text(x_zone, y_zone, group_name, fontsize=12, fontweight='bold',
ha='left', va='center', color='k', withdash=True,
bbox=zone_bbox_style)

# Fit Hawkes process on events
events = []
for g in range(n_groups_lats):
events += [group_timestamps[g]]
events[g].sort()

kernel_discretization = np.linspace(0, 10000 / 60, 6)
em = HawkesEM(kernel_discretization=kernel_discretization, tol=1e-9,
max_iter=1000, print_every=200)
em.fit(events)

plot_hawkes_kernel_norms(em, ax=ax_list[1], node_names=group_names)
ax_list[1].set_xticklabels(ax_list[1].get_xticklabels(), fontsize=11)
ax_list[1].set_yticklabels(ax_list[1].get_yticklabels(), fontsize=11)

fig.tight_layout()
plt.show()
13 changes: 11 additions & 2 deletions tick/plot/plot_hawkes.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@


def plot_hawkes_kernel_norms(kernel_object, show=True, pcolor_kwargs=None,
node_names=None, rotate_x_labels=0.):
node_names=None, rotate_x_labels=0., ax=None):
"""Generic function to plot Hawkes kernel norms.
Parameters
Expand Down Expand Up @@ -37,6 +37,10 @@ def plot_hawkes_kernel_norms(kernel_object, show=True, pcolor_kwargs=None,
Number of degrees to rotate the x-labels clockwise, to prevent
overlapping.
ax : `matplotlib.axes`, default=None
If not None, the figure will be plot on this axis and show will be
set to False. Used only with matplotlib
Notes
-----
Kernels are displayed such that it shows norm of column influence's
Expand All @@ -54,7 +58,6 @@ def plot_hawkes_kernel_norms(kernel_object, show=True, pcolor_kwargs=None,
column_labels = ['$\\rightarrow {}$'.format(i) for i in node_names]

norms = kernel_object.get_kernel_norms()
fig, ax = plt.subplots()

if rotate_x_labels != 0.:
# we want clockwise rotation because x-axis is on top
Expand All @@ -63,6 +66,12 @@ def plot_hawkes_kernel_norms(kernel_object, show=True, pcolor_kwargs=None,
else:
x_label_alignment = 'center'

if ax is None:
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
else:
fig = ax.figure
show = False

if pcolor_kwargs is None:
pcolor_kwargs = {}

Expand Down
14 changes: 14 additions & 0 deletions tools/travis/docker_run.sh
Original file line number Diff line number Diff line change
Expand Up @@ -6,11 +6,25 @@ set -e -x
cp -R /io src
cd src

apt-get update;
apt-get install -y libgeos*;

eval "$(pyenv init -)"

pyenv global ${PYVER}
pyenv local ${PYVER}

PYMAJ=$(python -c "import sys; print(sys.version_info[0])")
PYMIN=$(python -c "import sys; print(sys.version_info[1])")

if (( PYMAJ == 3 )) && (( PYMIN == 7 )); then
echo "Skipping installing only needed for doctest which are not run on Python 3.7 (see bellow)"
else
# needed for basemap see https://github.com/matplotlib/basemap/issues/414#issuecomment-436792915
python -m pip install https://github.com/jswhit/pyproj/archive/v1.9.5.1rel.zip
python -m pip install https://github.com/matplotlib/basemap/archive/v1.1.0.tar.gz
fi

## disabled for the moment - re-enable later
#python -m pip install yapf --upgrade
#python -m yapf --style tools/code_style/yapf.conf -i tick examples --recursive
Expand Down
10 changes: 10 additions & 0 deletions tools/travis/osx_install.sh
Original file line number Diff line number Diff line change
Expand Up @@ -45,5 +45,15 @@ python -m pip install -r requirements.txt
python -m pip install sphinx pillow
python -m pip install cpplint
[[ "${PYVER}" != "3.7.0" ]] && python -m pip install tensorflow # does not yet exist on python 3.7

# Check if not already installed
if brew ls --versions geos > /dev/null; then
echo "geos is already installed"
else
brew install geos
fi
# needed for basemap see https://github.com/matplotlib/basemap/issues/414#issuecomment-436792915
# python -m pip install https://github.com/jswhit/pyproj/archive/v1.9.5.1rel.zip
python -m pip install https://github.com/matplotlib/basemap/archive/v1.2.0rel.tar.gz
pyenv rehash

5 changes: 3 additions & 2 deletions tools/travis/osx_run.sh
Original file line number Diff line number Diff line change
Expand Up @@ -25,9 +25,10 @@ set -e
python setup.py build_ext --inplace cpptest pytest

export PYTHONPATH=${PYTHONPATH}:`pwd`
for f in $(find examples -maxdepth 1 -type f -name "*.py"); do
cd examples
for f in $(find . -maxdepth 1 -type f -name "*.py"); do
FILE=$(basename $f)
FILE="${FILE%.*}" # skipping it takes too long
[[ "plot_asynchronous_stochastic_solver" != "$FILE" ]] && \
DISPLAY="-1" python -c "import tick.base; import examples.$FILE"
DISPLAY="-1" python -c "import tick.base; import $FILE"
done

0 comments on commit 4b8cbef

Please sign in to comment.