From 473679e469f318c28ab692df2f9c8b8209803368 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Bj=C3=B6rn=20Dahlgren?= <bjodah@gmail.com>
Date: Tue, 19 Sep 2017 00:52:58 +0200
Subject: [PATCH 1/6] Support for order kwarg in Lambdify

---
 symengine/lib/symengine_wrapper.pyx | 253 ++++++++++++++--------------
 symengine/tests/test_lambdify.py    | 252 +++++++++++++++++++++++----
 symengine/tests/test_matrices.py    |  10 +-
 3 files changed, 352 insertions(+), 163 deletions(-)

diff --git a/symengine/lib/symengine_wrapper.pyx b/symengine/lib/symengine_wrapper.pyx
index 815a7c035..b99832698 100644
--- a/symengine/lib/symengine_wrapper.pyx
+++ b/symengine/lib/symengine_wrapper.pyx
@@ -9,6 +9,7 @@ from cpython cimport PyObject, Py_XINCREF, Py_XDECREF, \
 from libc.string cimport memcpy
 import cython
 import itertools
+import numbers
 from operator import mul
 from functools import reduce
 import collections
@@ -462,7 +463,8 @@ def sympy2symengine(a, raise_error=False):
             return conditionset(*(a.args))
 
     if raise_error:
-        raise SympifyError("sympy2symengine: Cannot convert '%r' to a symengine type." % a)
+        raise SympifyError(("sympy2symengine: Cannot convert '%r' (of type %s)"
+                            " to a symengine type.") % (a, type(a)))
 
 
 def sympify(a):
@@ -523,7 +525,7 @@ def _sympify(a, raise_error=True):
         return a
     elif isinstance(a, bool):
         return (true if a else false)
-    elif isinstance(a, (int, long)):
+    elif isinstance(a, numbers.Integral):
         return Integer(a)
     elif isinstance(a, float):
         return RealDouble(a)
@@ -543,7 +545,9 @@ def _sympify(a, raise_error=True):
         pass
 
     if raise_error:
-        raise SympifyError("sympify: Cannot convert '%r' to a symengine type." % a)
+        raise SympifyError(
+            "sympify: Cannot convert '%r' (of type %s) to a symengine type." % (
+                a, type(a)))
 
 funcs = {}
 
@@ -553,7 +557,7 @@ def get_function_class(function, module):
     return funcs[function]
 
 class Singleton(object):
-    
+
     __call__ = staticmethod(sympify)
 
     @property
@@ -980,7 +984,7 @@ cdef class Basic(object):
     @property
     def is_Not(self):
         return False
-    
+
     @property
     def is_Matrix(self):
         return False
@@ -1128,7 +1132,7 @@ class Symbol(Expr):
 
 
 class Dummy(Symbol):
-    
+
     def __init__(Basic self, name=None, *args, **kwargs):
         if name is None:
             self.thisptr = symengine.make_rcp_Dummy()
@@ -1163,7 +1167,6 @@ def symarray(prefix, shape, **kwargs):
     This function requires NumPy.
 
     """
-    import numpy as np
     arr = np.empty(shape, dtype=object)
     for index in np.ndindex(shape):
         arr[index] = Symbol('%s_%s' % (prefix, '_'.join(map(str, index))), **kwargs)
@@ -1273,13 +1276,13 @@ eulergamma = EulerGamma()
 
 
 cdef class Boolean(Expr):
-    
+
     def logical_not(self):
         return c2py(<RCP[const symengine.Basic]>(deref(symengine.rcp_static_cast_Boolean(self.thisptr)).logical_not()))
 
 
 cdef class BooleanAtom(Boolean):
-    
+
     @property
     def is_Boolean(self):
         return True
@@ -1364,7 +1367,7 @@ class Xor(Boolean):
 
 
 class Relational(Boolean):
-    
+
     @property
     def is_Relational(self):
         return True
@@ -2299,7 +2302,7 @@ class polygamma(Function):
         return sympy.polygamma(*self.args_as_sympy())
 
 class sign(OneArgFunction):
-    
+
     @property
     def is_complex(self):
         return True
@@ -2877,17 +2880,17 @@ class UniversalSet(Set):
 
 
 class FiniteSet(Set):
-    
+
     def __new__(self, *args):
         return finiteset(*args)
 
     def _sympy_(self):
         import sympy
-        return sympy.FiniteSet(*[arg._sympy_() for arg in self.args])    
+        return sympy.FiniteSet(*[arg._sympy_() for arg in self.args])
 
 
 class Contains(Boolean):
-    
+
     def __new__(self, expr, sset):
         return contains(expr, sset)
 
@@ -2897,7 +2900,7 @@ class Contains(Boolean):
 
 
 class Union(Set):
-    
+
     def __new__(self, *args):
         return set_union(*args)
 
@@ -2907,7 +2910,7 @@ class Union(Set):
 
 
 class Complement(Set):
-    
+
     def __new__(self, universe, container):
         return set_complement(universe, container)
 
@@ -2917,13 +2920,13 @@ class Complement(Set):
 
 
 class ConditionSet(Set):
-    
+
     def __new__(self, sym, condition):
         return conditionset(sym, condition)
 
 
 class ImageSet(Set):
-    
+
     def __new__(self, sym, expr, base):
         return imageset(sym, expr, base)
 
@@ -3540,7 +3543,11 @@ cdef class DenseMatrixBase(MatrixBase):
         return ImmutableDenseMatrix(self)
 
     def tolist(self):
-        return self[:]
+        return [[self[rowi, coli] for coli in range(self.ncols())]
+                for rowi in range(self.nrows())]
+
+    def __array__(self):
+        return np.array(self.tolist())
 
     def _mat(self):
         return self
@@ -4316,54 +4323,6 @@ def has_symbol(obj, symbol=None):
                 deref(symengine.rcp_static_cast_Symbol(s.thisptr)))
 
 
-cdef size_t _size(n):
-    try:
-        return n.size
-    except AttributeError:
-        return len(n)  # e.g. array.array
-
-
-def _get_shape_nested(ndarr):
-    # no checking of shape consistency is done
-    if isinstance(ndarr, (list, tuple)):
-        return (len(ndarr),) + _get_shape_nested(ndarr[0])
-    else:
-        return ()
-
-
-def get_shape(ndarr):
-    try:
-        return ndarr.shape
-    except AttributeError:
-        return _get_shape_nested(ndarr)
-
-
-def _nested_getitem(ndarr, indices):
-    if len(indices) == 0:
-        return ndarr
-    else:
-        return _nested_getitem(ndarr[indices[0]], indices[1:])
-
-
-def all_indices_from_shape(shape):
-    return itertools.product(*(range(dim) for dim in shape))
-
-
-def _ravel_nested(ndarr):
-    return [_nested_getitem(ndarr, indices) for indices in
-            all_indices_from_shape(get_shape(ndarr))]
-
-
-def ravel(ndarr):
-    try:
-        return ndarr.ravel()
-    except AttributeError:
-        try:
-            return _ravel_nested(ndarr.tolist())
-        except AttributeError:
-            return _ravel_nested(ndarr)
-
-
 cdef class _Lambdify(object):
     """
     Lambdify instances are callbacks that numerically evaluate their symbolic
@@ -4399,16 +4358,20 @@ cdef class _Lambdify(object):
     cdef list out_shapes
     cdef readonly bint real
     cdef readonly int n_exprs
+    cdef readonly str order
     cdef vector[int] accum_out_sizes
     cdef object numpy_dtype
 
-    def __cinit__(self, args, *exprs, cppbool real=True):
+    def __cinit__(self, args, *exprs, cppbool real=True, order='C'):
         cdef vector[int] out_sizes
         self.real = real
+        self.order = order
         self.numpy_dtype = np.float64 if self.real else np.complex128
-        self.out_shapes = [get_shape(expr) for expr in exprs]
+        self.out_shapes = [np.asarray(expr).shape for expr in exprs]
         self.n_exprs = len(exprs)
-        self.args_size = _size(args)
+        self.args_size = np.asarray(args).size
+        if self.args_size == 0:
+            raise NotImplementedError("Support for zero arguments not yet supported")
         self.tot_out_size = 0
         for idx, shape in enumerate(self.out_shapes):
             out_sizes.push_back(reduce(mul, shape or (1,)))
@@ -4418,7 +4381,7 @@ cdef class _Lambdify(object):
             for j in range(i):
                 self.accum_out_sizes[i] += out_sizes[j]
 
