[Numpy-discussion] Tensor contraction

[Alan Bromborsky]

Friedrich Romstedt wrote:
> 2010/6/13 Alan Bromborsky <abrombo at verizon.net>:
>   
>> Friedrich Romstedt wrote:
>>     
>>>> 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.
>>>       
>
> Now, after I read the Wikipedia article about the Christoffel symbols,
> I'm not so sure if I can help you with the subject.  It seems general
> relativity is a bit over my head, just a tiiiny bit ;-)
>
>   
>>>> 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.
>>>       
>
> Ok, this I understood now, nevertheless I'm completely lost with an
> over-all understanding ...
>
>   
>>>> 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.
>>>>         
>
> This is up to you, as I see it now.
>
>   
>>> 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
>>>       
>
> Btw, did this work for you?
>
> I like Sebastian's stride tricks, but I think the poor-man's
> .transpose() and contraction should also do it.  (and is maybe more
> readable.)
>
> Btw, how are you planning to implement the differentiation?  Hmm ...
> for the "normal" partial derivative ... it would be most elegant to
> create the \partial_{l} operator as an sympy object, is that possible?
>  If not, one can think about a wrapper class, yielding the
> differentiation upon multiplication with a normal sympy object.  If
> one has such operators, one can put them in an array, say p for
> "partial" and write:
>
> partial_T = p[numpy.newaxis, numpy.newaxis, numpy.newaxis] * T
>
> Or just write p as an instance of a class P with overloaded __mul__:
>
> def __mul__(self, T):
>     return self.partials[tuple([numpy.newaxis] * T.ndim)] * T
>
> Then partial differentiation is for all tensors of your
> multidimensional world just:
>
> partial_T = p * T
>
> . ? For the Christoffel symbols, it's I think not so easy.  But since
> I don't understand them, ...
>
> Friedrich
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>   
I must be stupid. 'numpy.trace' does exactly what I want!