Module pymake.util.math

Source code
# -*- coding: utf-8 -*-

import numpy as np
from numpy import ma
import scipy as sp
import networkx as nx
from .utils import nxG

import logging
lgg = logging.getLogger('root')

##########################
### Stochastic Process
##########################

def lognormalize(x):
    return np.exp(x - np.logaddexp.reduce(x))

def categorical(params):
    return np.where(np.random.multinomial(1, params) == 1)[0]

def bernoulli(param, size=1):
    return np.random.binomial(1, param, size=size)

### Power law distribution generator
def random_powerlaw(alpha, x_min, size=1):
    ### Discrete
    alpha = float(alpha)
    u = np.random.random(size)
    x = (x_min-0.5)*(1-u)**(-1/(alpha-1))+0.5
    return np.floor(x)

### A stick breakink process, truncated at K components.
def gem(gmma, K):
    sb = np.empty(K)
    cut = np.random.beta(1, gmma, size=K)
    for k in range(K):
        sb[k] = cut[k] * cut[0:k].prod()
    return sb

##########################
### Means and Norms
##########################

### Weighted  means
def wmean(a, w, mean='geometric'):
    if mean == 'geometric':
        kernel = lambda x : np.log(x)
        out = lambda x : np.exp(x)
    elif mean == 'arithmetic':
        kernel = lambda x : x
        out = lambda x : x
    elif mean == 'harmonic':
        num = np.sum(w)
        denom = np.sum(np.asarray(w) / np.asarray(a))
        return num / denom
    else:
        raise NotImplementedError('Mean Unknwow: %s' % mean)

    num = np.sum(np.asarray(w) * kernel(np.asarray(a)))
    denom = np.sum(np.asarray(w))
    return out(num / denom)

##########################
### Matrix/Image Operation
##########################
from scipy import ndimage

def draw_square(mat, value, topleft, l, L, w=0):
    tl = topleft

    # Vertical draw
    mat[tl[0]:tl[0]+l, tl[1]:tl[1]+w] = value
    mat[tl[0]:tl[0]+l, tl[1]+L-w:tl[1]+L] = value
    # Horizontal draw
    mat[tl[0]:tl[0]+w, tl[1]:tl[1]+L] = value
    mat[tl[0]+l-w:tl[0]+l, tl[1]:tl[1]+L] = value
    return mat

def dilate(y, size=1):
    dim = y.ndim
    mask = ndimage.generate_binary_structure(dim, dim)
    if size > 1:
        for i in range(1, size):
            mask = np.vstack((mask, mask[-1,:]))
            mask = np.column_stack((mask, mask[:, -1]))
    y_f = ndimage.binary_dilation(y, structure=mask).astype(y.dtype)
    return y_f


##########################
### Array routine  Operation
##########################
from collections import Counter

def sorted_perm(a, label=None, reverse=False):
    """ return sorted $a and the induced permutation.
        Inplace operation """
    # np.asarray applied this tuple lead to error, if label is string
    # because a should be used as elementwise comparison
    if label is None:
        label = np.arange(a.shape[0])
    hist, label = zip(*sorted(zip(a, label), reverse=reverse))
    hist = np.asarray(hist)
    label = np.asarray(label)
    return hist, label

def degree_hist_to_list(d, dc):
    degree = np.repeat(np.round(d).astype(int), np.round(dc).astype(int))
    return degree


def clusters_hist(clusters, labels=None, remove_empty=True):
    """ return non empty clusters histogramm sorted.

        parameters
        ---------
        clusters: np.array
            array of clusters membership of data.

        returns
        -------
        hist: np.array
            count of element by clusters (decrasing hist)
        label: np.array
            label of the cluster aligned with hist
    """
    block_hist = np.bincount(clusters)
    if labels is None:
        labels = range(len(block_hist))

    hist, labels = sorted_perm(block_hist, labels, reverse=True)

    if remove_empty is True:
        null_classes = (hist == 0).sum()
        if null_classes > 0:
            hist = hist[:-null_classes]; labels = labels[:-null_classes]

    return hist, labels

def adj_to_degree(y):
    # @debug: dont' call nxG or do a native integration !

    # To convert normalized degrees to raw degrees
    #ba_c = {k:int(v*(len(ba_g)-1)) for k,v in ba_c.iteritems()}
    G = nxG(y)
    #degree = sorted(dict(nx.degree(G)).values(), reverse=True)

    #ba_c = nx.degree_centrality(G)
    return dict(nx.degree(G))