-    def __init__(self, args, *exprs, cppbool real=True):
+    def __init__(self, args, *exprs, cppbool real=True, order='C'):
         cdef:
             Basic e_
             size_t ri, ci, nr, nc
@@ -4434,7 +4397,7 @@ cdef class _Lambdify(object):
                 for ci in range(nc):
                    args_.push_back(deref(mtx).get(ri, ci))
         else:
-            for arg in args:
+            for arg in np.ravel(args, order=order):
                 e_ = _sympify(arg)
                 args_.push_back(e_.thisptr)
 
@@ -4449,7 +4412,7 @@ cdef class _Lambdify(object):
                         b_ = deref(mtx).get(ri, ci)
                         outs_.push_back(b_)
             else:
-                for e in ravel(curr_expr):
+                for e in np.ravel(curr_expr, order=order):
                     e_ = _sympify(e)
                     outs_.push_back(e_.thisptr)
         self._init(args_, outs_)
@@ -4484,7 +4447,7 @@ cdef class _Lambdify(object):
             raise ValueError("Size of out incompatible with number of exprs.")
         self.unsafe_complex(inp, out)
 
-    def __call__(self, inp, *, out=None):
+    def __call__(self, inp, *, out=None, order=None):
         """
         Parameters
         ----------
@@ -4495,6 +4458,12 @@ cdef class _Lambdify(object):
             If ``None``: an output container will be allocated (NumPy ndarray).
             If ``len(exprs) > 0`` output is found in the corresponding
             order.
+        order : 'C' or 'F'
+            C- or Fortran-contiguous memory layout. Note that this affects
+            broadcasting: e.g. a (m, n) matrix taking 3 arguments and given a
+            (k, l, 3) (C-contiguous) input will give a (k, l, m, n) shaped output,
+            whereas a (3, k, l) (C-contiguous) input will give a (m, n, k, l) shaped
+            output. If ``None`` order is taken as ``self.order`` (from initialization).
 
         Returns
         -------
@@ -4507,36 +4476,73 @@ cdef class _Lambdify(object):
             tuple inp_shape
             double[::1] real_out, real_inp
             double complex[::1] cmplx_out, cmplx_inp
+        if order is None:
+            order = self.order
+        if order not in ('C', 'F'):
+            raise NotImplementedError("Only C & F order supported for now.")
+
         try:
-            inp = np.ascontiguousarray(inp, dtype=self.numpy_dtype)
+            inp = np.asanyarray(inp, dtype=self.numpy_dtype)
         except TypeError:
             inp = np.fromiter(inp, dtype=self.numpy_dtype)
-        inp_shape = inp.shape
 
         if self.real:
-            real_inp = inp.ravel()
+            real_inp = np.ascontiguousarray(inp.ravel(order=order))
         else:
-            cmplx_inp = inp.ravel()
+            cmplx_inp = np.ascontiguousarray(inp.ravel(order=order))
 
-        if inp.size % self.args_size != 0:
-            raise ValueError("Broadcasting failed")
+        if inp.size < self.args_size or inp.size % self.args_size != 0:
+            raise ValueError("Broadcasting failed (input/arg size mismatch)")
         nbroadcast = inp.size // self.args_size
-        if nbroadcast > 1 and self.args_size == 1 and inp.shape[-1] != 1:  # Implicit reshape
-            inp_shape = inp.shape + (1,)
-        else:
-            inp_shape = inp.shape
+
+        if inp.ndim > 1:
+            if self.args_size > 1:
+                if order == 'C':
+                    if inp.shape[inp.ndim-1] != self.args_size:
+                        raise ValueError(("C order implies last dim (%d) == len(args)"
+                                          " (%d)") % (inp.shape[inp.ndim-1], self.args_size))
+                    extra_dim = inp.shape[:inp.ndim-1]
+                elif order == 'F':
+                    if inp.shape[0] != self.args_size:
+                        raise ValueError("F order implies first dim (%d) == len(args) (%d)"
+                                         % (inp.shape[0], self.args_size))
+                    extra_dim = inp.shape[1:]
+            else:
+                extra_dim = inp.shape
+       else:
+            if nbroadcast > 1 and inp.ndim == 1:
+                extra_dim = (nbroadcast,)  # special case
+            else:
+                extra_dim = ()
+        extra_left = extra_dim if order == 'C' else ()
+        extra_right = () if order == 'C' else extra_dim
+        new_out_shapes = [extra_left + out_shape + extra_right
+                          for out_shape in self.out_shapes]
+
         new_tot_out_size = nbroadcast * self.tot_out_size
-        new_out_shapes = [inp_shape[:-1] + out_shape for out_shape in self.out_shapes]
         if out is None:
-            out = np.empty(new_tot_out_size, dtype=self.numpy_dtype)
+            out = np.empty(new_tot_out_size, dtype=self.numpy_dtype, order=order)
         else:
             if out.size < new_tot_out_size:
                 raise ValueError("Incompatible size of output argument")
-            if not (out.flags['C_CONTIGUOUS'] or out.flags['F_CONTIGUOUS']):
-                raise ValueError("Output argument needs to be C-contiguous")
+            if out.ndim > 1:
+                if order == 'C' and not out.flags['C_CONTIGUOUS']:
+                    raise ValueError("Output argument needs to be C-contiguous")
+                elif order == 'F' and not out.flags['F_CONTIGUOUS']:
+                    raise ValueError("Output argument needs to be F-contiguous")
+                if len(self.out_shapes) > 1:
+                    raise ValueError("output array with ndim > 1 assumes one output")
+                out_shape, = self.out_shapes
+                if order == 'C' and out.shape[-len(out_shape):] != tuple(out_shape):
+                    raise ValueError("shape mismatch for output array")
+                elif order == 'F' and out.shape[:len(out_shape)] != tuple(out_shape):
+                    raise ValueError("shape mismatch for output array")
+            else:
+                if not out.flags['F_CONTIGUOUS']:  # or C_CONTIGUOUS (ndim <= 1)
+                    raise ValueError("Output array need to be contiguous")
             if not out.flags['WRITEABLE']:
                 raise ValueError("Output argument needs to be writeable")
-            out = out.ravel()
+            out = out.ravel(order=order)
 
         if self.real:
             real_out = out
@@ -4552,10 +4558,18 @@ cdef class _Lambdify(object):
                 self.unsafe_complex(cmplx_inp, cmplx_out,
                                     idx*self.args_size, idx*self.tot_out_size)
 
-        out = out.reshape((nbroadcast, self.tot_out_size))
-        result = [out[:, self.accum_out_sizes[idx]:self.accum_out_sizes[idx+1]].reshape(
-            new_out_shapes[idx]) for idx in range(self.n_exprs)]
-
+        if order == 'C':
+            out = out.reshape((nbroadcast, self.tot_out_size), order='C')
+            result = [
+                out[:, self.accum_out_sizes[idx]:self.accum_out_sizes[idx+1]].reshape(
+                    new_out_shapes[idx], order='C') for idx in range(self.n_exprs)
+            ]
+        elif order == 'F':
+            out = out.reshape((self.tot_out_size, nbroadcast), order='F')
+            result = [
+                out[self.accum_out_sizes[idx]:self.accum_out_sizes[idx+1], :].reshape(
+                    new_out_shapes[idx], order='F') for idx in range(self.n_exprs)
+            ]
         if self.n_exprs == 1:
             return result[0]
         else:
@@ -4595,12 +4609,12 @@ IF HAVE_SYMENGINE_LLVM:
             self.lambda_double[0].call(&out[out_offset], &inp[inp_offset])
 
 
-def Lambdify(args, *exprs, cppbool real=True, backend=None):
+def Lambdify(args, *exprs, cppbool real=True, backend=None, order='C'):
     if backend is None:
         backend = os.getenv('SYMENGINE_LAMBDIFY_BACKEND', "lambda")
     if backend == "llvm":
         IF HAVE_SYMENGINE_LLVM:
-            return LLVMDouble(args, *exprs, real=real)
+            return LLVMDouble(args, *exprs, real=real, order=order)
         ELSE:
             raise ValueError("""llvm backend is chosen, but symengine is not compiled
                                 with llvm support.""")
@@ -4608,13 +4622,13 @@ def Lambdify(args, *exprs, cppbool real=True, backend=None):
         pass
     else:
         warnings.warn("Unknown SymEngine backend: %s\nUsing backend='lambda'" % backend)
-    return LambdaDouble(args, *exprs, real=real)
+    return LambdaDouble(args, *exprs, real=real, order=order)
 
 
-def LambdifyCSE(args, *exprs, cse=None, concatenate=None, **kwargs):
-    """
-    Analogous with Lambdify but performs common subexpression elimination
-    internally. See docstring of Lambdify.
+def LambdifyCSE(args, *exprs, cse=None, order='C', **kwargs):
+    """ Analogous with Lambdify but performs common subexpression elimination.
+
+    See docstring of Lambdify.
 
     Parameters
     ----------
@@ -4622,19 +4636,16 @@ def LambdifyCSE(args, *exprs, cse=None, concatenate=None, **kwargs):
     exprs: iterable of expressions (with symbols from args)
     cse: callback (default: None)
         defaults to sympy.cse (see SymPy documentation)
-    concatenate: callback (default: numpy.concatenate)
-        Examples when not using numpy:
-        ``lambda tup: tup[0]+list(tup[1])``
-        ``lambda tup: tup[0]+array.array('d', tup[1])``
-    \*\*kwargs: Keyword arguments passed onto Lambdify
+    order : str
+        order (passed to numpy.ravel and numpy.reshape).
+    \\*\\*kwargs: Keyword arguments passed onto Lambdify
 
     """
     if cse is None:
         from sympy import cse
