SINGA-347 Create a function that supports einsum 1. finish the einsum function in python but still needs some funciton to do in cpp 2. function TODO list to finish the function in cpp(just like numpy function):a.sum(A,axis = None) b.repeat(A,repeats) c.transpose(A,axes = None)
Project: http://git-wip-us.apache.org/repos/asf/incubator-singa/repo Commit: http://git-wip-us.apache.org/repos/asf/incubator-singa/commit/fa1b0237 Tree: http://git-wip-us.apache.org/repos/asf/incubator-singa/tree/fa1b0237 Diff: http://git-wip-us.apache.org/repos/asf/incubator-singa/diff/fa1b0237 Branch: refs/heads/master Commit: fa1b023717b029df6136579be4c3bf95790c914b Parents: 2ec1364 Author: sheyujian <[email protected]> Authored: Fri Apr 20 13:00:49 2018 +0800 Committer: sheyujian <[email protected]> Committed: Fri Apr 20 13:00:49 2018 +0800 ---------------------------------------------------------------------- python/singa/tensor.py | 262 +++----------------------------------------- 1 file changed, 17 insertions(+), 245 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/incubator-singa/blob/fa1b0237/python/singa/tensor.py ---------------------------------------------------------------------- diff --git a/python/singa/tensor.py b/python/singa/tensor.py index 38acf6d..caf5998 100644 --- a/python/singa/tensor.py +++ b/python/singa/tensor.py @@ -940,251 +940,16 @@ def mult(A, B, C=None, alpha=1.0, beta=0.0): return C -def tensordot(A, B, axes=2): - """Returns the tensor multiplication of two tensors along specified axes. - This is equivalent to compute dot product along the specified axes which - are treated as one axis by reshaping. - - Args: - A (Singa.Tensor): The first argument. - B (Singa.Tensor): The second argument. - axes: - - If it is an integer, then ''axes'' represent axes at the last of ''a`'' and - the first of ''b'' are used. - - If it is a pair of sequences of integers, then these two - sequences specify the list of axes for ''a'' and ''b''. The - corresponding axes are paired for sum-product. - - Return: - singa.tensor: The tensor product of ''a'' and ''b'' along the - axes specified by ''axes''. - - Thanks to numpy.tensordot. - the link is https://github.com/numpy/numpy/blob/v1.14.0/numpy/core/numeric.py#L1123-L1306 - """ - # when axes is an integer, axes_A and axes_B represent axes at the last of ''a`'' and - # the first of ''b''. For example, when axes is 1, we do the normal multiplication : - # if A is in shape(3,2,4), B is in shape(4,2,5), it will return a matrix in shape(3,2,2,5) - # when axes is 2 and A,B are in the shape (3,2,4) and (2,4,5), it will return a matrix in shape(3,5) - if type(axes) == int: - axes_A = list(range(-axes, 0)) - axes_B = list(range(0, axes)) - axes_B = axes_B - else: - axes_A,axes_B =axes - # when axes is a pair of sequences of integers.For example, A is in shape(3,2,4), - #B is in shape(4,2,5), we set axes as ([1,2],[1,0]), it will return a matrix in shape(3,5) - if isinstance(axes_A,list): - na = len(axes_A) - axes_A = list(axes_A) - else: - axes_A = [axes_A] - na = 1 - if isinstance(axes_B,list): - nb = len(axes_B) - axes_B = list(axes_B) - else: - axes_B = [axes_B] - nb = 1 - - # a_shape and b_shape are the shape of tensor A and B, while nda and ndb are the dim of A and B - a_shape = A.shape - nda = A.ndim - b_shape = B.shape - ndb = B.ndim - equal = True - # to check if the length of axe_A is equal to axes_B - if na != nb: - equal = False - else: - # to make the shape match - for k in range(na): - if a_shape[axes_A[k]] != b_shape[axes_B[k]]: - equal = False - break - if axes_A[k] < 0: - axes_A[k] += nda - if axes_B[k] < 0: - axes_B[k] += ndb - if not equal: - raise ValueError("shape-mismatch for sum") - '''start to do the calculation according to the axes''' - - notin = [k for k in range(nda) if k not in axes_A] - # nda is the dim of A, and axes_a is the axis for A, notin is the axis which is not in axes_A - newaxes_a = notin + axes_A - N2 = 1 - for axis in axes_A: - N2 *= a_shape[axis] - N1 = 1 - for ax in notin: - N1 *=a_shape[ax] - #newshape_a is the shape to do multiplication.For example, A is in shape(3,2,4), - #B is in shape(4,2,5), we set axes as ([1,2],[1,0]), then newshape_a should be (3,5) - #olda is the shape that will be shown in the result. - newshape_a = (N1,N2) - olda = [a_shape[axis] for axis in notin] - notin = [k for k in range(ndb) if k not in axes_B] - newaxes_b = axes_B + notin - N2 = 1 - for axis in axes_B: - N2 *= b_shape[axis] - N1 = 1 - for bx in notin: - N1 *= b_shape[bx] - newshape_b = (N2, N1) - oldb = [b_shape[axis] for axis in notin] - # do transpose and reshape to get the 2D matrix to do multiplication - at = A.transpose(newaxes_a).reshape(newshape_a) - bt = B.transpose(newaxes_b).reshape(newshape_b) - res = mult(at, bt) - #reshape the result - return res.reshape(olda + oldb) - -def einsum_(A,B,ops): - '''Do the matrix to matrix einsum calculation according to the operands - - Args: - A (Singa.Tensor): The first argument. - B (Singa.Tensor): The second argument. - ops(string): - the string specifies the subscripts for summation such as 'ki,kj->kij' - - Returns: Singa.Tensor - the output matirx of the einsum calculation - Thanks to nils-werner/sparse.einsum(), - the link is https://github.com/nils-werner/sparse/commit/449c75d21d3158630bc5be79c690a60cfc002578 +def einsum(A,B,ops): ''' - if len(ops) == 0: - raise ValueError("No input operands") - # to get the input and output ops - inputops, outputops = ops.split('->') - inputops = inputops.split(',') - - if A.ndim != len(inputops[0]) or B.ndim != len(inputops[1]): - raise ValueError("input dim doesn't match operands") - - # All indices that are in input AND in output are multiplies - # All indices that are in input BUT NOT in output are sum contractions - multiplies = sorted(list(set(inputops[0]) & set(inputops[1]) & set(outputops))) - sums = sorted(list((set(inputops[0]) & set(inputops[1])) - set(outputops))) - - # Map sums and indices to axis integers - multiplies = [[inops.find(x) for x in multiplies] for inops in inputops] - sums = [[inops.find(x) for x in sums] for inops in inputops] + function TODO list to finish the function in cpp(just like numpy function): + 1.sum(A,axis = None) + 2.repeat(A,repeats) + 3.transpose(A,axes = None) - - # Find output axes in input axes for final transpose - transpose = [''.join(inputops).find(x) for x in outputops] - #to make the transpose match to its rank - transpose = sorted(range(len(transpose)), key = transpose.__getitem__) - - return tensormult(A,B,sums,multiplies).transpose(transpose) - - -def tensormult(A, B, sum=None, multiply=None): - ''' - Args: - A (Singa.Tensor): The first input tensor from einsum. - B (Singa.Tensor): The second input tensor from einsum. - sum (list[list[int]]): - The axis to do the normal tensordot calculation - - multiply (list[list[int]]): - The axis to multiply - - Return: (Singa.Tensor) - The output of the tensormult calculation - - fix some error of source code in the nils-werner/sparse - ''' - if sum is None: - sum = [[], []] - else: - sum = list(sum) - - if multiply is None: - multiply = [[], []] - else: - multiply = list(multiply) - # For each multiply[0] we are adding one axis, thus we need to increment - # all following items by one: (0, 1, 2) -> (0, 2, 4) - idx = multipliessort(multiply[0]) - - post_multiply = multiply[0] - for i, v in enumerate(idx): - post_multiply[v] += i - - -# for i in post_multiply: -# A = diag(A,i) - for i in sorted(post_multiply): - - A = diag(A,i) - sum[0] = [x+1 if x >= i else x for x in sum[0]] - - sum[0] += post_multiply - sum[1] += multiply[1] - - return tensordot(A, B, axes=sum) - -def multipliessort(multiplies): - ''' - Returns the indices that would sort an array - Args: - multiplies(list[int]): - the list of the input multiplies - - Return(list[int]): - the indices list - ''' - - if multiplies is None: - multiplies = [] - else: - multiplies = list(multiplies) - - sort_multiplies = sorted(enumerate(multiplies), key=lambda x: x[1]) - idx = [x[0] for x in sort_multiplies] - return idx - -def diag(A,axis=-1): - ''' - Return the matrix that is diagonalize along the axis - Args: - A (Singa.Tensor):the input tensor - axis (int): - the axis which the matrix use to diagonalization - Return(Singa.Tensor): - the tensor which has been diagonalize along the given axis - - - It is like using numpy.diag() in the 1D condition, but this function can be used to do high dimension - matrix diagonalization - - ''' - # to get the shape of the diagonalize matirx - axis_diag = list(range(out.ndim))[axis] - shape = list(A.shape) - shape.insert(axis_diag,shape[axis_diag]) - - npA = to_numpy(A) - npA_diag = np.zeros(product(shape)) - index = np.argwhere(npA != 0) - for i in index: - i_diag = list(i) - i_diag.insert(axis, i[axis]) - index_npA = [[i[x]] for x in range(len(i))] - index_npA_diag = [[i_diag[x]] for x in range(len(i_diag))] - npA_diag[index_npA_diag] = npA[index_npA] - A_diag = from_numpy(npA_diag) - return A_diag - -def einsum2(A,B,ops): - ''' - Do the matrix to matrix einsum calculation according to the operands + Do the matrix to matrix einsum calculation according to the operands Args: A (Singa.Tensor): The first argument. @@ -1196,23 +961,27 @@ def einsum2(A,B,ops): the output matirx of the einsum calculation ''' + if len(ops) == 0: raise ValueError("No input operands") + # to get the input and output ops inputops, outputops = ops.split('->') inputops = inputops.split(',') if A.ndim != len(inputops[0]) or B.ndim != len(inputops[1]): raise ValueError("input dim doesn't match operands") + #to get the indices in input but not in output sums = sorted(list((set(inputops[0]) | set(inputops[1])) - set(outputops))) + #to get the indices that A and B use to broadcast to each other broadcast_A = sorted(list(set(inputops[1]) - set(inputops[0]))) broadcast_B = sorted(list(set(inputops[0]) - set(inputops[1]))) - - + # to get all the indices in input outputall = sorted(list(set(inputops[0]) | set(inputops[1]))) + ## Map indices to axis integers sums = [outputall.index(x) for x in sums] broadcast_idA = [inputops[1].find(x) for x in broadcast_A] broadcast_idB = [inputops[0].find(x) for x in broadcast_B] @@ -1220,10 +989,10 @@ def einsum2(A,B,ops): broadcast_a = [B.shape[x] for x in broadcast_idA] broadcast_b = [A.shape[x] for x in broadcast_idB] + # get the the transpose and reshape parameter used in the elementwise calculation transpose_A = [(list(inputops[0])+broadcast_A) .index(x) for x in outputall] transpose_B = [(list(inputops[1])+broadcast_B) .index(x) for x in outputall] - reshape_A = list(A.shape)+broadcast_a reshape_B = list(B.shape)+broadcast_b @@ -1231,9 +1000,10 @@ def einsum2(A,B,ops): mult_B = np.repeat(B, np.product(broadcast_b)).reshape(reshape_B).transpose(transpose_B) if mult_A.shape != mult_B.shape: - raise ValueError("error: matrix dimension mismatch") + raise ValueError("Error: matrix dimension mismatch") res = eltwise_mult(mult_A, mult_B) + # reduce the axis and find the final transpose for the output sum_R = sorted(sums, reverse=True) for i in sum_R: res = res.sum(axis=i) @@ -1243,6 +1013,8 @@ def einsum2(A,B,ops): + + def div(lhs, rhs, ret=None): '''Elementi-wise division.
