"""
A suite of functions that assist in carrying out the recurrence analysis
========================================================================
!! Please import as a module and use it as such -- not as a script !!
"""
# Created: Thu May 26, 2016 03:50PM
# Last modified: Tue Aug 15, 2017 01:21pm
# Copyright: Bedartha Goswami
import numpy as np
import igraph as ig
from utils import _printmsg
from utils import _progressbar_start
from utils import _progressbar_update
from utils import _progressbar_finish
def prob_recurrence_network(dist_list, var_span, ld, iqr,
thr_lims=(1E-3, 5E-1), tol=5E-4,
verbose=False, pbar=False):
"""
Return probability of recurrence network for given link density.
The returned recurrence network has real-valued edge weights from the
closed interval [0, 1] corresponding to a priori estimated
probabilities of recurrence between each pair of time points with the
help of upper and lower bounds on the difference distributions of those
two time points. Note that `time` here refers to the time-ordering of
the distribution list ``dist_list``.
Parameters
----------
dist_list : list
List of time-ordered cumulative probability distributions. Each item
of this list is of the same size as ``var_span`` and has float
entries corresponding to the cumulative probability evaluated at the
corresponding value in ``var_span``.
var_span : numpy 1D array
Contains the range of proxy values on which the cumulative
probabilities provided in each item of ``dist_list`` are evaluated.
For best results, it is advised to ensure that ``var_span`` should
cover a range of values wide enough so that the first and last
entries of all the CDFs are 0 and 1 respectively.
ld : number, float
Desired link density of the final recurrence network.
iqr : number, float
Inter-quartile range of the total probability distribution function
(PDF) of the variable when summed over all time instances.
thr_lims : tuple, float, optional
Tuple containing the lower and upper limit of recurrence thresholds
(as a fraction of ``iqr``) within which the bisection routine should
search for the appropriate network with desired link density. Note
that ``thr_lims = (thr_lo, thr_hi)``.
tol : number, float, optional
Tolerance level for the bisection routine. This number specifies the
error margin up to which the returned network's link density has to
be equal to the desired link density specified with ``ld``.
verbose: bool, optional
Specifies verbosity of the routine.
pbar: bool, optional
Specifies the display of a progress bar in terminal on execution.
Returns
-------
G : igraph.Graph
igraph.Graph object with edge weights corresponding to the estimated
probabilities of recurrence between each pair of nodes (i.e., time
instances of observation). The edge weights can be accessed by
``G.es["weight"]``.
See Also
--------
prob_recurrence_matrix
link_density
Notes
-----
The probability of recurrence matrix estimated here is obtained using
the function ``prob_recurrence_matrix`` with the only difference that
the diagonals are set to zero to avoid self-loops in the network. Also,
the method based on the upper and lower bounds on difference
distributions as outlined in Goswami (2015) does not allow us to specify
the link density beforehand. Thus, a bisection routine between
pre-specified upper and lower recurrence thresholds ``thr_lim`` is used.
This incurs an additional parameter, the tolerance level ``tol``, i.e.,
the error between the link density of the returned network and the
desired link density.
References
----------
.. [1] Goswami, Bedartha. "Uncertainties in climate data analysis".
Chapter 6. PhD Dissertation (2015).
https://publishup.uni-potsdam.de/frontdoor/index/index/docId/7831
.. [2] Williamson, Robert C., Downs, Tom. "Probabilistic arithmetic I.
Numerical methods for calculating convolutions and dependency bounds".
Int. J. Approx. Reasoning (1990). 4:89-158.
"""
_printmsg("Estimating network with link density = %.3f" % ld, verbose)
n = len(dist_list)
valid_lims, lds = _precnet_check_limits(dist_list, var_span, ld, iqr,
thr_lims, verbose, pbar)
if not valid_lims:
P = np.zeros((n, n))
G = _precmat_to_igraph(P)
else:
ld_lims = lds
G, _, _ = _precnet_bisection(dist_list, var_span, ld, ld_lims, iqr,
thr_lims, tol, verbose, pbar)
return G
def _precnet_check_limits(dist_list, var_span, ld, iqr,
thr_lims=(1E-3, 5E-1), verbose=False, pbar=False):
"""Checks if desired link density is within given threshold limits."""
_printmsg("Checking if threshold limits trap target density...",
verbose)
ld_lims = []
for thr in thr_lims:
e = thr * iqr
_printmsg("\tTHRESHOLD = %.2E OF INTER-QUARTILE RANGE"
% thr, verbose)
G = _precnet_igraph(dist_list, var_span, e, verbose, pbar)
ld_curr = link_density(G)
_printmsg("\tESTIMATED LINK DENSITY = %.4f" % ld_curr, verbose)
ld_lims.append(ld_curr)
cond = (ld >= ld_lims[0]) * (ld <= ld_lims[1])
if not cond:
str0 = "Target link density could not be bracketed!"
str1 = "Increase initial THR bracket or change target LD."
str2 = "LD = %.3f for THR = %.2E and %.3f for THR = %.2f" \
% (ld_lims[0], thr_lims[0], ld_lims[1], thr_lims[1])
print(str0 + "\n" + str1 + "\n" + str2)
return cond, ld_lims
def _precnet_bisection(dist_list, var_span, ld, ld_lims, iqr,
thr_lims=(1E-3, 5E-1), tol=5E-4,
verbose=False, pbar=False):
"""Bisection routine to get network of desired link density."""
_printmsg("Entering bisection routine...", verbose)
ld_lo, ld_hi = ld_lims[0], ld_lims[1]
thr_lo, thr_hi = thr_lims[0], thr_lims[1]
err = np.abs(ld_hi - ld_lo)
cond = err > tol
while cond:
thr_curr = 0.5 * (thr_lo + thr_hi)
_printmsg("Current values...", verbose)
_printmsg("\t...Thresholds = %.5f" % thr_curr, verbose)
_printmsg("\t...Thresholds (Lo, Hi) = (%.5f, %.5f)"
% (thr_lo, thr_hi), verbose)
_printmsg("\t...LDs (Lo, Hi) = (%.5f, %.5f)"
% (ld_lo, ld_hi), verbose)
_printmsg("\t...Error = %G" % err, verbose)
e = thr_curr * iqr
G = _precnet_igraph(dist_list, var_span, e, verbose, pbar)
ld_curr = link_density(G)
if ld_curr <= ld:
thr_lo = thr_curr
ld_lo = link_density(G)
else:
thr_hi = thr_curr
ld_hi = link_density(G)
err = np.abs(ld_hi - ld_lo)
cond = err > tol
return G, thr_curr, ld_curr
def _precnet_igraph(dist_list, var_span, e, verbose=False, pbar=False):
"""
Returns weighted igraph object by estimating prob. of rec. mat.
"""
P = prob_recurrence_matrix(dist_list, var_span, e, verbose, pbar)
np.fill_diagonal(P, 0.) # remove self-loops
G = _precmat_to_igraph(P)
return G
def _precmat_to_igraph(P):
"""
Returns igraph weighted graph for given prob. of. rec. matrix.
"""
G = ig.Graph.Weighted_Adjacency(P.tolist(), mode=ig.ADJ_UNDIRECTED)
return G
def link_density(G):
"""
Returns the link density of an igraph Graph object.
Since Python igraph does not return the correct link density for a
weighted graph taking into account the edge weights, this wrapper for
the igraph.Graph.density function is defined. It uses the
igraph.Graph.density function in case the graph is unweighted but if
not, it estimates the link density as the sum over all edge weights
divided by ``n * (n - 1) / 2``where ``n`` is the number of nodes in the
graph.
Parameters
----------
G : igraph Graph
An igraph Graph object that can be weighted or unweighted, but is
not directed, for which the link density is to be estimated.
Returns
-------
rho : scalar, float
The estimated link density for the un/weighted undirected graph G.
See Also
--------
igraph.Graph.density
"""
if G.is_weighted():
n = G.vcount()
sum_e = np.sum(G.es["weight"])
rho = 2. * sum_e / float(n * (n - 1))
else:
rho = G.density()
return rho
def prob_recurrence_matrix(dist_list, var_span, e=0.1,
verbose=False, pbar=False):
"""
Return probability of recurrence matrix for given recurrence threshold.
The estimated probability of recurrence matrix contains the a priori
probabilities of recurrence between each pair of time points with the
help of upper and lower bounds on the difference distributions of those
two time points. Note that `time` here refers to the time-ordering of
the distribution list ``dist_list``.
Parameters
----------
dist_list : list
List of time-ordered cumulative probability distributions. Each item
of this list is of the same size as ``var_span`` and has float
entries corresponding to the cumulative probability evaluated at the
corresponding value in ``var_span``.
var_span : numpy 1D array
Contains the range of proxy values on which the cumulative
probabilities provided in each item of ``dist_list`` are evaluated.
For best results, it is advised to ensure that ``var_span`` should
cover a range of values wide enough so that the first and last
entries of all the CDFs are 0 and 1 respectively.
e : number, float
Recurrence threshold for which to estimate probability of
recurrence. This is expressed as a fraction of the inter-quartile
range (see ``iqr``) and is thus a float between 0 and 1.
iqr : number, float
Inter-quartile range of the total probability distribution function
(PDF) of the variable when summed over all time instances.
verbose: bool, optional
Specifies verbosity of the routine.
pbar: bool, optional
Specifies the display of a progress bar in terminal on execution.
Returns
-------
P_rec : numpy 2D array
The probability of recurrence matrix containing float values falling
in the interval [0, 1]. The shape of this array is ``(n, n)`` where
``n``is the length of the dataset, i.e., ``n = len(distlist)``.
See Also
--------
prob_recurrence_network
Notes
-----
The probability of recurrence for each pair of time instances is
estimated by first estimating the upper and lower bounds on the
difference distribution (cf. Williamson & Down, 1990) and thereafter
using those bounds to derive upper and lower bounds on the probability
of recurrence itself. The expected probability of recurrence is then
proved to be the mean of the bounds obtained thus (cf. Goswami (2015)).
References
----------
.. [1] Goswami, Bedartha. "Uncertainties in climate data analysis".
Chapter 6. PhD Dissertation (2015).
https://publishup.uni-potsdam.de/frontdoor/index/index/docId/7831
.. [2] Williamson, Robert C., Downs, Tom. "Probabilistic arithmetic I.
Numerical methods for calculating convolutions and dependency bounds".
Int. J. Approx. Reasoning (1990). 4:89-158.
"""
u = var_span
n = len(dist_list)
P_rec = np.zeros((n, n))
_printmsg("\tPairwise recurrence probability bounds...", verbose)
f1, f2 = np.zeros((n, len(u))), np.zeros((n, len(u)))
for j in range(n):
f1[j, :] = np.interp(u - e, u, dist_list[j]) # for z = e
f2[j, :] = np.interp(u + e, u, dist_list[j]) # for z = -e
prog_bar = _progressbar_start(n, pbar)
for i in range(n):
fi = np.interp(u, u, dist_list[i]) # Xi.cdf(u)
plus_lo, plus_hi = _bounds_williamson(fi, f1)
mnus_lo, mnus_hi = _bounds_williamson(fi, f2)
diff = plus_lo - mnus_hi
Pij_l = np.fmax(diff, np.zeros(diff.shape))
diff = plus_hi - mnus_lo
Pij_u = np.fmin(diff, np.ones(diff.shape))
P_rec[i, :] = (Pij_l + Pij_u) / 2.
P_rec[i, i] = 1.
_progressbar_update(prog_bar, i)
_progressbar_finish(prog_bar)
return P_rec
def _bounds_williamson(fi, g):
"""
Returns bounds on the probability of diff. distribution.
"""
F = fi - g
m, M = np.max(F, axis=1), 1. + np.min(F, axis=1)
m[m < 0.] = 0.
M[M > 1.] = 1.
return m, M