Replacing subexpression by a symbol?

Question

I have a 3x3 matrix, of which I compute inverse. The inverse can be written legibly only when some subexpressions are replaced by new symbols, because they appear multiple times. Can I have sympy try hard to find those subexpressions and replace them? I tried the following, without success:

from sympy import *

Ex, Ez, nuxy, nuxz = symbols('E_x E_z nu_xy nu_xz')

# compliance matrix for cross-anisotropic material
compl = Matrix([[1/Ex, -nuxy/Ex, -nuxz/Ez],
                [-nuxy/Ex, 1/Ex, -nuxz/Ez],
                [-nuxz/Ex, -nuxz/Ex, 1/Ez]])

# stiffness matrix
stiff = compl.inv()

# symbols I want to introduce 
m, e = symbols('m e')  

meSubs = {Ex/Ez: e, (1 - nuxy - 2*e*nuxz**2): m}  # instead of these subexpressions

# stiff.simplify() returns None, is that a bug? that's why I apply simplify together with subs here:
stiff.applyfunc(lambda x: simplify(x.subs(meSubs)))
print stiff

Using sympy 0.6.7 (I could upgrade, if needed).

EDIT:

I upgraded to 0.7.1-git (cf9c01f8f9b4b749a7f59891f546646e4b38e580 to be precise), and run (thanks to @PreludeAndFugue for suggestion):

from sympy import *
Ex,Ez,nuxy,nuxz,m=symbols('E_x E_z nu_xy nu_xz m')
compl=Matrix([[1/Ex,-nuxy/Ex,-nuxz/Ez],[-nuxy/Ex,1/Ex,-nuxz/Ez],[-nuxz/Ex,-nuxz/Ex,1/Ez]])
stiff=compl.inv()
stiff.simplify()
stiff.subs({-nuxy-2*nuxz**2+1:m})    # tried other rearrangements of the equality, as well, same result.
stiff.applyfunc(lambda x: together(expand(x)))
pprint(stiff)

obtaining

⎡              ⎛    2    ⎞                         ⎛            2⎞                              ⎤
⎢           Eₓ⋅⎝ν_xz  - 1⎠                     -Eₓ⋅⎝-ν_xy - ν_xz ⎠                 Eₓ⋅ν_xz      ⎥
⎢ ──────────────────────────────────   ────────────────────────────────────  ───────────────────⎥
⎢     2              2         2             2              2         2                    2    ⎥
⎢ ν_xy  + 2⋅ν_xy⋅ν_xz  + 2⋅ν_xz  - 1   - ν_xy  - 2⋅ν_xy⋅ν_xz  - 2⋅ν_xz  + 1  -ν_xy - 2⋅ν_xz  + 1⎥
⎢                                                                                               ⎥
⎢            ⎛            2⎞                         ⎛    2    ⎞                                ⎥
⎢        -Eₓ⋅⎝-ν_xy - ν_xz ⎠                      Eₓ⋅⎝ν_xz  - 1⎠                   Eₓ⋅ν_xz      ⎥
⎢────────────────────────────────────   ──────────────────────────────────   ───────────────────⎥
⎢      2              2         2           2              2         2                     2    ⎥
⎢- ν_xy  - 2⋅ν_xy⋅ν_xz  - 2⋅ν_xz  + 1   ν_xy  + 2⋅ν_xy⋅ν_xz  + 2⋅ν_xz  - 1   -ν_xy - 2⋅ν_xz  + 1⎥
⎢                                                                                               ⎥
⎢              E_z⋅ν_xz                              E_z⋅ν_xz                  E_z⋅(ν_xy - 1)   ⎥
⎢        ───────────────────                   ───────────────────           ────────────────── ⎥
⎢                      2                                     2                            2     ⎥
⎣        -ν_xy - 2⋅ν_xz  + 1                   -ν_xy - 2⋅ν_xz  + 1           ν_xy + 2⋅ν_xz  - 1 ⎦

Hm, so why does not get "-ν_xy - 2⋅ν_xz² + 1" replaced with m?

Answer 1

Both of your substitutions fail because the patterns do not exist in your matrix. Let's look at them one at a time:

1. There is no Ex/Ez ratio; x.subs(x/y, z) is unchanged
2. You have an extra e that appears in 1 - nuxy - 2*e*nuxz**2 so that expression is never matched, either.

@asmeurer showed that the literal 1 - nuxy - 2*nuxz**2 can be replaced, but it appears, too, as a factor in some denominators. The more complicated replacement can be done by checking if the patterns evenly divides an expression.

Let's make a function that will do that replacement:

>>> from sympy import *
>>> t = 1 - nuxy - 2*e*nuxz**2
>>> def do(x):
...    w, r = div(x, t)
...    if not r:
...        return m*w
...    return x

Now we will apply this to each element of the matrix:

>>> stiff.applyfunc(lambda x: factor_terms(bottom_up(x, do)))
Matrix([
[   -E_x*(nu_xz**2 - 1)/(m*(nu_xy + 1)), E_x*(nu_xy + nu_xz**2)/(m*(nu_xy + 1)),        E_x*nu_xz/m],
[E_x*(nu_xy + nu_xz**2)/(m*(nu_xy + 1)),    -E_x*(nu_xz**2 - 1)/(m*(nu_xy + 1)),        E_x*nu_xz/m],
[                           E_z*nu_xz/m,                            E_z*nu_xz/m, -E_z*(nu_xy - 1)/m]])

With complicated expressions/matrices, it is sometimes good to get an overview of the structure by using cse:

>>> cse(stiff)
([(x0, nu_xz**2), (x1, E_x*x0), (x2, 2*x0), (x3, x2 - 1), (x4, 1/(nu_xy**2 + nu_xy*x2 + x3)), (x5, x4*(-E_x + x1)), (x6, x4*(-E_x*nu_xy - x1)), (x7, 1/(nu_xy + x3)), (x8, nu_xz*x7), (x9, -E_x*x8), (x10, -E_z*x8)], [Matrix([
[ x5,  x6,                   x9],
[ x6,  x5,                   x9],
[x10, x10, x7*(E_z*nu_xy - E_z)]])])