initial commit

README.md

+Example script for a synthetic example
+======================================
+
+This repo contains a few scripts that generate the synthetic example used in
+our paper [1] and thereafter uses the probability density function (PDF) time
+series to estimate the artifically imposed abrupt transitions in the data.
+
+In order to use this repo, it is advised that you have the 'uncertise' package
+in your Python path (it suffices to put the directory of the 'uncertise' [2]
+package in your current working directory, where the scripts are also located).
+Then, run the the following commands in the given order:
+
+    $ python pdfseries.py
+    $ python get_recurrencenetwork.py
+    $ python get_transitions.py
+    $ python plot_results.py
+
+At each step above, the scripts print messages ont he screen, telling you what
+is going on, and in the end, some kind of output file is saved to disk in the
+current working directory. In the last case, the output is a PNG figure,
+equivalent to the left column of Fig. 2 of the paper.
+
+[1] B. Goswami, N. Boers, A. Rheinwalt, N. Marwan, J. Heitzig, S. F. M.
+        Breitenbach, J. Kurths
+    Abrupt transitions in time series with uncertainties,
+    Nature Communications 9 (2018)
+    doi:10.1038/s41467-017-02456-6
+[2] https://gitlab.pik-potsdam.de/goswami/uncertise

get_pdfseries.py

+#! /usr/bin/env python
+"""
+Generate the time series for the synthetic example
+==================================================
+
+Reference
+=========
+
+B. Goswami, N. Boers, A. Rheinwalt, N. Marwan, J. Heitzig, S. F. M.
+Breitenbach, J. Kurths
+Abrupt transitions in time series with uncertainties,
+Nature Communications 9 (2018)
+doi:10.1038/s41467-017-02456-6
+
+"""
+# Created: Mon Feb 22, 2016  03:38PM
+# Last modified: Thu Jan 11, 2018  03:01PM
+# Copyright: Bedartha Goswami <goswami@pik-potsdam.de>
+
+
+import numpy as np
+from scipy.stats import norm
+from scipy.stats import gaussian_kde
+
+from uncertise.utils import _printmsg
+from uncertise.utils import _progressbar_start
+from uncertise.utils import _progressbar_update
+from uncertise.utils import _progressbar_finish
+
+from params import SF
+
+
+if __name__ == "__main__":
+    # set stdout parameters
+    VERBOSE = True
+    PROGBAR = True
+
+    # set parameters for analysis
+    _printmsg("Set up synthetic example...", VERBOSE)
+    T = 50                              # time period of sine curve
+    A = 1.                              # amplitude of sine curve
+    nt = 1000                           # total time range
+    tsplit = 675                        # time when the means split up
+    sample_res = SF                     # sampling frequency
+
+    # get basic dynamical variables
+    t = np.arange(nt)
+    x_ = A * np.sin((2 * np.pi / T) * t) 
+    I1 = np.r_[  # impulse causing jumps in dynamics
+               np.zeros(200),
+               5. * np.ones(200),
+               np.linspace(5., 0., 50),
+               np.zeros(550)
+               ]
+    # split to bimodality
+    x_ = x_ + I1
+    x1 = np.r_[
+                x_[:tsplit],
+                10.00 * x_[tsplit:]
+                ]
+    x2 = np.r_[
+                x_[:tsplit],
+                -10.00 * x_[tsplit:]
+                ]
+
+    # sample the time series with sampling noise
+    ts = t[::sample_res]              # sample every 'sample_res'th element
+    ns = len(ts)
+    es = 3.00 # dynamic noise
+    x1s = x1[::sample_res] + es * np.random.randn(ns)
+    x2s = x2[::sample_res] + es * np.random.randn(ns)
+
+    # create matrix to store series of PDFs
+    _printmsg("Creating probability distribution matrix...", VERBOSE)
+    lo1 = x1.min() - 5. * es
+    hi1 = x1.max() + 5. * es
+    lo2 = x2.min() - 5. * es
+    hi2 = x2.max() + 5. * es
+    lo = min(lo1, lo2)
+    hi = min(hi1, hi2)
+    var_res = 1000
+    var_span = np.linspace(lo, hi, var_res)
+    pdfmat = np.zeros((ns, var_res))
+    prog_bar = _progressbar_start(ns, PROGBAR)
+    ngrid = 1000
+    timeseries_grid = np.zeros((ngrid, nt))
+    ss = 1.50 # sampling noise
+    for i in range(ns):
+        pdf1 = norm(loc=x1s[i], scale=ss)
+        pdf2 = norm(loc=x2s[i], scale=ss)
+        samp = np.r_[
+                    pdf1.rvs(ngrid / 2),
+                    pdf2.rvs(ngrid / 2)
+                ]
+        timeseries_grid[:, i] = samp
+        kde = gaussian_kde(samp)
+        pdfmat[i] = kde(var_span)
+        _progressbar_update(prog_bar, i + 1)
+    _progressbar_finish(prog_bar)
+
+    # save output
+    FOUT = "pdfseries"
+    np.savez(FOUT,
+             pdfmat=pdfmat, ts=ts, var_span=var_span
+             )
+    _printmsg("saved to: %s.npz" % FOUT, VERBOSE)

