Inconsistent behavior of transformations.scale_matrix

After #959 I played around with lib/transformations trying to find ways to improve coverage, when I spotted the following inconsistency between the cython and python versions of scale_matrix if a non-default direction is used:

#Cython version:
MDAnalysis.lib.transformations.scale_matrix(factor=1.1, direction=[1.,2.,0])
array([[ 1.1,  0.2,  0. ,  0. ],
       [ 0.2,  1.4,  0. ,  0. ],
       [ 0. ,  0. ,  1. ,  0. ],
       [ 0. ,  0. ,  0. ,  1. ]])

#Python version:
MDAnalysis.lib.transformations._py_scale_matrix(factor=1.1, direction=[1.,2.,0])
array([[ 1.02,  0.04,  0.  ,  0.  ],
       [ 0.04,  1.08,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  1.  ,  0.  ],
       [ 0.  ,  0.  ,  0.  ,  1.  ]])

Things get worse when going for negative factor values:

#Cython version:
MDAnalysis.lib.transformations.scale_matrix(factor=-1.1, direction=[1.,2.,0])
array([[-1.1, -4.2, -0. ,  0. ],
       [-4.2, -7.4, -0. ,  0. ],
       [-0. , -0. ,  1. ,  0. ],
       [ 0. ,  0. ,  0. ,  1. ]])

#Python version:
MDAnalysis.lib.transformations._py_scale_matrix(factor=-1.1, direction=[1.,2.,0])
array([[ 0.58, -0.84,  0.  ,  0.  ],
       [-0.84, -0.68,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  1.  ,  0.  ],
       [ 0.  ,  0.  ,  0.  ,  1.  ]])

Not only are the cython and python results incompatible with one another, they are also incompatible with their negative-factor counterparts. In addition, these are scaling, not shearing operations. Why are we getting off-diagonal elements?

I think part of the problem lies on using both factor and direction, when a 3-tuple of factors would suffice for directional scaling (and it'd also make it easier to visualize the meaning of the scaling).

Finally, there is no usage of the scaling functions in our internal code, which is good, but I guess this must be dealt with or it'll bite us at some point in the future.

MDAnalysis / mdanalysis

Inconsistent behavior of transformations.scale_matrix #962