def degree_hist(_degree, filter_zeros=False):
    if isinstance(_degree, np.ndarray) and _degree.ndim == 2 :
        degree = list(dict(adj_to_degree(_degree)).values())
    elif isinstance(_degree, (list, np.ndarray)):
        degree = _degree
    else:
        # networkx
        degree = list(dict(_degree).values())

    max_c = np.max(degree)
    d = np.arange(max_c+1)
    dc = np.bincount(degree, minlength=max_c+1)

    if len(d) == 0:
        return [], []

    if dc[0] > 0:
        lgg.debug('%d unconnected vertex' % dc[0])
    d = d[1:]
    dc = dc[1:]


    if filter_zeros is True:
        #d, dc = zip(*filter(lambda x:x[1] != 0, zip(d, dc)))
        nzv = (dc != 0)
        d = d[nzv]
        dc = dc[nzv]

    return d, dc

def random_degree(Y, params=None):
    _X = []
    _Y = []
    N = Y[0].shape[0]
    nb_uniq_degree = []
    dc_list = []
    for y in Y:
        ba_c = adj_to_degree(y)
        d, dc = degree_hist(ba_c)
        nb_uniq_degree.append(len(dc))
        dc_list.append(dc)
    dc_mat = ma.array(np.empty((N, max(nb_uniq_degree))), mask=True)
    for i, degrees in enumerate(dc_list):
        size = nb_uniq_degree[i]
        dc_mat[i, :size] = degrees

    y = dc_mat.mean(0)
    yerr = dc_mat.std(0)
    # 0 are filtered out in degree_hist
    return np.arange(1, len(y)+1),  np.round(y), yerr

def reorder_mat(y, clusters, labels=False, reverse=True):
    """Reorder the  matrix according the clusters membership
        @Debug: square matrix
    """
    assert(y.shape[0] == y.shape[1] == len(clusters))
    if reverse is True:
        hist, label = clusters_hist(clusters)
        sorted_clusters = np.empty_like(clusters)
        for i, k in enumerate(label):
            if i != k:
                sorted_clusters[clusters == k] = i
    else:
        sorted_clusters = clusters

    N = y.shape[0]
    nodelist = [k[0] for k in sorted(zip(range(N), sorted_clusters),
                                     key=lambda k: k[1])]

    y_r = y[nodelist, :][:, nodelist]
    if labels is True:
        return y_r, nodelist
    else:
        return y_r

def shiftpos(arr, fr, to, axis=0):
    """ Move element In-Place, shifting backward (or forward) others """
    if fr == to: return
    x = arr.T if axis == 1 else arr
    tmp = x[fr].copy()
    if fr > to:
        x[to+1:fr+1] = x[to:fr]
    else:
        x[fr:to] = x[fr+1:to+1]
    x[to] = tmp


##########################
### Colors Operation
##########################
import math

def floatRgb(mag, cmin, cmax):
        """ Return a tuple of floats between 0 and 1 for the red, green and
                blue amplitudes.
        """

        try:
                # normalize to [0,1]
                x = float(mag-cmin)/float(cmax-cmin)
        except:
                # cmax = cmin
                x = 0.5
        blue = min((max((4*(0.75-x), 0.)), 1.))
        red  = min((max((4*(x-0.25), 0.)), 1.))
        green= min((max((4*math.fabs(x-0.5)-1., 0.)), 1.))
        return (red, green, blue)

def strRgb(mag, cmin, cmax):
   """ Return a tuple of strings to be used in Tk plots.
   """

   red, green, blue = floatRgb(mag, cmin, cmax)
   return "#%02x%02x%02x" % (red*255, green*255, blue*255)

def rgb(mag, cmin, cmax):
   """ Return a tuple of integers to be used in AWT/Java plots.
   """

   red, green, blue = floatRgb(mag, cmin, cmax)
   return (int(red*255), int(green*255), int(blue*255))

def htmlRgb(mag, cmin, cmax):
   """ Return a tuple of strings to be used in HTML documents.
   """
   return "#%02x%02x%02x"%rgb(mag, cmin, cmax)

Functions

def adj_to_degree(y)
Source code
def adj_to_degree(y):
    # @debug: dont' call nxG or do a native integration !

    # To convert normalized degrees to raw degrees
    #ba_c = {k:int(v*(len(ba_g)-1)) for k,v in ba_c.iteritems()}
    G = nxG(y)
    #degree = sorted(dict(nx.degree(G)).values(), reverse=True)

    #ba_c = nx.degree_centrality(G)
    return dict(nx.degree(G))
def bernoulli(param, size=1)
Source code
def bernoulli(param, size=1):
    return np.random.binomial(1, param, size=size)