get_recurrencenetwork.py

+#! /usr/bin/env python
+"""
+Synthetic example for abrupt transitions in time series with uncertainties
+==========================================================================
+
+Reference
+=========
+
+B. Goswami, N. Boers, A. Rheinwalt, N. Marwan, J. Heitzig, S. F. M.
+Breitenbach, J. Kurths
+Abrupt transitions in time series with uncertainties,
+Nature Communications 9 (2018)
+doi:10.1038/s41467-017-02456-6
+
+"""
+# Created: Mon Feb 22, 2016  03:38PM
+# Last modified: Thu Jan 11, 2018  02:24PM
+# Copyright: Bedartha Goswami <goswami@pik-potsdam.de>
+
+
+import numpy as np
+
+from uncertise.utils import _printmsg
+from uncertise import distributions
+from uncertise import networks
+
+from params import VERBOSE
+from params import PROGBAR
+from params import LD
+
+
+if __name__ == "__main__":
+    # load the data 
+    DAT = np.load("pdfseries.npz")
+    pdfmat = DAT["pdfmat"]
+    var_span = DAT["var_span"]
+
+    # get the series of CDF from the PDFs
+    cdfmat, iqr = distributions.cumulative(pdfmat, var_span, VERBOSE, PROGBAR)
+
+    # get the recurrence network as in igraph.Graph object
+    G = networks.prob_recurrence_network(cdfmat, var_span, LD, iqr, 
+                                         (1E-1, 5E-1), 5E-4,
+                                         VERBOSE, PROGBAR)
+    # note there will be a warning message if the required link density is not
+    # bracketed by the thresholds given by the function parameter 'thr_lims'.
+    # However, the code will not exit, but only return an empty graph G.
+    # Therefore, it is useful to have a check at this point
+    if G.density() == 0.:
+        import sys
+        print("Please increase the threshold limits!")
+        sys.exit()
+    
+    # save output
+    FOUT = "./recurrencenetwork.igraph"
+    G.write_pickle(FOUT)
+    _printmsg("saved to: %s" % FOUT, VERBOSE)

get_transitions.py

