After talking to Pascal we decided that having this kind of funcion can be
quite useful, such that not every optimization does this splitting over and
over again ( it is really easy to have bugs by either messing up the order
or the count).

I think this commit got reverted ( I'm not sure when and why). Without it
graphs produced with pydotprint might be extremely misleading.

Small refactoring at the same time.

The formula before was wrong, it does not work when inputs are matrices,
because it depends on which axis y is computed.  I believe this formulation
is much better.

Before only the gradient of inc_subtensor was ok ( and it also was the only
one checked)

numpy.all doesn't check dtype and shape

Rationale: its name shadowed the scipy package, and it was not successful in
simplifying other modules' code that had to do scipy version checks. (Which is
only two places anyway).