Module that defines indexed objects
The classes IndexedBase, Indexed and Idx would represent a matrix element M[i, j] as in the following graph:
1) The Indexed class represents the entire indexed object.

______
' '
M[i, j]
/ \__\______
 
 
 2) The Idx class represent indices and each Idx can
 optionally contain information about its range.

3) IndexedBase represents the `stem' of an indexed object, here `M'.
The stem used by itself is usually taken to represent the entire
array.
There can be any number of indices on an Indexed object. No transformation properties are implemented in these Base objects, but implicit contraction of repeated indices is supported.
Note that the support for complicated (i.e. nonatomic) integer expressions as indices is limited. (This should be improved in future releases.)
To express the above matrix element example you would write:
>>> from sympy.tensor import IndexedBase, Idx
>>> from sympy import symbols
>>> M = IndexedBase('M')
>>> i, j = map(Idx, ['i', 'j'])
>>> M[i, j]
M[i, j]
Repreated indices in a product implies a summation, so to express a matrixvector product in terms of Indexed objects:
>>> x = IndexedBase('x')
>>> M[i, j]*x[j]
M[i, j]*x[j]
If the indexed objects will be converted to component based arrays, e.g. with the code printers or the autowrap framework, you also need to provide (symbolic or numerical) dimensions. This can be done by passing an optional shape parameter to IndexedBase upon construction:
>>> dim1, dim2 = symbols('dim1 dim2', integer=True)
>>> A = IndexedBase('A', shape=(dim1, 2*dim1, dim2))
>>> A.shape
(dim1, 2*dim1, dim2)
>>> A[i, j, 3].shape
(dim1, 2*dim1, dim2)
If an IndexedBase object has no shape information, it is assumed that the array is as large as the ranges of it’s indices:
>>> n, m = symbols('n m', integer=True)
>>> i = Idx('i', m)
>>> j = Idx('j', n)
>>> M[i, j].shape
(m, n)
>>> M[i, j].ranges
[(0, m  1), (0, n  1)]
The above can be compared with the following:
>>> A[i, 2, j].shape
(dim1, 2*dim1, dim2)
>>> A[i, 2, j].ranges
[(0, m  1), None, (0, n  1)]
To analyze the structure of indexed expressions, you can use the methods get_indices() and get_contraction_structure():
>>> from sympy.tensor import get_indices, get_contraction_structure
>>> get_indices(A[i, j, j])
(set([i]), {})
>>> get_contraction_structure(A[i, j, j])
{(j,): set([A[i, j, j]])}
See the appropriate docstrings for a detailed explanation of the output.
Represent the base or stem of an indexed object
The IndexedBase class represent an array that contains elements. The main purpose of this class is to allow the convenient creation of objects of the Indexed class. The __getitem__ method of IndexedBase returns an instance of Indexed. Alone, without indices, the IndexedBase class can be used as a notation for e.g. matrix equations, resembling what you could do with the Symbol class. But, the IndexedBase class adds functionality that is not available for Symbol instances:
 An IndexedBase object can optionally store shape information. This can be used in to check array conformance and conditions for numpy broadcasting. (TODO)
 An IndexedBase object implements syntactic sugar that allows easy symbolic representation of array operations, using implicit summation of repeated indices.
 The IndexedBase object symbolizes a mathematical structure equivalent to arrays, and is recognized as such for code generation and automatic compilation and wrapping.
>>> from sympy.tensor import IndexedBase, Idx
>>> from sympy import symbols
>>> A = IndexedBase('A'); A
A
>>> type(A)
<class 'sympy.tensor.indexed.IndexedBase'>
When an IndexedBase object recieves indices, it returns an array with named axes, represented by an Indexed object:
>>> i, j = symbols('i j', integer=True)
>>> A[i, j, 2]
A[i, j, 2]
>>> type(A[i, j, 2])
<class 'sympy.tensor.indexed.Indexed'>
The IndexedBase constructor takes an optional shape argument. If given, it overrides any shape information in the indices. (But not the index ranges!)
>>> m, n, o, p = symbols('m n o p', integer=True)
>>> i = Idx('i', m)
>>> j = Idx('j', n)
>>> A[i, j].shape
(m, n)
>>> B = IndexedBase('B', shape=(o, p))
>>> B[i, j].shape
(o, p)
Represents an index, either symbolic or integer.
There are a number of ways to create an Idx object. The constructor takes two arguments:
Optionally you can specify a range as either
 Symbol or integer: This is interpreted as dimension. lower and upper ranges are set to 0 and range1
 tuple: This is interpreted as the lower and upper bounds in the range.