def categorical(params)
Source code
def categorical(params):
    return np.where(np.random.multinomial(1, params) == 1)[0]
def clusters_hist(clusters, labels=None, remove_empty=True)

return non empty clusters histogramm sorted.

parameters

clusters : np.array
array of clusters membership of data.

returns

hist : np.array
count of element by clusters (decrasing hist)
label : np.array
label of the cluster aligned with hist
Source code
def clusters_hist(clusters, labels=None, remove_empty=True):
    """ return non empty clusters histogramm sorted.

        parameters
        ---------
        clusters: np.array
            array of clusters membership of data.

        returns
        -------
        hist: np.array
            count of element by clusters (decrasing hist)
        label: np.array
            label of the cluster aligned with hist
    """
    block_hist = np.bincount(clusters)
    if labels is None:
        labels = range(len(block_hist))

    hist, labels = sorted_perm(block_hist, labels, reverse=True)

    if remove_empty is True:
        null_classes = (hist == 0).sum()
        if null_classes > 0:
            hist = hist[:-null_classes]; labels = labels[:-null_classes]

    return hist, labels
def degree_hist(_degree, filter_zeros=False)
Source code
def degree_hist(_degree, filter_zeros=False):
    if isinstance(_degree, np.ndarray) and _degree.ndim == 2 :
        degree = list(dict(adj_to_degree(_degree)).values())
    elif isinstance(_degree, (list, np.ndarray)):
        degree = _degree
    else:
        # networkx
        degree = list(dict(_degree).values())

    max_c = np.max(degree)
    d = np.arange(max_c+1)
    dc = np.bincount(degree, minlength=max_c+1)

    if len(d) == 0:
        return [], []

    if dc[0] > 0:
        lgg.debug('%d unconnected vertex' % dc[0])
    d = d[1:]
    dc = dc[1:]


    if filter_zeros is True:
        #d, dc = zip(*filter(lambda x:x[1] != 0, zip(d, dc)))
        nzv = (dc != 0)
        d = d[nzv]
        dc = dc[nzv]

    return d, dc
def degree_hist_to_list(d, dc)
Source code
def degree_hist_to_list(d, dc):
    degree = np.repeat(np.round(d).astype(int), np.round(dc).astype(int))
    return degree
def dilate(y, size=1)
Source code
def dilate(y, size=1):
    dim = y.ndim
    mask = ndimage.generate_binary_structure(dim, dim)
    if size > 1:
        for i in range(1, size):
            mask = np.vstack((mask, mask[-1,:]))
            mask = np.column_stack((mask, mask[:, -1]))
    y_f = ndimage.binary_dilation(y, structure=mask).astype(y.dtype)
    return y_f
def draw_square(mat, value, topleft, l, L, w=0)
Source code
def draw_square(mat, value, topleft, l, L, w=0):
    tl = topleft

    # Vertical draw
    mat[tl[0]:tl[0]+l, tl[1]:tl[1]+w] = value
    mat[tl[0]:tl[0]+l, tl[1]+L-w:tl[1]+L] = value
    # Horizontal draw
    mat[tl[0]:tl[0]+w, tl[1]:tl[1]+L] = value
    mat[tl[0]+l-w:tl[0]+l, tl[1]:tl[1]+L] = value
    return mat
def floatRgb(mag, cmin, cmax)

Return a tuple of floats between 0 and 1 for the red, green and blue amplitudes.

Source code
def floatRgb(mag, cmin, cmax):
        """ Return a tuple of floats between 0 and 1 for the red, green and
                blue amplitudes.
        """

        try:
                # normalize to [0,1]
                x = float(mag-cmin)/float(cmax-cmin)
        except:
                # cmax = cmin
                x = 0.5
        blue = min((max((4*(0.75-x), 0.)), 1.))
        red  = min((max((4*(x-0.25), 0.)), 1.))
        green= min((max((4*math.fabs(x-0.5)-1., 0.)), 1.))
        return (red, green, blue)
def gem(gmma, K)
Source code
def gem(gmma, K):
    sb = np.empty(K)
    cut = np.random.beta(1, gmma, size=K)
    for k in range(K):
        sb[k] = cut[k] * cut[0:k].prod()
    return sb
def htmlRgb(mag, cmin, cmax)

Return a tuple of strings to be used in HTML documents.

Source code
def htmlRgb(mag, cmin, cmax):
   """ Return a tuple of strings to be used in HTML documents.
   """
   return "#%02x%02x%02x"%rgb(mag, cmin, cmax)
