Commits

Michael Forbes committed 633be89

Added a SympyODE class:
- Base class works, but transformations not working yet.
Cleaned up some code to satisfy pep8 and other checkers (incomplete).
Docs:
- Use common acronyms in mmfgls.sty

  • Participants
  • Parent commits 3da7b68

Comments (0)

Files changed (5)

File doc/pweave/preamble.tex

 \accentedsymbol{\oV}{\op{V}}
 \accentedsymbol{\psisq}{\abs{\Psi}^2}
 
-% Acryonyms
-\newacronym{HFB}{hfb}{Hartree-Fock-Bogoliubov}
-\newacronym{BdG}{bdg}{Bogoliubov-de Gennes}
-\newacronym{QMC}{qmc}{Quantum Monte Carlo}
-\newacronym{DFT}{dft}{Density Functional Theory}
-\newacronym{KS}{ks}{Kohn-Sham}
-\newacronym{TDDFT}{tddft}{Time-Dependent Density Functional Theory}
-\newacronym{PDE}{dft}{Partial Differential Equation}
-\newacronym{SEQ}{seq}{Schr\"odinger Equation}
-\newacronym{GPE}{gpe}{Gross-Pitaevskii Equation}
-\newacronym{SLDA}{slda}{Superfluid Local Density Approximation}
-\newacronym{GP}{gp}{Gross-Pitaevskii}
-\newacronym{TF}{tf}{Thomas-Fermi}
-\newacronym{ETF}{etf}{Extended Thomas-Fermi}
-\newacronym{BEC}{bec}{Bose-Einstein Condensate}
-\newacronym{BCS}{bcs}{Bardeen-Cooper-Schrieffer}
-\newacronym{RHS}{rhs}{right-hand side}
-\newacronym{FiPy}{FiPy}{FiPy}
-\newacronym{FFT}{fft}{Fast Fourier Transform}
-\newacronym{LBFGS}{l-bfgs}{Limited-memory Broyden-Fletcher-Goldfarb-Shanno}
-
-\newacronym{LO}{lo}{Leading Order}
-\newacronym{NLO}{nlo}{Next-to Leading Order}
-\newacronym{NNLO}{nnlo}{Next-to-Next-to Leading Order}
+\RequirePackage{mmfgls}
 
 \DeclareMathOperator{\FFT}{FFT}
 \DeclareMathOperator{\IFFT}{IFFT}

File mmf/math/ode/deq.py

 ==========
  Overview
 ==========
+
 Here we provide a simplified interface to the `colnew` boundary value
 solver.  This can solve differential equations for a set of functions
 $Y_i(x)$ of the form
     Z = [&[Y_1, Y_1', \ldots, Y_1^{(m_1-2)}, Y_1^{(m_1-1)}],\\
         &[Y_2, Y_2', \ldots, Y_2^{(m_2-1)}],\\
         &[\ldots],\\
-        &[Y_n, Y_n', \ldots, Y_n^{(m_n-3)}, Y_n^{(m_n-2)}, Y_n^{(m_n-2)}]]
+        &[Y_n, Y_n', \ldots, Y_n^{(m_n-3)}, Y_n^{(m_n-2)}, Y_n^{(m_n-1)}]]
 
-is the set of derivatives $Y^{(m)} = \d{Y}(x)/\d{x}^m$. Thus, there are $n$
+is the set of derivatives $Y^{(m)} = \d^{m}{Y}(x)/\d{x}^m$. Thus, there are $n$
 differential equations of orders $m_1$, $m_2$, $\ldots$, $m_n`$.
 
-.. note:: The order $m$ may also be zero $m=0$ in which case the equation
-   represents a constraint to be solved.  Likewise, if one of the constraints is
-   that an integrated quantity has a fixed value (total particle number for
-   example), then you can introduce $N(x)$ through the first order 
+.. note:: Constants may be implemented using $m=1$ equations with $Y'= 0$ in
+   which case $Y$ is a constant to be found.  Likewise, if one of the
+   constraints is that an integrated quantity has a fixed value (total particle
+   number for example), then you can introduce $N(x)$ through the first order 
    equation $N'(x) = n(x)$ and give the boundary conditions $N(0)=0$
    and $N(x_{\text{max}}) = N$.  These two techniques are often used together to
    implement the constraint $N$ with a Lagrange multiplier $\mu$ that is entered
-   as a parameter with a zeroth-order equation.
-
-   A first-order equation with $N(0)=0$ can also be used to simply
-   integrate $N(x)$ using the adaptive mesh.
+   as a constant parameter.
 
 The `m_star` = `m_1` + `m_2` + ... `m_n` boundary conditions are specified by
 the conditions
 .. math::
    g_j(y_j, Z) = 0
 
+at boundary points $a_l \leq x_1 \leq x_2 \leq \cdots x_{m^*} \leq a_r$.  (In
+these formula, the abscissa $x_j$ is implicit.)  The solution will be tabulated
+on the range $x\in [a_l, a_r]$.
+
 To simplify the interface, it is assumed that the user will give names
 to the `n` variables `Y_n` and define the functions `f_(Y_n)`.
 
    :exc:`KeyboardInterrupt`... You might have to press Ctr-C many times while
    plotting...
 
+Transformations
+===============
+Any system can be reduced to a system of first-order equations by introducing
+auxiliary variables for the higher orders.  I.e. for $u''(x) = f(x, u, u'(x))$
+one can introduce $v(x) = u'(x)$ to express this as the system
+
+.. math::
+   v'(x) = f(x, u, v)\\
+   u'(x) = v.
+
+Let us combine all variables into a single vector $x_i$ so that the system
+becomes (the first equation is trivial $x_0 \equiv x$ so $x_0' = 1$)
+
+.. math::
+   \d{x_i} = f_i(\vect{x})\d{x}
+
+Now express these in terms of new variables as functions $x_i(y_j)$.  The
+resulting differentials satisfy
+
+.. math::
+   \d{x_i} = \sum_j \pdiff{x_{i}}{y_{j}}\d{y_j} = \sum_j J_{ij}y_j
+
+where $J_{ij}$ is the Jacobian of the transformation.  We may thus write
+
+.. math::
+   \d{y_i} = \sum_j J^{-1}_{ij} \d{x_j} = \sum_j J^{-1}_{ij} f_j\d{x}
+
+Dividing with $\d{y_0}$, we obtain
+
+.. math::
+   y'_a = \frac{\sum_j J^{-1}_{aj} f_j}{\sum_j J^{-1}_{0j} f_j}.
+
+Everything on the right-hand-side can be expressed in terms of the new
+variables. 
 """
 from __future__ import division
-__all__ = ['IDict', 'Solution', 'Transform', 'log_transform', 'ODE', 'ODE_1_2']
+__all__ = ['SympyODE', 'IDict', 'Solution', 'Transform', 'log_transform', 'ODE',
+           'ODE_1_2', 'bvp']
+
 import collections
 import time
+import inspect
 from warnings import warn
 
 import numpy as np
 from numpy import pi, sqrt, arctanh, tanh, array
 
 from mmf.objects import StateVars, Container, process_vars
-from mmf.objects import Required, Excluded, Computed, ClassVar
+from mmf.objects import Required, Delegate, Excluded, Fast, Computed, ClassVar
 
 import scikits.bvp1lg.jacobian
 import scikits.bvp1lg as bvp
 
 _tol = 1e-8
+
 class IDict(dict):
     def __init__(self, *v, **kw):
         dict.__init__(self, *v, **kw)
         for c in self.callbacks:
             c(self)
 
+
+class SympySolution(StateVars, bvp.colnew.Solution):
+    """Represents the solution to an ODE as returned by colnew.  The
+    default version simply ensures that copies are made of the
+    abscissa because the colnew solution will mutate these."""
+    _state_vars = [
+        ('ispace', Required,
+         "Integer info from colnew about solution."),
+        ('fspace', Required,
+         "Float info from colnew about solution."),
+        ('ncomp', Computed,
+         "Number of equations."),
+        ('mstar', Computed,
+         "Number of degrees of freedom."),
+        ('nmesh', Computed,
+         "Number of mesh points."),
+        ('ode', Required, "ODE object"), 
+        ('_names', Computed)
+        ]
+    process_vars()
+    def __init__(self, *varargin, **kwargs):
+        bvp.colnew.Solution.__init__(self,
+                                     ispace=self.ispace,
+                                     fspace=self.fspace)
+        self._names = set([_k[0] for _k in self.ode._dispatch])
+        
+    def __getattr__(self, name):
+        if name in self.__dict__.get('_names', []):
+            def f(x=None, d=0, _dispatch=self.ode._dispatch):
+                return self(x)[_dispatch[(name, d)]]
+            f.__name__ = name
+            return f
+        else:
+            raise StateVars.__getattr__(self, name)
+            
+    @property
+    def x(self):
+        return self.mesh
+        
+    def __call__(self, x=None):
+        """Evaluation the solution at given points.
+
+        :returns:
+            The solution vector, as::
+            
+                [u_1(x), u_1'(x), ..., u_1^{m_1 - 1}(x), u_2(x), ...
+                 u_{ncomp}^{m_{ncomp} - 1}]
+
+            broadcast to ``x``. Shape of the returned array
+            is x.shape + (mstar, ).
+
+        Also makes sure that x does not get mutated by making a
+        copy, unlike the default version.
+        """
+        if x is None:
+            x = self.x
+        return bvp.colnew.Solution.__call__(self, np.array(x)).T
+
+    def plot(self, *varagin, **kwargs):
+        """Plot the solution"""
+        x = self.x
+        import matplotlib.pylab as plt
+        plt.plot(x, self(x).T, *varagin, **kwargs)
+
+
+class SympyODE(StateVars):
+    r"""Represent an ODE.
+
+    Specify the original ODE here in terms of the original variables either as a
+    string or using :pkg:`sympy`.  The :pkg:`sympy` package will be used to
+    define the various functions required to satisfy the BVP interface.
+
+    Examples
+    --------
+    >>> class Example1(SympyODE):
+    ...     _state_vars = [
+    ...         ('vars', ['x', 'u', 'v']),
+    ...         ('eqs', 
+    ...          ['-u(x).diff(x)/x + ((nu/x)**2 -1)*u(x) - u(x).diff(x,x)',
+    ...           'x**(nu+1)*u(x) - v(x).diff(x)']),
+    ...         ('bcs', 
+    ...          [(1, 1e-5, 'u(x) - besselj(nu,x)'),
+    ...           (10, 0, 'u(x) - besselj(nu,x)'),
+    ...           (5, 1e-5, 'v(x) - x**(nu+1)*besselj(nu+1,x)')]),
+    ...         ('nu', 3.4123)]
+    ...     process_vars()
+    >>> e = Example1()
+    >>> s = e.solve(maximum_mesh_size=300)
+    """
+    _state_vars = [
+        ('vars', Required, \
+             r"""List of `sympy` vars `[x, Y_1, Y_2, ..., Y_n]`, the first of
+             which is the abscissa, and the rest correspond to the variables
+             `Y_i`."""),
+        ('eqs', Required, \
+             r"""List of equations defining the system. Each equation should be
+             zero -- i.e. the lhs of $Y^{(m)} - f(x, Z) = 0$ -- and linear in
+             the highest derivative."""),
+        ('bcs', Required, \
+             r"""List of tuples `[(x, g)]` where `xs = [x_1, x_2, ...]` is
+             the set of abscissa for the boundary conditions, and 
+             `gs = [g_1(x,Z), g_2(x,Z), ...]` is the list of equations.  Note:
+             the boundary points `xs` must evaluate to numbers (they cannot
+             contain parameters), the equations `gs` should explicitly contain
+             the abscissa (to allow later for transformations), but the
+             equations `gs` are not permitted to have the highest 
+             derivative.  The `xs` need not be sorted."""),
+        ('params', {}, \
+             r"""Dictionary of parameter values.  The equations can be expressed
+             in terms of additional constants.  One can specify values here, or
+             define them directly as state variables.  (This dictionary allows
+             you to avoid possible name-clashes.)  These values will be
+             preferred if a parameter appears both here and as a state var."""),
+        ]
+    process_vars()
+
+class TransformSympyODE(SympyODE):
+    _state_vars = [
+        ('ode', Required),
+        ('transform', Required,\
+             r"""List of expressions for the old variables in terms of the new
+             ones. The first element must be the old abscissa which must be a
+             simple function of the new abscissa, and the"""),
+        
+        ('params=ode.params'),        
+        ('vars', Computed),
+        ('eqs', Computed),
+        ('bcs', Computed),
+        ]        
+    process_vars()
+    def __init__(self, *v, **kw):
+        pass
+        
+class SympyBVP(StateVars):
+    _state_vars = [
+        ('ode', Required, "Instance of :class:`SympyODE`"),
+
+        ('left', Fast(None), "Left endpoint = `min(xs)`"),
+        ('right', Fast(None), "Right endpoint = `max(xs)`"),
+        ('tolerances', Fast(None), "Tolerances at each `xs`"),
+
+        # Optional class variables
+        ('tol', 1e-5, "Default tolerance."),
+        ('is_linear', ClassVar(False),\
+             "Set to True if the problem is linear"),
+        ('collocation_points', ClassVar(None),\
+             "List of explicit collocation points"),
+        ('extra_fixed_points', ClassVar(None),\
+             "List of explicit fixed points."),
+        ('Solution', SympySolution, "Should inherit from :cls:`SympySolution`"),
+        
+        # Instance variables
+        ('maximum_mesh_size', 1000, "Maximum size for dynamic mesh"),
+        ('problem_regularity', 0,\
+             "Set to 1 if the problem is particularly sensitive"),
+        ('verbosity', 0, ""),
+        ('coarsen_initial_guess_mesh', True,\
+             "If True, allow mesh to be coarsened."),
+
+        ('idict', Excluded(IDict()),\
+             "Dictionary in which to store intermediate results "\
+             "(I.e. for plotting.)"),
+        ('plot', False, "If True, then plot intermediate results."),
+        ('plot_pause', 0, "Number of seconds to wait after plotting."),
+        ('_debug', False,\
+             "If True, then some debugging information is printed during "\
+             "the iteration"),
+        ('_tols', Computed),
+        
+        ('bvp_solve_args', 
+         ClassVar(inspect.getargs(bvp.colnew.solve.func_code).args),
+         """Argument names for :func:`scikits.bvp1lg.colnew.solve`."""),
+
+        
+        # Computed parameters for bvp.colnew.solve
+        ('boundary_points', Computed),
+        ('degrees', Computed, r"$m_i$"),
+        ('vectorized', ClassVar(True), "We generate vectorized functions"),
+        
+        ('_f', Computed), ('_df', Computed), 
+        ('_g', Computed), ('_dg', Computed),
+        ('_dispatch', Computed, r"Dispatch table for solution"),
+        
+        # Excluded temporary variables.
+        ('m_star', Excluded, "Total number of 'degrees of freedom'.\n"\
+             "Computed in solve from vars and used in some methods"),
+        ('sol', Excluded, "Temporary storage for current solution."),
+        ('ispace_fspace', Excluded,\
+             "Temporary storage for the working arrays"),
+        ]
+    process_vars()
+    def __init__(self, *v, **kw):
+        ode = self.ode
+        env = dict(ode.__dict__, **ode.params)
+
+        x = sympy.S(ode.vars[0])
+        ys = map(sympy.S, ode.vars[1:])
+        ncomp = len(ys)
+        assert ncomp == len(ode.eqs)
+        ms = []
+        fs = []
+        zs = []
+        z_dict = {}
+        mstar = 0
+        for _i, _eq in enumerate(map(sympy.S, ode.eqs)):
+            _y = ys[_i]
+            m = sympy.ode.ode_order(_eq, _y(x))
+            _dy = _y(x).diff(x, m)
+            ms.append(m)
+            _f = sympy.solve(_eq, _dy)
+            if not len(_f) == 1:
+                raise ValueError(
+                    "Equation %i must be linear in '%s': got '%s'" %
+                    (_i, str(_dy), str(_eq)))
+            fs.append(_f[0])
+            for _j in xrange(m):
+                _z = sympy.symbols('_z%i' % (mstar,))
+                z_dict[_y(x).diff(x, _j)] = _z
+                zs.append(_z)
+                mstar += 1
+        
+        xs, gs = zip(*ode.bcs)
+        xs = [sympy.S(_x) for _x in xs]
+        gs = [sympy.S(_g) for _g in gs]
+
+        symbols = set()
+        for _vs in [xs, fs, gs]:
+            for _v in _vs:
+                symbols.update(_v.free_symbols)
+        symbols.difference_update([x])
+
+        if symbols.intersection(zs):
+            raise ValueError(
+                "Variable name class.  Do not use %s" %
+                str(tuple(symbols.intersection(zs))))
+
+        env = {}
+        for _v in map(str, symbols):
+            if _v in ode.params:
+                env[_v] = ode.params[_v]
+            elif hasattr(self, _v):
+                env[_v] = getattr(self, _v)
+            else:
+                raise ValueError(
+                    "Parameter '%s' not defined as an attribute or in 'params'"
+                    % (_v,))
+        
+        xs = [float(_x.subs(env)) for _x in xs]
+        _gs = gs
+        gs = [_g.subs(z_dict).subs(dict(x=_x)) for _x, _g in zip(xs, gs)]
+        xs, gs = zip(*sorted(zip(xs, gs)))
+
+        F = sympy.matrices.Matrix(fs).subs(env).subs(z_dict)
+        G = sympy.matrices.Matrix(gs).subs(env).subs(z_dict)
+        dF = F.jacobian(zs)
+        dG = G.jacobian(zs)
+        F = F.T.tolist()[0]
+        G = G.T.tolist()[0]
+        _zero = sympy.symbols('_zero')
+        for _i in xrange(ncomp):
+            if F[_i].is_number:
+                F[_i] += _zero
+            for _j in xrange(mstar):
+                if dF[_i, _j].is_number:
+                    dF[_i, _j] += _zero
+
+        for _i in xrange(mstar):
+            if G[_i].is_number:
+                G[_i] += _zero
+
+        self._f = sympy.lambdify([_zero, x] + zs, list(F))
+        self._df = sympy.lambdify([_zero, x] + zs, dF.tolist())
+        self._g = sympy.lambdify([_zero] + zs, G)
+        self._dg = _dg=sympy.lambdify(zs, dG)
+
+        dispatch = {}
+        _i = 0
+        for _m, _y in zip(ms, ys):
+            for _d in xrange(_m):
+                dispatch[(str(_y), _d)] = _i
+                _i +=1
+        self._dispatch = dispatch
+        self.boundary_points = xs
+        self.degrees = ms
+
+    def fsub(self, x, z):
+        return np.asarray(self._f(np.zeros(x.shape, x.dtype), x, *z))
+
+    def dfsub(self, x, z):
+        return np.asarray(self._df(np.zeros(x.shape, x.dtype), x, *z))
+
+    def gsub(self, z):
+        return np.diag(np.asarray(self._g(np.zeros(z.shape, z.dtype), *z)))
+
+    def dgsub(self, z): 
+        return np.asarray([np.asarray(self._dg(*_z))[:,_i] 
+                           for _i, _z in enumerate(z.T)]).T
+
+    def solve(self, **kw):
+        args = dict([(_n, getattr(self, _n)) for _n in self.bvp_solve_args
+                     if hasattr(self, _n)])
+
+        if self.tolerances is None:
+            tolerances = (self.tol,)*len(self.boundary_conditions)
+        else:
+            tolerances = self.tolerances
+
+        args.update(dict(tolerances=tolerances,
+                         left=self.left,
+                         right=self.right, **kw))
+        sol = bvp.colnew.solve(**args)
+        sol = self.check_sol(sol)
+        return SympySolution(ode=self, ispace=sol.ispace, fspace=sol.fspace)
+
+    def check_sol(self, sol):
+        """This function should check the solution sol, possibly
+        correct it and return the resulting solution."""
+        return sol
+
+class Example1_ODE(SympyODE):
+    _state_vars = [
+        ('vars', ['x', 'u', 'v']),
+        ('eqs', 
+         ['-u(x).diff(x)/x + ((nu/x)**2 -1)*u(x) - u(x).diff(x,x)',
+          'x**(nu+1)*u(x) - v(x).diff(x)']),
+        ('bcs', 
+         [(1, 'u(x) - besselj(nu,x)'),
+          (10, 'u(x) - besselj(nu,x)'),
+          (5, 'v(x) - x**(nu+1)*besselj(nu+1,x)')]),
+        ('nu', 3.4123), 
+        ]
+    process_vars()
+
+class Example1_BVP(SympyBVP):
+    _state_vars = [
+        ('ode', Delegate(Example1_ODE, ['nu'])),
+        ('tolerances', [1e-5, 0, 1e-5]),
+        ]
+    process_vars()
+
+
+#e = Example1_BVP()
+#s = e.solve(maximum_mesh_size=300)
+        
+try:
+    import sympy
+except:
+    del SympyODE
+
+
 class Solution(StateVars, bvp.colnew.Solution):
     """Represents the solution to an ODE as returned by colnew.  The
     default version simply ensures that copies are made of the

File mmf/objects/_objects.py

    does this, making sure that the delegates are copied while the
    other attributes are simply referenced.  Method
    :meth:`__deepcopy__` works as usual.
-        
+
 2) The method :meth:`items` returns a list of the data
    required to construct the object.  The object should be
    able constructed by passing all of these items to the
 __all__ = [
     'InitializationError', 'NameClashWarning', 'NameClashWarning1',
     'NameClashWarning2',
-    'Archivable', 'StateVars', 'Container', 
-    'process_vars', 
+    'Archivable', 'StateVars', 'Container',
+    'process_vars',
     'Attr', 'Required', 'Excluded', 'Internal', 'Fast', 'Computed',
     'Deleted', 'ClassVar', 'NoCopy', 'Delegate', 'Ref',
-    'attribute', 'parameter', 'calculate', 
-    'MemoizedContainer', 
-    'option', 'Options', 'class_NamedArray', 'recview', 
+    'attribute', 'parameter', 'calculate',
+    'MemoizedContainer',
+    'option', 'Options', 'class_NamedArray', 'recview',
     'RecView', 'Calculator']
 import sys
 import warnings
 
 try:
     from zope.interface.interface import Element as ZopeElement
-except ImportError:                     #pragma: no cover
+except ImportError:                     # pragma: no cover
     from types import NoneType as ZopeElement
 
 _ARCHIVE_CHECK = True
 __warningregistry__ = {}
 
+
 ###########################################################
 # Classes
 class Archivable(object):
     # interfaces.implements(interfaces.IArchivable)
     # Breaks help.  See zope Bug #181371.
     # https://bugs.launchpad.net/zope3/+bug/181371
-    
-    def items(self):        # pragma: no cover 
+
+    def items(self):        # pragma: no cover
         r"""Return a list `[(name, obj)]` of `(name, object)` pairs
         where the instance can be constructed as
         `ClassName(name1=obj1, name2=obj2, ...)`.
         b
         """
         return (k for (k, v) in self.items())
-    
+
     def archive_1(self, env=None):      # pylint: disable-msg=W0613
         r"""Return (rep, args, imports).
-        
+
         Defines a representation rep of the instance self where the
         instance can be reconstructed from the string rep evaluated in
         the context of dict args with the specified imports = list of
            o.a
         """
         arch = mmf.archive.Archive()
-        arch.insert(**{name:self})
+        arch.insert(**{name: self})
         return str(arch)
 
     def __repr__(self):
     def __str__(self):
         return self.__repr__()
 
+
 class InitializationError(Exception):
     r"""Initialization failed or incomplete."""
 
+
 class NameClashWarning(UserWarning):
     r"""Name clash.  Name already defined.
 
         __warningregistry__ = {}
         warnings.simplefilter(action, cls)
 
+
 class NameClashWarning1(NameClashWarning):
     r"""Name clash.  Name already defined.
 
     """
     default_action = 'ignore'
 
+
 class NameClashWarning2(NameClashWarning):
     r"""Name clash.  Name already defined.
 
 NameClashWarning1.simplefilter()
 NameClashWarning2.simplefilter()
 
+
 class _Required(object):                # pylint: disable-msg=R0903
     r"""Default value for required attributes.
-    
+
     Examples
     --------
     >>> Required
     def __init__(self, doc=NotImplemented):
         if doc is not NotImplemented:
             self.doc = doc
-            
+
     def __repr__(self):
         r"""Simple string for describing arguments etc.
         """
             return "Required(%r)" % (doc, )
         else:
             return "Required"
-        
+
     def __call__(self, doc):
         r"""Return a required object with documentation."""
         return _Required(doc=doc)
 
 Required = _Required()                  # pylint: disable-msg=C0103
 
+
 class _Deleted(Archivable):
     r"""Deleted variable.
 
     >>> 'a' in c._vars
     False
     """
-    
+
 Deleted = _Deleted()                    # pylint: disable-msg=C0103
 
+
 def _is(default, type_):
     r"""Return true if default is a subclass or instance of `type`."""
     try:
     except TypeError:
         return isinstance(default, type_)
 
+
 class Attr(Archivable):
     r"""Base class for all attribute types."""
 
         names are determined by inspecting the arguments to the
         constructor :meth:`__init__`."""
         items = []
-        func_code = self.__init__.im_func.func_code # pylint: disable-msg=E1101
+        func_code = self.__init__.im_func.func_code  # pylint: disable-msg=E1101
         nargs = func_code.co_argcount
 
         # pylint: disable-msg=E1101
             elif value != defaults[name - (nargs - 1 - len(defaults))]:
                 items.append((var, value))
         return items
-            
+
+
 class HasDefault(Attr):
     r"""Base class for types with default values.
 
     """
     value = NotImplemented              # Here so that one can always
                                         # access value, even on a class.
+
     def __init__(self, value):          # pylint: disable-msg=W0231
         self.value = value
 
+
 class Excluded(HasDefault):
     r"""Excluded variable.
 
     indicate that a particular method is executing, or for
     non-essential reporting of internal state/temporary data.
     """
+
+
 class Internal(Excluded):
     r"""Internal variable.
 
     Functionally equivalent to :class:`Excluded` attributes and are
     presently implemented as these (and reported as these).
     """
+
+
 class Fast(HasDefault):
     r"""Fast variable.
 
     Make sure that nothing depends on the value though, otherwise the object
     might not be properly updated."""
 
+
 class Computed(Attr):
     r"""Computed attributes.
 
         save = kw.pop('save', False)
         if kw:
             raise TypeError(
-                "Computed() got an unexpected keyword argument '%s'" 
+                "Computed() got an unexpected keyword argument '%s'"
                 % kw.keys()[0])
         elif v:
             raise TypeError(
     def items(self):
         r"Our constructor breaks the auto finding of attributes"
         return [('save', self.save)]
-    
+
+
 class ClassVar(HasDefault):
     r"""Class variables (not associated with instances).
 
     called.
     """
 
+
 class _Attr(HasDefault):
     r"""Denotes an attribute.
 
         self.doc = doc
         self.value = value
 
+
 class Delegate(HasDefault):
     r"""Delegate Class.
 
     --------
     This is the simplest delegation to `A`.  All variables in
     `A._state_vars` become members of `D1`:
-    
+
     >>> class A(StateVars):
     ...     _state_vars = [('a', 1), ('b', 2)]
     ...     process_vars()
 
     >>> d1
     D1(A_=A(a=3))
-    
+
     >>> d1.a=5; print d1
     a=5
     b=2
     >>> class D3(StateVars):
     ...     _state_vars = [('C_', Delegate(C, ['x','y']))]
     ...     process_vars()
-    
+
     >>> d3 = D3(x=10); d3
     D3(C_=C(y=2, x=10))
 
     >>> Delegate(None, vars=['x'], value=4)
     Delegate(value=4, vars=['x'], cls=None)
     """
+
     # pylint: disable-msg=W0231, W0622
-    def __init__(self, cls, vars=None, methods=None, 
+    def __init__(self, cls, vars=None, methods=None,
                  value=NotImplemented, cached=False):
         self.cls = cls
         self.vars = vars
         self.value = value
         self.cached = cached
 
+
 class Ref(Attr):
     r"""Ref Class.  Represents an attribute reference to another
     attribute.
         self.ref = ref
         self.cached = cached
         self.method = method
-    
+
+
 class NoCopy(object):                   # pylint: disable-msg=R0903
     r"""Default value for objects that should not be copied (or copied
     with a custom copier copy).  The default semantics are that
     NoCopy(3, copy=<function <lambda> at ...>)
     >>> n.copy                  # doctest: +ELLIPSIS
     <function <lambda> at ...>
-    
     """
-    def __init__(self, value, copy=False): # pylint: disable-msg=W0621
+
+    def __init__(self, value, copy=False):  # pylint: disable-msg=W0621
         self.value = value
         if copy is not False:
             self._copy = copy
+
     @property
     def copy(self):
         r"""Return the copier"""
             else:
                 return self._copy
         else:
-            return lambda x:x
+            return lambda x: x
+
     def __repr__(self):
         args = [repr(self.value)]
         if hasattr(self, '_copy'):
-            args.append('copy=%r'%self._copy)
+            args.append('copy=%r' % self._copy)
         return "NoCopy(%s)" % (", ".join(args),)
 
+
 def _classinitializer(proc):
     r"""Magic decorator to allow class decoration.
-    
+
     One added feature is that the return value of `proc` is used to
     set the docstring.  (This means that `proc` is called twice).
     """
         worms!)"""
         frame = sys._getframe(1)        # pylint: disable-msg=W0212
         if '__module__' in frame.f_locals and not \
-            '__module__' in frame.f_code.co_varnames: # we are in a class
+            '__module__' in frame.f_code.co_varnames:  # we are in a class
             thetype = frame.f_locals.get("__metaclass__", type)
+
             def makecls(name, bases, dic):
                 r"""Construct the new class and return it."""
                 # This dic is the source of a bug: See
                 # enclosing scope...
                 try:
                     cls = thetype(name, bases, dic)
-                except TypeError, err: # pragma: no cover
+                except TypeError, err:  # pragma: no cover
                     if "can't have only classic bases" in str(err):
                         cls = thetype(name, bases + (object, ), dic)
                     else:  # other strange errs, e.g. __slots__ conflicts
                 dic.update(res)
                 try:
                     cls = thetype(name, bases, dic)
-                except TypeError, err: # pragma: no cover
+                except TypeError, err:  # pragma: no cover
                     if "can't have only classic bases" in str(err):
                         cls = thetype(name, bases + (object, ), dic)
                     else:  # other strange errs, e.g. __slots__ conflicts
                 wrap__init__(cls)
                 # Need to do this after new class is created so
                 # _original__init__ has correct im_class.
-                
+
                 return cls
-            
+
             frame.f_locals["__metaclass__"] = makecls
         else:                   # pragma: no cover
             proc(*args, **kw)
     return newproc
 
 _NO_DESCRIPTION = "<no description>"
+
+
 def _normalize_state_vars(state_vars):
     r"""Return the normalized form of state_vars as a list of
     triplets.
     >>> _normalize_state_vars([('a'), # doctest: +NORMALIZE_WHITESPACE
     ...                        ('b', 1),
     ...                        ('c', 2, 'doc'),
-    ...                        ('d=c'), 
+    ...                        ('d=c'),
     ...                        ('e', Ref('c'), 'e doc'),
     ...                        ('f=a.x', 2, 'New f doc'),
     ...                        ('g=a.y', 2),
                 elif 3 > len(var_desc):
                     doc = _NO_DESCRIPTION
                 else:
-                    raise TypeError("Attr must be a str or tuple "+
+                    raise TypeError("Attr must be a str or tuple " +
                                     "of length 3 or less (got %s)" %\
                                         repr(var_desc))
-                
+
                 new_state_vars.append((var, Ref(ref), doc))
 
                 if 2 <= len(var_desc):
                                          (var_desc[1], var, ref, ref))
 
             elif (isinstance(var_desc, tuple)
-                  and 2 <= len(var_desc) 
+                  and 2 <= len(var_desc)
                   and isinstance(var_desc[1], Delegate)):
                 cls_state_vars = getattr(var_desc[1].cls, '_state_vars', [])
 
                     if not cls_state_vars:
                         # Add one element for construction zip(*) below
                         cls_state_vars = [(None, NotImplemented, '')]
-                    
+
                     extension.extend(
                         [(v, NotImplemented, _NO_DESCRIPTION)
                          for v in vars_
                     new_refs.append((_m,
                                      Ref(ref, method=True),
                                      doc))
-                        
+
                 for var, default, doc in extension:
                     ref = '%s.%s' % (var_desc[0], var)
                     # Insert new references at front so they don't hide
                     # explicit references.
                     new_refs.append((var,
-                                     Ref(ref,cached=var_desc[1].cached),
+                                     Ref(ref, cached=var_desc[1].cached),
                                      doc))
                     #if default is not NotImplemented:
                     #    state_vars.append((ref, default, ''))
                 elif 3 == len(var_desc):
                     var, default, doc = var_desc
                 else:
-                    raise TypeError("Attr must be a str or tuple "+
+                    raise TypeError("Attr must be a str or tuple " +
                                     "of length 3 or less (got %s)" %\
                                         repr(var_desc))
             else:
                 raise TypeError(
                     'state_vars must be a list of strings or tuples')
             if not isinstance(var, str):
-                raise TypeError('Name must be a string (got %s)'%\
+                raise TypeError('Name must be a string (got %s)' %\
                                     repr(var))
 
             # Now process references
             if '=' in var:
                 var, ref = var.split('=')
-                if default is not NotImplemented: # pragma: no cover
+                if default is not NotImplemented:  # pragma: no cover
                     raise ValueError("\n".join(
                             ["Cannot set default %r for ref '%s=%s'.",
                              "(Set default for %r instead)"]) %
             elif isinstance(key, tuple) and len(key) == 2:
                 state_vars_.append((key[0], state_vars[key], key[1]))
             else:
-                raise TypeError('Key must be name or (name, doc)'+
-                                ' (got %s)'%repr(key))
+                raise TypeError('Key must be name or (name, doc)' +
+                                ' (got %s)' % repr(key))
         return _normalize_state_vars(state_vars_)
     else:
         raise TypeError("state_vars must be a list or dict")
     return normalized_state_vars
 
+
 def _is_attr(attr):
     r"""Return `True` if `attr` is an attribute definition or a tuple of
     length 2 with the second element being an attribute definition.
             or _is(attr, Excluded)
             or _is(attr, Computed))
 
+
 def _preprocess(cls):
     r"""Preprocess `cls` converting all attributes in `cls.__dict__`
     that are attribute definitions as defined by :func:`_is_attr` to
             pass
         elif _is(attr, Computed):
             pass
-    
+
+
 def _gather_vars(cls):
     r"""Return `(vars, original_defaults, defaults, docs, excluded_vars,
     nodeps, computed_vars, refs, delegates, cached)` gathered from `cls`.
         to this set as the properties are stored in  `cls.__dict__`
         and we add the name to `cls._vars` so that the getters and
         setters will be called.
-    
+
     If `cls` has no :attr:`_state_vars`, then `vars` is `None`.
 
     >>> class A(object):
     >>> _gather_vars(B)         # doctest: +NORMALIZE_WHITESPACE
     (['a', 'b', 'n', 'c'],
      {'a': 3, 'c': 4, 'b': 2, 'n': NoCopy(3)},
-     {'a': 3, 'c': 4, 'b': 2, 'n': 3}, 
+     {'a': 3, 'c': 4, 'b': 2, 'n': 3},
      {'a': 'new_a', 'c': 'c', 'b': 'var b', 'n': 'var n'},
-     set([]), set([]), set([]), {}, set([]), {}, {}, {}, set([]), set([]), 
+     set([]), set([]), set([]), {}, set([]), {}, {}, {}, set([]), set([]),
      set(['a']))
 
     >>> class C(object):
     >>> class D(object):
     ...    _state_vars = [
     ...        ('c', Delegate(C)),
-    ...        ('y=c.x'), 
+    ...        ('y=c.x'),
     ...        ('c.x', 3),
     ...        ('h=c.a.n')]
     >>> _gather_vars(D)    # doctest: +NORMALIZE_WHITESPACE, +ELLIPSIS
         all_bases.append(Class)
     else:
         all_bases.append(cls)
-    
+
     # vars_ == None means no _state_vars which is different than empty
     # _state_vars.
     defaults = {}
                 if default is Deleted:
                     deleted.add(var)
                 elif var in vars_ and var not in deleted:
-                    if (defaults.get(var,None) is not Required # Fix issue 14 
+                    if (defaults.get(var,None) is not Required  # Fix issue 14
                         and var not in ignore_name_clash):
                         warnings.warn(NameClashWarning1(
                             "Redefining attr '%s' in class '%s'"
     # Delete all "Deleted" vars.
     for var in deleted:
         if var not in vars_:
-            raise ValueError("Attempt to delete non-existent state var '%s'" 
+            raise ValueError("Attempt to delete non-existent state var '%s'"
                              % (var,))
         else:
             while var in vars_:
                 del docs[var]
             while var in defaults:
                 del defaults[var]
-    
+
     # Now check for any data descriptors in the class dictionary.
     # If no corresponding var has been specified in _state_vars,
     # then a new entry will be added.  If there is no setter
             # Continue processing value.  This allows
             # NoCopy(Excluded) for example.
             default = default.value
-            
+
         if default is NotImplemented:
             del defaults[var]
         elif _is(default, ClassVar):
             raise ValueError(
                 "Required should not be called " +
                 "in _state_vars (got arg %r)" % (default.doc,))
-        
+
         if '.' in var and var in docs:
             # Removed documentation from overridden reference defaults.
             del docs[var]
             if inds:
                 # pylint: disable-msg=W0212
                 return cls._state_vars[inds[0]][-1]
-        except Exception:               #pragma: no cover
-            1+1 # pass, pylint: disable-msg=W0104
+        except Exception:               # pragma: no cover
+            1 + 1  # pass, pylint: disable-msg=W0104
 
         attr, _sep, tail = var.partition('.')
         #if sep and hasattr(cls, attr):
     refs = {}
     for var in refs_:
         ref_ = refs_[var]
-        ref = refs.get(ref_, ref_) # Dereference
+        ref = refs.get(ref_, ref_)  # Dereference
         if '.' in ref:
-            head, sep, tail  = ref.partition('.')
+            head, sep, tail = ref.partition('.')
             attr = refs.get(head, head)
             ref = "".join((head, sep, tail))
         else:
             refs, irefs, delegates, settable_properties, cached,
             ignore_name_clash)
 
+
 def _format_varstring(doc, width=70, indent=''):
     r"""Return the formatted varstring, processing it like a docstring
     and removing indentation.
     ...                        be preserved if it is indented:
     ...                           1) A series of points
     ...                           2) Here is a longer point
-    ...              
+    ...
     ...                        You can use paragraphs too!  They will
     ...                        also be wrapped for you.'''
     >>> print _format_varstring(a, width=66, indent='  ')
     """
     if with_docs is None:
         return default
-    
+
     if (not isinstance(with_docs, list)
         and not isinstance(with_docs, tuple)):
         with_docs = [with_docs]
             try:
                 doc = getattr(obj, name).__doc__
             except AttributeError:
-                1+1 # pass, pylint: disable-msg=W0104
+                1 + 1  # pass, pylint: disable-msg=W0104
         if doc:
             return doc
     return default
-            
+
+
 # pylint: disable-msg=W0621
 def _get_processed_vars(cls, copy, archive_check, with_docs=None):
     r"""Return (results, `__doc__`) where `results` is a dictionary of
 
     def wrap__init__(cls):
         r"""Wraps the original `cls.__init__` with a new version that
-        does processing."""        
+        does processing."""
         cls._original__init__ = cls.__init__
+
         def __init__(self, *varargin, **kwargs):
             r"""Wrapper for :meth:`__init__` to set _initializing flag."""
-            # pass, pylint: disable-msg=W0212            
+            # pass, pylint: disable-msg=W0212
             if '_no_init' in self.__dict__:
                 # Short circuit for copy construction:
                 del self.__dict__['_no_init']
             # Update all cached values.
             for name in self._cached:
                 self.__dict__[name] = getattr(self, self._refs[name])
-            
+ 
             try:
                 if self._call_base__init__:
                     for cl in cls.__bases__:
     doc = _get__doc__(cls, state_vars=results['_state_vars'])
     return (results, doc, wrap__init__)
 
+
 @_classinitializer
 def process_vars(cls, copy=copy.deepcopy, archive_check=None,
                  with_docs=None):

File mmf/utils/_utils.py

 __all__ = ['memory', 'print_mem',
            'get_args', 'all_bases', 'memoize', 'unique_list',
            'isequal', 'Timer', 'StatusBar',
-           'inline', 'vectorize']
+           'inline', 'vectorize', 'make_sure_path_exists']
 
 import sys
 import os
     print("%.2f%s" %((mem - since)/factors[unit], unit))
     return mem
 
+import os
+import errno
+
+def make_sure_path_exists(path):
+    r"""Make sure that all directories in `path` exist, creating them as
+    needed.  See http://stackoverflow.com/a/5032238/1088938"""
+    try:
+        os.makedirs(path)
+    except OSError as exception:
+        if exception.errno != errno.EEXIST:
+            raise
 
 class Timer(object):
     """An object that is True for a specified period of time, then

File mmf/utils/mmf_sympy.py

-r"""Tools for using sympy to generate code."""
-
+r"""Tools for using sympy to generate code.
+"""
 import numpy as np
 import math
 import sympy
 _ENVS = dict(math=dict(math.__dict__, math=math),
              numpy=dict(np.__dict__, numpy=np, np=np))
 
+
 def lambdify_cse(var, expr, modules='numpy', env={}):
     r"""Generates a function that evaluates 'expr' but uses cse to optimize
     common subexpressions."""
     env.update(env)
     exec code in env
     return env['f']
+
+
+class ODE(object):
+    r"""A set of functions deals with manipulating ODEs.  These are represented
+    as a set of equations in the form ``Derivative(y(x), x, m) == f(x, y(x),
+    ...)```.  A normalized form exists where each equation is first order and
+    the various derivatives are represented as dependent variables (by default
+    named ``y'(x)``, ``y''(x)`` etc.).  One can :meth:`compress` these equations
+    to obtain the form with higher orders but no intermediate variables.
+
+    Examples
+    --------
+    >>> from sympy import *
+    >>> x, u = symbols('x u')
+    >>> eq = ln(u(x).diff(x,x)) - x*u(x)**2
+    >>> ode = ODE([eq], [u])
+    >>> ode
+    [Derivative(u(x), x, x) == exp(x*u(x)**2)]
+    >>> ode.normalize()
+    [Derivative(u(x), x) == u'(x), Derivative(u'(x), x) == exp(x*u(x)**2)]
+    >>> ode.compress()
+    [Derivative(u(x), x, x) == exp(x*u(x)**2)]
+
+    The equations can be specified in any form that :func:`sympy.solve` can use
+    to isolate the dependent variables:
+    >>> eq = ln(u(x).diff(x,x)) - (x*u(x)).diff(x)
+    >>> ode = ODE([eq], [u(x)])
+    >>> ode
+    [Derivative(u(x), x, x) == exp(x*Derivative(u(x), x) + u(x))]
+    >>> ode.normalize()
+    {Derivative(u(x), x): u'(x)}
+    >>> ode
+    [Derivative(u(x), x) == u'(x), Derivative(u'(x), x) == exp(x*u'(x) + u(x))]
+    >>> ode.compress()
+    [Derivative(u(x), x, x) == exp(x*Derivative(u(x), x) + u(x))]
+
+    """
+
+    def __init__(self, eqs, ys=None):
+        r"""Accepts two types of inputs.
+
+        Parameters
+        ----------
+        eqs : list
+           If this is the sole parameter, then it should be a list of equations
+           with a lhs of the form `y(x).diff(x, order)` from which the
+           variables and order can be deduced.  Otherwise it is just a list of
+           equations.
+        ys : list
+           A list of dependent variables as functions of the abscissa:
+           ie. ``[y(x)]``"""
+        if ys is None:
+            ys = []
+            fs = []
+            orders = []
+            x = eqs[0].lhs.args[0].args[0]
+            for eq in eqs:
+                args = eq.lhs.args
+                y_x = args[0]
+                assert (1 == len(y_x.args)
+                        and x == y_x.args[0]
+                        and not [_x for _x in args[1:] if _x != x])
+                ys.append(y_x)
+                fs.append(eq.rhs)
+                orders.append(len(args[1:]))
+        else:
+            ys = map(sympy.S, ys)
+            eqs = map(sympy.S, eqs)
+            x = ys[0].args[0]
+
+            orders = []
+            dys = []
+            for y_x in ys:
+                assert (1 == len(y_x.args) and x == y_x.args[0])
+                orders.append(max(sympy.ode.ode_order(_eq, y_x) for _eq in eqs))
+                dys.append(y_x.diff(x, orders[-1]))
+            fs = sympy.solve(eqs, dys)
+            fs = [fs[_dy] for _dy in dys]
+            if not len(fs) == len(dys):
+                raise ValueError(
+                    "Could not define the ODE: must be able to isolate the " +
+                    "highest orders: %s" % (str(dys),))
+
+        self.x = x
+        self.ys = ys
+        self.fs = fs
+        self.orders = orders
+
+    def normalize(self):
+        r"""Convert to a series of first order equations by adding variables
+        such as `u'`, `u''`, etc."""
+        x = self.x
+        ys = []
+        fs = []
+        subs = {}
+        for y_x, f, o in zip(self.ys, self.fs, self.orders):
+            ys.append(y_x)
+            while 1 < o:
+                _y = sympy.Function(str(y_x.func) + "'")
+                assert _y(x) not in self.ys
+                ys.append(_y(x))
+                fs.append(ys[-1])
+                y_x = y_x.diff(x)
+                subs[y_x] = ys[-1]
+                o -= 1
+            fs.append(f)
+
+        # Replace derivatives with variables.
+        fs = [_f.subs(subs) for _f in fs]
+        self.ys = ys
+        self.fs = fs
+        self.orders = [1, ] * len(ys)
+        return subs
+        return self.as_expr()
+
+    def compress(self):
+        x = self.x
+        ys = self.ys
+        fs = self.fs
+        orders = self.orders
+
+        subs = {}
+        dummys = [y_x for y_x in ys if y_x in fs]
+        for _d in dummys:
+            i_y = ys.index(_d)
+            i_f = fs.index(y_x)
+            orders[i_f] += orders[i_y]
+            subs[ys[i_y]] = ys[i_f].diff(x)
+            fs[i_f] = fs[i_y]
+            del ys[i_y], fs[i_y], orders[i_y]
+        self.ys = ys
+        self.fs = [_f.subs(subs) for _f in fs]
+        self.orders = orders
+        return self.as_expr()
+
+    def as_expr(self):
+        x = self.x
+        return [sympy.Eq(_y_x.diff(x, _o), _f)
+                for _y_x, _f, _o in zip(self.ys, self.fs, self.orders)]
+
+    def __repr__(self):
+        return repr(self.as_expr())
+
+    def __str__(self):
+        return str(self.as_expr())
+
+
+def free_functions(expr, var=None):
+    r"""Return a set of all the free unknown functions in expr (similar to
+    `free_symbols`).
+
+    Examples
+    --------
+    >>> from sympy import S
+    >>> free_functions(S('u(x)*f(w)+sin(k(s(p)))'))
+    set([f(w), s(p), u(x)])
+    >>> free_functions(S('u(x)*f(w)+sin(k(s(p)))'), S('x'))
+    set([u(x)])
+    """
+    if expr.is_Atom:
+        return set()
+    elif (expr.is_Function
+          and isinstance(type(expr),
+                         sympy.function.UndefinedFunction)
+          and ((var is None and expr.args[0].is_Symbol)
+               or 
+               expr.args[0] == var)):
+        return set([expr])
+    else:
+        return set().union(*(free_functions(_a, var=var) for _a in expr.args))
+
+
+class Transform(object):
+    r"""Represents a change of variables.
+
+    Parameters
+    ----------
+    transforms : list
+       List of old variables defined in terms of the new variables:
+       `[Eq(x, f(y, v1(y), ...)), Eq(u1(x), ...), ...]`.
+    new_vars : list
+       List of new variables in the form `[y, v1(y), v2(y), v3(y), ...]`.  If
+       this is not provided, it will be deduced from the equations (but then all
+       function applications must be of the form `v1(y)`).
+
+    Examples
+    --------
+    >>> from sympy import S
+    >>> t = Transform(['x=2*y*v(y)', 'u(x)=y*v(y)*v(y).diff(y)'])
+    >>> t.transforms
+    [x == 2*y*v(y), u(x) == y*v(y)]
+    >>> t.old_vars
+    [x, u(x)]
+    >>> t.new_vars
+    [u, v(y)]
+
+    One can also include the derivatives.  These are not explicitly included as
+    new variables though
+    >>> t = Transform(['x=2*y*v(y)', 'u(x)=y*v(y)*v(y).diff(y)'])
+    >>> t.transforms
+    [x == 2*y*v(y), u(x) == y*v(y)*Derivative(v(y), y)]
+    >>> t.new_vars
+    [u, v(y)]
+
+    """
+    def __init__(self, transforms, new_vars=None):
+        transforms = map(self._S, transforms)
+        transforms = map(sympy.S, transforms)
+        if new_vars is None:
+            new_vars = list(set().union(*(free_functions(_f.rhs)
+                                          for _f in transforms)))
+            new_vars.insert(0, new_vars[0].args[0])
+        else:
+            new_vars = map(sympy.S, new_vars)
+
+        old_vars = [_t.lhs for _t in transforms]
+
+        assert np.all([1 == len(_v.args) for _v in old_vars[1:]])
+        assert np.all([1 == len(_v.args) for _v in new_vars[1:]])
+        assert np.all([old_vars[0] == _v.args[0] for _v in old_vars[1:]])
+        assert np.all([new_vars[0] == _v.args[0] for _v in new_vars[1:]])
+        assert len(new_vars) == len(old_vars)
+
+        self.transforms = transforms
+        self.new_vars = new_vars
+        self.old_vars = old_vars
+
+    @staticmethod
+    def _S(eq):
+        r"""Special sympify to allow for strings like 'u(x)=y + v(y)'."""
+        if isinstance(eq, sympy.Basic):
+            pass
+        elif '=' in eq:
+            eq = sympy.Eq(*map(sympy.S, eq.split('=')))
+        else:
+            eq = sympy.S(eq)
+        return eq
+
+    def as_dict(self, transforms=None):
+        if transforms is None:
+            transforms = self.transforms
+        x = self.old_vars[0]
+        y = self.new_vars[0]
+        return dict((_t.lhs, _t.rhs.subs({y: y(x)}))
+                    for _t in transforms)
+
+    def __repr__(self):
+        return repr(self.transforms)
+
+    def transform(self, ode):
+        ode.normalize()
+
+        x = ode.x
+        assert x == self.old_vars[0]
+        old_vars = [x]
+
+        y = self.new_vars[0]
+        t_dict = self.as_dict()
+        u_dict = self.as_dict()
+        del u_dict[x]
+        dy_dx = {y(x).diff(x):
+                     sympy.solve(t_dict[x].diff(x) - 1, y(x).diff(x))[0]}
+
+        # For simplicity, we assume that ys contains chains:
+        # [u(x), u'(x), u''(x), v(x), v'(x), ...] with primes
+        u = None
+        for _u in ode.ys:
+            _name = str(_u.func)
+            if _name[-1] == "'":
+                assert u is not None
+                u = u.diff(x)
+                t_dict[_u] = u.subs(u_dict).doit().subs(dy_dx)
+            else:
+                u = _u
+            old_vars.append(_u)
+
+            subs = {}
+        new_vars = [y]
+        for v_y in self.new_vars[1:]:
+            o = max(sympy.ode.ode_order(_t, v_y) for _t in t_dict.values())
+            _v_y = v_y
+            new_vars.append(_v_y)
+            while 0 < o:
+                _v = sympy.Function(str(_v_y.func) + "'")
+                _v_y = _v(y)
+                new_vars.append(_v_y)
+                v_y = v_y.diff(y)
+                subs[v_y] = _v_y
+                o -= 1
+        for _t in t_dict:
+            t_dict[_t] = t_dict[_t].subs(subs).doit()
+        return old_vars, new_vars, t_dict
+
+        # Sort transforms
+        old_vars = [ode.x] + ode.ys
+        transforms = ([self.transforms[self.old_vars.index(ode.x)]] +
+                      [self.transforms[self.old_vars.index(_v)]
+                       for _v in ode.ys])
+        t_dict = self.as_dict(transforms)
+
+        fs, xs, ys, subs = ode_normalize(
+            eqs=eqs, vars=vars, new_vars=new_vars)
+        trans = sympy.Matrix([[sympy.S(_t.rhs)
+                               for _t in transforms]]).subs(subs)
+        J = trans.jacobian(ys)
+        Jinv_fs = J.LUsolve(sympy.Matrix(fs))
+        y = ys[0]
+        new_eqs = [ys[_i](y).diff(y) - (Jinv_fs[_i] / Jinv_fs[0])
+                   for _i in xrange(len(fs) - 1)]
+        return new_eqs, ys
+
+
+def ode_transform(eqs, transforms, new_vars):
+    r"""Effects a change of variables.
+
+    Parameters
+    ----------
+    eqs : list
+       List of equations.  The i'th equation should have the highest order in
+      `u_i(x).diff(x, m_i)` and this must able to be solved for.
+    transforms : list
+       List of old variables defined in terms of the new variables:
+       `[Eq(x, f(y, v1, ...)), Eq(u1, ...), ...]`.
+    new_vars : list
+       List of new variables in the form `[y, v1, v2, v3, ...]`
+
+    Examples
+    --------
+    >>> from sympy import *
+    >>> x, u = symbols('x u')
+    >>> eq = ln(u(x).diff(x,x)) - x*u(x)**2
+    >>> ode_normalize([eq], [x, u])
+    ([1, u, exp(x*u**2)], [x, u', u''])
+    >>> eq = ln(u(x).diff(x,x)) - (x*u(x)).diff(x)
+    >>> ode_normalize([eq], [x, u])
+    ([1, u, exp(u + u'*x)], [x, u', u''])
+    >>> 
+    """
+    eqs = map(sympy.S, eqs)
+    vars = [sympy.S(_t.lhs) for _t in transforms]
+    new_vars = map(sympy.S, new_vars)
+    fs, xs, ys, subs = ode_normalize(
+        eqs=eqs, vars=vars, new_vars=new_vars)
+    trans = sympy.Matrix([[sympy.S(_t.rhs) for _t in transforms]]).subs(subs)
+    J = trans.jacobian(ys)
+    Jinv_fs = J.LUsolve(sympy.Matrix(fs))
+    y = ys[0]
+    new_eqs = [ys[_i](y).diff(y) - (Jinv_fs[_i]/Jinv_fs[0]) 
+               for _i in xrange(len(fs)-1)]
+    return new_eqs, ys
+
     
-        
+    
+    
+
+