-    if concatenate is None:
-        from numpy import concatenate
+    _exprs = [np.asanyarray(e) for e in exprs]
     from sympy import sympify as s_sympify
-    flat_exprs = list(itertools.chain(*map(ravel, exprs)))
+    flat_exprs = list(itertools.chain(*[np.ravel(e, order=order) for e in _exprs]))
     subs, flat_new_exprs = cse([s_sympify(expr) for expr in flat_exprs])
 
     explicit_subs = {}
@@ -4645,25 +4656,17 @@ def LambdifyCSE(args, *exprs, cse=None, concatenate=None, **kwargs):
         cse_symbs, cse_exprs = zip(*subs)
         new_exprs = []
         n_taken = 0
-        for expr in exprs:
-            shape = get_shape(expr) or (1,)
-            size = long(reduce(mul, shape))
-            if len(shape) == 1:
-                new_exprs.append(flat_new_exprs[n_taken:n_taken+size])
-            elif len(shape) == 2:
-                new_exprs.append(DenseMatrix(
-                    shape[0], shape[1], flat_new_exprs[n_taken:n_taken+size]))
-            else:
-                raise NotImplementedError("n-dimensional output not yet supported.")
-            n_taken += size
-        lmb = Lambdify(tuple(args) + cse_symbs, *new_exprs, **kwargs)
-        cse_lambda = Lambdify(args, [expr.xreplace(explicit_subs) for expr in cse_exprs], **kwargs)
+        for expr in _exprs:
+            new_exprs.append(np.reshape(flat_new_exprs[n_taken:n_taken+expr.size],
+                                        expr.shape, order=order))
+            n_taken += expr.size
+        lmb = Lambdify(tuple(args) + cse_symbs, *new_exprs, order=order, **kwargs)
+        cse_lambda = Lambdify(args, [expr.xreplace(explicit_subs) for expr in cse_exprs],
+                              **kwargs)
         def cb(inp, *, out=None, **kw):
             cse_vals = cse_lambda(inp, **kw)
-            print(inp, cse_vals) # DO-NOT-MERGE!
-            new_inp = concatenate((inp, cse_vals), axis=-1)
+            new_inp = np.concatenate((inp, cse_vals), axis=-1)
             return lmb(new_inp, out=out, **kw)
-
         return cb
     else:
         return Lambdify(args, *exprs, **kwargs)
@@ -4724,7 +4727,7 @@ def contains(expr, sset):
     cdef Set sset_ = sympify(sset)
     cdef RCP[const symengine.Set] s = symengine.rcp_static_cast_Set(sset_.thisptr)
     return c2py(<RCP[const symengine.Basic]>(symengine.contains(expr_.thisptr, s)))
-    
+
 
 def set_union(*args):
     cdef symengine.set_set s
diff --git a/symengine/tests/test_lambdify.py b/symengine/tests/test_lambdify.py
index a5cfade34..0392f5529 100644
--- a/symengine/tests/test_lambdify.py
+++ b/symengine/tests/test_lambdify.py
@@ -4,7 +4,9 @@
 
 import array
 import cmath
+from functools import reduce
 import itertools
+from operator import mul
 import math
 import sys
 
@@ -49,27 +51,13 @@ def allclose(vec1, vec2, rtol=1e-13, atol=1e-13):
             return False
     return True
 
-@unittest.skipUnless(have_numpy, "Numpy not installed")
-def test_get_shape():
-    get_shape = se.lib.symengine_wrapper.get_shape
-    assert get_shape([1]) == (1,)
-    assert get_shape([1, 1, 1]) == (3,)
-    assert get_shape([[1], [1], [1]]) == (3, 1)
-    assert get_shape([[1, 1, 1]]) == (1, 3)
-
-    x = se.symbols('x')
-    exprs = [x+1, x+2, x+3, 1/x, 1/(x*x), 1/(x**3.0)]
-    A = se.DenseMatrix(2, 3, exprs)
-    assert get_shape(A) == (2, 3)
-
 
 @unittest.skipUnless(have_numpy, "Numpy not installed")
 def test_ravel():
     x = se.symbols('x')
-    ravel = se.lib.symengine_wrapper.ravel
     exprs = [x+1, x+2, x+3, 1/x, 1/(x*x), 1/(x**3.0)]
     A = se.DenseMatrix(2, 3, exprs)
-    assert ravel(A) == exprs
+    assert np.all(np.ravel(A, order='C') == exprs)
 
 
 @unittest.skipUnless(have_numpy, "Numpy not installed")
@@ -141,6 +129,7 @@ def test_array():
 @unittest.skipUnless(have_numpy, "Numpy not installed")
 def test_numpy_array_out_exceptions():
     args, exprs, inp, check = _get_array()
+    assert len(args) == 3 and len(exprs) == 2
     lmb = se.Lambdify(args, exprs)
 
     all_right = np.empty(len(exprs))
@@ -156,9 +145,12 @@ def test_numpy_array_out_exceptions():
     read_only.flags['WRITEABLE'] = False
     raises(ValueError, lambda: (lmb(inp, out=read_only)))
 
-    all_right_broadcast = np.empty((4, len(exprs)))
+    all_right_broadcast_C = np.empty((4, len(exprs)), order='C')
     inp_bcast = [[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]]
-    lmb(np.array(inp_bcast), out=all_right_broadcast)
+    lmb(np.array(inp_bcast), out=all_right_broadcast_C, order='C')
+
+    all_right_broadcast_F = np.empty((len(exprs), 4), order='F')
+    lmb(np.array(np.array(inp_bcast).T), out=all_right_broadcast_F, order='F')
 
     noncontig_broadcast = np.empty((4, len(exprs), 3)).transpose((1, 2, 0))
     raises(ValueError, lambda: (lmb(inp_bcast, out=noncontig_broadcast)))
@@ -167,7 +159,8 @@ def test_numpy_array_out_exceptions():
 @unittest.skipUnless(have_numpy, "Numpy not installed")
 def test_broadcast():
     a = np.linspace(-np.pi, np.pi)
-    inp = np.vstack((np.cos(a), np.sin(a))).T  # 50 rows 2 cols
+    inp = np.ascontiguousarray(np.vstack((np.cos(a), np.sin(a))).T)  # 50 rows 2 cols
+    assert inp.flags['C_CONTIGUOUS']
     x, y = se.symbols('x y')
     distance = se.Lambdify([x, y], [se.sqrt(x**2 + y**2)])
     assert np.allclose(distance([inp[0, 0], inp[0, 1]]), [1])
@@ -183,7 +176,8 @@ def test_broadcast_multiple_extra_dimensions():
     cb = se.Lambdify([x], [x**2, x**3])
     assert np.allclose(cb([inp[0, 2]]), [4, 8])
     out = cb(inp)
-    assert out.shape == (4, 3, 2)
+    assert out.shape == (4, 3, 1, 2)
+    out = out.squeeze()
     assert abs(out[2, 1, 0] - 7**2) < 1e-14
     assert abs(out[2, 1, 1] - 7**3) < 1e-14
     assert abs(out[-1, -1, 0] - 11**2) < 1e-14
@@ -248,6 +242,7 @@ def test_cse_big():
 def test_broadcast_c():
     n = 3
     inp = np.arange(2*n).reshape((n, 2))
