def __getitem__(self, key):

        if isinstance(key, tuple):
            i, j = key
            try:
                i, j = self.key2ij(key)
                return self._smat.get((i, j), S.Zero)
            except (TypeError, IndexError):
                if isinstance(i, slice):
                    i = range(self.rows)[i]
                elif is_sequence(i):
                    pass
                else:
                    i = [i]
                if isinstance(j, slice):
                    j = range(self.cols)[j]
                elif is_sequence(j):
                    pass
                else:
                    j = [j]
                return self.extract(i, j)

        # check for single arg, like M[:] or M[3]
        if isinstance(key, slice):
            lo, hi = key.indices(len(self))[:2]
            L = []
            for i in range(lo, hi):
                m, n = divmod(i, self.cols)
                L.append(self._smat.get((m, n), S.Zero))
            return L

        i, j = divmod(a2idx(key, len(self)), self.cols)
        return self._smat.get((i, j), S.Zero)

    def __new__(cls, *args, **kw_args):
        """
        Constructor for the Permutation object.

        Examples
        ========

        >>> from sympy.combinatorics.permutations import Permutation
        >>> p = Permutation([0,1,2])
        >>> p
        Permutation([0, 1, 2])
        >>> q = Permutation([[0,1],[2]])
        >>> q
        Permutation([[0, 1], [2]])
        """
        if not args or not is_sequence(args[0]) or len(args) > 1 or \
           len(set(is_sequence(a) for a in args[0])) > 1:
            raise ValueError('Permutation argument must be a list of ints or a list of lists.')

        # 0, 1, ..., n-1 should all be present

        temp = [int(i) for i in flatten(args[0])]
        if set(range(len(temp))) != set(temp):
            raise ValueError("Integers 0 through %s must be present." % len(temp))

        cform = aform = None
        if args[0] and is_sequence(args[0][0]):
            cform = [list(a) for a in args[0]]
        else:
            aform = list(args[0])

        ret_obj = Basic.__new__(cls, (cform or aform), **kw_args)
        ret_obj._cyclic_form, ret_obj._array_form = cform, aform
        return ret_obj

    def _eval_subs(self, old, new):
        from sympy.geometry.point import Point

        if is_sequence(old) or is_sequence(new):
            old = Point(old)
            new = Point(new)
            return self._subs(old, new)

    def __new__(cls, periodical, limits=None):
        x, start, stop = None, None, None
        if limits is None:
            x, start, stop = Dummy("k"), 0, S.Infinity
        if is_sequence(limits, Tuple):
            if len(limits) == 3:
                x, start, stop = limits
            elif len(limits) == 2:
                x = Dummy("k")
                start, stop = limits

        if not isinstance(x, Symbol) or start is None or stop is None:
            raise ValueError("Invalid limits given: %s" % str(limits))

        if start is S.NegativeInfinity and stop is S.Infinity:
            raise ValueError("Both the start and end value" " cannot be unbounded")

        limits = sympify((x, start, stop))

        if is_sequence(periodical, Tuple):
            periodical = sympify(tuple(flatten(periodical)))
        else:
            raise ValueError("invalid period %s should be something " "like e.g (1, 2) " % periodical)

        if Interval(limits[1], limits[2]) is S.EmptySet:
            return S.EmptySequence

        return Basic.__new__(cls, periodical, limits)

 def __new__(cls, function, limits):
     fun = sympify(function)
     if not is_sequence(fun) or len(fun) != 2:
         raise ValueError("Function argument should be (x(t), y(t)) but got %s" % str(function))
     if not is_sequence(limits) or len(limits) != 3:
         raise ValueError("Limit argument should be (t, tmin, tmax) but got %s" % str(limits))
     return GeometryEntity.__new__(cls, tuple(sympify(fun)), tuple(sympify(limits)))

    def __getitem__(self, key):
        """Return portion of self defined by key. If the key involves a slice
        then a list will be returned (if key is a single slice) or a matrix
        (if key was a tuple involving a slice).

        Examples
        ========

        >>> from sympy import Matrix, I
        >>> m = Matrix([
        ... [1, 2 + I],
        ... [3, 4    ]])

        If the key is a tuple that doesn't involve a slice then that element
        is returned:

        >>> m[1, 0]
        3

        When a tuple key involves a slice, a matrix is returned. Here, the
        first column is selected (all rows, column 0):

        >>> m[:, 0]
        Matrix([
        [1],
        [3]])

        If the slice is not a tuple then it selects from the underlying
        list of elements that are arranged in row order and a list is
        returned if a slice is involved:

        >>> m[0]
        1
        >>> m[::2]
        [1, 3]
        """
        if isinstance(key, tuple):
            i, j = key
            try:
                i, j = self.key2ij(key)
                return self._mat[i*self.cols + j]
            except (TypeError, IndexError):
                if isinstance(i, slice):
                    i = range(self.rows)[i]
                elif is_sequence(i):
                    pass
                else:
                    i = [i]
                if isinstance(j, slice):
                    j = range(self.cols)[j]
                elif is_sequence(j):
                    pass
                else:
                    j = [j]
                return self.extract(i, j)
        else:
            # row-wise decomposition of matrix
            if isinstance(key, slice):
                return self._mat[key]
            return self._mat[a2idx(key)]