Note that the Idx constructor is rather pedantic, and will not accept noninteger symbols. The only exception is that you can use oo and oo to specify an unbounded range. For all other cases, both label and bounds must be declared as integers, in the sense that for a index label n, n.is_integer must return True.
For convenience, if the label is given as a string, it is automatically converted to an integer symbol. (Note that this conversion is not done for range or dimension arguments.)
Examples : 

>>> from sympy.tensor import IndexedBase, Idx
>>> from sympy import symbols, oo
>>> n, i, L, U = symbols('n i L U', integer=True)
0) Construction from a string. An integer symbol is created from the string.
>>> Idx('qwerty')
qwerty
>>> idx = Idx(i, (L, U)); idx
i
>>> idx.lower, idx.upper
(L, U)
>>> idx = Idx(i, n); idx.lower, idx.upper
(0, n  1)
>>> idx = Idx(i, 4); idx.lower, idx.upper
(0, 3)
>>> idx = Idx(i, oo); idx.lower, idx.upper
(0, oo)
>>> idx = Idx(i); idx.lower, idx.upper
(None, None)
4) for a literal integer instead of a symbolic label the bounds are still there:
>>> idx = Idx(2, n); idx.lower, idx.upper
(0, n  1)
Represents a mathematical object with indices.
>>> from sympy.tensor import Indexed, IndexedBase, Idx
>>> from sympy import symbols
>>> i, j = map(Idx, ['i', 'j'])
>>> Indexed('A', i, j)
A[i, j]
It is recommended that Indexed objects are created via IndexedBase:
>>> A = IndexedBase('A')
>>> Indexed('A', i, j) == A[i, j]
True
returns a list of tuples with lower and upper range of each index
If an index does not define the data members upper and lower, the corresponding slot in the list contains ``None’’ instead of a tuple.
returns a list with dimensions of each index.
Dimensions is a property of the array, not of the indices. Still, if the IndexedBase does not define a shape attribute, it is assumed that the ranges of the indices correspond to the shape of the array.
>>> from sympy.tensor.indexed import IndexedBase, Idx
>>> from sympy import symbols
>>> n, m = symbols('n m', integer=True)
>>> i = Idx('i', m)
>>> j = Idx('j', m)
>>> A = IndexedBase('A', shape=(n, n))
>>> B = IndexedBase('B')
>>> A[i, j].shape
(n, n)
>>> B[i, j].shape
(m, m)
Determine the outer indices of expression expr
By outer we mean indices that are not summation indices. Returns a set and a dict. The set contains outer indices and the dict contains information about index symmetries.
Examples : 

>>> from sympy.tensor.index_methods import get_indices
>>> from sympy import symbols
>>> from sympy.tensor import IndexedBase, Idx
>>> x, y, A = map(IndexedBase, ['x', 'y', 'A'])
>>> i, j, a, z = symbols('i j a z', integer=True)
The indices of the total expression is determined, Repeated indices imply a summation, for instance the trace of a matrix A:
>>> get_indices(A[i, i])
(set(), {})
In the case of many terms, the terms are required to have identical outer indices. Else an IndexConformanceException is raised.
>>> get_indices(x[i] + A[i, j]*y[j])
(set([i]), {})
Exceptions : 