+    assert inp.flags['C_CONTIGUOUS']
     lmb, check = _get_2_to_2by2()
     A = lmb(inp)
     assert A.shape == (3, 2, 2)
@@ -552,7 +547,6 @@ def _test_Lambdify_scalar_vector_matrix(Lambdify):
     f = Lambdify(args, x**y, vec, jac)
     assert f.n_exprs == 3
     s, v, m = f([2, 3])
-    print(s, v, m)
     assert s == 2**3
     assert np.allclose(v, [[2+3], [2*3]])
     assert np.allclose(m, [
@@ -560,22 +554,23 @@ def _test_Lambdify_scalar_vector_matrix(Lambdify):
         [3, 2]
     ])
 
-    s2, v2, m2 = f([[2, 3], [5, 7]])
-    assert np.allclose(s2, [2**3, 5**7])
-    assert np.allclose(v2, [
-        [[2+3], [2*3]],
-        [[5+7], [5*7]]
-    ])
-    assert np.allclose(m2, [
-        [
-            [1, 1],
-            [3, 2]
-        ],
-        [
-            [1, 1],
-            [7, 5]
-        ]
-    ])
+    for inp in [[2, 3, 5, 7], np.array([[2, 3], [5, 7]])]:
+        s2, v2, m2 = f(inp)
+        assert np.allclose(s2, [2**3, 5**7])
+        assert np.allclose(v2, [
+            [[2+3], [2*3]],
+            [[5+7], [5*7]]
+        ])
+        assert np.allclose(m2, [
+            [
+                [1, 1],
+                [3, 2]
+            ],
+            [
+                [1, 1],
+                [7, 5]
+            ]
+        ])
 
 
 def test_Lambdify_scalar_vector_matrix():
@@ -589,3 +584,188 @@ def test_Lambdify_scalar_vector_matrix_cse():
     _test_Lambdify_scalar_vector_matrix(lambda *args: se.LambdifyCSE(*args, backend='lambda'))
     if se.have_llvm:
         _test_Lambdify_scalar_vector_matrix(lambda *args: se.LambdifyCSE(*args, backend='llvm'))
+
+
+@unittest.skipUnless(have_numpy, "Numpy not installed")
+def test_Lambdify_gh174():
+    # Tests array broadcasting if the expressions form an N-dimensional array
+    # of say shape (k, l, m) and it contains 'n' arguments (x1, ... xn), then
+    # if the user provides a Fortran ordered (column-major) input array of shape
+    # (n, o, p, q), then the returned array will be of shape (k, l, m, o, p, q)
+    args = x, y = se.symbols('x y')
+    nargs = len(args)
+    vec1 = se.DenseMatrix([x, x**2, x**3])
+    assert vec1.shape == (3, 1)
+    assert np.asarray(vec1).shape == (3, 1)
+    lmb1 = se.Lambdify([x], vec1)
+    out1 = lmb1(3)
+    assert out1.shape == (3, 1)
+    assert np.all(out1 == [[3], [9], [27]])
+    out1a = lmb1([2, 3], order='F')  # another dimension
+    assert out1a.shape == (3, 1, 2)
+    ref1a_squeeze = [[2, 3],
+                     [4, 9],
+                     [8, 27]]
+    assert np.all(out1a.squeeze() == ref1a_squeeze)
+    assert out1a.flags['F_CONTIGUOUS']
+    assert not out1a.flags['C_CONTIGUOUS']
+
+    lmb2c = se.Lambdify(args, vec1, x+y, order='C')
+    lmb2f = se.Lambdify(args, vec1, x+y, order='F')
+    for out2a in [lmb2c([2, 3]), lmb2f([2, 3])]:
+        assert np.all(out2a[0] == [[2], [4], [8]])
+        assert out2a[0].ndim == 2
+        assert out2a[1] == 5
+        assert out2a[1].ndim == 0
+    inp2b = np.array([
+        [2.0, 3.0],
+        [1.0, 2.0],
+        [0.0, 6.0]
+    ])
+    raises(ValueError, lambda: (lmb2c(inp2b.T)))
+    out2c = lmb2c(inp2b)
+    out2f = lmb2f(np.asfortranarray(inp2b.T))
+    assert out2c[0].shape == (3, 3, 1)
+    assert out2f[0].shape == (3, 1, 3)
+    for idx, (_x, _y) in enumerate(inp2b):
+        assert np.all(out2c[0][idx, ...] == [[_x], [_x**2], [_x**3]])
+
+    assert np.all(out2c[1] == [5, 3, 6])
+    assert np.all(out2f[1] == [5, 3, 6])
+    assert out2c[1].shape == (3,)
+    assert out2f[1].shape == (3,)
+
+    def _mtx3(_x, _y):
+        return [[_x**row_idx + _y**col_idx for col_idx in range(3)]
+                for row_idx in range(4)]
+    mtx3c = np.array(_mtx3(x, y), order='C')
+    mtx3f = np.array(_mtx3(x, y), order='F')
+    lmb3c = se.Lambdify([x, y], x*y, mtx3c, vec1, order='C')
+    lmb3f = se.Lambdify([x, y], x*y, mtx3f, vec1, order='F')
+    inp3c = np.array([[2., 3], [3, 4], [5, 7], [6, 2], [3, 1]])
+    inp3f = np.asfortranarray(inp3c.T)
+    raises(ValueError, lambda: (lmb3c(inp3c.T)))
+    out3c = lmb3c(inp3c)
+    assert out3c[0].shape == (5,)
+    assert out3c[1].shape == (5, 4, 3)
+    assert out3c[2].shape == (5, 3, 1)  # user can apply numpy.squeeze if they want to.
+    for a, b in zip(out3c, lmb3c(np.ravel(inp3c, order='C'))):
+        assert np.all(a == b)
+
+    out3f = lmb3f(inp3f)
+    assert out3f[0].shape == (5,)
+    assert out3f[1].shape == (4, 3, 5)
+    assert out3f[2].shape == (3, 1, 5)  # user can apply numpy.squeeze if they want to.
+    for a, b in zip(out3f, lmb3f(np.ravel(inp3f, order='F'))):
+        assert np.all(a == b)
+
+    for idx, (_x, _y) in enumerate(inp3c):
+        assert out3c[0][idx] == _x*_y
+        assert out3f[0][idx] == _x*_y
+        assert np.all(out3c[1][idx, ...] == _mtx3(_x, _y))
+        assert np.all(out3f[1][..., idx] == _mtx3(_x, _y))
+        assert np.all(out3c[2][idx, ...] == [[_x],[_x**2],[_x**3]])
+        assert np.all(out3f[2][..., idx] == [[_x],[_x**2],[_x**3]])
+
+
+def _get_Ndim_args_exprs_funcs(order):
+    args = x, y = se.symbols('x y')
+
+    # Higher dimensional inputs
+    def f_a(index, _x, _y):
+        a, b, c, d = index
+        return _x**a + _y**b + (_x+_y)**-d
+
+    nd_exprs_a = np.zeros((3, 5, 1, 4), dtype=object, order=order)
+    for index in np.ndindex(*nd_exprs_a.shape):
+        nd_exprs_a[index] = f_a(index, x, y)
+
+    def f_b(index, _x, _y):
+        a, b, c = index
+        return b/(_x + _y)
+
+    nd_exprs_b = np.zeros((1, 7, 1), dtype=object, order=order)
+    for index in np.ndindex(*nd_exprs_b.shape):
+        nd_exprs_b[index] = f_b(index, x, y)
+    return args, nd_exprs_a, nd_exprs_b, f_a, f_b
+
+@unittest.skipUnless(have_numpy, "Numpy not installed")
+def test_Lambdify_Ndimensional_order_C():
+    args, nd_exprs_a, nd_exprs_b, f_a, f_b = _get_Ndim_args_exprs_funcs(order='C')
+    lmb4 = se.Lambdify(args, nd_exprs_a, nd_exprs_b, order='C')
+    nargs = len(args)
+
+    inp_extra_shape = (3, 5, 4)
+    inp_shape = inp_extra_shape + (nargs,)
+    inp4 = np.arange(reduce(mul, inp_shape)*1.0).reshape(inp_shape, order='C')
+    out4a, out4b = lmb4(inp4)
+    assert out4a.ndim == 7
+    assert out4a.shape == inp_extra_shape + nd_exprs_a.shape
+    assert out4b.ndim == 6
+    assert out4b.shape == inp_extra_shape + nd_exprs_b.shape
+    raises(ValueError, lambda: (lmb4(inp4.T)))
+    for b, c, d in np.ndindex(inp_extra_shape):
+        _x, _y = inp4[b, c, d, :]
+        for index in np.ndindex(*nd_exprs_a.shape):
+            assert np.isclose(out4a[(b, c, d) + index], f_a(index, _x, _y))
+        for index in np.ndindex(*nd_exprs_b.shape):
+            assert np.isclose(out4b[(b, c, d) + index], f_b(index, _x, _y))
+
+
+@unittest.skipUnless(have_numpy, "Numpy not installed")
+def test_Lambdify_Ndimensional_order_F():
+    args, nd_exprs_a, nd_exprs_b, f_a, f_b = _get_Ndim_args_exprs_funcs(order='F')
+    lmb4 = se.Lambdify(args, nd_exprs_a, nd_exprs_b, order='F')
+    nargs = len(args)
+
+    inp_extra_shape = (3, 5, 4)
+    inp_shape = (nargs,)+inp_extra_shape
+    inp4 = np.arange(reduce(mul, inp_shape)*1.0).reshape(inp_shape, order='F')
+    out4a, out4b = lmb4(inp4)
+    assert out4a.ndim == 7
+    assert out4a.shape == nd_exprs_a.shape + inp_extra_shape
+    assert out4b.ndim == 6
+    assert out4b.shape == nd_exprs_b.shape + inp_extra_shape
+    raises(ValueError, lambda: (lmb4(inp4.T)))
+    for b, c, d in np.ndindex(inp_extra_shape):
+        _x, _y = inp4[:, b, c, d]
+        for index in np.ndindex(*nd_exprs_a.shape):
+            assert np.isclose(out4a[index + (b, c, d)], f_a(index, _x, _y))
+        for index in np.ndindex(*nd_exprs_b.shape):
+            assert np.isclose(out4b[index + (b, c, d)], f_b(index, _x, _y))
+
+
+@unittest.skipUnless(have_numpy, "Numpy not installed")
+def test_Lambdify_inp_exceptions():
+    args = x, y = se.symbols('x y')
+    lmb1 = se.Lambdify([x], x**2)
+    raises(ValueError, lambda: (lmb1([])))
+    assert lmb1(4) == 16
+    assert np.all(lmb1([4, 2]) == [16, 4])
+
+    lmb2 = se.Lambdify(args, x**2+y**2)
+    assert lmb2([2, 3]) == 13
+    raises(ValueError, lambda: lmb2([]))
+    raises(ValueError, lambda: lmb2([2]))
+    raises(ValueError, lambda: lmb2([2, 3, 4]))
+    assert np.all(lmb2([2, 3, 4, 5]) == [13, 16+25])
+
+    def _mtx(_x, _y):
+        return [
+            [_x-_y, _y**2],
+            [_x+_y, _x**2],
+            [_x*_y, _x**_y]
+        ]
+
+    mtx = np.array(_mtx(x, y), order='F')
+    lmb3 = se.Lambdify(args, mtx, order='F')
+    inp3a = [2, 3]
+    assert np.all(lmb3(inp3a) == _mtx(*inp3a))
+    inp3b = np.array([2, 3, 4, 5, 3, 2, 1, 5])
+    for inp in [inp3b, inp3b.tolist(), inp3b.reshape((2, 4), order='F')]:
+        out3b = lmb3(inp)
+        assert out3b.shape == (3, 2, 4)
+        for i in range(4):
+            assert np.all(out3b[..., i] == _mtx(*inp3b[2*i:2*(i+1)]))
+    raises(ValueError, lambda: lmb3(inp3b.reshape((4, 2), order='F')))
+    raises(ValueError, lambda: lmb3(inp3b.reshape((2, 4)).T))
diff --git a/symengine/tests/test_matrices.py b/symengine/tests/test_matrices.py
index 2aada1dd6..7f4e55a06 100644
--- a/symengine/tests/test_matrices.py
+++ b/symengine/tests/test_matrices.py
@@ -37,6 +37,12 @@ def test_get():
     raises(IndexError, lambda: A.get(-3, 0))
 
 