def test_iterable_is_sequence():
    ordered = [list(), tuple(), Tuple(), Matrix([[]])]
    unordered = [set()]
    not_sympy_iterable = [{}, '', u('')]
    assert all(is_sequence(i) for i in ordered)
    assert all(not is_sequence(i) for i in unordered)
    assert all(iterable(i) for i in ordered + unordered)
    assert all(not iterable(i) for i in not_sympy_iterable)
    assert all(iterable(i, exclude=None) for i in not_sympy_iterable)

 def _eval_subs(self, old, new):
     from sympy.geometry.point import Point, Point3D
     if is_sequence(old) or is_sequence(new):
         if isinstance(self, Point3D):
             old = Point3D(old)
             new = Point3D(new)
         else:
             old = Point(old)
             new = Point(new)
         return  self._subs(old, new)

 def contains(self, o):
     """Return True if o is on this Line, or False otherwise."""
     if is_sequence(o):
         o = Point3D(o)
     if isinstance(o, Point3D):
         if self.arbitrary_point == 0:
             return True
         else:
             o = o.func(*[simplify(i) for i in o.args])
             eq = self.equation()
             a = []
             for i in range(3):
                 k = eq[i].subs(eq[i].free_symbols.pop(), o.args[i])
                 if k != nan:
                     a.append(k)
             if len(set(a)) == 1:
                 return True
             else:
                 return False
     elif not isinstance(o, LinearEntity3D):
         return False
     elif isinstance(o, Line3D):
         return self.__eq__(o)
     else:
         return o.p1 in self and o.p2 in self and o.p3 in self

def sequence(seq, limits=None):
    """Returns appropriate sequence object.
    If seq is a sympy sequence, returns SeqPer object
    otherwise returns SeqFormula object

    Examples
    ========

    >>> from sympy import sequence, SeqPer, SeqFormula
    >>> from sympy.abc import n

    >>> sequence(n**2, (n, 0, 5))
    SeqFormula(n**2, (n, 0, 5))

    >>> sequence((1, 2, 3), (n, 0, 5))
    SeqPer((1, 2, 3), (n, 0, 5))

    See Also
    ========

    sympy.series.sequences.SeqPer
    sympy.series.sequences.SeqFormula
    """
    seq = sympify(seq)

    if is_sequence(seq, Tuple):
        return SeqPer(seq, limits)
    else:
        return SeqFormula(seq, limits)

def mplot2d(f, var, show=True):
    """
    Plot a 2d function using matplotlib/Tk.
    """

    import warnings
    warnings.filterwarnings("ignore", "Could not match \S")
    
    p = import_module('pylab')
    if not p:
        sys.exit("Matplotlib is required to use mplot2d.")

    if not is_sequence(f):
        f = [f,]
    
    for f_i in f:
        x, y = sample2d(f_i, var)
        p.plot(x, y,'black')

    p.draw()
    
    p.ylabel("Transverse beam cordinate")
    p.xlabel("Propagation coordinate")
    p.grid(True)
    
    if show:
        p.show()

