Using sympy to sum over a range of symbols

Consider two sets i,j which both have m elements. Say we have an expression which describes a sum of terms. Each term can be described as a product of an element of i and j. Now, I would like to sum over each element of j, where each element has the range [i1,i2,...,im].

In the context of python & sympy, this is difficult since sympy's Sum describes the summation variable with (symbol,start,stop), which assume integer steps.

To demonstrate what I mean, consider the following code:

>>> from sympy import *
>>> i = symbols('i1,i2,i3,i4') # for the case m = 4
>>> j = symbols('j1,j2,j3,j4')

Here I use permutations to setup the expression:

>>> from itertools import permutations as perm
>>> c = list(perm(range(4),2))
>>> a,b = c[0]
>>> expr = i[a]*j[b]
>>> for a,b in c[1:]:
>>>     expr += i[a]*j[b]
>>> print(expr)
i1*j2 + i1*j3 + i1*j4 + i2*j1 + i2*j3 + i2*j4 + i3*j1 + i3*j2 + i3*j4 + i4*j1 + i4*j2 + i4*j3

Now, using Sum over each j with range of i. It would be ideal if I could write one of the following:

>>> s = Sum(expr,(j,i))
>>> s = Sum(expr,(j1,i),(j2,i),...,(jm,i))

But that's not canonical with the sympy documentation. Are there any other methods which can be used to solve this problem?

Edit:
In this post, I tried to isolate the problem by only using elements i,j in expr. The full context problem is where expr is a sum of Kronecker Delta functions of i,j and using a sum over index set j, where each element of j has range i. For example:

>>> from sympy import KroneckerDelta as KD
>>> expr = KD(i[0],j[1]) # Only doing j[1] to reduce clutter
>>> print(expr)
KroneckerDelta(i1,j2)
>>> s = Sum(expr,(j[1],i)).doit() 
>>> print(s)
# Desired output to look like:
1 + KroneckerDelta(i1,i2) + KroneckerDelta(i1,i3) + KroneckerDelta(i1,i4)

This is the reason for which I phrased my question as: summing over each element of j with range i.

IndexedBase can act as a symbolic, integer-indexed array:

>>> from sympy import *
>>> from sympy.abc import k,l
>>> i,j = map(IndexedBase,'ij')
>>> Sum(i[k]*j[l],(k,1,2),(l,1,2)).doit().expand()
i[1]*j[1] + i[1]*j[2] + i[2]*j[1] + i[2]*j[2]
>>> Sum(Piecewise((i[k]*j[l],Ne(k,l)),(0,True)),(k,1,2),(l,1,2)).doit()
i[1]*j[2] + i[2]*j[1]

It's not clear whether you want the cross terms when the indicices are the same, so both versions are shown.