An IndexConformanceException means that the terms ar not compatible, e.g.
>>> get_indices(x[i] + y[j])
(...)
IndexConformanceException: Indices are not consistent: x(i) + y(j)
Warning
The concept of outer indices applies recursively, starting on the deepest level. This implies that dummies inside parenthesis are assumed to be summed first, so that the following expression is handled gracefully:
>>> get_indices((x[i] + A[i, j]*y[j])*x[j])
(set([i, j]), {})
This is correct and may appear convenient, but you need to be careful with this as Sympy wil happily .expand() the product, if requested. The resulting expression would mix the outer j with the dummies inside the parenthesis, which makes it a different expression. To be on the safe side, it is best to avoid such ambiguities by using unique indices for all contractions that should be held separate.
Determine dummy indices of expr and describe it’s structure
By dummy we mean indices that are summation indices.
The stucture of the expression is determined and described as follows:
A conforming summation of Indexed objects is described with a dict where the keys are summation indices and the corresponding values are sets containing all terms for which the summation applies. All Add objects in the Sympy expression tree are described like this.
For all nodes in the Sympy expression tree that are not of type Add, the following applies:
If a node discovers contractions in one of it’s arguments, the node itself will be stored as a key in the dict. For that key, the corresponding value is a list of dicts, each of which is the result of a recursive call to get_contraction_structure(). The list contains only dicts for the nontrivial deeper contractions, ommitting dicts with None as the one and only key.
Note
The presence of expressions among the dictinary keys indicates multiple levels of index contractions. A nested dict displays nested contractions and may itself contain dicts from a deeper level. In practical calculations the summation in the deepest nested level must be calculated first so that the outer expression can access the resulting indexed object.
Examples : 

>>> from sympy.tensor.index_methods import get_contraction_structure
>>> from sympy import symbols
>>> from sympy.tensor import IndexedBase, Idx
>>> x, y, A = map(IndexedBase, ['x', 'y', 'A'])
>>> i, j, k, l = map(Idx, ['i', 'j', 'k', 'l'])
>>> get_contraction_structure(x[i]*y[i] + A[j, j])
{(i,): set([x[i]*y[i]]), (j,): set([A[j, j]])}
>>> get_contraction_structure(x[i]*y[j])
{None: set([x[i]*y[j]])}
A multiplication of contracted factors results in nested dicts representing the internal contractions.
>>> d = get_contraction_structure(x[i, i]*y[j, j])
>>> sorted(d.keys())
[None, x[i, i]*y[j, j]]
>>> d[None] # Note that the product has no contractions
set([x[i, i]*y[j, j]])
>>> sorted(d[x[i, i]*y[j, j]]) # factors are contracted ``first''
[{(i,): set([x[i, i]])}, {(j,): set([y[j, j]])}]
A parenthesized Add object is also returned as a nested dictionary. The term containing the parenthesis is a Mul with a contraction among the arguments, so it will be found as a key in the result. It stores the dictionary resulting from a recursive call on the Add expression.
>>> d = get_contraction_structure(x[i]*(y[i] + A[i, j]*x[j]))
>>> sorted(d.keys())
[(i,), x[i]*(y[i] + A[i, j]*x[j])]
>>> d[(i,)]
set([x[i]*(y[i] + A[i, j]*x[j])])
>>> d[x[i]*(A[i, j]*x[j] + y[i])]
[{None: set([y[i]]), (j,): set([A[i, j]*x[j]])}]
Powers with contractions in either base or exponent will also be found as keys in the dictionary, mapping to a list of results from recursive calls:
>>> d = get_contraction_structure(A[j, j]**A[i, i])
>>> d[None]
set([A[j, j]**A[i, i]])
>>> nested_contractions = d[A[j, j]**A[i, i]]
>>> nested_contractions[0]
{(j,): set([A[j, j]])}
>>> nested_contractions[1]
{(i,): set([A[i, i]])}
The description of the contraction structure may appear complicated when represented with a string in the above examples, but it is easy to iterate over:
>>> from sympy import Expr
>>> for key in d:
... if isinstance(key, Expr):
... continue
... for term in d[key]:
... if term in d:
... # treat deepest contraction first
... pass
... # treat outermost contactions here