I’d like to re-arrange a complex expression, a simplified form can be expressed as:

i, j, N = symbols("i j N", integers=True, finite=True, nonnegative=True)
a, b, c, d, e = symbols("a b c d e", cls=IndexedBase)

expr = (Sum(a[i] * Sum(b[i] * d[j], (j, 1, N)) * c[i], (i, 1, N)) 
        + Sum(a[i] * Sum(b[i] * e[j], (j, 1, N)) * c[i], (i, 1, N)))
expr

Ideally, I’d like to convert to something like the following:

I’ve not used replace much before but it seemed like it might be useful for this type of manipulation. Unfortunately, my attempt doesn’t properly capture the dummy variables:

ws1, ws2, ws3, wd1, wd2, wl1 = symbols("ws1 ws2 ws3 wd1 wd2 wl1", cls=Wild)

expr.replace(
    Sum(ws1 * Sum(ws3, (wd1, 1, wl1)), (wd2, 1, wl1)) * ws2,  
    Sum(ws1 * ws2 * ws3, (wd1, 1, wl1), (wd2, 1, wl1))
)

Wondering if anyone could offer any suggestions?