def flatten(iterable, levels=None, cls=None):
    """
    Recursively denest iterable containers.

    >>> from sympy.utilities.iterables import flatten

    >>> flatten([1, 2, 3])
    [1, 2, 3]
    >>> flatten([1, 2, [3]])
    [1, 2, 3]
    >>> flatten([1, [2, 3], [4, 5]])
    [1, 2, 3, 4, 5]
    >>> flatten([1.0, 2, (1, None)])
    [1.0, 2, 1, None]

    If you want to denest only a specified number of levels of
    nested containers, then set ``levels`` flag to the desired
    number of levels::

    >>> ls = [[(-2, -1), (1, 2)], [(0, 0)]]

    >>> flatten(ls, levels=1)
    [(-2, -1), (1, 2), (0, 0)]

    If cls argument is specified, it will only flatten instances of that
    class, for example:

    >>> from sympy.core import Basic
    >>> class MyOp(Basic):
    ...     pass
    ...
    >>> flatten([MyOp(1, MyOp(2, 3))], cls=MyOp)
    [1, 2, 3]

    adapted from http://kogs-www.informatik.uni-hamburg.de/~meine/python_tricks
    """
    if levels is not None:
        if not levels:
            return iterable
        elif levels > 0:
            levels -= 1
        else:
            raise ValueError("expected non-negative number of levels, got %s" % levels)

    if cls is None:
        reducible = lambda x: is_sequence(x, set)
    else:
        reducible = lambda x: isinstance(x, cls)

    result = []

    for el in iterable:
        if reducible(el):
            if hasattr(el, 'args'):
                el = el.args
            result.extend(flatten(el, levels=levels, cls=cls))
        else:
            result.append(el)

    return result

def hessian(f, varlist):
    """Compute Hessian matrix for a function f

    see: http://en.wikipedia.org/wiki/Hessian_matrix

    See Also
    ========

    sympy.matrices.mutable.Matrix.jacobian
    wronskian
    """
    # f is the expression representing a function f, return regular matrix
    if is_sequence(varlist):
        m = len(varlist)
        if not m:
            raise ShapeError("`len(varlist)` must not be zero.")
    elif isinstance(varlist, MatrixBase):
        m = varlist.cols
        if not m:
            raise ShapeError("`varlist.cols` must not be zero.")
        if varlist.rows != 1:
            raise ShapeError("`varlist` must be a row vector.")
    else:
        raise ValueError("Improper variable list in hessian function")
    if not getattr(f, 'diff'):
        # check differentiability
        raise ValueError("Function `f` (%s) is not differentiable" % f)
    out = zeros(m)
    for i in range(m):
        for j in range(i, m):
            out[i, j] = f.diff(varlist[i]).diff(varlist[j])
    for i in range(m):
        for j in range(i):
            out[i, j] = out[j, i]
    return out

    def distance(self, o):
        """
        Finds the shortest distance between the ray and a point.

        Raises
        ======

        NotImplementedError is raised if o is not a Point

        Examples
        ========

        >>> from sympy import Point, Ray
        >>> p1, p2 = Point(0, 0), Point(1, 1)
        >>> s = Ray(p1, p2)
        >>> s.distance(Point(-1, -1))
        sqrt(2)
        >>> s.distance((-1, 2))
        3*sqrt(2)/2
        """
        if not isinstance(o, Point):
            if is_sequence(o):
                o = Point(o)
        s = self.perpendicular_segment(o)
        if isinstance(s, Point):
            if self.contains(s):
                return S.Zero
        else:
            # since arg-order is arbitrary, find the non-o point
            non_o = s.p1 if s.p1 != o else s.p2
            if self.contains(non_o):
                return Line(self).distance(o)  # = s.length but simpler
        # the following applies when neither of the above apply
        return self.source.distance(o)

    def copyin_list(self, key, value):
        """Copy in elements from a list.

        Parameters
        ==========

        key : slice
            The section of this matrix to replace.
        value : iterable
            The iterable to copy values from.

        Examples
        ========

        >>> from sympy.matrices import eye
        >>> I = eye(3)
        >>> I[:2, 0] = [1, 2] # col
        >>> I
        [1, 0, 0]
        [2, 1, 0]
        [0, 0, 1]
        >>> I[1, :2] = [[3, 4]]
        >>> I
        [1, 0, 0]
        [3, 4, 0]
        [0, 0, 1]

        See Also
        ========

        copyin_matrix
        """
        if not is_sequence(value):
            raise TypeError("`value` must be an ordered iterable, not %s." % type(value))
        return self.copyin_matrix(key, Matrix(value))

    def contains(self, other):
        """
        Is the other GeometryEntity contained within this Segment?

        Examples
        ========

        >>> from sympy import Point3D, Segment3D
        >>> p1, p2 = Point3D(0, 1, 1), Point3D(3, 4, 5)
        >>> s = Segment3D(p1, p2)
        >>> s2 = Segment3D(p2, p1)
        >>> s.contains(s2)
        True
        """
        if is_sequence(other):
            other = Point3D(other)
        if isinstance(other, Segment3D):
            return other.p1 in self and other.p2 in self
        elif isinstance(other, Point3D):
            if Point3D.are_collinear(self.p1, self.p2, other):
                if other.distance(self.p1) + other.distance(self.p2) == self.length:
                    return True
                else:
                    return False
        return False

    def extract_leading_order(self, symbols, point=None):
        """
        Returns the leading term and it's order.

        Examples
        ========

        >>> from sympy.abc import x
        >>> (x + 1 + 1/x**5).extract_leading_order(x)
        ((x**(-5), O(x**(-5))),)
        >>> (1 + x).extract_leading_order(x)
        ((1, O(1)),)
        >>> (x + x**2).extract_leading_order(x)
        ((x, O(x)),)

        """
        lst = []
        symbols = list(symbols if is_sequence(symbols) else [symbols])
        if not point:
            point = [0]*len(symbols)
        seq = [(f, C.Order(f, *zip(symbols, point))) for f in self.args]
        for ef, of in seq:
            for e, o in lst:
                if o.contains(of) and o != of:
                    of = None
                    break
            if of is None:
                continue
            new_lst = [(ef, of)]
            for e, o in lst:
                if of.contains(o) and o != of:
                    continue
                new_lst.append((e, o))
            lst = new_lst
        return tuple(lst)

    def __new__(cls, label, shape=None, **kw_args):
        from sympy import MatrixBase, NDimArray

        if isinstance(label, string_types):
            label = Symbol(label)
        elif isinstance(label, Symbol):
            pass
        elif isinstance(label, (MatrixBase, NDimArray)):
            return label
        elif isinstance(label, Iterable):
            return _sympify(label)
        else:
            label = _sympify(label)

        if is_sequence(shape):
            shape = Tuple(*shape)
        elif shape is not None:
            shape = Tuple(shape)

        offset = kw_args.pop('offset', S.Zero)
        strides = kw_args.pop('strides', None)

        if shape is not None:
            obj = Expr.__new__(cls, label, shape)
        else:
            obj = Expr.__new__(cls, label)
        obj._shape = shape
        obj._offset = offset
        obj._strides = strides
        obj._name = str(label)
        return obj

def idiff(eq, y, x, n=1):
    """Return ``dy/dx`` assuming that ``eq == 0``.

    Parameters
    ==========

    y : the dependent variable or a list of dependent variables (with y first)
    x : the variable that the derivative is being taken with respect to
    n : the order of the derivative (default is 1)

    Examples
    ========

    >>> from sympy.abc import x, y, a
    >>> from sympy.geometry.util import idiff

    >>> circ = x**2 + y**2 - 4
    >>> idiff(circ, y, x)
    -x/y
    >>> idiff(circ, y, x, 2).simplify()
    -(x**2 + y**2)/y**3

    Here, ``a`` is assumed to be independent of ``x``:

    >>> idiff(x + a + y, y, x)
    -1

    Now the x-dependence of ``a`` is made explicit by listing ``a`` after
    ``y`` in a list.

    >>> idiff(x + a + y, [y, a], x)
    -Derivative(a, x) - 1

    See Also
    ========

    sympy.core.function.Derivative: represents unevaluated derivatives
    sympy.core.function.diff: explicitly differentiates wrt symbols

    """
    if is_sequence(y):
        dep = set(y)
        y = y[0]
    elif isinstance(y, Symbol):
        dep = set([y])
    else:
        raise ValueError("expecting x-dependent symbol(s) but got: %s" % y)

    f = dict([(s, Function(
        s.name)(x)) for s in eq.atoms(Symbol) if s != x and s in dep])
    dydx = Function(y.name)(x).diff(x)
    eq = eq.subs(f)
    derivs = {}
    for i in range(n):
        yp = solve(eq.diff(x), dydx)[0].subs(derivs)
        if i == n - 1:
            return yp.subs([(v, k) for k, v in f.items()])
        derivs[dydx] = yp
        eq = dydx - yp
        dydx = dydx.diff(x)

    def __new__(cls, term, *symbols, **assumptions):
        term = sympify(term)

        if term is S.NaN:
            return S.NaN

        if len(symbols) == 1:
            symbol = symbols[0]

            if isinstance(symbol, C.Equality):
                k = symbol.lhs
                a = symbol.rhs.start
                n = symbol.rhs.end
            elif is_sequence(symbol):
                k, a, n = symbol
            else:
                raise ValueError("Invalid arguments")

            k, a, n = map(sympify, (k, a, n))

        else:
            raise NotImplementedError

        obj = Expr.__new__(cls, **assumptions)
        obj._args = (term, Tuple(k, a, n))

        return obj