+def test_tolist():
+    A = DenseMatrix([2, 3])
+    assert A.shape == (2, 1)
+    assert A.tolist() == [[2], [3]]
+
+
 def test_get_item():
     A = DenseMatrix(3, 3, [1, 2, 3, 4, 5, 6, 7, 8, 9])
 
@@ -432,14 +438,14 @@ def test_col_insert():
 def test_rowmul():
     M = ones(3)
     assert M.rowmul(2, 2) == DenseMatrix([[1, 1, 1],
-                                          [1, 1, 1], 
+                                          [1, 1, 1],
                                           [2, 2, 2]])
 
 
 def test_rowadd():
     M = ones(3)
     assert M.rowadd(2, 1, 1) == DenseMatrix([[1, 1, 1],
-                                             [1, 1, 1], 
+                                             [1, 1, 1],
                                              [2, 2, 2]])
 
 

From e949b3f5363bf6fc4dd87f5b06eb078b590905db Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Bj=C3=B6rn=20Dahlgren?= <bjodah@gmail.com>
Date: Tue, 19 Sep 2017 15:09:59 +0200
Subject: [PATCH 2/6] Simplify Lambdify's __init__

---
 symengine/lib/symengine_wrapper.pyx | 43 ++++++++++++-----------------
 1 file changed, 18 insertions(+), 25 deletions(-)

diff --git a/symengine/lib/symengine_wrapper.pyx b/symengine/lib/symengine_wrapper.pyx
index b99832698..8ff1cbf73 100644
--- a/symengine/lib/symengine_wrapper.pyx
+++ b/symengine/lib/symengine_wrapper.pyx
@@ -3264,8 +3264,15 @@ cdef class DenseMatrixBase(MatrixBase):
         def __get__(self):
             return self.nrows()*self.ncols()
 
-    def ravel(self):
-        return [self._get(i, j) for i in range(self.nrows()) for j in range(self.ncols())]
+    def ravel(self, order='C'):
+        if order == 'C':
+            return [self._get(i, j) for i in range(self.nrows())
+                    for j in range(self.ncols())]
+        elif order == 'F':
+            return [self._get(i, j) for j in range(self.ncols())
+                    for i in range(self.nrows())]
+        else:
+            raise NotImplementedError("Unknown order '%s'" % order)
 
     def reshape(self, rows, cols):
         if len(self) != rows*cols:
@@ -4389,28 +4396,14 @@ cdef class _Lambdify(object):
             RCP[const symengine.Basic] b_
             symengine.vec_basic args_, outs_
 
-        if isinstance(args, DenseMatrixBase):
-            nr = args.nrows()
-            nc = args.ncols()
-            mtx = (<DenseMatrixBase>args).thisptr
-            for ri in range(nr):
-                for ci in range(nc):
-                   args_.push_back(deref(mtx).get(ri, ci))
-        else:
-            for arg in np.ravel(args, order=order):
-                e_ = _sympify(arg)
-                args_.push_back(e_.thisptr)
-
-
-        for curr_expr in exprs:
-            if isinstance(curr_expr, DenseMatrixBase):
-                nr = curr_expr.nrows()
-                nc = curr_expr.ncols()
-                mtx = (<DenseMatrixBase>curr_expr).thisptr
-                for ri in range(nr):
-                    for ci in range(nc):
-                        b_ = deref(mtx).get(ri, ci)
-                        outs_.push_back(b_)
+        for arg in np.ravel(args, order=order):
+            e_ = _sympify(arg)
+            args_.push_back(e_.thisptr)
+
+        for curr_expr in map(np.array, exprs):
+            if curr_expr.ndim == 0:
+                e_ = _sympify(curr_expr.item())
+                outs_.push_back(e_.thisptr)
             else:
                 for e in np.ravel(curr_expr, order=order):
                     e_ = _sympify(e)
@@ -4509,7 +4502,7 @@ cdef class _Lambdify(object):
                     extra_dim = inp.shape[1:]
             else:
                 extra_dim = inp.shape
-       else:
+        else:
             if nbroadcast > 1 and inp.ndim == 1:
                 extra_dim = (nbroadcast,)  # special case
             else:

From 93d268299ce82aba919123c78fdfc6ebbe5c43bc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Bj=C3=B6rn=20Dahlgren?= <bjodah@gmail.com>
Date: Tue, 19 Sep 2017 18:57:23 +0200
Subject: [PATCH 3/6] Support broadcasting also for CSE

