Thanks for the link. I've been wondering how to lift N-D operations into the type system. If I write with NumPy

  x = r_[:n].reshape((3, -1, n/15))[:, :2].sum(axis=-1)

a compiler could figure out that the shape of x is (3, 2) and allow

  y = x + randn(m, 1, 2)

while forbidding

  z = y - r_[:7]

With NumPy you have to wait for a runtime error, and only if the shapes can't be broadcast. We're not even talking correct use of dimensions, etc.

I suppose this may have been tackled in Haskell in the Repa library, but I'd need to read papers like

http://benl.ouroborus.net/papers/repa/repa-icfp2010.pdf

to know. Hence my question if we could lift this into the type system in a conventional language, but I don't think it's possible in a general way.