+#! /usr/bin/env python
+"""
+Synthetic example for abrupt transitions in time series with uncertainties
+==========================================================================
+
+Reference
+=========
+
+B. Goswami, N. Boers, A. Rheinwalt, N. Marwan, J. Heitzig, S. F. M.
+Breitenbach, J. Kurths
+Abrupt transitions in time series with uncertainties,
+Nature Communications 9 (2018)
+doi:10.1038/s41467-017-02456-6
+
+"""
+# Created: Mon Feb 22, 2016  03:38PM
+# Last modified: Thu Jan 11, 2018  02:27PM
+# Copyright: Bedartha Goswami <goswami@pik-potsdam.de>
+
+
+import numpy as np
+import igraph as ig
+
+from uncertise.utils import _printmsg
+from uncertise import events
+
+from params import VERBOSE
+from params import PROGBAR
+from params import WS
+from params import ST
+from params import NS
+from params import ALPHA
+
+if __name__ == "__main__":
+    _printmsg("load data ...", VERBOSE)
+    # load the recurrence network stored as an igraph.Graph object
+    FNAME = "recurrencenetwork.igraph"
+    G = ig.Graph.Read_Pickle(FNAME)
+
+    # load time scale data
+    FNAME = "pdfseries.npz"
+    D = np.load(FNAME)
+    t = D["ts"]
+
+    # get the community strength values for the observed data as well as the
+    # degree configuration model in order to identify non-random transitions
+    Q, tmid = events.community_strength_data(G, t, WS, ST, VERBOSE, PROGBAR)
+    res = events.community_strength_random_model(G, t, WS, ST, NS,
+                                                 VERBOSE, PROGBAR)
+    Qsurr = res[0]
+    tmid_surr = res[1]
+
+    # get p-values 
+    _printmsg("get p-values ...", VERBOSE)
+    pvals = events.pvalues(Q, Qsurr, len(tmid), VERBOSE, PROGBAR)
+
+    # get statistically significant windows using Holm's method
+    sig_idx = events.holm(pvals, alpha=ALPHA, corr_type="dunn")
+
+    # save output
+    _printmsg("saving output ...", VERBOSE)
+    FOUT = "pvalues"
+    np.savez(FOUT,
+             Q=Q, tmid=tmid, Qsurr=Qsurr, tmid_surr=tmid_surr,
+             pvals=pvals, sig_idx=sig_idx
+             )
+    _printmsg("saved to: %s.npz" % FOUT, VERBOSE)

params/__init__.py

+#! /usr/bin/env python
+"""
+Contains the various fixed parameters to use for the analysis
+=============================================================
+
+
+"""
+# Created: Mon Feb 22, 2016  03:38PM
+# Last modified: Thu Jan 11, 2018  03:00PM
+# Copyright: Bedartha Goswami <goswami@pik-potsdam.de>
+
+# initialize parameters for stdout
+VERBOSE = True
+PROGBAR = True
+
+# initalize parameters for the analysis
+SF = 2              # sampling frequency
+LD = 0.30           # link density
+WS = 100            # window size for transition analysis
+ST = 10             # step size by which the window is moved
+NS = 1000           # number of suggogate networks to generate
+ALPHA = 0.01        # significance level of statistical test
+
+
+
+

plot_results.py