---
 symengine/lib/symengine_wrapper.pyx | 27 +++++++++++++++++----------
 symengine/tests/test_lambdify.py    | 12 ++++++++++++
 2 files changed, 29 insertions(+), 10 deletions(-)

diff --git a/symengine/lib/symengine_wrapper.pyx b/symengine/lib/symengine_wrapper.pyx
index 8ff1cbf73..18f6f2b18 100644
--- a/symengine/lib/symengine_wrapper.pyx
+++ b/symengine/lib/symengine_wrapper.pyx
@@ -4637,15 +4637,16 @@ def LambdifyCSE(args, *exprs, cse=None, order='C', **kwargs):
     if cse is None:
         from sympy import cse
     _exprs = [np.asanyarray(e) for e in exprs]
+    _args = np.ravel(args, order=order)
     from sympy import sympify as s_sympify
     flat_exprs = list(itertools.chain(*[np.ravel(e, order=order) for e in _exprs]))
     subs, flat_new_exprs = cse([s_sympify(expr) for expr in flat_exprs])
 
-    explicit_subs = {}
-    for k, v in subs:
-        explicit_subs[k] = v.xreplace(explicit_subs)
-
     if subs:
+        explicit_subs = {}
+        for k, v in subs:
+            explicit_subs[k] = v.xreplace(explicit_subs)
+
         cse_symbs, cse_exprs = zip(*subs)
         new_exprs = []
         n_taken = 0
@@ -4653,13 +4654,19 @@ def LambdifyCSE(args, *exprs, cse=None, order='C', **kwargs):
             new_exprs.append(np.reshape(flat_new_exprs[n_taken:n_taken+expr.size],
                                         expr.shape, order=order))
             n_taken += expr.size
-        lmb = Lambdify(tuple(args) + cse_symbs, *new_exprs, order=order, **kwargs)
-        cse_lambda = Lambdify(args, [expr.xreplace(explicit_subs) for expr in cse_exprs],
-                              **kwargs)
+        new_lmb = Lambdify(tuple(_args) + cse_symbs, *new_exprs, order=order, **kwargs)
+        cse_lambda = Lambdify(_args, [ce.xreplace(explicit_subs) for ce in cse_exprs], **kwargs)
         def cb(inp, *, out=None, **kw):
-            cse_vals = cse_lambda(inp, **kw)
-            new_inp = np.concatenate((inp, cse_vals), axis=-1)
-            return lmb(new_inp, out=out, **kw)
+            _order = kw.pop('order', order)
+            _inp = np.asanyarray(inp)
+            cse_vals = cse_lambda(_inp, order=_order, **kw)
+            if order == 'C':
+                new_inp = np.concatenate((_inp[(Ellipsis,) + (np.newaxis,)*(cse_vals.ndim - _inp.ndim)],
+                                          cse_vals), axis=-1)
+            else:
+                new_inp = np.concatenate((_inp[(np.newaxis,)*(cse_vals.ndim - _inp.ndim) + (Ellipsis,)],
+                                          cse_vals), axis=0)
+            return new_lmb(new_inp, out=out, order=_order, **kw)
         return cb
     else:
         return Lambdify(args, *exprs, **kwargs)
diff --git a/symengine/tests/test_lambdify.py b/symengine/tests/test_lambdify.py
index 0392f5529..0b9b7fc32 100644
--- a/symengine/tests/test_lambdify.py
+++ b/symengine/tests/test_lambdify.py
@@ -202,6 +202,18 @@ def test_cse():
     out = lmb(inp)
     assert allclose(out, ref)
 
+
+@unittest.skipUnless(have_numpy, "Numpy not installed")
+@unittest.skipUnless(have_sympy, "SymPy not installed")
+def test_cse_gh174():
+    x = se.symbols('x')
+    funcs = [se.cos(x)**i for i in range(5)]
+    f_lmb = se.Lambdify([x], funcs)
+    f_cse = se.LambdifyCSE([x], funcs)
+    a = np.array([1, 2, 3])
+    assert np.allclose(f_lmb(a), f_cse(a))
+
+
 def _get_cse_exprs_big():
     # this is essentially a performance test (can be replaced by a benchmark)
     x, p = se.symarray('x', 14), se.symarray('p', 14)

From 8376e3a6b3a3749750dc58efc0b673bcccfa3103 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Bj=C3=B6rn=20Dahlgren?= <bjodah@gmail.com>
Date: Thu, 21 Sep 2017 01:09:54 +0200
Subject: [PATCH 4/6] Remove order as kwarg from Lambdify.__call__

---
 symengine/lib/symengine_wrapper.pyx | 51 ++++++++++++++---------------
 symengine/tests/test_lambdify.py    | 18 ++++++----
 2 files changed, 36 insertions(+), 33 deletions(-)

diff --git a/symengine/lib/symengine_wrapper.pyx b/symengine/lib/symengine_wrapper.pyx
index 18f6f2b18..6e62e343e 100644
--- a/symengine/lib/symengine_wrapper.pyx
+++ b/symengine/lib/symengine_wrapper.pyx
@@ -4365,7 +4365,7 @@ cdef class _Lambdify(object):
     cdef list out_shapes
     cdef readonly bint real
     cdef readonly int n_exprs
-    cdef readonly str order
+    cdef public str order
     cdef vector[int] accum_out_sizes
     cdef object numpy_dtype
 
@@ -4440,7 +4440,7 @@ cdef class _Lambdify(object):
             raise ValueError("Size of out incompatible with number of exprs.")
         self.unsafe_complex(inp, out)
 
