<gh_stars>10-100
from __future__ import division
from builtins import zip
from ektelo import util
from ektelo.matrix import EkteloMatrix
from functools import reduce
import math
import numpy as np
from scipy import sparse
from ektelo import workload

def cantor_pairing(a, b):
    """
    A function returning a unique positive integer for every pair (a,b) of positive integers
    """
    return (a+b)*(a+b+1)/2 + b
    
def _replace(vector, new_values):
    for i in range(len(vector)):
        vector[i] = new_values[ vector[i] ]
    return vector

def get_partition_vec(rank,n,cluster,closeRange=False):
    """ get the partition vector from clusters returned by partition algorithms
        rank: If the bins are sorted, give the rank of each item in the input list.
            Used by AHP partition. Set rank = None if not used.
        n: Length of vector in original domain
        cluster: Cluster/partition returned by partition algorithms
        closeRange: if set to True, ranges in clusters are close range. (DAWA partition)
            i.e. [a,b] indicates [a,a+1,...b-1,b]
            if set to False, ranges in clusters are default python representation. (AHP partition)
            i.e. [a,b] indicates [a,a+1,...b-1]
    """
    partition_vec_sorted = np.empty(n,int)
    assert cluster[0][0] == 0,"First bin of partition must start with 0"
    # assign groupID to elements in sorted list.
    for i in range(len(cluster)):
        if closeRange:
            assert cluster[-1][1] == n-1, " Last bin of partition must end with length of original data"
            partition_vec_sorted[cluster[i][0]:cluster[i][1]+1] = i
        else:
            assert cluster[-1][1] == n, " Last bin of partition must end with length of original data"
            partition_vec_sorted[cluster[i][0]:cluster[i][1]] = i
    # get index in sorted list for elements in original domain, then get groupID.
    if rank is None:
        partition_vec = partition_vec_sorted
    else:
        partition_vec = np.array([partition_vec_sorted[rank[i]] for i in range(n)] )

    return partition_vec


def update_corners(corner, groupID, row, start, end):
    ''' helper function for get_subdomain
        update corners coordinates for a certain group.
        return False if the domain is not rectangular
    '''
    # if it is the first ocurrence of the group
    # update the upper left and upper right corner
    if groupID not in corner:
        corner[groupID] = {'ul':(row, start),'ur':(row,end), 'll':(row, start),'lr':(row,end)}
    else:
        temp = corner[groupID]
        if row == temp['ll'][0]: # incontinous group on the upper line
            return False 

        # update the lower corners
        # make sure the columns match and rows are continous.
        if temp['ll'][1] == start and temp['lr'][1] == end and temp['ll'][0] == row-1:
            # move the lower corners one line lower
            corner[groupID]['ll'] = (temp['ll'][0]+1, temp['ll'][1])
            corner[groupID]['lr'] = (temp['lr'][0]+1, temp['lr'][1])
        else:
            return False

    return True



def get_subdomain_grid(mapping, domain_shape):
    '''
    Given a mapping, return the domain size of all the subdomain when it is 
    used by the SplitByPartition operator.
    The original domain needs to be 2D and the mapping should split the domain 
    to smaller grids. Non-rectangular subdomain shapes are not supported,
    None will be returned.

    '''
    assert len(domain_shape) == 2 , 'Only works for 2D domains'
    m, n = domain_shape
    # unflatten the mapping vector
    mapping = mapping.reshape(domain_shape)
    corners = {}
    # record corners of each group in one pass of the mapping vector
    for i in range(m):
        start = 0
        for j in range(n):
            if j+1 >= n or mapping[i][j] != mapping[i][j+1]:
                groupID = mapping[i][start]
                status = update_corners(corners, groupID, i, start, j)
                start = j+1
                if status == False:
                    return None

    # calculate subdomains from corners
    sub_domains = {}
    for g in corners:
        temp = corners[g]
        sub_domains[g] = (temp['ll'][0] - temp['ul'][0] + 1, temp['ur'][1] - temp['ul'][1] + 1)

    return sub_domains


def canonical_ordering(mapping):
    """ remap according to the canonical order.
     if bins are noncontiguous, use position of first occurrence.
     e.g. [3,4,1,1] => [1,2,3,3]; [3,4,1,1,0,1]=>[0,1,2,2,3,2]
    """
    unique, indices, inverse, counts = mapping_statistics(mapping)

    uniqueInverse, indexInverse = np.unique(inverse,return_index =True)

    indexInverse.sort()
    newIndex = inverse[indexInverse]
    tups = list(zip(uniqueInverse, newIndex)) 
    tups.sort(key=lambda x: x[1])
    u = np.array( [u for (u,i) in tups] )
    mapping = u[inverse].reshape(mapping.shape)

    return mapping