+#! /usr/bin/env python
+"""
+Fig. 2: Synthetic example and schematic of the analysis.
+========================================================
+
+Note: This script plots only the left column of the Fig. 2 from the paper!
+
+Reference
+=========
+
+B. Goswami, N. Boers, A. Rheinwalt, N. Marwan, J. Heitzig, S. F. M.
+Breitenbach, J. Kurths
+Abrupt transitions in time series with uncertainties,
+Nature Communications 9 (2018)
+doi:10.1038/s41467-017-02456-6
+
+"""
+# Created: Mon Jun 06, 2016  09:56PM
+# Last modified: Thu Jan 11, 2018  03:02PM
+# Copyright: Bedartha Goswami <goswami@pik-potsdam.de>
+
+import sys
+import numpy as np
+import igraph as ig
+import datetime as dt
+
+import matplotlib.pyplot as pl
+import matplotlib.dates as mdates
+from matplotlib.collections import PolyCollection
+from matplotlib.colors import Normalize
+from matplotlib.colorbar import ColorbarBase
+from matplotlib.ticker import MultipleLocator, FixedLocator
+
+from uncertise.utils import _printmsg
+
+from params import VERBOSE
+from params import WS
+from params import SF as RES
+
+# set matplotlib rcParams for Sans Serif Latex math font
+pl.rcParams["mathtext.fontset"] = "custom"
+pl.rcParams["text.usetex"] = True
+pl.rcParams["text.latex.preamble"] = r"""
+                                         \usepackage{helvet}
+                                         \renewcommand{\familydefault}{\sfdefault}
+                                         \usepackage{sfmath}
+                                     """
+
+
+if __name__ == "__main__":
+    _printmsg("Loading data...", True)
+    
+    # load the synthetic dataset
+    f = "pdfseries.npz"
+    d = np.load(f)
+    t, p, v = d["ts"], d["pdfmat"], d["var_span"]
+    # res = d["sample_res"]
+    tlo, thi = t.min(), t.max()
+    vlo, vhi = v.min(), v.max()
+
+    # load the probability of recurrence network using igraph
+    f = "recurrencenetwork.igraph"
+    G = ig.Graph.Read_Pickle(f)
+
+    # transform igraph.Graph to numpy 2d array
+    print("convert to numpy 2D array ...")
+    n = len(G.vs)
+    P = np.zeros((n, n), dtype="float")
+    for edge in G.es:
+        i, j = edge.source, edge.target
+        P[i, j] = edge["weight"]
+    P = P + P.T
+
+    # load the p-value data for event analysis
+    f = "pvalues.npz"
+    d = np.load(f)
+    te, pv = d["tmid"], d["pvals"]
+    idx_sig = d["sig_idx"]
+
+    _printmsg("Setting up figure...", True)
+    # set up figure
+    fig = pl.figure(figsize=[6.28, 10.24])
+    lm1, lm2 = 0.130, 0.800
+    bm1, bm2, bm3 = 0.715, 0.300, 0.090
+    wd1, wd2 = 0.65, 0.04
+    ht1, ht2 = 0.195, 0.400
+    axlabfs, tklabfs = 12, 12
+
+    # set up axes
+    # for analysis of PDFs
+    ax_a = fig.add_axes([lm1, bm1, wd1, ht1])
+    ax_b = ax_a.twinx()
+    ax_b.set_position([lm1, bm2, wd1, ht2])
+    ax_c = fig.add_axes([lm1, bm3, wd1, ht1])
+    cx_a = fig.add_axes([lm2, bm1, wd2, ht1])
+    cx_b = fig.add_axes([lm2, bm2, wd2, ht2])
+
+    _printmsg("Plotting...", True)
+    # a. sequence of PDFs of the synthetic data
+    xx, yy = np.meshgrid(t, v)
+    zz = p.T
+    im_a = ax_a.pcolormesh(xx, yy, zz, cmap="viridis_r")
+
+    # b. probability of recurrence matrix estimated from the PDFs
+    im_b = ax_b.imshow(P,
+                       origin="lower", cmap="gist_heat_r",
+                       extent=(tlo, thi, tlo, thi),
+                       interpolation="nearest")
+
+    # c. p-values of the event detection / community-strength analysis
+    ax_c.fill_between(te, 1., pv,
+                      edgecolor="none", facecolor="SkyBlue",
+                      )
+    ax_c.plot(te[idx_sig], 1E-2 * np.ones(len(idx_sig)), "+",
+              mec="Maroon", mfc="Maroon", mew=1.50, ms=8,
+              zorder=5,)
+
+    # set up colorbars
+    im_a.set_clim(0., 0.3)
+    cb_a = pl.colorbar(im_a, cax=cx_a, extend="max")
+    im_b.set_clim(0., 0.4)
+    cb_b = pl.colorbar(im_b, cax=cx_b, extend="max")
+
+    # set/format axes ticks
+    for ax in [ax_a]:
+        ax.tick_params(which="both",
+                       labelbottom="off", labeltop="on",
+                       labelsize=tklabfs,
+                       )
+        YminorLocator = FixedLocator(np.arange(-20, 21, 2.5))
+        ax.yaxis.set_minor_locator(YminorLocator)
+        YmajorLocator = FixedLocator(np.arange(-20, 21, 10))
+        ax.yaxis.set_major_locator(YmajorLocator)
+    XminorLocator = MultipleLocator(50)
+    XmajorLocator = MultipleLocator(200)
+    for ax in [ax_a, ax_b, ax_c]:
+        ax.tick_params(which="major", axis="both", size=6)
+        ax.tick_params(which="minor", axis="both", size=4)
+        ax.xaxis.set_minor_locator(XminorLocator)
+        ax.xaxis.set_major_locator(XmajorLocator)
+        if ax in [ax_b]:
+            ax.yaxis.set_minor_locator(XminorLocator)
+            ax.yaxis.set_major_locator(XmajorLocator)
+    for ax in [ax_a, ax_b, ax_c]:
+        ax.tick_params(which="both",
+                       left="on", right="off",
+                       labelleft="on", labelright="off",
+                       labelsize=tklabfs,
+                       )
+    # colorbar axes ticks
+    cb_a.set_ticks(np.arange(0., 0.31, 0.10))
+    cb_b.set_ticks(np.arange(0., 0.41, 0.10))
+
+    # set axis limits
+    for ax in [ax_a]:
+        ax.set_xlim(0., 1001)
+        ax.set_ylim(-20., 20.)
+    ax_c.set_xlim(ax_a.get_xlim())
+    ax_c.set_yscale("log", nonposy='clip')
+    ax_c.set_ylim(1E-4, 1E0)
+
+    # set_axis labels
+    for ax in [ax_a]:
+        ax.set_xlabel("Time",
+                      fontsize=axlabfs)
+        ax.xaxis.set_label_position("top")
+        ax.set_ylabel("Signal",
+                      fontsize=axlabfs)
+    for ax in [ax_b]:
+        ax.set_ylabel("Time",
+                      fontsize=axlabfs)
+    ax_b.yaxis.set_label_position("left")
+    for ax in [ax_c]:
+        ax.set_xlabel("Time",
+                      fontsize=axlabfs)
+        ax.set_ylabel(r"$p$-value",
+                      fontsize=axlabfs)
+    cb_a.set_label("Prob. density",
+                   fontsize=axlabfs)
+    cb_b.set_label("Prob. of recurrence",
+                   fontsize=axlabfs)
+
+    # hatched rectangles in subplots c and f where there are no windows
+    for ax in [ax_c]:
+        ylo, yhi = ax.get_ylim()
+        xlim = ax.get_xlim()
+        time = te
+        nh, nv = 5, 5E-2
+        hlw = 0.5
+        # ###### left-hand side rectangle
+        xlo = xlim[0]
+        xhi = time[0]
+        N = int( (xhi - xlo) / nh )
+        x = np.linspace(xlo, xhi, N)
+        for i in range(N):
+            ax.vlines(x[i], ylo, yhi,
+                      linestyles="-", linewidths=hlw, colors="0.15")
+        N = int( (yhi - ylo) / nv )
+        y = np.logspace(np.log10(ylo), np.log10(yhi), N)
+        for i in range(N):
+            ax.hlines(y[i], xlo, xhi,
+                      linestyles="-", linewidths=hlw, colors="0.15")
+        # ###### right-hand side rectangle
+        xlo = time[-1]
+        xhi = xlim[1]
+        N = int( (xhi - xlo) / nh )
+        x = np.linspace(xlo, xhi, N)
+        for i in range(N):
+            ax.vlines(x[i], ylo, yhi,
+                      linestyles="-", linewidths=hlw, colors="0.15")
+        N = int( (yhi - ylo) / nv )
+        y = np.logspace(np.log10(ylo), np.log10(yhi), N)
+        for i in range(N):
+            ax.hlines(y[i], xlo, xhi,
+                      linestyles="-", linewidths=hlw, colors="0.15")
+
+
+    # add legends for the bars showing signifcant event windows
+    rect = pl.Rectangle
+    ax_c.plot(1080, 2E-1, "+",
+              mec="Maroon", mfc="Maroon", mew=1.50, ms=8,
+              zorder=5, clip_on=False)
+    ax_c.text(1060, 2E-4,
+              r"Sig. at $\alpha$ =0.05",
+              ha="left", va="bottom", rotation=90,
+              fontsize=tklabfs - 1, zorder=8, clip_on=False)
+    rbbox = rect(xy=(1045, 1E-4),
+                 width=70, height=5E-1,
+                 color="k", fill=None, linewidth=1., clip_on=False,
+                 alpha=1., zorder=15)
+    ax_c.add_artist(rbbox)
+
+    # add "no data" text to hatched portion of subplots c and f.
+    for ax in [ax_c]:
+        ax.text(x=50, y=1E-2,
+                s="No data",
+                ha="center", va="center",
+                fontsize=tklabfs, rotation=90, zorder=10,
+                bbox={'facecolor': 'w',
+                      'edgecolor': 'none',
+                      'alpha': 1.0,
+                      'pad': 1}
+                )
+        ax.text(x=950., y=1E-2,
+                s="No data",
+                ha="center", va="center",
+                fontsize=tklabfs, rotation=90, zorder=10,
+                bbox={'facecolor': 'w',
+                      'edgecolor': 'none',
+                      'alpha': 1.0,
+                      'pad': 1}
+                )
+
+    # show top and bottom ticks
+    for ax in [ax_a, ax_b, ax_c]:
+        ax.tick_params(which="both", top="on", bottom="on")
+    for ax in [ax_a, ax_c]:
+        ax.grid(which="minor", axis="x")
+
+    # add rectangle in subplot b to show sliding window
+    t0 = 150
+    rslid = rect(xy=(t0, t0),
+                 width=WS * RES, height=WS * RES,
+                 color="w", linewidth=2.5, fill=None, zorder=10)
+    ax_b.add_artist(rslid)
+
+    # add arrows to denote the movement of the sliding window
+    tof = 100
+    t01 = t0 + WS * RES
+    ax_b.annotate("",
+                  xy=(t01 + tof, t01 + tof),
+                  xytext=(t01 + tof / 4., t01 + tof / 4.),
+                  xycoords="data",
+                  arrowprops={
+                      "facecolor": "w",
+                      "edgecolor": "none",
+                      "shrink": 0.05,
+                  },
+                  zorder=15)
+
+    # annotate the three shifts at T = 200, 400-450, and 675.
+    for ax in [ax_a]:
+        if ax == ax_a:
+            YY = -17.
+        else:
+            YY = -15.
+        ax.text(200., YY,
+                "(1)",
+                fontsize=tklabfs, ha="center", va="center",
+                bbox={'facecolor': 'w',
+                      'edgecolor': 'none',
+                      'alpha': 0.0,
+                      'pad': 2},
+                )
+        ax.text(425., YY,
+                "(2)",
+                fontsize=tklabfs, ha="center", va="center",
+                bbox={'facecolor': 'w',
+                      'edgecolor': 'none',
+                      'alpha': 0.0,
+                      'pad': 2},
+                )
+        ax.text(675., YY,
+                "(3)",
+                fontsize=tklabfs, ha="center", va="center",
+                bbox={'facecolor': 'w',
+                      'edgecolor': 'none',
+                      'alpha': 0.0,
+                      'pad': 2},
+                )
+
+    # save figure
+    _printmsg("save figure ...", VERBOSE)
+    FOUT = "results.png"
+    pl.savefig(FOUT)
+    _printmsg("saved to: %s" % FOUT, VERBOSE)
+
+