-    def __call__(self, inp, *, out=None, order=None):
+    def __call__(self, inp, *, out=None):
         """
         Parameters
         ----------
@@ -4469,9 +4469,7 @@ cdef class _Lambdify(object):
             tuple inp_shape
             double[::1] real_out, real_inp
             double complex[::1] cmplx_out, cmplx_inp
-        if order is None:
-            order = self.order
-        if order not in ('C', 'F'):
+        if self.order not in ('C', 'F'):
             raise NotImplementedError("Only C & F order supported for now.")
 
         try:
@@ -4480,9 +4478,9 @@ cdef class _Lambdify(object):
             inp = np.fromiter(inp, dtype=self.numpy_dtype)
 
         if self.real:
-            real_inp = np.ascontiguousarray(inp.ravel(order=order))
+            real_inp = np.ascontiguousarray(inp.ravel(order=self.order))
         else:
-            cmplx_inp = np.ascontiguousarray(inp.ravel(order=order))
+            cmplx_inp = np.ascontiguousarray(inp.ravel(order=self.order))
 
         if inp.size < self.args_size or inp.size % self.args_size != 0:
             raise ValueError("Broadcasting failed (input/arg size mismatch)")
@@ -4490,12 +4488,12 @@ cdef class _Lambdify(object):
 
         if inp.ndim > 1:
             if self.args_size > 1:
-                if order == 'C':
+                if self.order == 'C':
                     if inp.shape[inp.ndim-1] != self.args_size:
                         raise ValueError(("C order implies last dim (%d) == len(args)"
                                           " (%d)") % (inp.shape[inp.ndim-1], self.args_size))
                     extra_dim = inp.shape[:inp.ndim-1]
-                elif order == 'F':
+                elif self.order == 'F':
                     if inp.shape[0] != self.args_size:
                         raise ValueError("F order implies first dim (%d) == len(args) (%d)"
                                          % (inp.shape[0], self.args_size))
@@ -4507,35 +4505,37 @@ cdef class _Lambdify(object):
                 extra_dim = (nbroadcast,)  # special case
             else:
                 extra_dim = ()
-        extra_left = extra_dim if order == 'C' else ()
-        extra_right = () if order == 'C' else extra_dim
+        extra_left = extra_dim if self.order == 'C' else ()
+        extra_right = () if self.order == 'C' else extra_dim
         new_out_shapes = [extra_left + out_shape + extra_right
                           for out_shape in self.out_shapes]
 
         new_tot_out_size = nbroadcast * self.tot_out_size
         if out is None:
-            out = np.empty(new_tot_out_size, dtype=self.numpy_dtype, order=order)
+            out = np.empty(new_tot_out_size, dtype=self.numpy_dtype, order=self.order)
         else:
             if out.size < new_tot_out_size:
                 raise ValueError("Incompatible size of output argument")
             if out.ndim > 1:
-                if order == 'C' and not out.flags['C_CONTIGUOUS']:
-                    raise ValueError("Output argument needs to be C-contiguous")
-                elif order == 'F' and not out.flags['F_CONTIGUOUS']:
-                    raise ValueError("Output argument needs to be F-contiguous")
                 if len(self.out_shapes) > 1:
                     raise ValueError("output array with ndim > 1 assumes one output")
                 out_shape, = self.out_shapes
-                if order == 'C' and out.shape[-len(out_shape):] != tuple(out_shape):
-                    raise ValueError("shape mismatch for output array")
-                elif order == 'F' and out.shape[:len(out_shape)] != tuple(out_shape):
-                    raise ValueError("shape mismatch for output array")
+                if self.order == 'C':
+                    if not out.flags['C_CONTIGUOUS']:
+                        raise ValueError("Output argument needs to be C-contiguous")
+                    if out.shape[-len(out_shape):] != tuple(out_shape):
+                        raise ValueError("shape mismatch for output array")
+                elif self.order == 'F':
+                    if not out.flags['F_CONTIGUOUS']:
+                        raise ValueError("Output argument needs to be F-contiguous")
+                    if out.shape[:len(out_shape)] != tuple(out_shape):
+                        raise ValueError("shape mismatch for output array")
             else:
                 if not out.flags['F_CONTIGUOUS']:  # or C_CONTIGUOUS (ndim <= 1)
                     raise ValueError("Output array need to be contiguous")
             if not out.flags['WRITEABLE']:
                 raise ValueError("Output argument needs to be writeable")
-            out = out.ravel(order=order)
+            out = out.ravel(order=self.order)
 
         if self.real:
             real_out = out
@@ -4551,13 +4551,13 @@ cdef class _Lambdify(object):
                 self.unsafe_complex(cmplx_inp, cmplx_out,
                                     idx*self.args_size, idx*self.tot_out_size)
 
-        if order == 'C':
+        if self.order == 'C':
             out = out.reshape((nbroadcast, self.tot_out_size), order='C')
             result = [
                 out[:, self.accum_out_sizes[idx]:self.accum_out_sizes[idx+1]].reshape(
                     new_out_shapes[idx], order='C') for idx in range(self.n_exprs)
             ]
-        elif order == 'F':
+        elif self.order == 'F':
             out = out.reshape((self.tot_out_size, nbroadcast), order='F')
             result = [
                 out[self.accum_out_sizes[idx]:self.accum_out_sizes[idx+1], :].reshape(
@@ -4657,16 +4657,15 @@ def LambdifyCSE(args, *exprs, cse=None, order='C', **kwargs):
         new_lmb = Lambdify(tuple(_args) + cse_symbs, *new_exprs, order=order, **kwargs)
         cse_lambda = Lambdify(_args, [ce.xreplace(explicit_subs) for ce in cse_exprs], **kwargs)
         def cb(inp, *, out=None, **kw):
-            _order = kw.pop('order', order)
             _inp = np.asanyarray(inp)
-            cse_vals = cse_lambda(_inp, order=_order, **kw)
+            cse_vals = cse_lambda(_inp, **kw)
             if order == 'C':
                 new_inp = np.concatenate((_inp[(Ellipsis,) + (np.newaxis,)*(cse_vals.ndim - _inp.ndim)],
                                           cse_vals), axis=-1)
             else:
                 new_inp = np.concatenate((_inp[(np.newaxis,)*(cse_vals.ndim - _inp.ndim) + (Ellipsis,)],
                                           cse_vals), axis=0)
-            return new_lmb(new_inp, out=out, order=_order, **kw)
+            return new_lmb(new_inp, out=out, **kw)
         return cb
     else:
         return Lambdify(args, *exprs, **kwargs)
diff --git a/symengine/tests/test_lambdify.py b/symengine/tests/test_lambdify.py
index 0b9b7fc32..bd379cfd6 100644
--- a/symengine/tests/test_lambdify.py
+++ b/symengine/tests/test_lambdify.py
@@ -147,14 +147,16 @@ def test_numpy_array_out_exceptions():
 
     all_right_broadcast_C = np.empty((4, len(exprs)), order='C')
     inp_bcast = [[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]]
-    lmb(np.array(inp_bcast), out=all_right_broadcast_C, order='C')
-
-    all_right_broadcast_F = np.empty((len(exprs), 4), order='F')
-    lmb(np.array(np.array(inp_bcast).T), out=all_right_broadcast_F, order='F')
+    lmb(np.array(inp_bcast), out=all_right_broadcast_C)
 
     noncontig_broadcast = np.empty((4, len(exprs), 3)).transpose((1, 2, 0))
     raises(ValueError, lambda: (lmb(inp_bcast, out=noncontig_broadcast)))
 
+    all_right_broadcast_F = np.empty((len(exprs), 4), order='F')
+    lmb.order = 'F'
+    lmb(np.array(np.array(inp_bcast).T), out=all_right_broadcast_F)
+
+
 
 @unittest.skipUnless(have_numpy, "Numpy not installed")
 def test_broadcast():
@@ -613,7 +615,9 @@ def test_Lambdify_gh174():
     out1 = lmb1(3)
     assert out1.shape == (3, 1)
     assert np.all(out1 == [[3], [9], [27]])
-    out1a = lmb1([2, 3], order='F')  # another dimension
+    assert lmb1([2, 3]).shape == (2, 3, 1)
+    lmb1.order = 'F'  # change order
+    out1a = lmb1([2, 3])
     assert out1a.shape == (3, 1, 2)
     ref1a_squeeze = [[2, 3],
                      [4, 9],
@@ -661,7 +665,7 @@ def _mtx3(_x, _y):
     assert out3c[0].shape == (5,)
     assert out3c[1].shape == (5, 4, 3)
     assert out3c[2].shape == (5, 3, 1)  # user can apply numpy.squeeze if they want to.
-    for a, b in zip(out3c, lmb3c(np.ravel(inp3c, order='C'))):
+    for a, b in zip(out3c, lmb3c(np.ravel(inp3c))):
         assert np.all(a == b)
 
     out3f = lmb3f(inp3f)
@@ -779,5 +783,5 @@ def _mtx(_x, _y):
         assert out3b.shape == (3, 2, 4)
         for i in range(4):
             assert np.all(out3b[..., i] == _mtx(*inp3b[2*i:2*(i+1)]))
-    raises(ValueError, lambda: lmb3(inp3b.reshape((4, 2), order='F')))
+    raises(ValueError, lambda: lmb3(inp3b.reshape((4, 2))))
     raises(ValueError, lambda: lmb3(inp3b.reshape((2, 4)).T))

From 248ac70db2a873298f7324339283058b1b1e34fd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Bj=C3=B6rn=20Dahlgren?= <bjodah@gmail.com>
Date: Thu, 21 Sep 2017 10:00:36 +0200
Subject: [PATCH 5/6] Make sure np.asanyarray is applied least number of times

---
 symengine/lib/symengine_wrapper.pyx | 101 ++++++++++++++--------------
 1 file changed, 52 insertions(+), 49 deletions(-)

diff --git a/symengine/lib/symengine_wrapper.pyx b/symengine/lib/symengine_wrapper.pyx
index 6e62e343e..ab6159d94 100644
--- a/symengine/lib/symengine_wrapper.pyx
+++ b/symengine/lib/symengine_wrapper.pyx
@@ -4331,36 +4331,6 @@ def has_symbol(obj, symbol=None):
 
 
 cdef class _Lambdify(object):
-    """
-    Lambdify instances are callbacks that numerically evaluate their symbolic
-    expressions from user provided input (real or complex) into (possibly user
-    provided) output buffers (real or complex). Multidimensional data are
-    processed in their most cache-friendly way (i.e. "ravelled").
-
-    Parameters
-    ----------
-    args: iterable of Symbols
-    \*exprs: array_like of expressions
-        the shape of exprs is preserved
-    real : bool
-        Whether datatype is ``double`` (``double complex`` otherwise).
-
-    Returns
-    -------
-    callback instance with signature f(inp, out=None)
-
-    Examples
-    --------
-    >>> from symengine import var, Lambdify
-    >>> var('x y z')
-    >>> f = Lambdify([x, y, z], [x+y+z, x*y*z])
-    >>> f([2, 3, 4])
-    [ 9., 24.]
-    >>> out = np.array(2)
-    >>> f(x, out); out
-    [ 9., 24.]
-
-    """
     cdef size_t args_size, tot_out_size
     cdef list out_shapes
     cdef readonly bint real
@@ -4368,15 +4338,19 @@ cdef class _Lambdify(object):
     cdef public str order
     cdef vector[int] accum_out_sizes
     cdef object numpy_dtype
+    cdef args
+    cdef tuple exprs
 
     def __cinit__(self, args, *exprs, cppbool real=True, order='C'):
         cdef vector[int] out_sizes
+        self.args = np.asanyarray(args)
+        self.args_size = self.args.size
+        self.exprs = tuple(np.asanyarray(expr) for expr in exprs)
+        self.out_shapes = [expr.shape for expr in self.exprs]
+        self.n_exprs = len(self.exprs)
         self.real = real
         self.order = order
         self.numpy_dtype = np.float64 if self.real else np.complex128
-        self.out_shapes = [np.asarray(expr).shape for expr in exprs]
-        self.n_exprs = len(exprs)
-        self.args_size = np.asarray(args).size
         if self.args_size == 0:
             raise NotImplementedError("Support for zero arguments not yet supported")
         self.tot_out_size = 0
@@ -4388,7 +4362,7 @@ cdef class _Lambdify(object):
             for j in range(i):
                 self.accum_out_sizes[i] += out_sizes[j]
 
-    def __init__(self, args, *exprs, cppbool real=True, order='C'):
+    def __init__(self, *args, **kwargs):
         cdef:
             Basic e_
             size_t ri, ci, nr, nc
@@ -4396,16 +4370,16 @@ cdef class _Lambdify(object):
             RCP[const symengine.Basic] b_
             symengine.vec_basic args_, outs_
 
-        for arg in np.ravel(args, order=order):
+        for arg in np.ravel(self.args, order=self.order):
             e_ = _sympify(arg)
             args_.push_back(e_.thisptr)
 
-        for curr_expr in map(np.array, exprs):
+        for curr_expr in self.exprs:
             if curr_expr.ndim == 0:
                 e_ = _sympify(curr_expr.item())
                 outs_.push_back(e_.thisptr)
             else:
-                for e in np.ravel(curr_expr, order=order):
+                for e in np.ravel(curr_expr, order=self.order):
                     e_ = _sympify(e)
                     outs_.push_back(e_.thisptr)
         self._init(args_, outs_)
@@ -4422,18 +4396,14 @@ cdef class _Lambdify(object):
                          int inp_offset=0, int out_offset=0):
         raise ValueError("Not supported")
 
-    cpdef eval_real(self,
-                    inp,
-                    out):
+    cpdef eval_real(self, inp, out):
         if inp.size != self.args_size:
             raise ValueError("Size of inp incompatible with number of args.")
         if out.size != self.tot_out_size:
             raise ValueError("Size of out incompatible with number of exprs.")
         self.unsafe_real(inp, out)
 
-    cpdef eval_complex(self,
-                       inp,
-                       out):
+    cpdef eval_complex(self, inp, out):
         if inp.size != self.args_size:
             raise ValueError("Size of inp incompatible with number of args.")
         if out.size != self.tot_out_size:
@@ -4451,12 +4421,6 @@ cdef class _Lambdify(object):
             If ``None``: an output container will be allocated (NumPy ndarray).
             If ``len(exprs) > 0`` output is found in the corresponding
             order.
-        order : 'C' or 'F'
-            C- or Fortran-contiguous memory layout. Note that this affects
-            broadcasting: e.g. a (m, n) matrix taking 3 arguments and given a
-            (k, l, 3) (C-contiguous) input will give a (k, l, m, n) shaped output,
-            whereas a (3, k, l) (C-contiguous) input will give a (m, n, k, l) shaped
-            output. If ``None`` order is taken as ``self.order`` (from initialization).
 
         Returns
         -------
@@ -4603,6 +4567,45 @@ IF HAVE_SYMENGINE_LLVM:
 
 
 def Lambdify(args, *exprs, cppbool real=True, backend=None, order='C'):
+    """
+    Lambdify instances are callbacks that numerically evaluate their symbolic
+    expressions from user provided input (real or complex) into (possibly user
+    provided) output buffers (real or complex). Multidimensional data are
+    processed in their most cache-friendly way (i.e. "ravelled").
+
+    Parameters
+    ----------
+    args: iterable of Symbols
+    \*exprs: array_like of expressions
+        the shape of exprs is preserved
+    real : bool
+        Whether datatype is ``double`` (``double complex`` otherwise).
+    backend : str
+        'llvm' or 'lambda'. When ``None`` the environment variable
+        'SYMENGINE_LAMBDIFY_BACKEND' is used (taken as 'lambda' if unset).
+    order : 'C' or 'F'
+        C- or Fortran-contiguous memory layout. Note that this affects
+        broadcasting: e.g. a (m, n) matrix taking 3 arguments and given a
+        (k, l, 3) (C-contiguous) input will give a (k, l, m, n) shaped output,
+        whereas a (3, k, l) (C-contiguous) input will give a (m, n, k, l) shaped
+        output. If ``None`` order is taken as ``self.order`` (from initialization).
+
+    Returns
+    -------
+    callback instance with signature f(inp, out=None)
+
+    Examples
+    --------
+    >>> from symengine import var, Lambdify
+    >>> var('x y z')
+    >>> f = Lambdify([x, y, z], [x+y+z, x*y*z])
+    >>> f([2, 3, 4])
+    [ 9., 24.]
+    >>> out = np.array(2)
+    >>> f(x, out); out
+    [ 9., 24.]
+
+    """
     if backend is None:
         backend = os.getenv('SYMENGINE_LAMBDIFY_BACKEND', "lambda")
     if backend == "llvm":

From 21f950d4c6a6e954f9ae273825eca8992e629fcd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Bj=C3=B6rn=20Dahlgren?= <bjodah@gmail.com>
Date: Fri, 22 Sep 2017 18:33:03 +0200
Subject: [PATCH 6/6] Get rid of __cinit__ in lambdify

---
 symengine/lib/symengine_wrapper.pyx | 37 +++++++++++++----------------
 1 file changed, 17 insertions(+), 20 deletions(-)

diff --git a/symengine/lib/symengine_wrapper.pyx b/symengine/lib/symengine_wrapper.pyx
index ab6159d94..3a97c58c9 100644
--- a/symengine/lib/symengine_wrapper.pyx
+++ b/symengine/lib/symengine_wrapper.pyx
@@ -4338,16 +4338,21 @@ cdef class _Lambdify(object):
     cdef public str order
     cdef vector[int] accum_out_sizes
     cdef object numpy_dtype
-    cdef args
-    cdef tuple exprs
-
-    def __cinit__(self, args, *exprs, cppbool real=True, order='C'):
-        cdef vector[int] out_sizes
-        self.args = np.asanyarray(args)
-        self.args_size = self.args.size
-        self.exprs = tuple(np.asanyarray(expr) for expr in exprs)
-        self.out_shapes = [expr.shape for expr in self.exprs]
-        self.n_exprs = len(self.exprs)
+
+    def __init__(self, args, *exprs, cppbool real=True, order='C'):
+        cdef:
+            Basic e_
+            size_t ri, ci, nr, nc
+            symengine.MatrixBase *mtx
+            RCP[const symengine.Basic] b_
+            symengine.vec_basic args_, outs_
+            vector[int] out_sizes
+
+        args = np.asanyarray(args)
+        self.args_size = args.size
+        exprs = tuple(np.asanyarray(expr) for expr in exprs)
+        self.out_shapes = [expr.shape for expr in exprs]
+        self.n_exprs = len(exprs)
         self.real = real
         self.order = order
         self.numpy_dtype = np.float64 if self.real else np.complex128
@@ -4362,19 +4367,11 @@ cdef class _Lambdify(object):
             for j in range(i):
                 self.accum_out_sizes[i] += out_sizes[j]
 
-    def __init__(self, *args, **kwargs):
-        cdef:
-            Basic e_
-            size_t ri, ci, nr, nc
-            symengine.MatrixBase *mtx
-            RCP[const symengine.Basic] b_
-            symengine.vec_basic args_, outs_
-
-        for arg in np.ravel(self.args, order=self.order):
+        for arg in np.ravel(args, order=self.order):
             e_ = _sympify(arg)
             args_.push_back(e_.thisptr)
 
-        for curr_expr in self.exprs:
+        for curr_expr in exprs:
             if curr_expr.ndim == 0:
                 e_ = _sympify(curr_expr.item())
                 outs_.push_back(e_.thisptr)