def mapping_statistics(mapping):
    return np.unique(mapping, return_index=True, return_inverse=True, return_counts=True)   




def reduction_matrix(mapping, canonical_order=False):
    """ Returns an m x n matrix R where n is the dimension of 
        the original data and m is the dimension of the reduced data.

        Reduces data vector x with R x
        Expands workload matrix W with W' R
    """
    assert mapping.ndim == 1, "Can only handle 1-dimesional mappings for now, domain should be flattened"

    unique, indices, inverse, counts = mapping_statistics(mapping)

    if canonical_order:
        mapping = canonical_ordering(mapping)

    n = mapping.size
    m = unique.size
    data = np.ones(n)
    cols = np.arange(n)
    rows = inverse

    return EkteloMatrix(sparse.csr_matrix((data, (rows, cols)), shape=(m, n), dtype=int))

def expansion_matrix(mapping, canonical_order=False):
    """ Returns an n x m matrix E where n is the dimension of 
        the original data and m is the dimension of the reduced data.

        Expands data vector x with E x'
        Reduces workload matrix W with W E
    """
    assert mapping.ndim == 1, "Can only handle 1-dimesional mappings for now, domain should be flattened"

    unique, indices, inverse, counts = mapping_statistics(mapping)

    if canonical_order:
        mapping = canonical_ordering(mapping)

    n = mapping.size
    m = unique.size
    data = np.ones(n)
    cols = np.arange(n)
    rows = inverse

    R = sparse.csr_matrix((data, (rows, cols)), shape=(m, n), dtype=int)
    scale = sparse.spdiags(1.0 /counts, 0, m, m)

    return EkteloMatrix(R.T * scale)

def projection_matrix(mapping, idx):
    """ Returns m x n matrix P where n is the dimension of the 
        original data and m is the number of occurence of idx
        in mapping.

        :param mapping: vector with indices representing groups
        :param idx: index of group from which to create projection

        Projects vector x with P x and matrix W with W P^T
        Unprojects vector x with P^T x and matrix W with W P
    """
    mask = np.ma.masked_where(mapping!=idx, mapping).mask

    if np.all(~mask): # when all entries are False, a single False will be returned
        mask = np.array([False]*len(mapping))

    cols = np.where(~mask)[0]
    rows = np.arange(cols.size)
    vals = np.ones_like(rows)
    P = sparse.csr_matrix((vals, (rows, cols)), (rows.size, mask.size))

    return EkteloMatrix(P)

def combine(p1, p2):
    """ Returns p3, an (n+m) dimensional array of integers such that
        p3[i,j] = p3[i', j'] iff p1[i] = p1[i'] and p2[j] = p2[j']

        :param p1:  an n dimensional array of integers
        :param p2:  an m dimensional array of integers
    """
    _, i1 = np.unique(p1.flatten(), return_inverse=True)
    _, i2 = np.unique(p2.flatten(), return_inverse=True)
    lookup = np.arange(i1.size * i2.size).reshape(i1.size, i2.size)
    # note: cartesian product, but order is very important
    # this order works when flattening/reshaping is done in row-major form
    pairs = np.dstack(np.meshgrid(i1, i2, indexing='ij')).reshape(-1,2)
    flat = lookup[pairs[:,0], pairs[:,1]]

    return flat.reshape(p1.shape + p2.shape)


def combine_all(mappings):
    """ Returns an ndarray with each dimension corresponding to one
        of mapping.
    """
    return reduce(combine, mappings, np.ones((), dtype=int))


def extract_M(W):
    assert type(W) is sparse.csr_matrix, 'W must by csr_sparse'

    return W.getrow(W.nonzero()[0][0])


def complement(A, grid_size=None):
    '''return the queries on the complementary domain
    :param grid_size: The griding size of the new queris, if None, return total on the complementary domain
    Currently complementary domain are those indices with column norm(L1) 0.
    '''
    comp = []
    if isinstance(A, np.ndarray) is False:
        A = A.toarray()
    norm = np.linalg.norm(A,ord = 1,axis = 0)

    compl_size = len(norm) - np.count_nonzero(norm)
    grid_size = compl_size if grid_size is None else grid_size
    grid_num = int(math.ceil(compl_size/float(grid_size)))
    if grid_num==0:
        return None

    ind = 0
    for g in range(grid_num):
        q = np.zeros(len(norm))
        remain_in_group = grid_size
        while (remain_in_group>0) and (ind<len(norm)):
            if np.isclose(norm[ind],0.0):
                q[ind]=1
                remain_in_group-=1
            ind +=1
        comp.append(q)

    return sparse.csr_matrix(comp)
