sympify and matrices
Robert Cimrman
In [16]: sympy.__version__
Out[16]: '0.5.15-hg'
In [19]: K = sympy.Matrix(nm.array(1, ndmin = 2))
In [20]: sympy.sympify(K)
...
NotImplementedError: matrix support
Is there a reason for this? Why a Matrix instance is not just passed
along, since it is already a sympy object?
Then, what is the preferred way of dealing with symbolic matrices w.r.t.
their usability in sympy functions:
- use regular lists (numpy arrays would do too?) and call sympify() with
sympify_lists=True at the beginning of a function accepting a
matrix-like argument
- use Matrix
Ideally, both should work the same way(?)
r.
Ondrej Certik
>
> Hi,
>
> In [16]: sympy.__version__
> Out[16]: '0.5.15-hg'
> In [19]: K = sympy.Matrix(nm.array(1, ndmin = 2))
> In [20]: sympy.sympify(K)
> ...
> NotImplementedError: matrix support
>
> Is there a reason for this? Why a Matrix instance is not just passed
> along, since it is already a sympy object?
The patch fixing this is trivial:
$ hg di
diff --git a/sympy/core/sympify.py b/sympy/core/sympify.py
--- a/sympy/core/sympify.py
+++ b/sympy/core/sympify.py
@@ -107,7 +107,7 @@ def sympify(a, sympify_lists=False, loca
if isinstance(a, Polynomial):
return a
if isinstance(a, Matrix):
- raise NotImplementedError('matrix support')
+ return a
if not isinstance(a, str):
# At this point we were given an arbitrary expression
Is anyone against committing it? All tests pass.
>
> Then, what is the preferred way of dealing with symbolic matrices w.r.t.
> their usability in sympy functions:
> - use regular lists (numpy arrays would do too?) and call sympify() with
> sympify_lists=True at the beginning of a function accepting a
> matrix-like argument
> - use Matrix
> Ideally, both should work the same way(?)
It depends what you want to do.
- for simple stuff I suggest to just use lists
- if you need linear algebra (which I think you almost always need),
I'd suggest Matrices
- in some special (numerical) cases numpy arrays could be used too
I just discovered we probably need to think about this:
In [1]: a = Matrix((1, x), (2, y))
In [2]: from numpy import array
In [3]: a
Out[3]:
⎡1 x⎤
⎣2 y⎦
In [4]: Matrix(sin(array(a)))
Out[4]:
⎡sin(1) sin(x)⎤
⎣sin(2) sin(y)⎦
In [5]: sin(a)
[...]
AttributeError: 'Matrix' object has no attribute 'is_Number'
What should [5] return? The same as [4] or the sin of a Matrix, which
is defined as sin(A) = Q*sin(D)*Q^T if A = Q*D*Q^T and D is
diagonal, so sin(D) is just the D with sin applied to the diagonal
terms.
Ondrej
Ondrej Certik
numpy must have exactly the same problem btw. So let's do the same as
numpy to be compatible, e.g. I am for returning [4]. And let's provide
some methods for returning [5] if needed.
Ondrej
Robert Cimrman
> I just discovered we probably need to think about this:
>
> In [1]: a = Matrix((1, x), (2, y))
>
> In [2]: from numpy import array
>
> In [3]: a
> Out[3]:
> ⎡1 x⎤
> ⎣2 y⎦
>
> In [4]: Matrix(sin(array(a)))
> Out[4]:
> ⎡sin(1) sin(x)⎤
> ⎣sin(2) sin(y)⎦
>
> In [5]: sin(a)
> [...]
>
> AttributeError: 'Matrix' object has no attribute 'is_Number'
>
>
> What should [5] return? The same as [4] or the sin of a Matrix, which
> is defined as sin(A) = Q*sin(D)*Q^T if A = Q*D*Q^T and D is
> diagonal, so sin(D) is just the D with sin applied to the diagonal
> terms.
This stuff generates neverending discussions on numpy lists. Numpy
currently solves this by having array (n-dimensional container) and
matrix (special 2d array w.r.t. multiplication and some functions). As
Matrix in SymPy is a linear algebra matrix, it should IMHO behave as
such, i.e. sin(A) = Q*sin(D)*Q^T. Then maybe an Array would be usefull,
something like sympy version of list of lists [of lists...].
r.
Ondrej Certik
>
> Ondrej Certik wrote:
>> I just discovered we probably need to think about this:
>>
>> In [1]: a = Matrix((1, x), (2, y))
>>
>> In [2]: from numpy import array
>>
>> In [3]: a
>> Out[3]:
>> ⎡1 x⎤
>> ⎣2 y⎦
>>
>> In [4]: Matrix(sin(array(a)))
>> Out[4]:
>> ⎡sin(1) sin(x)⎤
>> ⎣sin(2) sin(y)⎦
>>
>> In [5]: sin(a)
>> [...]
>>
>> AttributeError: 'Matrix' object has no attribute 'is_Number'
>>
>>
>> What should [5] return? The same as [4] or the sin of a Matrix, which
>> is defined as sin(A) = Q*sin(D)*Q^T if A = Q*D*Q^T and D is
>> diagonal, so sin(D) is just the D with sin applied to the diagonal
>> terms.
>
> This stuff generates neverending discussions on numpy lists. Numpy
> currently solves this by having array (n-dimensional container) and
> matrix (special 2d array w.r.t. multiplication and some functions). As
Well, not for sin:
In [8]: from numpy import matrix, array, sin
In [11]: a = matrix([[1,2],[3,0]])
In [12]: sin(a)
Out[12]:
matrix([[ 0.84147098, 0.90929743],
[ 0.14112001, 0. ]])
In [13]: sin(array(a))
Out[13]:
array([[ 0.84147098, 0.90929743],
[ 0.14112001, 0. ]])
So that's not how I would like this to be. It should either work for
all functions or for no functions.
> Matrix in SymPy is a linear algebra matrix, it should IMHO behave as
> such, i.e. sin(A) = Q*sin(D)*Q^T. Then maybe an Array would be usefull,
> something like sympy version of list of lists [of lists...].
Yes, I think so too. Plus we could create some pure python arrray
implementation, fully compatible with numpy.
Ondrej
Robert Cimrman
I see. I do not use numpy matrices anyway, just arrays, as they are not
really the linear algebra matrices, only arrays with some peculiar
modifications - they have no extra value that arrays have not. SymPy
Matrix could differ in this respect, as otherwise the Matrix is just a
misnomed Array. But my voting power is zero here :)
>> Matrix in SymPy is a linear algebra matrix, it should IMHO behave as
>> such, i.e. sin(A) = Q*sin(D)*Q^T. Then maybe an Array would be usefull,
>> something like sympy version of list of lists [of lists...].
>
> Yes, I think so too. Plus we could create some pure python arrray
> implementation, fully compatible with numpy.
But internally, SymPy should stick to one, otherwise checking all the
argument combinations would be a nightmare, no?
r.
Ondrej Certik
>> all functions or for no functions.
>
> I see. I do not use numpy matrices anyway, just arrays, as they are not
> really the linear algebra matrices, only arrays with some peculiar
> modifications - they have no extra value that arrays have not. SymPy
> Matrix could differ in this respect, as otherwise the Matrix is just a
> misnomed Array. But my voting power is zero here :)
We listen to our users and actually we'd like our users to also become
developers. :)
>
>>> Matrix in SymPy is a linear algebra matrix, it should IMHO behave as
>>> such, i.e. sin(A) = Q*sin(D)*Q^T. Then maybe an Array would be usefull,
>>> something like sympy version of list of lists [of lists...].
>>
>> Yes, I think so too. Plus we could create some pure python arrray
>> implementation, fully compatible with numpy.
>
> But internally, SymPy should stick to one, otherwise checking all the
> argument combinations would be a nightmare, no?
I don't think it will be particularly difficult, especially with a
good test suite. It's just a matter of what we want.
Ondrej
Ondrej Certik
http://projects.scipy.org/pipermail/numpy-discussion/2008-June/034801.html
The ultimate reason that Matrix is not a subclass of a Basic is that
it is mutable, while SymPy objects need to be immutable, so that when
you have an expression, like (1+x)**2, you can be sure that all the
instances (like "x") doesn't change in the middle of a calculation.
Q: Why is a Matrix mutable?
A: Because we want to play with matrices using the natural syntax:
A[0,0] = x**2, etc.
Q: Do we want to use Matrices in expressions?
A: I think it could be very useful. Currently we can use Symbol("A",
commutative=False), but it'd be nice to have an option to "attach"
some matrix for it, that represents it. Also it will make the whole
SymPy cleaner, if we adapt for example this approach:
Make Matrix immutable, while introducing an Array class, that will
behave like numpy array, be mutable etc., and a 2D array could be
converted to Matrix (and back of course), thus Matrices would become
regular SymPy objects.
Disadvantage -- if one needs to set entries of a Matrix using the
syntax "A[1, 2] = y" not only at the beginning of the calculaton, but
also in the middle, then he would have to do: A = Matrix(Array(A)[1,2]
= y). Of course A[1,2] = y would create an exception saying: "use the
A = Matrix(Array(A)[1,2] = y)" syntax.
What do you think?
Ondrej
Ondrej Certik
> syntax "A[1, 2] = y" not only at the beginning of the calculaton, but
> also in the middle, then he would have to do: A = Matrix(Array(A)[1,2]
> = y). Of course A[1,2] = y would create an exception saying: "use the
> A = Matrix(Array(A)[1,2] = y)" syntax.
Advantage -- the vectorization will use Arrays/lists, while Matrix
will be reserved for strictly mathematical thing.
BTW, it's exactly like the tuple vs list in Python.
Ondrej
Robert Cimrman
> Let's clean this up. I sent an email to the numpy-list asking for an advice:
>
> http://projects.scipy.org/pipermail/numpy-discussion/2008-June/034801.html
I agree with Robert Kern's view that sympy.Matrix should not be tied to
numpy.matrix - the latter in numpy has a limited use basically in cases
one is lazy to write dot( A.T, B ) instead of A * B (:-}). Sympy is
about math, so sympy.Matrix could become a real matrix in the
mathematical sense.
> The ultimate reason that Matrix is not a subclass of a Basic is that
> it is mutable, while SymPy objects need to be immutable, so that when
> you have an expression, like (1+x)**2, you can be sure that all the
> instances (like "x") doesn't change in the middle of a calculation.
>
> Q: Why is a Matrix mutable?
> A: Because we want to play with matrices using the natural syntax:
> A[0,0] = x**2, etc.
>
> Q: Do we want to use Matrices in expressions?
> A: I think it could be very useful. Currently we can use Symbol("A",
> commutative=False), but it'd be nice to have an option to "attach"
> some matrix for it, that represents it. Also it will make the whole
> SymPy cleaner, if we adapt for example this approach:
>
> Make Matrix immutable, while introducing an Array class, that will
> behave like numpy array, be mutable etc., and a 2D array could be
> converted to Matrix (and back of course), thus Matrices would become
> regular SymPy objects.
>
> Disadvantage -- if one needs to set entries of a Matrix using the
> syntax "A[1, 2] = y" not only at the beginning of the calculaton, but
> also in the middle, then he would have to do: A = Matrix(Array(A)[1,2]
> = y). Of course A[1,2] = y would create an exception saying: "use the
> A = Matrix(Array(A)[1,2] = y)" syntax.
By the "beginning of the calculation" you mean that the assignement
A[i,j] = value will work only once after an empty A is created?
> What do you think?
Fine!
r.
Pearu Peterson
Just a quick note that could be relevant to the discussion:
check also out
http://code.google.com/p/sympycore/wiki/MatrixSupportIdeas
that contains ideas how to deal with mutable matrices in
operations where they should be immutable, all in a very efficient
way.
The idea is based on using views of matrices.
Pearu
On Jun 7, 7:52 pm, "Ondrej Certik" <ond...@certik.cz> wrote:
> Let's clean this up. I sent an email to the numpy-list asking for an advice:
>
Ondrej Certik
<pearu.p...@gmail.com> wrote:
>
> Hi,
>
> Just a quick note that could be relevant to the discussion:
> check also out
>
> http://code.google.com/p/sympycore/wiki/MatrixSupportIdeas
>
> that contains ideas how to deal with mutable matrices in
> operations where they should be immutable, all in a very efficient
> way.
> The idea is based on using views of matrices.
Thanks for the tip. So what you do is that Matrix is writeable until
it is used in an expression? And you use this flag:
"
Matrices are mutable until .is_writable==True.
"
You seem you implemented some interesting ideas in sympycore -- would
you be interested in more collaboration and/or even merging the two
codebases? At least you write on sympycore webpage that "it is created
to fix SymPy performance and robustness issues", so let's point out
the robustness issues (I think we fixed all of them already, but if
there are some more left, let's fix that too) and as to speed, I'd be
interested in your opinions how to move forward. I'd be interested in
coming together for example before the EuroScipy sprint and work on it
together. But I think it'd be efficient to prepare some plan how to
merge it, or what to do.
Anyone interested in this, please look how things are done in
sympycore and in sympy and let's start a discussion which way is good
to do things and implement it. And if there are two good ways (which
is very well possible), then let's find a solution to have both ways,
but with the same interface.
Ondrej
Robert Cimrman
>> Disadvantage -- if one needs to set entries of a Matrix using the
>> syntax "A[1, 2] = y" not only at the beginning of the calculaton, but
>> also in the middle, then he would have to do: A = Matrix(Array(A)[1,2]
>> = y). Of course A[1,2] = y would create an exception saying: "use the
>> A = Matrix(Array(A)[1,2] = y)" syntax.
>
> Advantage -- the vectorization will use Arrays/lists, while Matrix
> will be reserved for strictly mathematical thing.
I really like the idea of sympy.Array.
We should also resolve the somewhat ad-hoc naming/call signature of
functions that create special matrices in the process, consider:
numpy: zeros(), ones(), accepting arbitrary shape tuple
sympy: zeronm( n, m ), one( n )
This is not good, and easy to fix even now, at least for the 2D matrices
(or should I say arrays?). I do not propose removing zeronm and one, of
course, as it would probably break a lot of code, but implementing also
zeros and ones. The former may print some future deprecation warning, as
usual.
r.
Pearu Peterson
On Jun 9, 9:18 am, "Ondrej Certik" <ond...@certik.cz> wrote:
> On Mon, Jun 9, 2008 at 10:40 AM, Pearu Peterson
>
>
> > Hi,
>
> > Just a quick note that could be relevant to the discussion:
> > check also out
>
> >http://code.google.com/p/sympycore/wiki/MatrixSupportIdeas
>
> > that contains ideas how to deal with mutable matrices in
> > operations where they should be immutable, all in a very efficient
> > way.
> > The idea is based on using views of matrices.
>
> Thanks for the tip. So what you do is that Matrix is writeable until
> it is used in an expression? And you use this flag:
>
> "
> Matrices are mutable until .is_writable==True.
> "
computed, that is, when matrix instance is used as a dictionary
key. After that, the matrix will become automatically immutable.
The is_writable flag can be used by inplace methods that
then know whether matrices need to be copied or it will
be safe to modify matrices inplace.
> You seem you implemented some interesting ideas in sympycore -- would
> you be interested in more collaboration and/or even merging the two
> codebases? At least you write on sympycore webpage that "it is created
> to fix SymPy performance and robustness issues", so let's point out
> the robustness issues (I think we fixed all of them already, but if
> there are some more left, let's fix that too)
that
robust assumption model requires an appropiate data model and I am
not sure that sympy has it.
> and as to speed, I'd be
> interested in your opinions how to move forward. I'd be interested in
> coming together for example before the EuroScipy sprint and work on it
> together. But I think it'd be efficient to prepare some plan how to
> merge it, or what to do.
start:
try to use sympycore as an internal implementation of symbolic model
in sympy and if it will be succesful, gradually remove unneeded
layers.
One can start with taking, say, the Add class and use sympycore data
structures to represent its internal content and do manipulations with
sympycore methods.
To the rest of sympy, the Add would look like nothing had changed.
Such a solution would not have the same speed as sympycore but
if it will improve sympy speed in general, then it is already a
positive sign
and improving the speed further becomes just a matter of moving
unneeded interface layers.
> Anyone interested in this, please look how things are done in
> sympycore and in sympy and let's start a discussion which way is good
> to do things and implement it. And if there are two good ways (which
> is very well possible), then let's find a solution to have both ways,
> but with the same interface.
results:
http://code.google.com/p/sympycore/wiki/PerformanceHistory
And interesting result is that sympycore is 3.5x slower than
singular in this test. Not too bad for a Python solution compared
to fastest symbolic library around. It would be interesting to
try also the new sage.calculus code that Gary is working on but
I was not able to find his code in sage sources..
Pearu
Ondrej Certik
I got your email just now, don't know what's wrong. Anyway, I created
a new issue for that:
http://code.google.com/p/sympy/issues/detail?id=890
Ondrej