Friedrich Romstedt wrote: > 2010/6/13 Alan Bromborsky <abro...@verizon.net>: > >> I am writing symbolic tensor package for general relativity. In making >> symbolic tensors concrete >> I generate numpy arrays stuffed with sympy functions and symbols. >> > > That sound's interesting. > > >> The >> operations are tensor product >> (numpy.multiply.outer), permutation of indices (swapaxes), partial and >> covariant (both vector operators that >> increase array dimensions by one) differentiation, and contraction. >> > > I would like to know more precisely what this differentiations do, and > how it comes that they add an index to the tensor. > > >> I think I need to do the contraction last >> to make sure everything comes out correctly. Thus in many cases I would >> be performing multiple contractions >> on the tensor resulting from all the other operations. >> > > Hm, ok, so I guess I shall give my 1 cent now. > > Ok. > > # First attempt (FYI, failed): > > # The general procedure is, to extract a multi-dimensional diagonal array. > # The sum \sum_{ij = 0}^{M} \sum_{kl = 0}^{N} is actually the sum over a > # 2D array with indices I \equiv i \equiv j and K \equiv k \equiv > l. Meaning: > # \sum_{(I, K) = (0, 0)}^{(M, N)}. > # Thus, if we extract the indices with 2D arrays [[0], [1], ..., > [N - 1]] for I and > # [[0, 1, ..., M - 1]] on the other side for K, then numpy's broadcasting > # mechanism will broadcast them to the same shape, yielding (N, M) arrays. > # Then finally we sum over this X last dimensions when there were X > # contractions, and we're done. > > # Hmmm, when looking closer at the problem, it seems that this isn't > # adequate. Because we would have to insert open slices, but cannot > # create them outside of the [] operator ... > > # So now follows second attemt: > > def contract(arr, *contractions): > """*CONTRACTIONS is e.g.: > (0, 1), (2, 3) > meaning two contractions, one of 0 & 1, and one of 2 & 2, > but also: > (0, 1, 2), > is allowed, meaning contract 0 & 1 & 2.""" > > # First, we check if we can contract using the *contractions* given ... > > for contraction in contractions: > # Extract the dimensions used. > dimensions = numpy.asarray(arr.shape)[list(contraction)] > > # Check if they are all the same. > dimensionsdiff = dimensions - dimensions[0] > if numpy.abs(dimensionsdiff).sum() != 0: > raise ValueError('Contracted indices must be of same dimension.') > > # So now, we can contract. > # > # First, pull the contracted dimensions all to the front ... > > # The names of the indices. > names = range(arr.ndim) > > # Pull all of the contractions. > names_pulled = [] > for contraction in contractions: > names_pulled = names_pulled + list(contraction) > # Remove the pulled index names from the pool: > for used_index in contraction: > # Some more sanity check > if used_index not in names: > raise ValueError('Each index can only be used in one > contraction.') > names.remove(used_index) > > # Concatenate the pulled indices and the left-over indices. > names_final = names_pulled + names > > # Perform the swap. > arr = arr.transpose(names_final) > > # Perform the contractions ... > > for contraction in contractions: > # The respective indices are now, since we pulled them, the > frontmost indices: > ncontraction = len(contraction) > # The index array: > # shape[0] = shape[1] = ... = shape[ncontraction - 1] > I = numpy.arange(0, arr.shape[0]) > # Perform contraction: > index = [I] * ncontraction > arr = arr[tuple(index)].sum(axis=0) > > # If I made no mistake, we should be done now. > return arr > > Ok, it didn't get much shorter than Pauli's solution, so you decide ... > > >> One question to >> ask would be considering that I am stuffing >> the arrays with symbolic objects and all the operations on the objects >> would be done using the sympy modules, >> would using numpy operations to perform the contractions really save any >> time over just doing the contraction in >> python code with a numpy array. >> > > I don't know anything about sympy. I think there's some typo around: > I guess you mean creating some /sympy/ array and doing the operations > using that instead of using a numpy array having sympy dtype=object > content? > > Friedrich > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion > > No I stuff a numpy array with sympy objects. Except for contraction and differentiation numpy does everything I need to do with the tensors. Namely -
addition subtraction multiplication by a scalar (sympy object) tensor product (multiply.outer) index permutation (swapaxes) Plus - contraction differentiation There are two types of vector derivatives (think gradient) partial and covariant. In index notation the partial derivative of an array T_{ijk} is the array \partial_{l}T_{ijk} = \frac{\partial T_{ijk}}{\partial x^{l}} The covariant derivative is D_{l}T_{ijk} = \partial_{l}T_{ijk} +T_{uvw}\Gamma_{i}^{u}_{l}\Gamma_{j}^{v}_{l}\Gamma_{k}^{w}_{l} Where the \Gamma_{i}^{j}_{k} are the Christoffle symbols with appropriately raised or lowered indices. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion