Repository: sympy/sympy
File: sympy/tensor/array/expressions/arrayexpr_derivatives.py (Python)

Code before:
    diag_indices = [(b + i, b + c + i) for i in range(c)]
    return _array_diagonal(tp, *diag_indices)


@array_derive.register(MatrixExpr)
def _(expr: MatrixExpr, x: Expr):
    cg = convert_matrix_to_array(expr)
    if cg == expr:
        # Avoid infinite looping:
        raise NotImplementedError()
    return array_derive(cg, x)


@array_derive.register(MatPow)
def _(expr: MatPow, x: Expr):
    base = expr.base
    exponent = expr.exp
    dbase = array_derive(base, x)
    dexponent = array_derive(exponent, x)
    if not isinstance(dexponent, ZeroArray) or (dexponent == 0) == True:
        raise NotImplementedError()
    d = Dummy("d")
    tp = _array_tensor_product(base**d, dbase, base**(exponent-d-1))
    tc = _array_contraction(tp, (1, 4), (5, 6))
    pd = _permute_dims(tc, [1, 2, 0, 3])
    return ArraySum(pd,(d, 1, exponent))


@array_derive.register(HadamardProduct)
def _(expr: HadamardProduct, x: Expr):
    raise NotImplementedError()


@array_derive.register(ArrayContraction)
def _(expr: ArrayContraction, x: Expr):
    fd = array_derive(expr.expr, x)
    rank_x = len(get_shape(x))
    contraction_indices = expr.contraction_indices
    new_contraction_indices = [tuple(j + rank_x for j in i) for i in contraction_indices]
    return _array_contraction(fd, *new_contraction_indices)


@array_derive.register(ArrayDiagonal)
def _(expr: ArrayDiagonal, x: Expr):
    dsubexpr = array_derive(expr.expr, x)
    rank_x = len(get_shape(x))
    diag_indices = [[j + rank_x for j in i] for i in expr.diagonal_indices]
    return _array_diagonal(dsubexpr, *diag_indices)


@array_derive.register(ArrayAdd)

Review comment (suggestion):
This is still causing issues if `exponent = 0`, creating a summation from 1 to... zero!

As a consequence, this code is currently not working:
```python
expr = Determinant((X**j).subs(j, 0))
assert expr.diff(X) == ZeroMatrix(k, k)
```
adding a check for `exponent == 0` would be an easy fix, but this will probably hide the issue, as `exponent` could be a symbolic expression equivalent to zero.

Code after:
    diag_indices = [(b + i, b + c + i) for i in range(c)]
    return _array_diagonal(tp, *diag_indices)


@array_derive.register(MatrixExpr)
def _(expr: MatrixExpr, x: Expr):
    cg = convert_matrix_to_array(expr)
    if cg == expr:
        # Avoid infinite looping:
        raise NotImplementedError()
    return array_derive(cg, x)


@array_derive.register(MatPow)
def _(expr: MatPow, x: Expr):
    base = expr.base
    exponent = expr.exp
    dbase = array_derive(base, x)
    dexponent = array_derive(exponent, x)
    if not isinstance(dexponent, ZeroArray) or (dexponent == 0) == True:
        raise NotImplementedError()
    d = Dummy("d")
    tp = _array_tensor_product(base**d, dbase, base**(exponent-d-1))
    tc = _array_contraction(tp, (1, 4), (5, 6))
    pd = _permute_dims(tc, [1, 2, 0, 3])
    return ArraySum(pd,(d, 0, exponent-1))


@array_derive.register(HadamardProduct)
def _(expr: HadamardProduct, x: Expr):
    raise NotImplementedError()


@array_derive.register(ArrayContraction)
def _(expr: ArrayContraction, x: Expr):
    fd = array_derive(expr.expr, x)
    rank_x = len(get_shape(x))
    contraction_indices = expr.contraction_indices
    new_contraction_indices = [tuple(j + rank_x for j in i) for i in contraction_indices]
    return _array_contraction(fd, *new_contraction_indices)


@array_derive.register(ArrayDiagonal)
def _(expr: ArrayDiagonal, x: Expr):
    dsubexpr = array_derive(expr.expr, x)
    rank_x = len(get_shape(x))
    diag_indices = [[j + rank_x for j in i] for i in expr.diagonal_indices]
    return _array_diagonal(dsubexpr, *diag_indices)


@array_derive.register(ArrayAdd)