""" 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