def lognormalize(x)
Source code
def lognormalize(x):
    return np.exp(x - np.logaddexp.reduce(x))
def random_degree(Y, params=None)
Source code
def random_degree(Y, params=None):
    _X = []
    _Y = []
    N = Y[0].shape[0]
    nb_uniq_degree = []
    dc_list = []
    for y in Y:
        ba_c = adj_to_degree(y)
        d, dc = degree_hist(ba_c)
        nb_uniq_degree.append(len(dc))
        dc_list.append(dc)
    dc_mat = ma.array(np.empty((N, max(nb_uniq_degree))), mask=True)
    for i, degrees in enumerate(dc_list):
        size = nb_uniq_degree[i]
        dc_mat[i, :size] = degrees

    y = dc_mat.mean(0)
    yerr = dc_mat.std(0)
    # 0 are filtered out in degree_hist
    return np.arange(1, len(y)+1),  np.round(y), yerr
def random_powerlaw(alpha, x_min, size=1)
Source code
def random_powerlaw(alpha, x_min, size=1):
    ### Discrete
    alpha = float(alpha)
    u = np.random.random(size)
    x = (x_min-0.5)*(1-u)**(-1/(alpha-1))+0.5
    return np.floor(x)
def reorder_mat(y, clusters, labels=False, reverse=True)

Reorder the matrix according the clusters membership @Debug: square matrix

Source code
def reorder_mat(y, clusters, labels=False, reverse=True):
    """Reorder the  matrix according the clusters membership
        @Debug: square matrix
    """
    assert(y.shape[0] == y.shape[1] == len(clusters))
    if reverse is True:
        hist, label = clusters_hist(clusters)
        sorted_clusters = np.empty_like(clusters)
        for i, k in enumerate(label):
            if i != k:
                sorted_clusters[clusters == k] = i
    else:
        sorted_clusters = clusters

    N = y.shape[0]
    nodelist = [k[0] for k in sorted(zip(range(N), sorted_clusters),
                                     key=lambda k: k[1])]

    y_r = y[nodelist, :][:, nodelist]
    if labels is True:
        return y_r, nodelist
    else:
        return y_r
def rgb(mag, cmin, cmax)

Return a tuple of integers to be used in AWT/Java plots.

Source code
def rgb(mag, cmin, cmax):
   """ Return a tuple of integers to be used in AWT/Java plots.
   """

   red, green, blue = floatRgb(mag, cmin, cmax)
   return (int(red*255), int(green*255), int(blue*255))
def shiftpos(arr, fr, to, axis=0)

Move element In-Place, shifting backward (or forward) others

Source code
def shiftpos(arr, fr, to, axis=0):
    """ Move element In-Place, shifting backward (or forward) others """
    if fr == to: return
    x = arr.T if axis == 1 else arr
    tmp = x[fr].copy()
    if fr > to:
        x[to+1:fr+1] = x[to:fr]
    else:
        x[fr:to] = x[fr+1:to+1]
    x[to] = tmp
def sorted_perm(a, label=None, reverse=False)

return sorted $a and the induced permutation. Inplace operation

Source code
def sorted_perm(a, label=None, reverse=False):
    """ return sorted $a and the induced permutation.
        Inplace operation """
    # np.asarray applied this tuple lead to error, if label is string
    # because a should be used as elementwise comparison
    if label is None:
        label = np.arange(a.shape[0])
    hist, label = zip(*sorted(zip(a, label), reverse=reverse))
    hist = np.asarray(hist)
    label = np.asarray(label)
    return hist, label
def strRgb(mag, cmin, cmax)

Return a tuple of strings to be used in Tk plots.

Source code
def strRgb(mag, cmin, cmax):
   """ Return a tuple of strings to be used in Tk plots.
   """

   red, green, blue = floatRgb(mag, cmin, cmax)
   return "#%02x%02x%02x" % (red*255, green*255, blue*255)
def wmean(a, w, mean='geometric')
Source code
def wmean(a, w, mean='geometric'):
    if mean == 'geometric':
        kernel = lambda x : np.log(x)
        out = lambda x : np.exp(x)
    elif mean == 'arithmetic':
        kernel = lambda x : x
        out = lambda x : x
    elif mean == 'harmonic':
        num = np.sum(w)
        denom = np.sum(np.asarray(w) / np.asarray(a))
        return num / denom
    else:
        raise NotImplementedError('Mean Unknwow: %s' % mean)

    num = np.sum(np.asarray(w) * kernel(np.asarray(a)))
    denom = np.sum(np.asarray(w))
    return out(num / denom)