Planet SymPy

July 23, 2014

Sushant Hiray

GSoC Week 9: Complexions in Complex Module

Week 9 for GSoC ended a couple of days back. I’m still working on the Complex Module.

Progress

I started of the week continuing the work on pull 248.

Initially I started working on the add and mul functions and integrating the complex module with them. This part was pretty easy as Complex addition and Complex multiplication was completely implemented. Intially I just added special cases for the case of Complex in pow. But this quickly ran into issues as is_canonical failed in some testcases.

Next I simplified the complex power for the following cases:

  • When exp is 1 or -1
  • When the base is purely imaginary and exp is integer.

Following this I made appropriate changes to is_canonical for pow

Meanwhile @certik sent a couple of PR’s sushant-hiray/csympy/#1 and sushant-hiray/csympy/#2 for implementing the Python wrappers for Complex Module. I spent some time trying to sync the wrappers and make sure they work as intended.

One important fix needed was to improve the printing. Over the entire PR, the printing has been changed multiple times to match the SymPy output. With this PR I’ve updated the current existing printing for add, mul as well as pow

The expand function is also updated to handle the case of Complex Numbers. Currently it doesn’t expand terms like: (2 + 3I)^n. Apart from this all the issues which were found have been fixed.

The Week Ahead

Merge the current PR and send a new PR for handling the cases like (2 + 3I)^n in expand


Thats all for now :) Will get back next week!

July 23, 2014 05:30 PM

July 22, 2014

Avichal Dayal

This week I converted all the operations on FPS to work with generators.
For e.g.:-
Addition for two series is something like
 def add(s1, s2):  
yield s1[0] + s2[0]
s1, s2 = s1.tail, s2.tail
for t in add(s1, s2):
yield t

Almost all the operations use recursive generators. In this case argument to the add function is different each time.
However in some cases the same function with same argument is called recursively.  For e.g:- In case of division
 def div(s1, s2):  
a, b = s1[0], s2[0]
def gen():
yield a/b
iterator = s2 * gen()
for t in iterator:
yield -t/b 
   for t in gen():
      yield t 
gen() is called recursively with the same arguments (no arguments in this case).
I'll try to optimize recursive generators if possible over the next few days.

Other problem I'm facing is args of Stream class. It is necessary that args of all Basic instances should be of type Basic. But the Stream class takes a generator object as an argument.

Next week:-
- Optimize recursive generators if possible
- Generator object as Basic instance? (Find solution for this)
- Implement all necessary ``eval`` routines for FPS like _eval_derivative, _eval_leadterm etc.

by Avichal Dayal (noreply@blogger.com) at July 22, 2014 01:11 PM

Akshay Narasimha

Gsoc Week-9

Sorry for the late post as I had to travel a lot this week. Other than that I finished the implementation of the Hyperbola class including the tests. Here is the link to the PR. Apart from that I fixed a few bugs in the Plane and Line3D class.

This week I plan to implement the Parabola class and hopefully will get the PR on Hyperbola merged.

Until then cheers!

by Akshay Narasimha (noreply@blogger.com) at July 22, 2014 06:05 AM

July 20, 2014

Sudhanshu Mishra

GSoC'14 Progress: Finished refraction at planar surface

I'm sorry for this delayed post. I couldn't work last week due to some other issues so had nothing much to report.
This week I completed refraction function for planar surface, fixed some bugs and wrote tests for it. I also added a function to calculate angle of deviation. Details can be found on the following link
I'm waiting for Sean to start reviewing this.
My next priority is to complete a pending PR on mirrors.

by Sudhanshu Mishra (noreply@blogger.com) at July 20, 2014 02:32 PM

Kundan Kumar

Week 9: Checksysodesol

This week a lot of time is spent debugging solutions. The PR solves the issue #7736 and #7723 which was arising due to wrong calculation of func method while solving non-linear system of ODEs of type5. I sent another PR which contains Checksysodesol, a method for checking the solution of system of ODEs. The PR remains to be updated with the test cases.

The PR of non-linear system of ODEs of three equation for 1st order still remains to be merged.

Thats all for this week. On monday I will be going to my college, hope to work harder there. The most haunting thing at home is electricity, having tough time with it sometime getting meagerly for 3-4 hrs a day. Well thats home, I suppose, getting a delicious meal of mom’s hand.


by Kundan at July 20, 2014 10:13 AM

Thilina Rathnayake

[GSoC] Week 9: Matrix Inverse and Sparse Matrices

Hi All, Sorry for a late blog post. I was kind of busy during last week, preparing for a competition.

During the last week, I mainly did two things, implementing matrix inverse and starting the implementation of sparse matrices.

Implementing Matrix Inverse

I implemented matrix inverse in two different methods, using Gauss Jordan elimination and fraction free LU decomposition. I had only implemented gauss Jordan elimination to solve a system Ax = b where b is a column matrix. I had to enhance the algorithm so that it can be used to solve several systems at once, i.e. b can be a collection of column matrix.

This was not a problem in fraction free LU method because L and U factors can be used to solve column vectors e_{1}, e_{2}, . . . e_{n} after calculating L and U once.

Implementing Sparse Matrices

We decided to adapt the implementation of SciPy sparse matrices. For the time being I have implemented CSR form of a sparse matrix. CSR is an acronym for “Compress Sparse Row”. You can learn more about it in the following links.

Wiki Article

Netlib Article

You can find information about scipy implementation of CSR matrices here.


by Thilina Rathnayake at July 20, 2014 08:32 AM

July 19, 2014

Harsh Gupta

week 9

This week I moved back to college and my classes have restarted. This week I worked on a PR to allow infinitely indexed Range. See https://github.com/sympy/sympy/pull/7741. While doing this PR I discovered that you cannot monkey patch object to assign some attribute. And I want to mention that Sergey(one of my mentors) is damn good reviewer.

by Harsh Gupta at July 19, 2014 06:50 PM

Sachin Joglekar

GSoC Week 9: Dyadics done

Not much to report this week, except that the code for dyadics is now done, with a all-tests-passing PR being pushed. A happy moment was the basic vector framework code finally getting pushed. Phew :-). Now the next PR in line is the dyadics one itself.
I also spent time mapping out the API for the Rotator classes, which I will proceed to discuss with Jason during this week’s meeting. I am still a little doubtful about how useful these classes may be, API wise. Lets see about this one.
A problem that I am facing currently is the confusion about the Del operator. I started out trying to write it as a SymPy function, but then I realised that the _outputs_ from the methods of Del should be the ones being defined as unevaluated functions- like Gradient, Divergence, etc. Will have to take this up with Jason too.

Anyways, I hope the dyadics code goes in soon, and I plan to send one or two PRs with more functionality soon. Till then, have a great week!


by srjoglekar246 at July 19, 2014 06:02 PM

July 18, 2014

Jim Crist

GSoC Week 9: Docs!

This week I spent time on all sorts of little things:

  • Finished up the refactoring of KanesMethod
  • Little fixes to my current PR. Just waiting on my mentors to review this, I want to get it merged soon-ish.
  • Documentation.

Writing documentation is the worst1. After taking time to implement all sorts of new interesting things, the last thing I want to do is go back and write about them in detail. Which is why it's so important to do early on. Good documentation needs to accomplish three things:

  1. Provide motivation for why your software is necessary/better/useful.
  2. Describe the user interface, showing how to use each function or class.
  3. Provide real world examples showing how to tie everything together.

Python's documentation is interesting in that there are varying ways to do it. Some of Sympy's documentation is just a nicely rendered form of the docstrings for all the methods. Other modules have a more prose-y explanation of their functionality. mechanics is one of those modules.

In my opinion the prose documentation approach is the better way. Having good docstrings is important, but they aren't the end-all of documentation2. Of course, if I have a question the first thing I'm going to do is read the docstrings (IPython makes this trivially easy). Only if I still have questions afterwards will I turn to the online documentation. However, it'd be extremely off-putting if the online documentation was just the docstrings again.

With the various changes I've made so far I needed to:

  1. Update the LagrangesMethod documentation to reflect the interface change.
  2. Create a documentation page all about the linearization methods.
  3. Update all the examples to reflect the new functionality.

All of these are "done". I still need to go through and proofread, but overall I'd say that the current state of the documentation is acceptable. I would like to take some time to reorganize the layout of the whole mechanics documentation at some point. The current layout isn't the easiest to navigate for what you're looking for.

With this out of the way, the linearization portion of my project is tentatively done. I say tentatively because I'm still waiting on my PRs to get merged, and am also still playing around with solving the nan issue that I've been writing about these last couple weeks.

With that done, I hope to move on to code generation. I've read the current code generation code and documentation, as well as this pydy wiki page on Jason's ideas about code generation. I'm still a little iffy about the intention of this functionality, so I'm waiting until we can all meet to discuss what needs to be done. That was supposed to have happened this week, but fell through. Hopefully we can set some time aside next week, and I can finally get to work on it.


  1. Not actually the worst. 

  2. This article by Steve Losh is a really good read on this. 

by Jim Crist at July 18, 2014 08:00 PM

Soumya Dipta Biswas

GSoC 2014: Week 8

Hello Everyone,

I have majorly been working on fixes and optimizations this week. There was an interesting error that was coming up because of the way quantifiers had been designed. I have made some significant changes to the Quantifier class. An interesting problem regarding constants has cropped up, which I am yet to fix. Please find the discussion regarding it at FOL Constants and Unification. For the rest or the week, I have majorly been writing tests for the FOL class and making small changes. Additionally I have spent some time writing tutorials for the logic module of SymPy. The documentation available currently only talks about how to do something but not what all can be done with it. The tutorials are far from complete but I will definitely try to finish the entire thing before the summer is over.

That’s all for now.

Farewell

by SD at July 18, 2014 01:12 PM

July 15, 2014

Akshay Narasimha

Gsoc Week - 8

This week was great in terms of the work done. My PR on the Plane class got merged finally, so this completes the  implementation of 3-D geometry. Though my code had less test coverage to which I have sent a PR which improves it .

Apart from that I have been working on the Hyperbola class and have added most of the methods, tests remain though. Here is the link to the PR.

This week I will add  tests for the Hyperbola class and will start implementing the Parabola class once this is done.

Until then cheers!

by Akshay Narasimha (noreply@blogger.com) at July 15, 2014 03:19 AM

July 13, 2014

Avichal Dayal

As it turned out, I had to replace most of my last two week's code. I used lambda functions to implement infinite structures. Lambda function stored the delayed part of the structure.

I scraped that out and now instead use the more Python-ic generators for infinite structures.

Earlier FPS was just for series expansion around x0=0. I added support for asymptotic expansion and series expansion around arbitrary points. This just included the substitutions as done in the series() method.

However I faced an irritating bug regarding this.
The Lazyseries class has arguments: <x>, <generator object>
So if ``l = Lazyseries(x, generator)``
When  ``l.subs(x, x+1)`` is done, it turned out that <generator object> was being iterated. I noticed this after sometime, fixing it by writing a custom _eval_subs()

Regarding the operation on series, I had to rewrite all of them using generators.
Adding series is just term-wise addition which can be accomplished by a simple one-liner i.e. ``imap(operator.add, self, other)``.
But it turned out that this was actually wrong. This didn't give sorted output when adding series.
E.g.:-
s1 = x + x**2 + x**3 + ...
s2 = x + x**3 + x**5 + ...
Now if s = s1 + s2,
s = 2x + (x**2 + x**3) + (x**3 + x**5) + ...
If I want series expansion of s with order = 4, it terminates at the third term after getting the x**5 term. Whereas, it should go further as there is a x**4 term too.


I wrote a custom add function merging two iterators similar to the mergesort algorithm. There is a built-in function called heapq.merge that does the same thing but it does not provide a key to compare values. It turns out that it is one of Cpython's ideas floating around to provide a key for heapq.

I hope to get this done in a few days and get it ready for the merge.


by Avichal Dayal (noreply@blogger.com) at July 13, 2014 06:41 PM

Aaron Meurer

SciPy 2014

I just finished SciPy 2014, a week-long conference in Austin, TX for scientific computing with Python.

This is my third SciPy (I have been to 2011 and 2013). This year, the conference was noticeably larger. Last year there were ~350 people, this year, there were ~450 people. Aside from there being a lot more people, and the main keynotes taking place in a larger room, the most noticeable consequence of this is that there were three days of talks this year, and three concurrent tracks of talks all three days (last year there were two of each). The conference consisted of two days of tutorials, three days of talks, and two days of sprints, running from July 5 to July 12.

Tutorials

The conference started on Sunday with tutorials. I gave a tutorial on SymPy with Matthew Rocklin and Jason Moore. The videos are on YouTube (parts one, two, three, and four). I gave tutorials for SymPy the previous two times I was at SciPy, although with different people (with Mateusz Paprocki in 2011 and Ondřej Čertík in 2013). I really enjoy seeing new people learn about SymPy, and working with Matthew Rocklin, who is a very good speaker and teacher.

I also attended the tutorial on PyDy by Jason Moore, Gilbert Gede, and Obinna Nwanna (parts one and two). This tutorial was also well done, and I highly recommend it if you are interested in Newtonian mechanics.

I unfortunately was unable to attend any of the other tutorials, but I heard good things about them, especially the Julia tutorial.

Talks

From Tuesday to Thursday were talks. The quality of talks this year was very high. The SciPy talks have always been high quality talks, but this year I felt that they were particularly good. I don't think I saw a bad talk.

Thus, I can't really recommend the good talks that I saw without recommending all of them. You should go to YouTube and the SciPy schedule and watch any talk that looks interesting.

I therefore am going to recommend here the very best talks. Two talks in particular stood out to me as the best.

First is Greg Wilson's Thursday keynote, which is among the best talks I've ever seen from any conference.

<iframe allowfullscreen="allowfullscreen" frameborder="0" height="315" src="http://www.youtube.com/embed/1e26rp6qPbA" width="560"></iframe>

Greg mentions a lot of ideas, quite a few of which are controversial, which I think always makes for an interesting talk (it also means that I don't agree with everything he said, although I do agree with most of it). Most of the talk is about pedagogy, especially regarding his experiences at Software Carpentry. Some things he posited:

  • There is actually good research about what methods work and don't work in teaching. He referenced this presentation, which lists just about every possible pedagogical method, and the net difference that it has on students, based on over 50,000 studies. For example, individualized instruction has a very small positive effect, whereas teacher feedback has a very large positive effect. Since each takes time and resources, we should focus on those effects that have the highest impact. Greg pointed out that web-based learning has very little positive effect, and hence is a waste of time and money. The most effective change is removing disruptive students.

    In particular, I liked the quote, "if you want more computing in high school, you have to tell me what to take out." People like to go on that schools need to teach more of this or more of that, and computing and programming tends to be high on that list these days, but anyone who does not discuss what things should be removed from the curriculum, which is already quite full, is not being honest about the discussion.

  • The other big point Greg made is that we need more incremental massive collaboration in teaching. This is the same model that has built open source and Wikipedia, but is mostly absent from teaching. Incremental change is important here, as well. It is more useful for someone to contribute fixes to existing lesson plans, so that they become better for the students, but in his experience, people are much more willing to add new lessons. Greg calls for a "culture of patching". If incremental change could be adopted in teaching, teachers could aggregate methods and lesson plans, removing the massive duplication, and most importantly, making teaching materials that actually work for students to learn. Greg Wilson asks why open source and Wikipedia are able to thrive on massive incremental change, but teaching is not, a question he hasn't found the answer to.

    My thought on the matter is that unlike writing software or collecting and presenting facts, pedagogy is very difficult. If I contribute a patch to an open source project that fixes a bug, I can run the tests to see if my fix is "correct". If I fix an incorrect fact on Wikipedia, it is less easy, but I can still cite and check references to make sure it is correct. But for teaching, it is very difficult to know what methods work and what don't, and as Greg pointed out at the beginning of his talk, the effects of different methods can be very counterintuitive.

The second talk that I recommend is Jake VanderPlas's talk about Frequentism and Bayesianism.

<iframe allowfullscreen="allowfullscreen" frameborder="0" height="315" src="http://www.youtube.com/embed/KhAUfqhLakw" width="560"></iframe>

I won't summarize this talk, as Jake has done a much better job in his blog (parts one, two, three, and four). The best thing is to just watch the talk. I will just point out that before the talk, I did not really understand the difference, not being a statistician or someone who works with statistics regularly, and having seen the talk, I now feel that I do. It's a controversial topic, and if you care about the matter, you should know that Jake is a strong Bayesian, although I felt that he gave both sides a fair exposition.

Again, all talks I saw at the conference were good. But those two I felt were the best. I should also mention here that I myself gave a talk on Conda (more on that later).

The Conference

Of course, the talks are only a part of any conference. The best part of SciPy is the gathering of the community. Each year I meet more new people, as well as talk with people I already know, but don't get to see outside of SciPy.

For me, the biggest part of the interactions this year were on Conda and packaging. The background is that I have been working full time for Continuum since January, and I had interned last summer, working primarily on the Conda package manager and Anaconda, the Python distribution. This year, some of the biggest buzz at the conference was about Conda. I'm obviously selection biased, because people came to me specifically to talk about Conda, but I also overheard it in other people's conversations, in several of the presentations, and frankly, the people who did talk to me about Conda were very excited about it. Just like everyone was talking about the IPython Notebook last year and how it has solved the fundamental problems of sharing and presenting data analysis, this year, everyone thanked me for my work on Conda and how it has basically solved the packaging problem, the ubiquitous problem in Python since people started using it.

Conda: The Packaging Problem Solved

Here is the talk I gave on Conda:

<iframe allowfullscreen="allowfullscreen" frameborder="0" height="315" src="http://www.youtube.com/embed/UaIvrDWrIWM" width="560"></iframe>

I made the claim in my talk that Conda has solved the packaging problem, and the general feel from people I talked to who are using Conda is that it has.

I think this slide from my presentation summarizes why Conda solves the packaging problem.

One of the most amazing things about the scientific Python community, and one of the things that I think really sets it apart from other Python communities, is the use of Python alongside other languages, such as C, C++, Fortran, R, or Julia. No one language is enough to get the job done for serious scientific work. The fundamental brokenness of Python packaging has been that it has focused too much on Python specific tools and processes. The distutils/setuptools/pip/virtualenv stack works great if your code begins and ends with Python. Where it falls over is when you want to link against a C library, compile some Fortran or Cython code, and communicate with other languages like R and Julia. By being a system level package manager, which is fundamentally Python agnostic, Conda is able to deal with all packages equally, whether that package be a Python package, a C extension which other packages link against, or Python itself.

By being truly cross-platform and user installable, Conda is able to reach the maximal number of users, especially those who have historically been hit by the packaging problem the hardest: those who are on Windows or those who do not have admin rights to install necessary tools to install the packages they need.

Finally, Conda installs binaries, not source packages, and its metadata is entirely static (you do not need to execute arbitrary Python code to capture the metadata of a package). These two things remove two of the largest sources of issues with the existing Python packaging tools, such as compiler errors, and nonuniformity in metadata standards (there seem to be as many different ways of writing setup.py as there are packages on PyPI), by removing arbitrary code execution from package installation.

Conda opens up its ecosystem to anybody by making it easy for people to build their own Conda packages using reproducible Conda recipes. And Binstar makes it easy to share those packages. I'm very excited about Binstar, as I think it does for packaging what GitHub has done for open source, i.e., distributes and democratizes it. There are challenges on how to deal with this, of course. As with any distributed democratized system, Binstar can be a wild west of packages. Continuum is thinking about ways to manage this complexity, while still reaping the benefits it provides. If you have any thoughts on things that can be done, let me know in the comments section below.

Of course, solving the packaging problem and removing it are different things. Conda does not make it easier to compile difficult packages. It only makes it so that fewer people have to do it. And there is still work to be done before Conda really takes over the world.

Sprints

The conference ended with two days of sprints. I mainly helped people with Conda packaging. One key thing that happened is that I worked with Aron Ahmadia so that HashDist can generate Conda packages. HashDist is a package compiling framework that makes it easy to have completely reproducible builds by hashing all the information that was used to compile a package, and recompiling when any of that information changes. You can learn more about HashDist by watching Aron's talk from the conference:

<iframe allowfullscreen="allowfullscreen" frameborder="0" height="315" src="http://www.youtube.com/embed/wviHkzk0AkY" width="560"></iframe>

I am now convinced that HashDist is a good solution for people who still want the control of compiling their own packages. Once HashDist is able to produce Conda packages, then you can gain the benefits of both worlds: Conda's powerful package management and environment notion, with HashDist's modular and reproducible package building framework.

Other thoughts

The organizers of SciPy did an excellent job this year. The video crew did something which I have not seen before, which is that they uploaded the videos of the talks on the same day that the talks were held. My talk, which was held right before lunch, was uploaded before the last talk of the day. Something that I saw come out of this is that people not attending the conference were able to watch the talks and take part of the conversation with the conference attendees, via Twitter and other social media, or by joining the sprints after the conference.

The extended three days of talks really took their toll on me. The conference starts early enough in the morning and the social events after go so late in the evening that each day of the conference I become a little more sleep deprived. Usually by two days of tutorials and two days of talks I have hit my limit, and this year, I really had a hard time making it through that fifth day. Fortunately for the sprints I was able sleep in a little bit, as it's not a big deal if you miss the beginning.

This year the conference organizers made a push for diversity, and it shows. There were noticeably more women at the conference this year, and not just insomuch as there were more people at all.

Finally, I leave you with the greatest lightening talk. Ever.

<iframe allowfullscreen="allowfullscreen" frameborder="0" height="315" src="http://www.youtube.com/embed/ln4nE_EVDCg?start=3254" width="560"></iframe>

by Aaron Meurer at July 13, 2014 04:52 PM

Kundan Kumar

Week 7 and Week 8 : Non-linear system of ODEs

In week 7, there was negligible work done, had health issues. Because of that I was unable to update things here. During week 8, I got in non-linear system of ODEs for 3 equations of first order (PR). But most of them is not responding in my system. It takes very long time to respond but still I am not getting the final solutions (before that my system crashes).

Though I have written method for these types, test cases for two types have been included and for other 3 I am not getting any solution. I will needing help from Tim and Sean for making any advances in this.

The one PR of non-linear system of ODEs for 2 equation of first order has been merged and the PR for docs is all ready to go.

The next week I plan to finish the PR of non-linear one and start new for checksysodesol which I have just started and start working on non-linear system of ODEs for 2 equation of second order.


by Kundan at July 13, 2014 09:18 AM

Sushant Hiray

This Week in CSymPy: #8

Week 8 for GSoC just ended. I’m particularly pleased that I’ve managed to get a basic version of Complex Module merged into master.

Progress

I hadn’t made much progress last week as I was traveling back home. So I had only managed to implement some of the functions. This week I implemented all the virtual functions which were necessary. Updated Complex::from_mpq to return the appropriate instance of Complex and Rational. Pull 223 implemented the basic version of complex module.

I also noticed that some of the files I had commented out earlier had a mismatched tab configuration. So I fixed the mismatch identation for all the files in the project. Pull 242 fixed this issue.

While writing tests for the Complex Module, I noticed that there was no explicit divide by zero check in the code and as a result, on division by zero, we were getting a coredumped. This issue was filed as #239 which was subsequently fixed via Pull 241.

After some fixes and adding divide by zero checks in the complex module, the Pull 223 was merged. With this in place, we can perform all basic operations on complex integers as well as rationals. The only exception being pow. It was decided that complex pow shouldn’t evaluated by default. Also since the current definition of pow returns a RCP<const Number> so it needs to be refactored to account the current design.

The Week Ahead

Make sure complex module works with all other submodules of CSymPy and add tests accordingly. Pull 248 has been started to deal with this.


Thats all for now :) Will get back next week!

July 13, 2014 05:45 AM

July 12, 2014

Harsh Gupta

week 7

Week 7, 8

Fu Simplification

In the early part of the week 7 was thinking and working on design decisions. I want the code to be very modular so that it could be easily extended by people other than me. This came up in a meeting with Matthew and I want the solvers to be like the Fu simplification explained by Matthew in this Scipy Talk. The idea was that we can see solving an equation as series of transformations. If we have a lot of small transformations such that the input type is same as output type, and some notion of what makes a "better" output we can search though the list of transformations running one on top of other. I also posted about it on the mailing list which brought out some flaws in the preliminary design. The idea is pretty crude in the current stage and I'll have to look deeper into it, but not now.

I also discussed about the implementing a pattern dispached based solver suggested by F.B on the mailing list. But we decided that it will be better if we finish the equation solver by the current technique first.

Intersection with S.Reals

I decribed in the last post that one way to solve trigonometric equation is rewriting them in terms of \( exp \). But that is \( exp \) in the complex domain and the solution of \(exp(x) = a \) is \( \left\{i \left(2 \pi n + \arg{\left (a \right )}\right) + \log{\left (\left\lvert{a}\right\rvert \right )}\; |\; n \in \mathbb{Z}\right\} \). Hence we have to filter out real solutions from the obtained solutions. The filering is equivalent to the intersection of the solutions with the \( \mathbb{R} \) set. Suppose \( g(x) \) and \( h(x) \) are real valued functions and we have to perform $$ \mathbb{R} \cap \left\{g{\left (n \right )} + i h{\left (n \right )}\; |\; n \in \mathbb{Z}\right\} $$ then the answer will be simply $$ \left\{g{\left (n \right )}\; |\; n \in \left\{h{\left (n \right )} = 0\; |\; n \in \mathbb{Z}\right\}\right\} $$

Separate the real and imaginary parts and equate the imaginary to zero but the problem was with the assumptions on the symbols. For example while separating real and imaginary parts of the equation.

In[]: (n + I*n).as_imag_real()
Out[]:(re(n) - im(n), re(n) + im(n))

That is because n is by default complex, even in the Lambda(..).expr. I wrote some code to decide the assumption on the variable of imageset from the baseset. See PR 7694. There was another issue that needs to be resolved S.Integers.intersect(S.Reals) doesn't evaluate to S.Reals.

LambertW and Multivariate Solvers

The method to solve equation containing exp and log function is using the LambertW function. LambertW function is the inverse of \( x \exp(x) \). The function is multivariate function both for the real and complex domains and Sympy has only one branch implemented. This also leads us to loss of solutions. Aaron gave an example here. But I'm pretty unfamiliar with solving by LambertW and LambertW itself and it will take me some time to build an understanding of them. As an ad hoc solution I'm using the code in the _tsolve in the solvers.solvers module to do at least what the current solvers can do.

When the importing of _tsolve method was done. I started working on the multivariate solvers. Here's how the current multivariate solvers work:

Solving single multivariate equation

  1. count the number of free symbols in f - no of symbols for equation. If the equation has exactly one symbol which is not asked for then use solve_undetermined_coeffs, the solve_undetermined_coeffs find the values of the coefficient in a univariate polynomial such that it always equates to zero.

  2. Then for each symbol solve_linear is tried which tries to find a solution of that symbol in terms of constants or other symbols, the docstring says

    No simplification is done to f other than and mul=True expansion,
    so the solution will correspond strictly to a unique solution.
    

So we don't have to worry about loosing a solution. For every symbol it is checked if doesn't depend on previously solved symbols, if it does that solution is discarded.

  1. For the symbols for which the above method failed, the _solve function is called for the equation for that variable and as above if the solution contains a variable already solved then that solution is discarded.

System of equations in multiple variables

  • Try to convert the system of equations into a system of polynomial equation in variables

  • If all the equations are linear solve then using solve_linear_system, check the result and return it. If asked for particular solution solve using minsolve_linear_system

  • If the number of symbols is same as the size of the system solve the polynomial system using solve_poly_system. In case the system is over-determined All the free symbols intersection the variables asked for are calculated. Then for every subset of such symbols of length equal to that of the system, an attempt to solve the equations by solve_poly_system is made. Here if any of the solution depends on previously solved system the solution is discarded.

  • In the case there are failed equations:

    • For every know result:
    • Substitute every thing into the failed equation and see if the equation turns to zero. if it does accept the result otherwise put it in the bad_results group.
    • Then try to solve try to solve the failed equation using solve for each symbol.
    • If that solution depends on any other previously solved symbols discard it.
    • If it doesn't satisfy other equations, discard it.
    • Check if the solution doesn't set any denominator to zero, if it does discard that solution.
    • If it satisfies the above conditions substitute this value in know solutions and add it as a new result.

by Harsh Gupta at July 12, 2014 02:00 PM

Sachin Joglekar

GSoC Weeks 7,8: Beginning with Dyadics

The last two weeks have been slow, hence a common blog post for both of them. Last week was spent in improving the current PR even more, especially with dedicated rotator methods for each type of rotation (Body/Space/Axis/Quaternion) of Coordinate Systems. Its not much, just methods like ‘orient_new_axis’, ‘orient_new_body’, etc. as add-ons on top of the all-powerful ‘orient_new’ method. The coding involved wasn’t much, since the new methods essentially call the ‘orient_new’ method with the relevant parameters. However, the rotation API is now much simpler than using the ‘orient_new’ method for all kinds of orientations – mainly because the extra flexibility provided by the ‘orient_new’ method implies a more awkward API for the same.
Jason, Mathew and Aaron were together at SciPy, during our last week’s GSoC meeting. There, the common consensus was to try and have Rotator classes to make the orientation framework even more user-friendly. Discussing the interface for the same is on this week’s meeting agenda.
A big change that we thought of doing in the timeline was to let go of spherical/cylindrical systems for a while, and instead focus on dyadic tensors – since the code is mostly done (my work during last year’s GSoC), they will be a useful and powerful addition to the vector module.
So my next few (much shorter than the current one) PRs would be-
1) The code for dyadic tensors
2) Documentation for all the basic, non-time-dependent framework
3) Rotator classes(?)
4) Implementation of the Del operator as an SymPy Function
5) Just some basic helper methods to convert vectors from one form (rect/spherical/cylindrical) to another
The code for Dyadics is mostly done, the docs are on-going, and I am thinking of having a common super-class for Vector and Dyadic- called Tensor(inheriting from Expr) – since many of the SymPy manipulation procedures work in the exact same manner for Vectors as well as Dyadics. Will discuss this with Jason soon.

Anyways, thats all for now, have a great week!


by srjoglekar246 at July 12, 2014 07:37 AM

July 11, 2014

Jim Crist

GSoC Week 8: Try, try, try again...

I'm still struggling to solve the nan and oo issue I've discussed in my post a couple weeks ago. Last week I showed off a custom written subs function for use inside sympy.physics.mechanics that helped with speed considerably, and attempted to solve the nan problem. This worked great for small-medium expressions, but failed on large ones. Or did it? I'm not sure anymore.

This pull request brought up something that I had witnessed, but never really thought about as a potential source of my issues. To summarize, Sympy's current (hopefully soon to be old) caching system never clears. Ever. For interactive work, or short running sessions this is fine. However, for the huge expressions generated in mechanics, this can be a source of memory issues, as the cache grows to hold all sub-expressions that were cached.

It turns out that simplify is one of those functions that is cached. This may explain why when I tried to use msubs with smart_subs=True (which crawls the expression tree and does selective simplification) this resulted in all of my RAM being used up (4 GB!!!). I haven't had a chance to pull in this PR into my repo and test it out, but it sounds like it should fix the problem. Instead of growing infinitely, the cache uses a least recently used (LRU) algorithm to determine what stays and what is removed. The cache size can be set by the user, so those that prefer speed over memory use can still cache everything. Per his benchmarks it seems to be only 10% slower, which shouldn't be much of a problem. Overall, I'm really psyched to start using this. Perhaps with this the smart_subs I wrote up will work, even if it takes a while. If not, I'm kind-of out of ideas.

I spent some time this week trying out a few other ways of solving this problem. So far none of them have worked.

1. Using cse, and applying simplify selectively to the sub-expressions.

The basic idea here was to apply cse on the expression, and then evaluate each sub-expression. If it evaluated to nan, simplify it, then evaluate it again.

This seemed like a good idea at first, but upon closer examination it falls apart. The issue is that the expressions that could cancel/simplify out are often broken into separate sub-expressions. This means that they are evaluated numerically separately, and only once combined will they result in a nan, at which point they can't be simplified anyway.

2. Taking the limit of the bad sub-expressions.

This was another idea that seemed good until I tried it. Similar to the smart_subs I talked about last week, except this time it's taking the limit of the bad sub-expressions as they approach the operating point. The thought being that it may be computationaly cheaper to find the limit than to apply simplify and then evaluate.

There were several problems iwth this design. The first being that Sympy has no functionality for finding multivariable limits. These can't be calculated iteratively either (by that I mean find the limit for x, then the limit for y, then the limit for z, etc...), as the part that could "simplify out" could already be gone.

The second, and more serious issue, is that there was no way to tell if the limit at that point was equal to the value the expression should actually evaluate too, or if it is just the value of the limit at that point. For example:

>>> expr = (a - 1)/(a**2 - 1)
>>> op_point = {a: 1}
>>> expr.subs(op_point)
nan
>>> limit(expr, a, 1, '+')
1/2
>>> limit(expr, a, 1, '-')
1/2

Using the method described above, it would seem that the expression should just evaluate to 1/2. However, if you actually plot this expression, you'll find that there is a discontinuity at a = 1. From either side it approaches 1/2, but at 1 it is actually nan.

3. Numerical perturbation about the setpoint to find the limit of the bad sub-expressions.

The idea here was to calculate the limit of the sub-expressions through numerical evaluation and perturbation. This fails for all the reasons described above, as well as the fact that Sympy is a symbolic computation library, and we should be able to do this symbolically.


Unfortunately those were all the ideas I had to solve this problem. If the algorithm described last week doesn't end up working using the new cacheing system, I'm kind of stumped. Back on the struggle bus...


Meanwhile...

As another potential solution, I've set about refactoring the KanesMethod class in the hope that I'll find some way of generating expressions that are smaller than they currently are. The first step was rewriting to make it readable, more modular, and remove the dead code that had built up over the years. This is done. In it's current state it passes all tests, and runs them in half the time that it had before!!! Still no major reduction in expression size, but I'll hopefully find some magical place in the code that could be made more efficient. We'll see.

I'm also working on the documentation for the linearization stuff that's already done, as well as waiting on someone to finally review my PR for LagrangesMethod support. I hope to get that merged soon so that I can get started on the code generation portion of this project.

by Jim Crist at July 11, 2014 10:00 PM

Thilina Rathnayake

[GSoC] Week 8: Solving Ax = b and Determinant

Hi All,

Sorry, I couldn’t write a blog post last week. During past two weeks, I contributed in solving some of the bugs in CSymPy and the determinant of a Matrix. Also, I implemented the solutions of the system Ax = b by using various decomposition methods.

Determinant of a Matrix

I implemented determinant using two different methods, Bareiss’s method and Berkowitz algorithm.

Bareiss’s algorithm

Bareiss’s algorithm[1] can be used to find the row echelon form or determinant of a Matrix. This is not a division free algorithm, but guarantees that all the divisions are exact, i.e. there is no remainder [2]. The algorithm is based on Sylvester’s identity and a transformation that yields the determinant after successive application. This can be also used to solve a system of equations. You can read more about the algorithm in [2].

Berkowitz Algorithm

Berkowitz algorithm can also be used in calculating the determinant. This algorithms has various other applications as well, like calculating characteristic polynomial of a matrix, principal minors and Eigen values. I am yet to implement the calculation of characteristic polynomial and Eigen values using this algorithm but that won’t be a difficult thing. I wish to do it over the weekend.

Ax = b

I used various matrix decompositions implemented earlier to solve the system Ax = b. Only the fraction free LU decomposition was used earlier, but now we can solve linear systems using LU decomposition and LDL decomposition. QR and Cholesky decomposition can be enabled after figuring out a good way to do expression simplification in CSymPy.

I hope to work on Sparse matrices in upcoming weeks and do some benchmarks of CSymPy’s algorithms with GiNaC.

References

[1] Erwin H. Bareiss, Sylvester’s Identity and Multi step Integer-Preserving Gaussian Elimination

[2] Wikipedia article, http://en.wikipedia.org/wiki/Bareiss_algorithm

[3] Wolfram Mathworld article, http://mathworld.wolfram.com/SylvestersDeterminantIdentity.html

 


by Thilina Rathnayake at July 11, 2014 03:51 PM

July 10, 2014

Soumya Dipta Biswas

GSoC 2014: Week 7

Hello everyone,

I have been extremely busy this week. So I haven’t managed to get the time to write out a proper post. However since it is already very late in the week it is better to create a small outline than wait for the actual product. Firstly, I worked on standardization of FOL expression. However the bulk of the week was consumed in writing code for the most general unifier and resolution. While these functions mark almost the end of my FOL proposal, the entire module is yet not ready and still needs some work. I will explain all these concepts in detail (preferably some time next week, once I have more time) along with examples and how to do it in SymPy.

Cheers!!!

by SD at July 10, 2014 07:13 PM

July 08, 2014

Sudhanshu Mishra

GSoC'14 Progress: Working with geometry

This week I completed refraction_angle. Now it also works with Ray3D and Plane. This function calculates transmitted vector after refraction. medium1 and medium2 can be Mediumor any sympifiable object. Ifincident is an object of Ray3Dnormal also has to be an instance of Ray3D in order to get the output as a Ray3D. If plane of separation is not provided and normal is an instance of Ray3D, normal will be assumed to be intersecting incident ray at the plane of separation. This will not be the case when normal is a Matrix or any other sequence. 
If incident is an instance of Ray3D and plane has not been provided and normal is notRay3D, output will be a Matrix. It is dependent on Plane so I haven’t added tests for it.
Here’s the link to the PR https://github.com/sympy/sympy/pull/7626
This week I also worked on making spherical mirrors in SymPy. There are few issues that I am facing. Currently I’m working(stuck) on locating mirrors in space. It’s an amalgamation of 2D and 3D geometry. I still have to subclass Ray of geometry and make it useful for this. I’ve sent a WIP PR.
I’m waiting for Sean to merge this long waited PR for moving Gaussian optics module.
That's all for now.

by Sudhanshu Mishra (noreply@blogger.com) at July 08, 2014 05:24 PM

Akshay Narasimha

Gsoc Week-7

This has been a decent week in terms of progress. One of my long standing PR's has been merged (Line3D class). The other one on Plane has been reviewed by a lot of people and luckily I got a lot of suggestions as to how I can improve the code. I have made the required changes and am waiting for a final go on this PR as that would speed up the work on the optics module as it has a dependency on this class.

Meanwhile this week I also worked on the Hyperbola class. I have added a few more methods. Here is the link to the PR. I plan to implement more methods this week along with the Rectangular Hyperbola class.

Until then cheers!

by Akshay Narasimha (noreply@blogger.com) at July 08, 2014 04:42 AM

July 06, 2014

Jim Crist

GSoC Week 7: Expression Trees and Substitution

Late post this week due to the celebration of everything that is American. This week I finally got my first PR merged. The general linearization code is now part of the Sympy codebase. I currently have a PR for Lagrange support waiting to be reviewed. After that I just need to write the ReST docs and the first part of the project is "complete". The rest of the week was spent on more optimization work. I'm getting closer to being able to solve the bicycle example in a reasonable amount of time!

Last week's post showed the issues with the large expression size and subs (it takes forever to run). I took some time this week to look into how expressions work in sympy, and wrote a specialized subs function for use in sympy.physics.mechanics. The rest of this post will show give an overview and some benchmarks of this code.

Expression Trees

In sympy, expressions are stored as trees. Each node is an object that has an attribute args that contains a list of it's child nodes. This is best shown by an example:

In [1]:
from sympy import *
from sympy.printing.dot import dotprint
In [2]:
a, b, c = symbols('a, b, c')
In [3]:
test = a*cos(a + b) + a**2/b
test
Out[3]:
a**2/b + a*cos(a + b)
In [4]:
with open('test.dot', 'w') as fil:
    fil.write(dotprint(test))
In [5]:
%%bash
dot -Tpng test.dot -o test.png
In [6]:
from IPython.core.display import Image 
Image(filename='test.png')
Out[6]:

The root node of this tree is an Add object, as the outermost operation is adding two smaller expressions together. Going down the left side first, we see that \(a^2/b\) is stored as multiplying \(a^2\) times \(b^{-1}\). Sympy doesn't have Div objects, fractions are all expressed using Mul, with the denominator wrapped in a Pow(den, -1). Traversing the right side \(a \cos(a + b)\) is stored as multiplying a and a cos object together. The cos object itself contains one child node - an Add - which holds a and b.

Crawling the Tree

The design of sympy expression objects uses two key features: args, and func. As mentioned above, object.args holds a tuple of all the child nodes for that object. Likewise, object.func is a class method that takes in arguments, and returns an instance of that class. For most objects in sympy, running

object.func(*object.args) == object

will be a true expression. I say most because not all objects adhere to this (leaf nodes). There has been a lot of discussion about this here, if you're interested.

One of the great things about this design is that it makes it very easy to write operations that modify these trees using recursion. Normally, recursion in python is frowned upon because the function calls add overhead that could be removed if rewritten as a loop. There is also a maximum recursion depth (default of 1000) to prevent stackoverflow conditions. However, a sympy expression that has 1000 nested nodes is highly unlikely (even in mechanics), and the recursion makes the code much more readable. As the zen of python says, "Readability counts".

A simple crawler to print every node is written below. It consists of two functions.

The first one is a generic function crawl, that crawls the expression tree, calls func on each node, and returns the result if there is one. Otherwise it recurses down a level into the child nodes, forming a list of new_args, and then calls expr.func to rebuild the expression from those args.

The second one is a printer function. As we don't want to modify the expression at all, we'll just print the node, and then return the expression if it doesn't have args (it's a leaf node). Note that there are more efficient ways to traverse the tree and print all the nodes - this is mostly to demonstrate crawl, as it will be used later.

Using these two functions a function that crawls the tree, prints every node, and returns the original expression can be composed using a simple lambda statement.

In [7]:
def crawl(expr, func, *args, **kwargs):
    """Crawl the expression tree, and apply func to every node."""
    val = func(expr, *args, **kwargs)
    if val is not None:
        return val
    new_args = (crawl(arg, func, *args, **kwargs) for arg in expr.args)
    return expr.func(*new_args)

def printer(expr):
    """Print out every node"""
    print(expr)
    if not expr.args:
        return expr

print_expr = lambda expr: crawl(expr, printer)

Let's test this function on our expression from above:

In [8]:
temp = print_expr(test)
# Checking that the expression was unchanged (temp == test)
assert temp == test
a**2/b + a*cos(a + b)
a**2/b
a**2
a
2
1/b
b
-1
a*cos(a + b)
a
cos(a + b)
a + b
b
a

Comparing the printed results with the tree diagram from before, one can see how each big expression can be decomposed into smaller expressions. Further, the rebuilt expression after traversing was identical to the input expression. In the next section we'll write another function that changes the expression tree using crawl.

A Custom subs Function

In sympy.physics.mechanics, we deal with symbols, dynamic symbols (which are of type AppliedUndef), and derivatives of these dynamicsymbols. Unfortunately, the provided subs function traverses inside the Derivative terms, giving underdesired results:

In [9]:
from sympy.physics.mechanics import dynamicsymbols
x, y = dynamicsymbols('x, y')
test = a*(x + x.diff())
test
Out[9]:
a*(x(t) + Derivative(x(t), t))
In [10]:
# Subbing in b for x. Desired result is a*(b + x.diff())
test.subs({x: b})
Out[10]:
a*(b + Derivative(b, t))

To get around this problem, we've been using a custom function _subs_keep_derivs. This function creates two substitution dictionaries - one with Derivative, and one without. Four substitutions then take place:

  1. Perform subs with the terms in the derivative dictionary
  2. Substitute in Dummy symbols for all Derivative terms in the resulting expression
  3. Perform subs with the terms in the non-derivative dictionary
  4. Substitute back the original Derivative terms from the Dummy symbols

This is slow due to the need for four calls to expr.subs. Also, subs applies substitutions sequentially (i.e. each term in the substitution dict requires its own tree traversal). For our purposes in mechanics, this is unecessary. Using the already written crawl function, we can compose our own subs that ignores terms inside derivative objects fairly easily:

In [11]:
def sub_func(expr, sub_dict):
    """Perform expression subsitution, ignoring derivatives."""
    if expr in sub_dict:
        return sub_dict[expr]
    elif not expr.args or expr.is_Derivative:
        return expr
    
new_subs = lambda expr, sub_dict: crawl(expr, sub_func, sub_dict)

That's it. Due to the composable nature of crawl, the code needed to perform this operation is incredibly simple. Let's test it to make sure it on our previous expression:

In [12]:
# Simple example
new_subs(test, {x: b})
Out[12]:
a*(b + Derivative(x(t), t))

So it leaves terms inside derivaties alone, exactly as desired. We can see how this compares to the previous implementation:

In [13]:
# Old way of doing things, taken from sympy.physics.mechanics.functions
def _subs_keep_derivs(expr, sub_dict):
    """Performs subs exactly as subs normally would be,
    but doesn't sub in expressions inside Derivatives."""

    ds = expr.atoms(Derivative)
    gs = [Dummy() for d in ds]
    items = sub_dict.items()
    deriv_dict = dict((i, j) for (i, j) in items if i.is_Derivative)
    sub_dict = dict((i, j) for (i, j) in items if not i.is_Derivative)
    dict_to = dict(zip(ds, gs))
    dict_from = dict(zip(gs, ds))
    return expr.subs(deriv_dict).subs(dict_to).subs(sub_dict).subs(dict_from)
In [14]:
# Benchmark substitution
print("Using _subs_keep_derivs")
%timeit _subs_keep_derivs(test, {x: b})
print("Using new_subs")
%timeit new_subs(test, {x: b})
Using _subs_keep_derivs
100 loops, best of 3: 3 ms per loop
Using new_subs
10000 loops, best of 3: 39.9 µs per loop

So it's significantly faster. For this small benchmark, approximately 75x faster. Also, as it works in only one traversal of the tree it has a smaller computational complexity - meaning that for larger expressions this speed increase will be even higher. For kicks, let's see how it compares to normal subs for an expression without derivatives:

In [15]:
# For kicks, see how it compares to subs for expr without derivative:
test = a*cos(x + b) + (a - x)**2
print("Using subs")
%timeit test.subs(test, {a: 1, b: 2})
print("Using new subs")
%timeit new_subs(test, {a: 1, b: 2})
Using subs
1000 loops, best of 3: 676 µs per loop
Using new subs
1000 loops, best of 3: 247 µs per loop

So our implementation of subs is faster than sympy's. However, this micro benchmark isn't all that meaningful. Normal subs works in multiple passes so that expressions inside the sub_dict are also affected. subs also incorporates math knowledge about the expressions, while ours just does a naive direct match and replace. For our purposes though, this is sufficient. Also, as with all micro-benchmarks, this should be taken with a grain of salt. The main point is that our custom subs function is significantly faster than the old method of _subs_keep_derivs.

Custom Simplification

The other major issue I discussed last week was that some expressions result in nan or oo (infinity) when not simplified, but after simplification result in a realizable expression. For example:

In [16]:
test = sin(a)/tan(a)
print("Before simplification:", test.subs(a, 0))
print("After simplification:", simplify(test).subs(a, 0))
Before simplification: nan
After simplification: 1

For small expressions, doing simplify before evaluation is acceptable, but for larger ones simplification is way too slow. However, these divide by zero errors are caused by subexpressions, not the whole expression. Using some knowledge of the types of expressions present in the mechanics module, we can come up with some simple heuristics for what can result in nan, oo and zoo:

  1. tan(pi/2)
  2. Fractions with 0 in the denominator

In reality, these are the same thing, because

\[ \tan(\pi/2) = \frac{\sin(\pi/2)}{\cos(\pi/2)} = \frac{1}{0} \]

Using this knowledge, we can come up with a simple algorithm for performing subs and catching these issues at the same time:

  1. Replace all tan(*) with sin(*)/cos(*). This will allow us to only have to check for denominator = 0 conditions.
  2. For nodes that are fractions, check if the denominator evaluates to 0. If so, apply simplify to the fraction, and then carry on as normal.

This may not catch every instance that results in would result in a nan without simplification, but it should catch most of them. Also, the algorithm is so simple, that it can be implemented in only a few lines. First, the tan replacement. This requires almost no new code, as it can be composed using the already written crawl function:

In [17]:
def tan_repl_func(expr):
    """Replace tan with sin/cos"""
    if isinstance(expr, tan):
        return sin(*expr.args)/cos(*expr.args)
    elif not expr.args or expr.is_Derivative:
        return expr

tan_repl = lambda expr: crawl(expr, tan_repl_func)

Testing it:

In [18]:
tan_repl(tan(a) + tan(a)/tan(b))
Out[18]:
sin(a)/cos(a) + sin(a)*cos(b)/(sin(b)*cos(a))

So that works as expected. Now for the second pass; the subs with denominator checking and selective simplification. This takes a little bit more code than before, but I've heavily commented it so you should be able to see what's going on:

In [19]:
def smart_subs_func(expr, sub_dict):
        # Decompose the expression into num, den
        num, den = fraction_decomp(expr)
        if den != 1:
            # If there is a non trivial denominator, we need to handle it
            denom_subbed = smart_subs_func(den, sub_dict)
            if denom_subbed.evalf() == 0:
                # If denom is 0 after this, attempt to simplify the bad expr
                expr = simplify(expr)
            else:
                # Expression won't result in nan, find numerator
                num_subbed = smart_subs_func(num, sub_dict)
                return num_subbed/denom_subbed
        # We have to crawl the tree manually, because `expr` may have been
        # modified in the simplify step. First, perform subs as normal:
        val = sub_func(expr, sub_dict)
        if val is not None:
            return val
        new_args = (smart_subs_func(arg, sub_dict) for arg in expr.args)
        return expr.func(*new_args)

def fraction_decomp(expr):
    """Return num, den such that expr = num/den"""
    if not isinstance(expr, Mul):
        return expr, 1
    num = []
    den = []
    for a in expr.args:
        if a.is_Pow and a.args[1] == -1:
            den.append(1/a)
    else:
        num.append(a)
    if not den:
        return expr, 1
    num = Mul(*num)
    den = Mul(*den)
    return num, den

Finally, we can put everything from above inside a nice wrapper function:

In [20]:
def smart_subs(expr, sub_dict):
    """Performs subs, checking for conditions that may result in `nan` or 
    `oo`, and attempts to simplify them out.

    The expression tree is traversed twice, and the following steps are
    performed on each expression node:
    - First traverse: 
        Replace all `tan` with `sin/cos`.
    - Second traverse:
        If node is a fraction, check if the denominator evaluates to 0.
        If so, attempt to simplify it out. Then if node is in sub_dict,
        sub in the corresponding value."""
    expr = crawl(expr, tan_repl_func)
    return smart_subs_func(expr, sub_dict)

Let's see if it works as expected. Using a simple test case:

In [21]:
test = tan(a)*cos(a)
sub_dict = {a: pi/2}
print("Without `smart_subs`:")
print(new_subs(test, sub_dict))
print("With `smart_subs`:")
print(smart_subs(test, sub_dict))
Without `smart_subs`:
nan
With `smart_subs`:
1

And some timings:

In [22]:
print("Using `smart_subs`:")
%timeit smart_subs(test, sub_dict)
print("Using simplification, then normal subs")
%timeit simplify(test).subs(sub_dict)
print("Using trigsimp, then normal subs")
%timeit trigsimp(test).subs(sub_dict)
Using `smart_subs`:
10000 loops, best of 3: 92 µs per loop
Using simplification, then normal subs
10 loops, best of 3: 42.9 ms per loop
Using trigsimp, then normal subs
10 loops, best of 3: 30.8 ms per loop

Using selective simplification, the same results can be obtained for a fraction of the cost. 360 times faster for this small test.

Let's see what the overhead of smart_subs is for an expression that doesn't need it:

In [23]:
test = sin(x + y)/cos(x + y)* (a + b) + a/(a + x)
sub_dict = {x: a, y: 1}
print("new_subs")
%timeit new_subs(test, sub_dict)
print("smart_subs")
%timeit smart_subs(test, sub_dict)
print("Using simplification, then normal subs")
%timeit simplify(test).subs(sub_dict)
print("Using trigsimp, then normal subs")
%timeit trigsimp(test).subs(sub_dict)
new_subs
1000 loops, best of 3: 222 µs per loop
smart_subs
1000 loops, best of 3: 2.34 ms per loop
Using simplification, then normal subs
1 loops, best of 3: 1.11 s per loop
Using trigsimp, then normal subs
1 loops, best of 3: 815 ms per loop

So there's considerable overhead, which was expected. Still, it's much faster than using simplify first, and then running subs.

The best use method would be to first try it with new_subs, and if you get a nan or oo, then try using smart_subs. To aid in this, we can write a nice wrapper function msubs that contains both methods:

In [24]:
def msubs(expr, sub_dict, smart=False):
    """A custom subs for use on expressions derived in physics.mechanics.

    Traverses the expression tree once, performing the subs found in sub_dict.
    Terms inside `Derivative` expressions are ignored:

    >>> x = dynamicsymbols('x')
    >>> msubs(x.diff() + x, {x: 1})
    Derivative(x, t) + 1

    If smart=True, also checks for conditions that may result in `nan`, but
    if simplified would yield a valid expression. For example:

    >>> (sin(a)/tan(a)).subs(a, 0)
    nan
    >>> msubs(sin(a)/tan(a), {a: 0}, smart=True)
    1

    It does this by first replacing all `tan` with `sin/cos`. Then each node
    is traversed. If the node is a fraction, subs is first evaluated on the
    denominator. If this results in 0, simplification of the entire fraction
    is attempted. Using this selective simplification, only subexpressions
    that result in 1/0 are targeted, resulting in faster performance."""

    if smart:
        func = smart_subs
    else:
        func = new_subs
    if isinstance(expr, Matrix):
        # For matrices, func is applied to each element using `applyfunc`
        return expr.applyfunc(lambda x: func(x, sub_dict))
    else:
        return func(expr, sub_dict)
In [25]:
#Check that the results are the same:
msubs(test, sub_dict) == new_subs(test, sub_dict)
msubs(test, sub_dict, smart=True) == smart_subs(test, sub_dict)
Out[25]:
True

This code has been included in my optimization branch, and is performing admirably against all the included tests (so much speed up!!!). The big test though is the bicycle example.

The Big Test

To really put this new code to the test I applied it to the bicycle example, which has operation counts ranging from 400,000 to several million depending on which expression you're working with. How did it fair? Mixed results...

Using msubs with smart=False in the formation of the linearizer object resulted in a huge speed increase. Previously using _subs_keep_derivs there resulted in a run time of several hours. Now it runs in 14 seconds!

The second test is in the application of the operating point to the M, A, and B matrices. Previously evaluating these had resulted in nan and oo, and with operation counts exceeding 400,000 for M pre-simplification is out of the question. I tried using msubs with smart=True on this M, and let it run for 6 hours. It ended up consuming 98% of my RAM (4 GB worth), and it still wasn't done :( So for reasonable sized expressions the smart_subs we've implemented is acceptable, but it still doesn't work for the huge expressions. I'll have to keep working on optimizations to make this faster/reduce the initial expression size.

Still, all is not lost. Going from formulation to M, A, and B without the operating point substitution now only takes 20 minutes - down from the several hours before. This is actually faster than the previous linearize method, but unforunately results in expressions that evaluate to nan.

Future Work

Per my original timeline I was supposed to be working on "pre-linearization" routines for generating the equations of motion directly in linearized form. I think this will be dropped in favor of cleaning up and speeding up the existing codebase. There's some dead code that needs removing, and I have some other ideas for speed-ups that I'd like to try.

After that, I hope to get to the original part 3 of the project, which is adding support for matrices to the code-generation code. Not sure if this would live in sympy or pydy, but it would be extremely useful to myself and others. For now I plan on taking it week-by-week.


If you have thoughts on how to deal with the nan and oo issue, please let me know in the comments. Thanks!

by Jim Crist at July 06, 2014 10:25 PM

July 05, 2014

Sushant Hiray

This Week in CSymPy: #6-#7

Week 7 for GSoC just ended. I successfully cleared the mid-term review, thanks @certik

Progress

I started the week 6 by implementing refactoring the various instances of subs by adding a method create to construct classes with canonicalization. #210 was the relevant issue which was raised. I opened pull 213 to address this issue. This was a good addition as it reduced the line of codes by 117.



After successfully resolving #210 I started working on the Complex Module. By the end of the sixth week, I had discussed with Ondrej regarding the basic structure for the Complex Class. We had decided to use the mpq_class itself for the real and the imaginary part, as opposed to creating an alternate real or complex part.

Pull 223 implements the complex module. The class design structure could be found in the PR.

In Week 7 I was travelling back from University to home, so I couldn’t contribute much code for a couple of days. By the end of Week 7, I had managed to implement some part of the basic API. All the virtual functions are yet to be implemented, so right now I’m not returning an object in from_mpq

I will speed up during the next week and hope to merge the Complex Module.

Week Highlights!

Somewhere around this period, I crossed 100 commits into CSymPy master. Its been a long journey since I first started contributing to CSymPy.

The Week Ahead

Complete the basic API for Complex Module and merge it into master.


Thats all for now :) Will get back next week!

July 05, 2014 01:45 PM

Avichal Dayal

This week was unproductive coding-wise. I spent most of the time designing and writing pseudo-code.

For the Formal Power Series part of my project, last week I implemented the algorithm to find general term for series expansion of a function. Also various operations on series like addition, multiplication, composition were implemented.
I store the coefficients of series in a Stream class which works on lazy evaluation scheme.

However the Stream class currently is not a very generic one. It can only be used to represent series expansion of a function with parameters like "x", "general_term" and methods for series arithmetics.

To make it more generic, I'll make the Stream class act only as a container to store elements (infinitely many in lazy scheme).
To represent infinite series, there will be another class - Lazyseries. It will store its coefficients using the Stream class and will contain the methods for series arithmetics.
This way, the stream class can also be used to store other infinite structures. E.g.:- Fibonacci series can be represented in the following way using "Stream"

 def fib(x, y):  
head = (x, y)
tail = lambda: fib(y, x+y)
return Stream(head, tail)


But I'm a bit unsure if this is the best way for lazy evaluation. Currently I store the first element and the rest of the elements as a lambda function. Lambda function delays the evaluation of the tail. I don't know if this is the most efficient method in terms of memory and speed.

As an example:- Multiplication of two series is defined as:-
head = head1 * head2
tail = head1 * tail2 + head2 * tail1 + tail1 * tail2
Now when tail is evaluated continuously, many stream objects are created in the memory resulting in memory overflow. I'm working on this.

by Avichal Dayal (noreply@blogger.com) at July 05, 2014 08:53 AM

July 04, 2014

Tarun Gaba

GSoC 14: Midterm Evaluations have ended!

{% include JB/setup %} [ <-Back to posts ](/gsoc14) It has been a week after the midterm evaluations are over, and I am back to work after a small break(with permission from my mentor, off course!). I have been working on writing a test suite for the Dynamics Visualizer. This is the wrapping up part of the visualizer for this gsoc. [Here](/blog/visualization/index.html?load=samples/scene_desc.json) is a visualization of a rolling disc(it is slightly buggy though), that i prepared. To view the animation, allow the visualization to load in the webpage(it shall load automatically), and then hit the `Play Animation` button. After writing some tests for visualizer, I am going to start fleshing out API for the module, to provide IPython support to the visualizer. The main aim of writing this module is to make visualizer interactive, in the sense, that a user should be able to change all the variables from the GUI(which is rendered inside notebook's output cell) and then rerun the simulations without having to write any code, or execute any of the code manually. The data of the new simulations will be automatically fed into visualizer, and then it can be viewed as an animation. This whole workflow will be very convenient for the existing PyDy users, as well as the new ones. It will be particularly convenient for those who want to just play around with the existing systems, by changing the system variables, and view how it affects the resulting animations. With the development of this module, as well as ongoing improvements in the other PyDy modules(by my fellow GSoC'ers from PyDy), we should be able to perform lightening fast simulations for a system, as well as view them on a canvas. I will keep posting the new work I will be doing, with better details(once I actually start implementing new stuff!). [ <-Back to posts ](/gsoc14)

July 04, 2014 05:40 PM

July 03, 2014

Soumya Dipta Biswas

GSoC 2014: Week 6

Welcome back everyone,

This week I set about to fix one of my long pending PR associated with Propositional Logic. Following that I worked on conversion to Prenex Normal Form and Skolem Normal Form. Firstly, let me talk about the PR associated with Propositional Logic. It basically involved adding small functions like validity, entailment and some small modifications to the Propositional Knowledge Base. I assume most of the readers would be familiar with the concept of validity and entailment but I will briefly mention the algorithm to compute this and the intuition behind the same. Validity is given by not satisfiable(Not(expr))</expr>. Why does this work? If an expr is valid then it is True under all interpretations. Therefore Not(expr) is always False i.e. it cannot be satisfied by any model. Hence if the negation of the expr is unsatisfiable then the original formula is valid. Moving on to entailment, the conventional definition stands that a set of formulas formula_set entails a formula expr iff expr is True whenever formula_set is True. That basically boils down to And(*formula_set) >> expr. However there is another method, to accomplish the same, namely not satisfiable(And(~expr, *formula_set)). Ok, firstly why does this non-intuitive expression even work? Well whenever formula set is True, expr must be True (for the entailment to hold) hence ~expr must be False thereby making the entire conjunction False. If the formula_set is False then the conjunction is clearly False. So if the entailment is to hold then the conjunction must be unsatisfiable, which is exactly what was implied earlier. Let us look at one special case before moving on. What if the formula_set is empty. While the solution is a matter of convention to a good extent, SymPy uses ideas from the second formula to resolve the ambiguity which yields not satisfiable(And(~expr)) which is exactly the same as validity of expr. Theoretically entailment means “Given this what can we infer”. So, given nothing we can only infer that, which is always the truth. Except for these changes I also made changes to PropKB which is probably best not discussed at the moment.

Coming now to the major interest of the post. The major time was spent in the functionality for conversion to Prenex Normal Form and Skolem Normal Form. Without going into a lot of details let me introduce these to you. A formula in PNF consists of a sequence of Quantifiers (containing both Universal and Existential quantifiers) followed by a Quatifier-less expression. The series of Quantifiers is called a prefix while the expression is called the matrix. Now, every formula in FOL can be converted to an equivalent formula in PNF. The matrix can be compared to a formula in propositional logic (if we think of the predicates as literals of PL) and most of the operations of PL become relevant. However we cannot completely ignore the prefix and hence we try to perform some operations such that we can look at the matrix without any regard to the prefix. This brings us to the Skolem Normal Form. A formula is converted to SNF by first converting it to PNF then dropping the Existential quantifiers (according to some rule). Now the prefix contains only Universal quantifiers, which are assumed to be implicitly implied and hence simply ignored while looking at the matrix. I will be updating the various rules for conversion and examples soon. Till then stay tuned.

Arrivederci

by SD at July 03, 2014 06:13 PM

July 02, 2014

Harsh Gupta

week 6

Solving Trigonometric Function (part II)

There is another technique to solve trigonometric function. Just as every trigonometric function can be written in term of \( \tan \) it can also be written in terms of \( \exp \).

$$ sin(x) = - \frac{i}{2} \left(e^{i x} - e^{- i x}\right) $$ $$ cos(x) = \frac{e^{i x}}{2} + \frac{1}{2} e^{- i x} $$ $$ tan(x) = \frac{i \left(- e^{i x} + e^{- i x}\right)}{e^{i x} + e^{- i x}} $$ $$ cot(x) = \frac{i \left(e^{i x} + e^{- i x}\right)}{e^{i x} - e^{- i x}} $$

So, solving a trigonometric equation is equivalent to solving a rational function in \( \exp \). Note: here the \( \exp \) is in complex domain and equation \( exp(x) = y \) has solution \( \left\{i \left(2 \pi n + \arg{\left (y \right )}\right) + \log{\left (\left\lvert{y}\right\rvert \right )}\; |\; n \in \mathbb{Z}\right\} \) when solved for \( x \).

by Harsh Gupta at July 02, 2014 08:43 AM

Sudhanshu Mishra

GSoC'14 progress: Week 6

It was a busy week though I managed to do some work.
I've been working on following things:
  • Integrating TWave with refraction_angle(7626)
  • Interference class
Some blocking 3D geometry code got merged yesterday which gives me a lot of new things to work on. :-)
I also tried to take a look on very old patch for gaussopt.
This week I wrote a script to update the development documentation when a PR gets merged into the master. This is up and running.
I passed mid term evaluation that took place last week and I would like to thank SeanOndrej,Aaron and Chris for their constant support.
That's all for now. Cheers!

by Sudhanshu Mishra (noreply@blogger.com) at July 02, 2014 03:27 AM

July 01, 2014

Akshay Narasimha

Gsoc Week-6

Sorry for the late post ,had an issue with my data card . First of all I couldn't get my PR's (Line3D, Plane) merged as Stefan was busy, though I addressed all of his issues and hopefully it is time before they get merged.

I started working on the Hyperbola class, you can find the code here.
This week I plan to continue the work on the Hyperbola class and hopefully will get my previous PR's merged.

Until then cheers!

by Akshay Narasimha (noreply@blogger.com) at July 01, 2014 04:22 AM

June 29, 2014

Kundan Kumar

GSoC Week 6: Docs of system of ODEs

This week has been hell of a week of work and enjoyment with marriage preparation of my cousin. This lead to very less progress in my project work. I was unable to do much with just completion of docs of all implemented methods of system of ODEs. The non-linear system of ODEs which was left on last week has not been much progressed since then.

I have confusion too on representation of few things which I posted on mailing list. I need community opinion on these things before I proceed on these facts. This will be all for this week. Rest I intend to implement non-linear system of first order for 3 equations during next week.


by Kundan at June 29, 2014 04:52 PM

Sachin Joglekar

GSoC Week 6: Completing coordinate systems

Phew. This was a busy week. Initially, the plan was to push a WIP PR (without tests maybe) with the code for coordinate systems and point classes. Midway through the week, Jason and I decided to push this new code to the first PR itself, and merge the entire fundamental framework together – in a single PR. Initially I did get a little worked up looking at the amount of work that needed to be done.

However, things moved much faster than expected – though with a lot of issues along the way. There were problems with the args for CoordSysRect, then some issues with the Tree-algorithm for inter-point distance calculation, and then some more with the code that implements orientation of systems wrt each other…you get the drift. But the code is finally done and polished, along with the unit tests and doctests. Thats why the late post this week – I ‘bravely’ decided to get the code fully working, without glitches, and then only do the blogpost for the week. Thankfully, the plan was a success :-D.

Some minor things still remain, like renaming of some attributes to follow conventions  – mostly cosmetic changes, nothing that will affect how the code works.

The next immediate steps now would be -

1. Finish off the PR and get it merged.

2. Study caching mechanisms and implement them in a new ‘patch’ PR.

The first PR is here.

See you later then, next week :-).


by srjoglekar246 at June 29, 2014 11:28 AM

June 28, 2014

Avichal Dayal

<script src="http://latex.codecogs.com/latexit.php?p&amp;li&amp;div" type="text/javascript"></script>
A bit late for this week's blog.

Implementing the solveRE() method took longer than I expected. I spent most of this week doing that. This method solves recurrence equation to get the general term for series expansion of a function.
rsolve() can solve some equations but not recurrences of hypergeometric type.

Hypergeometric recurrence equation is of type:-
Q(k) * r(k + m) = P(k) * r(k)
where P(k) and Q(k) are rational functions in k.

For the case when m = 1, there exists a simple formula to get r(k).
When m > 1, we can represent it as sum of m-shifted m-fold symmetric functions.
For e.g.:- Let the RE be,
r(k + m) = R(k) . r(k)
We replace k with mk and r(mk) with c(k) we get,
r(m(k+1)) = R(mk) . r(mk) which is equivalent to
c(k+1) = R(mk) . c(k)
Now we can use the same formula for the case when m = 1.

Similarly by replacing m with km + 1, km + 2 ...., we can solve for r(k) as sum of m-shifted m-fold symmetric functions.

Here are some results:-


The general terms are not in their simplest form. RisingFactorial must be re-written in terms of factorial to bring it to the usual form.

Now, I'm able to find FPS or infinite series for a wide range of functions. Next week, I'll take care of efficiency and making the code more in line with SymPy's style.

by Avichal Dayal (noreply@blogger.com) at June 28, 2014 07:21 AM

Jim Crist

GSoC Week 6: Just the little things

I was rather busy with my research this week, so no time for a long-winded post like some of my previous updates. There's not much interesting to say anyway. This week was mostly spent on little fixes to get my current pull request merged.

Topping the list of things that are better than they were last week is speed. The profiling I did last week showed that the current function sympy.physics.mechanics uses to solve a system of linear equations (_mat_inv_mul) is sloooooooooow. The underlying reason is because subs is slow - more on that later. I spent some time swapping out all forms of solving ($A x = B$) for LUsolve, the clear winner of last weeks benchmarks. This resulted in a 10x speedup of the formulation of equations for the bicycle model example.

This bicycle example has become the bane of my existence for the last couple weeks. It's a super slow test that I'd never actual gotten to run before. But with the speed improvements made, it actual finishes in a reasonable time. Except it still doesn't work. I'm able to run all the way up to

M, A, B = KM.linearize()

But when I go to sub in values for symbols in these matrices, things get hairy. There are two issues:

Issue 1: Get nan when not simplified

M.subs(val_dict) results in nan and oo upon after subs. But doesn't if it's simplified before the subs. An example of this behavior would be:

>>> M = sin(q1)/tan(q1)
>>> M.subs({q1: 0}
nan

Note that if this is simplified, this results in something completely different:

>>> M = sin(q1)/tan(q1)
>>> M = M.trigsimp()
>>> M
cos(q1)
>>> M.subs({q1: 0})
1

However, for the bicycle case M has over 19 thousand operations. This doesn't simplify quickly. Also, by default we don't simplify before subs in Linearizer (you can opt in to simplify, but it's done right before the return, so it won't affect the subbed result at all). Right now I'm looking through ways to make the resulting expressions smaller after the formulation, as this will result in speedups for all operations. This could be extremely helpful for issue 2...

Issue 2: subs is slow

because A has over 38 million operations!!! In this case subs doesn't even return. Ever. I left it running on my computer for 4 hours and came back and it was still whirring along, fans on high, eating up all my ram. No idea how to solve this. One possible solution is csympy, a fast core written in C++. Once this matures, subs, trigsimp, and other time consuming operations used heavily in sympy.physics.mechanics could rely on the equivalent, faster, C++ versions. I filed an issue with an example expression generated from the bicycle example (this one only had 147,841 operations, not nearly as bad). Hopefully Ondrej and the team can use this as a benchmark problem to help improve subs in csympy.

If you have thoughts on how to overcome these issues, please let me know. I'm kind of stumped right now.

The Good News

I didn't want to end this post on a bad note, so I'll close with the remainder of the things I did last week that actually worked:

  1. Improved documentation! Docstrings that are worth reading, and a start on the sphinx documentation.

  2. Added a deprecation warning for KanesMethod.linearize to warn people about the method change.

  3. Major interface changes. Now all operating points are specified as a single dictionary, or an iterable of dictionaries. This is to aid in consistency across different system implementations. Referring to a dictionary as u_op in LagrangesMethod doesn't really make any sense, as Lagrange's method only uses $q$, $\dot{q}$, and $\ddot{q}$. Also added a kwarg to make simplification of the results optional.

  4. Added a method to the LagrangesMethod class to calculate the value of the multipliers at different points. This is useful for multiple reasons. The multipliers have meaning, so knowing what the solution is symbolically is nice for calculating the constraint forces. Also, when linearizing with Lagrange's method, the multipliers have operating points as well, and these need to be calculated based on the operating point for the other states ($q$, $\dot{q}$, etc...). Now a user can go:

    op_point = dict_or_iterable_of_dicts
    lam_op = LM.solve_multipliers(op_point)
    op_point.append(lam_op)     # Or op_point.update if op_point is a dict, not a list of dicts
    M, A, B = LM.linearize(q_ind=q_ind, qd_ind=qd_ind, op_point=op_point)
    

Hopefully in the next week I can get my PR merged, so the Lagrange stuff can finally be submitted.

<script type="text/javascript"> if (!document.getElementById('mathjaxscript_pelican_#%@#$@#')) { var mathjaxscript = document.createElement('script'); mathjaxscript.id = 'mathjaxscript_pelican_#%@#$@#'; mathjaxscript.type = 'text/javascript'; mathjaxscript.src = 'https:' == document.location.protocol ? 'https://c328740.ssl.cf1.rackcdn.com/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' : 'http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'; mathjaxscript[(window.opera ? "innerHTML" : "text")] = "MathJax.Hub.Config({" + " config: ['MMLorHTML.js']," + " TeX: { extensions: ['AMSmath.js','AMSsymbols.js','noErrors.js','noUndefined.js'], equationNumbers: { autoNumber: 'AMS' } }," + " jax: ['input/TeX','input/MathML','output/HTML-CSS']," + " extensions: ['tex2jax.js','mml2jax.js','MathMenu.js','MathZoom.js']," + " displayAlign: 'center'," + " displayIndent: '0em'," + " showMathMenu: true," + " tex2jax: { " + " inlineMath: [ ['$','$'] ], " + " displayMath: [ ['$$','$$'] ]," + " processEscapes: true," + " preview: 'TeX'," + " }, " + " 'HTML-CSS': { " + " styles: { '.MathJax_Display, .MathJax .mo, .MathJax .mi, .MathJax .mn': {color: 'black ! important'} }" + " } " + "}); "; (document.body || document.getElementsByTagName('head')[0]).appendChild(mathjaxscript); } </script>

by Jim Crist at June 28, 2014 02:21 AM

June 27, 2014

Thilina Rathnayake

[GSoC] Week 6: Cholesky and LDL Algorithms

This week I implemented Cholesky decomposition and LDL decomposition. In addition to that I fixed two errors In CSymPy. I was also bale to finish work with LU decomposition and merge it to master. Also, I could solve the system Ax = b using fraction free LU factorization.

Cholesky Decomposition

Cholesky decomposition can be applied to a Hermitian positive definite matrix. Hermitian matrix is a matrix with complex entries that is equal to it’s conjugate transpose [1]. Hence a symmetric matrix with real entries can be considered as a Hermitian matrix and can be decomposed using Cholesky decomposition if it’s positive definite. A Symmetric matrix A is positive definite if z^TAz is greater than zero for every non zero column matrix z [2]. If the above conditions are satisfied, Cholesky decomposition for a matrix A can be written as A = LL^* where L is an lower triangular Matrix. This is equal to A = LL^T when L is a real matrix. This factorization can be used for fast solution of the system Ax = b. I am yet to use this decomposition in solving above mentioned system.

LDL Factorization

LDL decomposition is closely related to Cholesky decomposition. As the name implies, in LDL decomposition of a matrix A can be written as A = LDL^* where L is a lower triangular matrix and D is a diagonal matrix [3]. This decomposition can be used for some matrices which don’t have a Cholesky decomposition.

CSymPy printing error and simplification errror

I also worked on a printing error of CSymPy and a simplification error in Mul class which is used to represent multiplication types in CSymPy. There is still some work to be done to fix simplification error completely. The most important thing was that we were able to introduce a fix which doesn’t create a considerable slow down in speed after being applied.

References

[1] Hermitian Matrix, Wikipedia Article: http://en.wikipedia.org/wiki/Hermitian_matrix

[2] Positive definite Matrix, Wikipedia Article: http://en.wikipedia.org/wiki/Positive-definite_matrix

[3] LDL Decomposition, Wikipedia Article: http://en.wikipedia.org/wiki/Cholesky_decomposition#LDL_decomposition


by Thilina Rathnayake at June 27, 2014 03:58 PM

June 26, 2014

Harsh Gupta

week 5

Solving Trigonometric Function (part I)

This week I spend time on making trigonometric solvers work. Every trigonometric function can be written in terms of tan.

$$ sin(x) = \frac{2*tan(x/2)}{tan^{2}(x/2)} $$

$$ cos(x) = \frac{-tan^{2}(x/2) + 1}{tan^{2}(x/2) + 1} $$

$$ cot(x) = \frac{1}{tan(x)} $$

A basic technique to solve trigonometric equations can be rewriting the equation in terms of tan. And if the equation is made by addition, multiplication or quotient of trigonometric functions then the transformed equation is a equivalent to a rational function in tan. That equation can be solved by the usual polynomial solving techniques.

Taking the example from the doc \( cos(x) + sin(x) \) gets converted to \( \frac{-tan^{2}(x/2) + 2*tan(x/2) + 1}{tan^{2}(x/2) + 1} \)

The solution of this equations is \( tan(x/2) = 1 +- sqrt(2) \). Since the inverse of tan is \( \left\{2 \pi n + \operatorname{atan}{\left (y \right )}\; |\; n \in \mathbb{Z}\right\} \) the solution of the given equation is $$ \left\{2 \pi n - \frac{\pi}{8}\; |\; n \in \mathbb{Z}\right\} \cup \left\{2 \pi n + \frac{3 \pi}{8}\; |\; n \in \mathbb{Z}\right\} $$

Though it appears this technique should work universally for trigonometric equation it fails for even \( sin(x) = 0 \). From the table above \( sin(x) = \frac{2*tan(x/2)}{tan^{2}(x/2)} \) So, the \( sin(x) = 0 \) occurs at \( tan(x/2) = 0 \) which has solution \( \left\{2 \pi n\; |\; n \in \mathbb{Z}\right\} \) But the solution is \( \left\{ \pi n\; |\; n \in \mathbb{Z}\right\} \) . Why are we missing some solutions? The reason is \( sin(x) = 0 \) also occurs when denominator tends to \( \infty \), i.e., the values where \( tan^{2}(x/2) + 1 \) tends to \( \infty \). We had encountered a similar problem for the solution of $$ \frac{1}{\left(\frac{x}{x + 1} + 3\right)^{2}} $$

here \( x = -1 \) is not a point in the domain of the of the equation. The solver simplifies the equation to

$$ \frac{\left(x + 1\right)^{2}}{\left(4 x + 3\right)^{2}} $$

which extends the domain to include the point \( x = -1 \) which is also the solution to the transformed equation. There we wrote a sub procedure domain_check to verify if the returned solution is part of the domain of the original equation. The problem here is slightly different in the sense that transforming the equation decreases the domain of the solutions and not increase it.

To find such solution we have allow \( \infty \) to be solution to equations, we will be working on extended reals instead of just reals. I think this change will simplify a lot of things.

Another thing which should be taken care off is that we cannot naively search for the values for which the denominator tends to infinity as for the same value numerator might also attain infinitely large value, we will have to conceder the limiting value of the equation.

by Harsh Gupta at June 26, 2014 11:43 AM

June 25, 2014

Soumya Dipta Biswas

GSoC 2014: Week 5

Hello Folks,

This week was probably not as productive as I would have liked it to be. Firstly, it seems like the algorithm for faster conversion to CNF/DNF I talked about before works fast conditionally. If a particular formula is already quite near to a particular normal form, then the conversion takes a longer time than the recursive algorithm. The overhead incurred for the conversion in this case starts to dominate the actual conversion. So, I am now parallelly working on fixing this problem and the FOL module (as scheduled). Hence, I majorly managed to do only 2 things this week namely Interpretation and conversion to Prenex Normal Form.

Interpretation
Akin to pl_true(expr, model) present in the propositional logic module, I have implemented a fol_true(expr, model) for the FOL module. Before jumping into the concept of interpretation for First Order Logic, let us see what it means for Propositional Logic. If you are already familiar with the same, feel free to jump to the next paragraph. Now, the interpretation of an expression is the value of expression under an assignment. So given the formula A | B such that {A: True, B: False} the formula is True. So, this is essentially the evaluated result of a formula given an assignment. This is quite similar to a + b > 0 which has no value inherently, but if a = 1 and b = 0, then the expression has a value (True). In some cases only a partial assignment is sufficient to determine the value of the expr. For e.g. A & B with A = False is clearly False without any regard for the value of B.

In Propositional Logic, the most basic elements are symbols which can have one of True or False as its value. FOL on the other hand has many elements which can take on any value in the domain. While the concept of interpretation remains the same there, the number of things to be evaluated change. Interpretation in FOL is a rather complex job which will become quite apparent very soon. Let us set about to find the value of each of the basic elements.

  • Constants are the easiest to evaluate, well because they stand for themselves.
  • Variables are similar to variables in any mathematical expression and it is sufficient to simply give them a value. However, one may also give a variable a domain. This domain is the set of all values that the corresponding variable may assume.
  • Next come the functions. For functions, one needs to provide a mapping of every possible set of arguments to its corresponding value. For e.g. for f(X, Y) such that X = {1, 2, 3} and Y = {10, 20, 30}, one needs to provide a value for every feasible possible combination of X and Y.
  • Predicates are almost the same as functions with the exception that the values they yield must be boolean. The examples below will probably make things much clearer.

Let’s look at an example that will give better perspective to the idea of Interpretation.

>>> from sympy.abc import X, T
>>> from sympy.logic.FOL import Predicate, ForAll, Exists, fol_true
>>> Person = Predicate('Person')
>>> Time = Predicate('Time')
>>> CanFool = Predicate('CanFool')
>>> X_domain = ['John', 'Jack']
>>> T_domain = [1, 2, 3]
>>> def person(X): return X in X_domain
>>> def time(T): return T in T_domain
>>> CanFoolMap = {('John',2):False, ('John',3):False, 'default':True}
>>> model = {X:X_domain, T:T_domain, Person:person, Time:time, CanFool:CanFoolMap}

# You can fool some of the people all of the time
>>> expr = Exists(X, ForAll(T, (Person(X) & Time(T)) >> CanFool(X, T)))
>>> fol_true(expr, model)
True

# You can fool all of the people some of the time
>>> expr = ForAll(X, Exists(T, (Person(X) & Time(T)) >> CanFool(X, T)))
>>> fol_true(expr, model)
True

# You can fool all of the people all of the time
>>> expr = ForAll(X, ForAll(T, (Person(X) & Time(T)) >> CanFool(X, T)))
>>> fol_true(expr, model)
False

I hope that makes the idea of interpreation in general and how to use them in SymPy a little clearer. I will update the post giving more examples soon.

Adios!!!

by SD at June 25, 2014 08:44 PM

June 24, 2014

Sushant Hiray

This Week in CSymPy: #5

Week 5 for GSoC just ended and this week I completed implementing the Exponential Module.

Progress

I started the week by implementing LambertW function. It wasn’t a part of the original proposal but it is a nice addition to have so I implemented it.

If you’re unaware, LambertW function is defined as the inverse function of x*exp(x). This function represents the principal branch of this inverse function, which is multivalued. For more information, see: wiki. Pull 200 implemented the LambertW class!

After this, I worked on Hyperbolic module, which was the main target for this week. Pull 203 implemented the hyperbolic module. In addition to the normal functions which were implemented as a part of the TrigFunction API, we implemented expand_as_exp in the Hyperbolic Module.

With this in place, we can now do things like sinh(x)->expand_as_exp() and it will return the exponential form: (e^x - e^{-x})/2 using a combination of exp, adds and divs!

Follow up to this PR we have openened a couple of issues:

  • #207 is a bug which I encountered while writing test cases. Essentially same symbolic expressions written in form of add and mul are not equated as equal. This is a crucial issue and we are yet to see how to proceed further on this. SymPy also faced the same issue sympy/sympy#4596 where it was decided to remove such automatic redistribution. As @certik suggested, it is perhaps the best to compare underneath via expand.

  • #210 is basically a more generic way to cleanup the existing Trignometric and Hyperbolic Module. This for instance can be used to define the subs method in the base class rather than writing almost the same code in each subsequent inherited class. The only reason why this was being defined in all the classes was that we needed to access the method to create the class via canonicalization. I’ve proposed a couple of methods in which we can approach the problem. Once decided we can fix this issue.

Discussions

As usual most of the actual discussions happened on PR, we were supposed to have a meetup on gitter but it was postponed.

Week Highlights!

  • CSymPy just crossed the 1000 commits mark! Yay! It is 1036 at the time of writing.
  • The issue serial number has crossed the 200 mark!

Quite an interesting week this has been as far as milestones are concerned!

The Week Ahead

  • Implement the long awaited Complex Module.


Thats all for now :) Will get back next week!

June 24, 2014 01:45 PM

Akshay Narasimha

Gsoc14 Week-5

I did some good amount of work this week. I finished the implementation of the Plane class. Here is the link to the code. I haven't sent a PR though as I am waiting for the merging of the Line3D class, which will happen soon enough as it had a few styling issues which I happen to fix in the latest commit.

Here is the implementation of the interaction between 2D and 3D entities with the use of Planes.

>>> a = Point(1, 2)
>>> b = Point3D(1, 3, 5)
>>> c = Plane(Point3D(x, y, z), normal_vector=[0, 0 ,1]) # Parallel to xy plane
>>> d = c.projection(c) # A 3D point
>>> d
Point3D(1, 2, 0)
>>> d.distance(b)
sqrt(26)

Here are the things that I have implemented till now
  • Point3D
  • Line3D
  • Ray3D
  • Segment3D
  • Plane
Among these only Point3D class has been merged but I will get most of the other classes merged this week, so that would complete the implementation of 3D Geometry as a whole.So I will spend most of this week trying to get these PR's merged.

Until then cheers!

by Akshay Narasimha (noreply@blogger.com) at June 24, 2014 04:16 AM

June 23, 2014

Sudhanshu Mishra

GSoC'14 progress, Week 5

I did good amount of work this week. I sent a small part of it related to ray tracing for review.
I got some great suggestions from Sean to make it compatible directly with TWave class. It is a WIP as I'm waiting for few 3D geometry classes. I'm also working on inteference of light waves and I hope that I'll send it for review in next couple of days. My last PR on medium is still unmerged and I really need it to be in master to work on more implementations.
I'm also trying to move physics/gaussopt to physics/optics/gaussopt as it'll be good to keep things related to optics at the same place. I'm waiting for opinion of the community about it. Here's the link to that PR:
Here is some other minor work that I've done this week outside of optics:
That's all for now. Cheers!

by Sudhanshu Mishra (noreply@blogger.com) at June 23, 2014 05:59 PM

June 22, 2014

Tarun Gaba

GSoC 14: Second Week!

{% include JB/setup %} [ <-Back to posts ](/gsoc14) Second week of GSoC '14 has ended. It has been a very hectic week. A good part of my time was spent in studying and tweaking the MGView code(whose fork we had decided to work on). The new _generic-visualizer_ has to be compatible with both PyDy as well as MG(MotionGenesis), so great care is being taken to maintain the backwards compatibility of the package, as well to add the support for PyDy visualization data. I have also worked on Python side, to modify the JSON data generated for the visualizations(from PyDy side). ###Accomplishments: There have been three major accomplishments this week: - **Studying MGView code:** This has to be the most time consuming part! I had to go over a module spanning multiple files to study the workflow, and how each methods are linked together. Some conversations with Adam Leeper(who wrote the original package) have been very useful in this regard. - **Modifying JSON generated from PyDy:** Certain modifications have been done on the JSON generated from the PyDy backend, so that they are compatible with MotionView. A major change is that the generated JSON data is split over two files, one `scene_desc.json` shall contain all the parameters required to create a static scene, and another file `simulation_data.json` which shall contain all the simulation data(in the form of 4x4 matrices). It would be very helpful in un-cluttering the data. The API changes in Python side are totally backward compatible, which means no previous code shall break once these changes are merged. - **Tweaking MGView code:** I have started working on MGView code, to provide the support for PyDy generated JSON data. This part will span the coming weeks too. ###Objectives: The objectives of the upcoming week are to be able to produce some visualizations on the frontend side. For that I will need to tweak into the MGView code to allow support for PyDy JSONs. ###Issues: The main issues I am facing is the inhomogenities in the output produced by PyDy and MG. Since the visualizer (MOtionView) must be able to visualize data generated from both of these packages, without any bias it has to be designed like so. Since It is already built to specifically visualize data generated from MotionGenesis software, There are certain parts of code which need to be totally rewritten, and it has to be done without breaking the workflow of the visualizer. So it can get somewhat messy to do so. I am working on to find a clean way to do so. The upcoming two weeks are somewhat crucial to the project timeline. I have to work really hard in the coming weeks to keep up with my deliverables before the midterm evaluations. I think I shall get back to coding now! :) [ <-Back to posts ](/gsoc14)

June 22, 2014 07:00 PM

Kundan Kumar

GSoC Week 5: Non-linear system of ODEs

Considering the 5th week has passed, PR for linear system of 3 equations of first order still remains to be merged. Also, PR for non-linear system of ODEs for 2 equations of first order is on the way. The non-linear system of ODEs which I have included are these -

>x’ = x**n*F(x, y)

   y’ = g(y)*F(x, y)

>x’ = exp(λx)*F(x, y)

   y’ = g(y)*F(x, y)

>x’ = F(x, y)

   y’ = G(x, y)

>x’ = f1(x)*g1(y)*Φ(x, y, t)

   y’ = f2(x)*g2(y)*Φ(x, y, t)

>x = t*x’ + F(x’, y’)

   y = t*y’ + G(x’, y’)

Though the representation of output of these equation are due, because of `Eq` representation not working with a list. So the representation of test cases is not complete.

There are some equations which may be of above format but is not solvable due to problem of dsolve with certain single equation ODEs, whose solutions, being entwined in terms of function and variable, are inseparable. This is also due to inefficiency of `solve` function that `solve` is not able seperate function and variable in a solution returned by dsolve.

This week’s I spent a lot of time correcting conflicts in git rebase which created a lot of mess. I tried to solve those but wasn’t very succesful so I deleted the previous branches and created new one.

This week I intend to complete the documentation part which I have left for implemented differential equations and complete the system of non-linear ODEs for 1st order.


by Kundan at June 22, 2014 10:28 AM

June 21, 2014

Sachin Joglekar

GSoC Week 5: PR reviews and starting CoordSysRect

This week had the dual focus of -

1. Polishing off the Vector framework based on PR reviews

2. Starting the implementation of CoordSysRect+Point

I had thought last week, that the vector classes I implemented in the first PR were the polished version. However, with the reviews given, I have come to realise that a lot of improvements had/have to be made. Francesco Bonazzi, one of the SymPy contributors, was a huge help in improving the code and making it more efficient than what it initially was.

The biggest suggestion for improvement that was provided was to not implement custom hashing and equality checking methods for the Vector classes. Initially, I had overridden the inherited versions of these methods in the new classes, based on tuples generated with respect to the vectorial components of any Vector instance. However, SymPy already has a strong framework for these functionalities, based on the args of a Basic object. Moreover, the SymPy core is apparently being re-written in C++. This would make the SymPy-core versions of the aforementioned methods much, much more efficient than my Python versions.

Luckily, I had implemented the ‘__new__’ methods of the Vector classes in such a way that ‘equivalent’ Vector instances had exactly the same args. This made the usage of the inherited methods much, much easier – the only thing that remained to do was to pass the appropriate args to the superclass __new__ methods. Hence, all Vector classes rely on the SymPy core now, for as much as they can. Phew. Still a problem is the BaseScalar class, whose superclass is Dummy. The __new__ method for Dummy makes it not-possible to pass coordinate-system and index-related information to it – making the appropriate hashing impossible. If I instead inherit from AtomicExpr, I will have to rewrite many methods – not sure if its a good idea. Lets see what the community has to say about this.

About the implementation of the new stuff, I have already finished implementing the Point class. Now I am currently working on implementing CoordSysRect – I am still writing the HUGE __new__ method for this class. For now, I am referring to the tests of sympy.physics.vector while doing the coding. Hope I am successful in writing the classes soon, and integrating it with a polishing and ‘perfected’ version of the vector core.

Thats all for this week. Have a great week!


by srjoglekar246 at June 21, 2014 11:45 AM

Thilina Rathnayake

[GSoC] Week 5: Matrix Decompositions

Hi All,

Week 5 of the GSoC is over and I have made some progress in Matrix factorizations. I implemented LU factorization and QR factorization. Still I have to implement Cholesky factorization. I hope to do it during the Weekend.

LU factorization

LU factorization was invented by Alan Turing and this is very important when solving a system of equations. In LU decomposition we try to factorize a given square matrix into a product of two matrices, L and U where L is a lower triangular matrix and U is a upper triangular matrix. The factorization can be a true factorization i.e. A = LU or it may be the process of finding two matrices L, U which can be used in solving the system Ax = B. In the latter case we don’t have A = LU. I implemented three different algorithms for LU decomposition.

Fraction Free LU decomposition in [1]

Fraction free version described in [1] results in two matrices, a lower triangular matrix L and a upper triangular matrix U. But this is not a true factorization i.e. A != LU in this case. The algorithm is more focused in a factorization that helps to solve a system Ax = B. So we are only interested in finding two matrices which can be used in backward / forward substitution to solve the system.

Normal and Fraction Free version in SymPy

I also implemented the normal and fraction free version of the algorithm in SymPy which gives an exact factorization A = LU. Normal algorithm results in fractions and not good for numerical calculations if you are not using arbitrary precision arithmetic. In contrast fraction free versions do not generate these kinds of fraction as a result of the decomposition. These algorithms can be found in [2] and [3] respectively.

QR decomposition

In QR decomposition we find two matrices, an [orthogonal matrix](http://en.wikipedia.org/wiki/Orthogonal_matrixand  a upper triangular matrix R such that A = QR. Here A doesn’t have to be a square matrix.There are several method to calculate QR factorization but SymPy uses Gram–Schmidt process. The Algorithm can be found in [4].

Issues

QR decomposition is not a `symbolic friendly` algorithm as it involves finding square root. I searched thoroughly for an algorithm that is suited for symbolic entries but couldn’t find one.

References

[1] Algorithm 3, page 14, Nakos, G. C., Turner, P. R., Williams, R. M. (1997). Fraction-free algorithms for linear and polynomial equations. ACM SIGSAM Bulletin, 31(3), 11–19. doi:10.1145/271130.271133.

by Thilina Rathnayake at June 21, 2014 06:20 AM

June 20, 2014

Jim Crist

GSoC Week 5: Adventures in Profiling

This week I've been working on little fixes. Improving the docstrings, a couple interface improvements, finalizing some design decisions. The plan is to finish up this work next week, and then move on to the code generation portion of the project. In it's current state, the linearizer can handle KanesMethod and LagrangesMethod with ease. I'll probably spend some time hacking away at the matrix_to_linearizer function as well, but that's not priority.

The rest of the week was spent profiling and doing some basic optimizations. Being a mechanical engineer, I haven't spent much time learning about numerical methods (although, I am taking this MIT OCW course this summer). As such, my natural way of solving the system

\[ A x = B \]

where \(A\), \(B\), and \(x\) are matrices, is to take the inverse of \(A\), and multiply it by \(B\):

\[ x = A^{-1} B \]

Turns out, this is horribly inefficient. But what is the best way? I had to do some profiling.

Creating Benchmark Matrices

For most systems of the form \(A x = B\) in sympy.physics.mechanics, \(A\) is a symmetric matrix, and \(B\) is a column vector. To do proper benchmarks, I'd need to create random matrices of this form. Looking at some example equations of motion, it can be seen that these mostly consisted of products and sums of terms composed of symbols, trigonometric functions, power functions, sqrt, and inv. After some time, I was able to create a couple functions to create these matrices in a composable way:

In [1]:
# Below are some functions to create random matrices of varying sizes.
# These will be used in the benchmarks.

from random import sample, randint
from sympy import Matrix, symbols, sin, cos, tan, sqrt, zeros, \
    Symbol, diag, prod

# Define some operations
pow2 = lambda x: pow(x, 2)
pow3 = lambda x: pow(x, 3)
pow4 = lambda x: pow(x, 4)
inv = lambda x: 1/x

# Defaults
# OPS is a list of common operations in sympy.physics.mechanics
OPS = [sin, cos, tan, pow2, pow3, pow4, sqrt, inv]
# SYMS is a list of symbols that could be found in a matrix
SYMS = symbols('a, b, c, d, e, f, g, h, i, j')

def sum_or_prod(vec):
    """ Return either the sum or product of a vector """
    func = sample([sum, prod], 1)[0]
    return func(vec)

def randterms(n, max_terms, ops=OPS, syms=SYMS):
    """ Creates a list of random terms of size n. Each cell is composed of
    0 to max_terms, composed of randomly sampled functions from ops, and
    symbols from syms """
    ntermlist = [randint(1, max_terms) for i in range(n)]
    pick = lambda vec, nlist: [[vec[i] for i in 
            sample(range(len(vec)), n)] for n in nlist]
    rops = pick(ops, ntermlist)
    rsyms = pick(syms, ntermlist)
    terms = [zip(*opsym) for opsym in zip(rops, rsyms)]
    return [sum_or_prod(op(s) for (op, s) in term) for term in terms]

def randmat(m, n, max_terms, ops=OPS, syms=SYMS):
    """ Creates a random matrix of size (m,n). Each cell is composed of
    1 to max_terms, composed of randomly sampled functions from ops, and
    symbols from syms """
    return Matrix(randterms(m*n, max_terms, ops, syms)).reshape(m, n)

def randuppertriag(n, max_terms, ops=OPS, syms=SYMS):
    """ Creates a random upper triangular matrix of size (n, n). Each cell is
    composed of 1 to max_terms, composed of randomly sampled functions from
    ops, and symbols from syms """
    nupper = sum(range(n))
    terms = randterms(nupper, max_terms, ops, syms)
    t = 0
    rows = []
    for i in range(1, n):
        rows.append(zeros(1, i).row_join(Matrix(terms[t:t+n-i]).T))
        t += n-i
    rows.append(zeros(1, n))
    return Matrix(rows)

def symrandmat(n, max_terms, ops=OPS, syms=SYMS):
    """ Creates a random symmetric matrix of size (n, n). Each cell is
    composed of 1 to max_terms, composed of randomly sampled functions from
    ops, and symbols from syms """
    upper = randuppertriag(n, max_terms, ops, syms)
    lower = upper.T
    D = diag(*randterms(n, max_terms, ops, syms))
    return upper + D + lower
In [2]:
# 3x3 Random Upper Triangular Matrix with max number of terms 2
randuppertriag(3, 2)
Out[2]:
Matrix([
[0, cos(e),      tan(g)],
[0,      0, f**4*cos(j)],
[0,      0,           0]])
In [3]:
# 3x3 Random Matrix with max number of terms 2
randmat(3, 3, 2)
Out[3]:
Matrix([
[   b**4, j**2 + tan(c),           d**4],
[    1/d,        cos(f), c**2 + sqrt(d)],
[sqrt(a),        cos(g),           j**4]])
In [4]:
# 3x3 Symmetric matrix with max number of terms 2
symrandmat(3, 2)
Out[4]:
Matrix([
[          cos(a), sqrt(g) + sin(b),        tan(f)],
[sqrt(g) + sin(b),   sqrt(i)*cos(j),           1/b],
[          tan(f),              1/b, cos(f)*tan(c)]])

Solution Methods

There are 4 different methods we'll be testing:

  1. LUsolve: Solve the problem with LU decomposition
  2. LDLsolve: For symmetric matrices, solve with LDL decomposition
  3. cholesky_solve: For symmetric matrices, solve with cholesky decomposition
  4. _mat_inv_mul: Solve using LDL decomposition, and an intermediate substitution dictionary. This is what mechanics is currently using.

After doing some reading up on these, LDLsolve and cholesky_solve should be the fastest, as the \(A\) matrix is symmetric. _mat_inv_mul also uses LDL decomposition, but it has the overhead of substitution. However, for larger matrices this may yield benefits, as the LDL decomposition is performed only on a matrix of single symbols.

In [5]:
# This is what sympy.physics.mechanics is currently using.
# Copied into this file, to show what it's doing.
def _mat_inv_mul(A, B):
    """
    Computes A^-1 * B symbolically w/ substitution, where B is not
    necessarily a vector, but can be a matrix.

    """

    r1, c1 = A.shape
    r2, c2 = B.shape
    temp1 = Matrix(r1, c1, lambda i, j: Symbol('x' + str(j) + str(r1 * i)))
    temp2 = Matrix(r2, c2, lambda i, j: Symbol('y' + str(j) + str(r2 * i)))
    for i in range(len(temp1)):
        if A[i] == 0:
            temp1[i] = 0
    for i in range(len(temp2)):
        if B[i] == 0:
            temp2[i] = 0
    temp3 = []
    for i in range(c2):
        temp3.append(temp1.LDLsolve(temp2[:, i]))
    temp3 = Matrix([i.T for i in temp3]).T
    return temp3.subs(dict(list(zip(temp1, A)))).subs(dict(list(zip(temp2, B))))

Benchmark: Constant Matrix Complexity, Varying n

Here we check increasing matrix dimensions, with a constant complexity of each cell (max_terms = 3). Only running up to n=6, as computation got excessive after that.

In [6]:
from timeit import repeat

# Run the benchmark for varying values of n, with max_terms = 3
timestmt = lambda stmt, setup, n, r: min(repeat(stmt, setup=setup, repeat=r, number=n))/n
ns = range(2, 7)
lu_t = []
ldl_t = []
chol_t = []
matinv_t = []
for n in ns:
    A = symrandmat(n, 3)
    B = randmat(n, 1, 3)
    lu_t.append(timestmt('A.LUsolve(B)', 'from __main__ import A, B', 1, 3))
    ldl_t.append(timestmt('A.LDLsolve(B)', 'from __main__ import A, B', 1, 3))
    chol_t.append(timestmt('A.cholesky_solve(B)', 'from __main__ import A, B', 1, 3))
    matinv_t.append(timestmt('_mat_inv_mul(A, B)',
        'from __main__ import A, B, _mat_inv_mul', 1, 1))
In [7]:
# Plot the results
import matplotlib.pyplot as plt
% matplotlib inline

plt.figure()
plt.plot(ns, lu_t, label='A.LUsolve(B)')
plt.plot(ns, ldl_t, label='A.LDLsolve(B)')
plt.plot(ns, chol_t, label='A.cholesky_solve(B)')
plt.plot(ns, matinv_t, label='mat_inv_mul(A, B)')
plt.xlabel('Matrix Dimension')
plt.ylabel('Time (s)')
plt.title('Timings for Solution to A*x = b, Varying Dimensions')
plt.legend(loc=2)
plt.figure()
plt.plot(ns, lu_t, label='A.LUsolve(B)')
plt.plot(ns, ldl_t, label='A.LDLsolve(B)')
plt.plot(ns, chol_t, label='A.cholesky_solve(B)')
plt.xlabel('Matrix Dimension')
plt.ylabel('Time (s)')
plt.title('Timings for Solution to A*x = b,  Varying Dimensions')
plt.legend(loc=2)
plt.show()

Observing the above plots it can be seen that _mat_inv_mul is several orders of magnitude slower than the other three methods, and increases in time at a faster rate. As the underlying algorithm is the same as LDLsolve, this indicates that the subs operations are slow, and of larger complexity than LDLsolve. While the other three methods are all close, LUsolve is the fastest. This is interesting, because for purely numeric computation (i.e. floating point only), LDLsolve and cholesky_solve should be faster for symmetric matrices.

Benchmark: Varying Matrix Complexity, Constant n

Here we vary the complexity of the matrix (how many terms in each cell), but keep the matrix size constant. I picked n = 4, because it was big enough to take some time, but not too big to take forever. I'd expect to see the computation time increase with complexity, but not as rapidly as it did with matrix dimension.

In [8]:
# Run the benchmark for varying values of max_terms, with n = 4
ms = range(2, 9)
lu_t = []
ldl_t = []
chol_t = []
matinv_t = []
for m in ms:
    A = symrandmat(4, m)
    B = randmat(4, 1, m)
    lu_t.append(timestmt('A.LUsolve(B)', 'from __main__ import A, B', 1, 3))
    ldl_t.append(timestmt('A.LDLsolve(B)', 'from __main__ import A, B', 1, 3))
    chol_t.append(timestmt('A.cholesky_solve(B)', 'from __main__ import A, B', 1, 3))
    matinv_t.append(timestmt('_mat_inv_mul(A, B)',
        'from __main__ import A, B, _mat_inv_mul', 1, 1))
In [9]:
# Plot the results
plt.figure()
plt.plot(ms, lu_t, label='A.LUsolve(B)')
plt.plot(ms, ldl_t, label='A.LDLsolve(B)')
plt.plot(ms, chol_t, label='A.cholesky_solve(B)')
plt.plot(ms, matinv_t, label='mat_inv_mul(A, B)')
plt.xlabel('Matrix Complexity')
plt.ylabel('Time (s)')
plt.title('Timings for Solution to A*x = b, Varying Complexity')
plt.legend(loc=2)
plt.figure()
plt.plot(ms, lu_t, label='A.LUsolve(B)')
plt.plot(ms, ldl_t, label='A.LDLsolve(B)')
plt.plot(ms, chol_t, label='A.cholesky_solve(B)')
plt.xlabel('Matrix Complexity')
plt.ylabel('Time (s)')
plt.title('Timings for Solution to A*x = b, Varying Complexity')
plt.legend(loc=2)
plt.show()

Interestingly, it seems that LUsolve, LDLsolve, and cholesky_solve don't increase in time with complexity, while _mat_inv_mul does. This must be due to the intermediate substitution step.

Benchmark: Turning Matrix Term Complexity up to 11

Using the matrices composed by the above methods, all cells are either a sum of terms, or a product. This is close to how the matrices we need to solve look, but it's not a perfect representation. Let's try with a really ugly matrix. We can use the fact that for any matrix \(A\), \(A A^T\) is a symmetric matrix.

In [10]:
temp = prod([randmat(4, 4, 2) for i in range(2)])
A = temp*temp.T
B = randmat(4,1,6)
In [11]:
A
Out[11]:
Matrix([
[                                                                                                                                                                                                                                                                                                   (h**(9/2) + (sqrt(a) + cos(i))*(c**4 + tan(f)) + sin(i)*cos(c) + sin(f)*sin(g)/c)**2 + (c*sin(g) + (sqrt(a) + cos(i))*tan(g) + sin(i)*cos(c) + sqrt(h)*tan(a)/e)**2 + (sqrt(h)*(b**4 + sqrt(h)) + sqrt(i)*(sqrt(a) + cos(i)) + sin(a)*cos(c) + i**2*sin(g)*sin(j)/c)**2 + (b**2*e**4*cos(c) + d**3*(sqrt(a) + cos(i)) + g**2*sqrt(h) + (g**4 + cos(f))*sin(g)/c)**2, (h**(9/2) + (sqrt(a) + cos(i))*(c**4 + tan(f)) + sin(i)*cos(c) + sin(f)*sin(g)/c)*(d**3*sin(i) + j**3*sin(f) + h**4/f + h**2*(c**4 + tan(f))/f) + (c*sin(g) + (sqrt(a) + cos(i))*tan(g) + sin(i)*cos(c) + sqrt(h)*tan(a)/e)*(c**2*j**3 + d**3*sin(i) + h**2*tan(g)/f + tan(a)/(e*f)) + (d**3*sin(a) + i**2*j**3*sin(j) + h**2*sqrt(i)/f + (b**4 + sqrt(h))/f)*(sqrt(h)*(b**4 + sqrt(h)) + sqrt(i)*(sqrt(a) + cos(i)) + sin(a)*cos(c) + i**2*sin(g)*sin(j)/c) + (b**2*d**3*e**4 + d**3*h**2/f + j**3*(g**4 + cos(f)) + g**2/f)*(b**2*e**4*cos(c) + d**3*(sqrt(a) + cos(i)) + g**2*sqrt(h) + (g**4 + cos(f))*sin(g)/c), (h**(9/2) + (sqrt(a) + cos(i))*(c**4 + tan(f)) + sin(i)*cos(c) + sin(f)*sin(g)/c)*(c**2*sin(f)/h + h**4*(sin(a) + cos(b)) + (c**4 + tan(f))/c + sin(i)/c) + (c*sin(g) + (sqrt(a) + cos(i))*tan(g) + sin(i)*cos(c) + sqrt(h)*tan(a)/e)*(c**4/h + (sin(a) + cos(b))*tan(a)/e + sin(i)/c + tan(g)/c) + (sqrt(h)*(b**4 + sqrt(h)) + sqrt(i)*(sqrt(a) + cos(i)) + sin(a)*cos(c) + i**2*sin(g)*sin(j)/c)*(c**2*i**2*sin(j)/h + (b**4 + sqrt(h))*(sin(a) + cos(b)) + sqrt(i)/c + sin(a)/c) + (b**2*e**4/c + c**2*(g**4 + cos(f))/h + g**2*(sin(a) + cos(b)) + d**3/c)*(b**2*e**4*cos(c) + d**3*(sqrt(a) + cos(i)) + g**2*sqrt(h) + (g**4 + cos(f))*sin(g)/c), (h**(9/2) + (sqrt(a) + cos(i))*(c**4 + tan(f)) + sin(i)*cos(c) + sin(f)*sin(g)/c)*(sqrt(b)*g**2*sin(f) + sqrt(e)*sin(i) + f**2*h**4 + (c**4 + tan(f))*tan(i)) + (c*sin(g) + (sqrt(a) + cos(i))*tan(g) + sin(i)*cos(c) + sqrt(h)*tan(a)/e)*(sqrt(b)*c**2*g**2 + sqrt(e)*sin(i) + tan(g)*tan(i) + f**2*tan(a)/e) + (sqrt(h)*(b**4 + sqrt(h)) + sqrt(i)*(sqrt(a) + cos(i)) + sin(a)*cos(c) + i**2*sin(g)*sin(j)/c)*(sqrt(b)*g**2*i**2*sin(j) + sqrt(e)*sin(a) + f**2*(b**4 + sqrt(h)) + sqrt(i)*tan(i)) + (sqrt(b)*g**2*(g**4 + cos(f)) + b**2*e**(9/2) + d**3*tan(i) + f**2*g**2)*(b**2*e**4*cos(c) + d**3*(sqrt(a) + cos(i)) + g**2*sqrt(h) + (g**4 + cos(f))*sin(g)/c)],
[                                                  (h**(9/2) + (sqrt(a) + cos(i))*(c**4 + tan(f)) + sin(i)*cos(c) + sin(f)*sin(g)/c)*(d**3*sin(i) + j**3*sin(f) + h**4/f + h**2*(c**4 + tan(f))/f) + (c*sin(g) + (sqrt(a) + cos(i))*tan(g) + sin(i)*cos(c) + sqrt(h)*tan(a)/e)*(c**2*j**3 + d**3*sin(i) + h**2*tan(g)/f + tan(a)/(e*f)) + (d**3*sin(a) + i**2*j**3*sin(j) + h**2*sqrt(i)/f + (b**4 + sqrt(h))/f)*(sqrt(h)*(b**4 + sqrt(h)) + sqrt(i)*(sqrt(a) + cos(i)) + sin(a)*cos(c) + i**2*sin(g)*sin(j)/c) + (b**2*d**3*e**4 + d**3*h**2/f + j**3*(g**4 + cos(f)) + g**2/f)*(b**2*e**4*cos(c) + d**3*(sqrt(a) + cos(i)) + g**2*sqrt(h) + (g**4 + cos(f))*sin(g)/c),                                                                                                                                                                                                                                                                                                                                       (c**2*j**3 + d**3*sin(i) + h**2*tan(g)/f + tan(a)/(e*f))**2 + (d**3*sin(a) + i**2*j**3*sin(j) + h**2*sqrt(i)/f + (b**4 + sqrt(h))/f)**2 + (d**3*sin(i) + j**3*sin(f) + h**4/f + h**2*(c**4 + tan(f))/f)**2 + (b**2*d**3*e**4 + d**3*h**2/f + j**3*(g**4 + cos(f)) + g**2/f)**2,                                                                                      (c**2*j**3 + d**3*sin(i) + h**2*tan(g)/f + tan(a)/(e*f))*(c**4/h + (sin(a) + cos(b))*tan(a)/e + sin(i)/c + tan(g)/c) + (d**3*sin(a) + i**2*j**3*sin(j) + h**2*sqrt(i)/f + (b**4 + sqrt(h))/f)*(c**2*i**2*sin(j)/h + (b**4 + sqrt(h))*(sin(a) + cos(b)) + sqrt(i)/c + sin(a)/c) + (d**3*sin(i) + j**3*sin(f) + h**4/f + h**2*(c**4 + tan(f))/f)*(c**2*sin(f)/h + h**4*(sin(a) + cos(b)) + (c**4 + tan(f))/c + sin(i)/c) + (b**2*e**4/c + c**2*(g**4 + cos(f))/h + g**2*(sin(a) + cos(b)) + d**3/c)*(b**2*d**3*e**4 + d**3*h**2/f + j**3*(g**4 + cos(f)) + g**2/f),                                                                                      (c**2*j**3 + d**3*sin(i) + h**2*tan(g)/f + tan(a)/(e*f))*(sqrt(b)*c**2*g**2 + sqrt(e)*sin(i) + tan(g)*tan(i) + f**2*tan(a)/e) + (d**3*sin(a) + i**2*j**3*sin(j) + h**2*sqrt(i)/f + (b**4 + sqrt(h))/f)*(sqrt(b)*g**2*i**2*sin(j) + sqrt(e)*sin(a) + f**2*(b**4 + sqrt(h)) + sqrt(i)*tan(i)) + (d**3*sin(i) + j**3*sin(f) + h**4/f + h**2*(c**4 + tan(f))/f)*(sqrt(b)*g**2*sin(f) + sqrt(e)*sin(i) + f**2*h**4 + (c**4 + tan(f))*tan(i)) + (sqrt(b)*g**2*(g**4 + cos(f)) + b**2*e**(9/2) + d**3*tan(i) + f**2*g**2)*(b**2*d**3*e**4 + d**3*h**2/f + j**3*(g**4 + cos(f)) + g**2/f)],
[                 (h**(9/2) + (sqrt(a) + cos(i))*(c**4 + tan(f)) + sin(i)*cos(c) + sin(f)*sin(g)/c)*(c**2*sin(f)/h + h**4*(sin(a) + cos(b)) + (c**4 + tan(f))/c + sin(i)/c) + (c*sin(g) + (sqrt(a) + cos(i))*tan(g) + sin(i)*cos(c) + sqrt(h)*tan(a)/e)*(c**4/h + (sin(a) + cos(b))*tan(a)/e + sin(i)/c + tan(g)/c) + (sqrt(h)*(b**4 + sqrt(h)) + sqrt(i)*(sqrt(a) + cos(i)) + sin(a)*cos(c) + i**2*sin(g)*sin(j)/c)*(c**2*i**2*sin(j)/h + (b**4 + sqrt(h))*(sin(a) + cos(b)) + sqrt(i)/c + sin(a)/c) + (b**2*e**4/c + c**2*(g**4 + cos(f))/h + g**2*(sin(a) + cos(b)) + d**3/c)*(b**2*e**4*cos(c) + d**3*(sqrt(a) + cos(i)) + g**2*sqrt(h) + (g**4 + cos(f))*sin(g)/c),                                                     (c**2*j**3 + d**3*sin(i) + h**2*tan(g)/f + tan(a)/(e*f))*(c**4/h + (sin(a) + cos(b))*tan(a)/e + sin(i)/c + tan(g)/c) + (d**3*sin(a) + i**2*j**3*sin(j) + h**2*sqrt(i)/f + (b**4 + sqrt(h))/f)*(c**2*i**2*sin(j)/h + (b**4 + sqrt(h))*(sin(a) + cos(b)) + sqrt(i)/c + sin(a)/c) + (d**3*sin(i) + j**3*sin(f) + h**4/f + h**2*(c**4 + tan(f))/f)*(c**2*sin(f)/h + h**4*(sin(a) + cos(b)) + (c**4 + tan(f))/c + sin(i)/c) + (b**2*e**4/c + c**2*(g**4 + cos(f))/h + g**2*(sin(a) + cos(b)) + d**3/c)*(b**2*d**3*e**4 + d**3*h**2/f + j**3*(g**4 + cos(f)) + g**2/f),                                                                                                                                                                                                                                                                                                                                       (c**4/h + (sin(a) + cos(b))*tan(a)/e + sin(i)/c + tan(g)/c)**2 + (b**2*e**4/c + c**2*(g**4 + cos(f))/h + g**2*(sin(a) + cos(b)) + d**3/c)**2 + (c**2*sin(f)/h + h**4*(sin(a) + cos(b)) + (c**4 + tan(f))/c + sin(i)/c)**2 + (c**2*i**2*sin(j)/h + (b**4 + sqrt(h))*(sin(a) + cos(b)) + sqrt(i)/c + sin(a)/c)**2,                                                     (c**4/h + (sin(a) + cos(b))*tan(a)/e + sin(i)/c + tan(g)/c)*(sqrt(b)*c**2*g**2 + sqrt(e)*sin(i) + tan(g)*tan(i) + f**2*tan(a)/e) + (sqrt(b)*g**2*(g**4 + cos(f)) + b**2*e**(9/2) + d**3*tan(i) + f**2*g**2)*(b**2*e**4/c + c**2*(g**4 + cos(f))/h + g**2*(sin(a) + cos(b)) + d**3/c) + (sqrt(b)*g**2*sin(f) + sqrt(e)*sin(i) + f**2*h**4 + (c**4 + tan(f))*tan(i))*(c**2*sin(f)/h + h**4*(sin(a) + cos(b)) + (c**4 + tan(f))/c + sin(i)/c) + (sqrt(b)*g**2*i**2*sin(j) + sqrt(e)*sin(a) + f**2*(b**4 + sqrt(h)) + sqrt(i)*tan(i))*(c**2*i**2*sin(j)/h + (b**4 + sqrt(h))*(sin(a) + cos(b)) + sqrt(i)/c + sin(a)/c)],
[(h**(9/2) + (sqrt(a) + cos(i))*(c**4 + tan(f)) + sin(i)*cos(c) + sin(f)*sin(g)/c)*(sqrt(b)*g**2*sin(f) + sqrt(e)*sin(i) + f**2*h**4 + (c**4 + tan(f))*tan(i)) + (c*sin(g) + (sqrt(a) + cos(i))*tan(g) + sin(i)*cos(c) + sqrt(h)*tan(a)/e)*(sqrt(b)*c**2*g**2 + sqrt(e)*sin(i) + tan(g)*tan(i) + f**2*tan(a)/e) + (sqrt(h)*(b**4 + sqrt(h)) + sqrt(i)*(sqrt(a) + cos(i)) + sin(a)*cos(c) + i**2*sin(g)*sin(j)/c)*(sqrt(b)*g**2*i**2*sin(j) + sqrt(e)*sin(a) + f**2*(b**4 + sqrt(h)) + sqrt(i)*tan(i)) + (sqrt(b)*g**2*(g**4 + cos(f)) + b**2*e**(9/2) + d**3*tan(i) + f**2*g**2)*(b**2*e**4*cos(c) + d**3*(sqrt(a) + cos(i)) + g**2*sqrt(h) + (g**4 + cos(f))*sin(g)/c),                                    (c**2*j**3 + d**3*sin(i) + h**2*tan(g)/f + tan(a)/(e*f))*(sqrt(b)*c**2*g**2 + sqrt(e)*sin(i) + tan(g)*tan(i) + f**2*tan(a)/e) + (d**3*sin(a) + i**2*j**3*sin(j) + h**2*sqrt(i)/f + (b**4 + sqrt(h))/f)*(sqrt(b)*g**2*i**2*sin(j) + sqrt(e)*sin(a) + f**2*(b**4 + sqrt(h)) + sqrt(i)*tan(i)) + (d**3*sin(i) + j**3*sin(f) + h**4/f + h**2*(c**4 + tan(f))/f)*(sqrt(b)*g**2*sin(f) + sqrt(e)*sin(i) + f**2*h**4 + (c**4 + tan(f))*tan(i)) + (sqrt(b)*g**2*(g**4 + cos(f)) + b**2*e**(9/2) + d**3*tan(i) + f**2*g**2)*(b**2*d**3*e**4 + d**3*h**2/f + j**3*(g**4 + cos(f)) + g**2/f),                                    (c**4/h + (sin(a) + cos(b))*tan(a)/e + sin(i)/c + tan(g)/c)*(sqrt(b)*c**2*g**2 + sqrt(e)*sin(i) + tan(g)*tan(i) + f**2*tan(a)/e) + (sqrt(b)*g**2*(g**4 + cos(f)) + b**2*e**(9/2) + d**3*tan(i) + f**2*g**2)*(b**2*e**4/c + c**2*(g**4 + cos(f))/h + g**2*(sin(a) + cos(b)) + d**3/c) + (sqrt(b)*g**2*sin(f) + sqrt(e)*sin(i) + f**2*h**4 + (c**4 + tan(f))*tan(i))*(c**2*sin(f)/h + h**4*(sin(a) + cos(b)) + (c**4 + tan(f))/c + sin(i)/c) + (sqrt(b)*g**2*i**2*sin(j) + sqrt(e)*sin(a) + f**2*(b**4 + sqrt(h)) + sqrt(i)*tan(i))*(c**2*i**2*sin(j)/h + (b**4 + sqrt(h))*(sin(a) + cos(b)) + sqrt(i)/c + sin(a)/c),                                                                                                                                                                                                                                                                                                                                       (sqrt(b)*c**2*g**2 + sqrt(e)*sin(i) + tan(g)*tan(i) + f**2*tan(a)/e)**2 + (sqrt(b)*g**2*(g**4 + cos(f)) + b**2*e**(9/2) + d**3*tan(i) + f**2*g**2)**2 + (sqrt(b)*g**2*sin(f) + sqrt(e)*sin(i) + f**2*h**4 + (c**4 + tan(f))*tan(i))**2 + (sqrt(b)*g**2*i**2*sin(j) + sqrt(e)*sin(a) + f**2*(b**4 + sqrt(h)) + sqrt(i)*tan(i))**2]])
In [12]:
B
Out[12]:
Matrix([
[                 a**3 + d**4 + g**2],
[a**4 + f**2 + sin(g) + tan(e) + 1/d],
[                         f**2 + 1/g],
[                             sin(d)]])

Yikes, that's some ugly stuff. Let's see how the various methods fair against that:

In [13]:
# Run the benchmarks
print("Time for LUsolve: ", timestmt('A.LUsolve(B)', 'from __main__ import A, B', 1, 3))
print("Time for LDLsolve: ", timestmt('A.LDLsolve(B)', 'from __main__ import A, B', 1, 3))
print("Time for cholesky_solve: ", timestmt('A.cholesky_solve(B)', 'from __main__ import A, B', 1, 3))
print("Time for _mat_inv_mul: ", timestmt('_mat_inv_mul(A, B)', 'from __main__ import A, B, _mat_inv_mul', 1, 1))
Time for LUsolve:  0.0034651060013857204
Time for LDLsolve:  0.009068403000128455
Time for cholesky_solve:  0.006806592999055283
Time for _mat_inv_mul:  144.967967898001

From this it can be seen that _mat_inv_mul is far slower than the other three methods, even for extremely complex matrices. There doesn't seem to be any benefit to the intermediate subs step.

Benchmark: Expression Size

This isn't a time benchmark. Rather, it's a measurment of the resulting expression's readability. If LUsolve is faster, but results in a less compact expression than LDLsolve, then it may be more beneficial to go with the latter. To test this we'll measure the number of operations using the count_ops method. This is a rough metric of the compactness of the solution.

In [14]:
from sympy import count_ops
# Run the benchmark for varying values of n, with max_terms = 2
ns = range(2, 6)
lu_nops = []
ldl_nops = []
chol_nops = []
matinv_nops = []
for n in ns:
    A = symrandmat(n, 2)
    B = randmat(n, 1, 2)
    # Solve the system
    sol_LU = A.LUsolve(B)
    sol_LDL = A.LDLsolve(B)
    sol_chol = A.cholesky_solve(B)
    sol_matinv = _mat_inv_mul(A, B)
    # Count the number of ops
    lu_nops.append(count_ops(sol_LU))
    ldl_nops.append(count_ops(sol_LDL))
    chol_nops.append(count_ops(sol_chol))
    matinv_nops.append(count_ops(sol_matinv))
In [15]:
# Plot the results
plt.figure()
plt.plot(ns, lu_nops, label='A.LUsolve(B)')
plt.plot(ns, ldl_nops, label='A.LDLsolve(B)')
plt.plot(ns, chol_nops, label='A.cholesky_solve(B)')
plt.plot(ns, matinv_nops, label='mat_inv_mul(A, B)')
plt.xlabel('Matrix Dimensions')
plt.ylabel('Number of Operations')
plt.title('Expression Size for Solution to A*x = b')
plt.legend(loc=2)
plt.figure()
plt.plot(ns, lu_nops, label='A.LUsolve(B)')
plt.plot(ns, ldl_nops, label='A.LDLsolve(B)')
plt.plot(ns, chol_nops, label='A.cholesky_solve(B)')
plt.xlabel('Matrix Dimensions')
plt.ylabel('Number of Operations')
plt.title('Expression Size for Solution to A*x = b')
plt.legend(loc=2)
plt.show()

Observing the above plots, once again _mat_inv_mul comes out as the worst (by far). The remaining three barely differ, but LUsolve results in the most compact expression.

Discussion of Results

Overall, _mat_inv_mul comes out the worst for every benchmark, and LUsolve comes out the best. Unless I did something wrong with my benchmarks LUsolve should replace every call to _mat_inv_mul in mechanics. It results in a more compact form, and has several order of magnitude faster running speed.

This surprised me. I would have thought that a more complicated expression would take longer to solve, but as seen in the second benchmark matrix complexity had no effect on running speed for the three solution algorithms (although it did affect the substitution speed for _mat_inv_mul). I suppose that's why you're always told to profile before you optimize. Often your intuition is wrong.

Something else that surprised me was that LDLsolve and cholesky_solve had a slower running time than LUsolve. For numerical symmetric matrices, this shouldn't be the case. I wasn't able to find anything about symbolic calculation of these decompositions, but I assume it should be about the same. Either way, the more general LU decomposition seems to be the fastest.


Did I do something wrong? Disagree with these benchmarks? Let me know in the comments below!

by Jim Crist at June 20, 2014 07:00 PM

June 19, 2014

Fabian Pedregosa

Surrogate Loss Functions in Machine Learning

TL; DR These are some notes on calibration of surrogate loss functions in the context of machine learning. But mostly it is an excuse to post some images I made.

In the binary-class classification setting we are given $n$ training samples $\{(X_1, Y_1), \ldots, (X_n, Y_n)\}$, where $X_i$ belongs to some sample space $\mathcal{X}$, usually $\mathbb{R}^p$ but for the purpose of this post we can keep i abstract, and $y_i \in \{-1, 1\}$ is an integer representing the class label.

We are also given a loss function $\ell: \{-1, 1\} \times \{-1, 1\} \to \mathbb{R}$ that measures the error of a given prediction. The value of the loss function $\ell$ at an arbitrary point $(y, \hat{y})$ is interpreted as the cost incurred by predicting $\hat{y}$ when the true label is $y$. In classification this function is often the zero-one loss, that is, $\ell(y, \hat{y})$ is zero when $y = \hat{y}$ and one otherwise.

The goal is to find a function $h: \mathcal{X} \to [k]$, the classifier, with the smallest expected loss on a new sample. In other words, we seek to find a function $h$ that minimizes the expected $\ell$-risk, given by $$ \mathcal{R}_{\ell}(h) = \mathbb{E}_{XY}[\ell(Y, h(X))] $$

In theory, we could directly minimize the $\ell$-risk and we would have the optimal classifier, also known as Bayes predictor. However, there are several problems associated with this approach. One is that the probability distribution of $XY$ is unknown, thus computing the exact expected value is not feasible. It must be approximated by the empirical risk. Another issue is that this quantity is difficult to optimize because the function $\ell$ is discontinuous. Take for example a problem in which $\mathcal{X} = \mathbb{R}^2, k=2$, and we seek to find the linear function $f(X) = \text{sign}(X w), w \in \mathbb{R}^2$ and that minimizes the $\ell$-risk. As a function of the parameter $w$ this function looks something like

loss as function of w

This function is discontinuous with large, flat regions and is thus extremely hard to optimize using gradient-based methods. For this reason it is usual to consider a proxy to the loss called a surrogate loss function. For computational reasons this is usually convex function $\Psi: \mathbb{R} \to \mathbb{R}_+$. An example of such surrogate loss functions is the hinge loss, $\Psi(t) = \max(1-t, 0)$, which is the loss used by Support Vector Machines (SVMs). Another example is the logistic loss, $\Psi(t) = 1/(1 + \exp(-t))$, used by the logistic regression model. If we consider the logistic loss, minimizing the $\Psi$-risk, given by $\mathbb{E}_{XY}[\Psi(Y, f(X))]$, of the function $f(X) = X w$ becomes a much more more tractable optimization problem:

In short, we have replaced the $\ell$-risk which is computationally difficult to optimize with the $\Psi$-risk which has more advantageous properties. A natural questions to ask is how much have we lost by this change. The property of whether minimizing the $\Psi$-risk leads to a function that also minimizes the $\ell$-risk is often referred to as consistency or calibration. For a more formal definition see [1] and [2]. This property will depend on the surrogate function $\Psi$: for some functions $\Psi$ it will be verified the consistency property and for some not. One of the most useful characterizations was given in [1] and states that if $\Psi$ is convex then it is consistent if and only if it is differentiable at zero and $\Psi'(0) < 0$. This includes most of the commonly used surrogate loss functions, including hinge, logistic regression and Huber loss functions.


  1. P. L. Bartlett, M. I. Jordan, and J. D. McAuliffe, “Convexity , Classification , and Risk Bounds,” J. Am. Stat. Assoc., pp. 1–36, 2003. 

  2. A. Tewari and P. L. Bartlett, “On the Consistency of Multiclass Classification Methods,” J. Mach. Learn. Res., vol. 8, pp. 1007–1025, 2007. 

by Fabian Pedregosa at June 19, 2014 10:00 PM

Avichal Dayal

In this week I worked on Formal Power Series and getting my two PRs merged.

While getting my PRs related to asymptotic expansion merged, I learned about how important documentation is and the little details that matter. I made some very silly mistakes causing delay in getting the PRs merged.

Formal Power Series:
I made some good progress regarding FPS this week. As I said in the earlier post, I decided to go with lazy recursion to implement infinite power series.
In this a series is represented as a tuple of head and tail where head is the first term of the series while tail is rest of the infinite series.
Algorithms to do addition, subtraction, multiplication are very simple if we use recursion.
E.g.:-
series1 = (head1, tail1)
series2 = (head2, tail2)
Then,
- series1 + series2 = (head1+head2, tail1+tail2)
- series1 * series2 = (head1*head2, head1*tail2 + head2*tail1 + tail1*tail2)

Similarly other operations like division, inversion, composition can be written.

Along with this, I also finished implementing simpleDE, DEtoRE and somewhat solveRE methods. These methods are used to find the generator for a given function in x.
- simpleDE() converts function into a simple differential equation. This involved finding the number of rationally independent terms among a set of expressions. I coded a very inefficient solution with repeatedly dividing each term by other and checking if there quotient is a rational function. This is O(n**2) approach. There might be a better approach.
- DEtoRE() converts differential equation to a recurrence equation. This is a straightforward algorithm with simple substitutions.

Next week:
Since the operations on infinite series are a recursive process, I guess there is a problem of memory management with my code. Operations like division, multiplication are highly recursive creating Stream class numerous times. Upon giving high value of "n" while printing terms, it gives a "maximum recursion depth" error.
I will work on this and add more operations to Stream class.

by Avichal Dayal (noreply@blogger.com) at June 19, 2014 04:40 PM

June 18, 2014

Harsh Gupta

week 4

Intersections of Infinitely indexed sets

This week I implemented a method to do intersection of imagesets using the solutions of Diophantine equations at PR.

Say you have to find the intersection of sets 2*n| n in Integers and 3*m| m in Integers. The intersection of these sets is the set of the common values in the two sets, which in this case is equivalent to the values of n for which the equation 2*n - 3*m has some integral solution in m. Or the values of m for which the 2*n - 3*m has some integral solution in n. Diophantine equations are equations for which only integral solutions are searched for. The Diophantine module was written by Thilina as his GSoC project last year. It gives the parametric solution for such equation.

In [17]: diophantine(2*n - 3*m)
Out[17]: {(-2*t, -3*t)}

The Solution is sorted according to alphabetic order of the variables involved. So the value of LHS (2*n) for which the equation is 2*(-3*t) that is -6*t and it is the intersection of the sets described above -6*t| t in Integers. Since -6*t| t in Integers is same as 6*t| t in Integers I also wrote some simplification rules for the imagesets with Integers as baseset.

Sets for Invert Function

The sets module turned out to be better than I expected. I had a perception that substitutions doesn't work properly with sets and I have even opened an issue for that but it turned out I hadn't looked closely enough. It worked well for the free variables and it didn't worked for the things it shouldn't work i.e., the bound variables in the imagesets.

Using sets simplified the code. All the list comprehensions like this [i.subs(symbol, symbol/g) for i in _invert(h, symbol)] were converted to simple substitutions for sets and other sets operations. _invert(h, symbol).subs(symbol, symbol/g)

Just by changing the output of invert to sets, then by adding the inverse of trigonometric function and writing the code to rewrite then as tan I was able to return all the solutions of the equations like cos(x) + sin(x) == 0 it turned to out to easier than I thought. Using sets as output makes thinking about the mathematics of the solvers much more easier and the code comes out to be pretty natural. Now when we can see the results I can surely say there can be no better output for solvers than sets.

This week I'll study LambertW function and then code additional techniques to solve real equations. I'll also try to figure out techniques to perform Union on infinitely indexed sets.

by Harsh Gupta at June 18, 2014 05:59 AM

June 17, 2014

Soumya Dipta Biswas

GSoC 2014: Week 4

Hello Everyone,

The majority of work this week was finishing up my previous work and creating the architecture for the First Order Logic module. Consequently, this post is concerned with First Order Logic (FoL) and is intended mostly for those who are not familiar with the concepts. If you find yourself comfortable with writing formulas in FoL then you can safely skip this post.

Before we go on to the discussion about First Order Logic, one very interesting topic that came up this week was the idea of an Atom. What is an atom for propositional logic? The classic answer amounts to anything that cannot be broken down further to yield simpler terms. Analogously, let us define a literal to mean any atom or its negation. So, A is a literal (and so is ~A). How about A > B. Without bordering on any ideas of higher order logic, let us tweak the definition of atom to mean anything that cannot be broken down using propositional logic to yield simpler terms. So, how does this definition change things? It means we can treat any expression that propositional logic can’t break down as a black box with the assumption that it will yield a boolean value once evaluated. Analogously, a literal is now any expression which cannot be simplified (in propositional logic) or its negation. This definition gives a lot more freedom and power to the idea of formulas. Ofcourse, one needs to keep in mind that it is compulsory for these literals to yield a boolean once evaluated (typically before propositional formula evaluation methods are called on the formula).

Ok, with that, we finish the first part our Propositional Logic post series. The future posts will be associated with First Order Logic before coming back to SAT (in propositional formula). For the remainder of the series the words First Order Logic and Predicate Logic will be used interchangeably (while there is some difference between the two, we will ignore them. Finally beginning the actual discussion about FoL, let us start with some definitions.

  • Constant: Any object of the universe whose value remains the same throughout the expression. e.g. john, 1, computer, socrates
  • Variable: An object which can take up any value from the domain. e.g. X, Y
  • Function: A mapping of n-ary terms (see below) to a term. e.g. F(X, G(a, Y)), teacher(teacher(aristotle))
  • Term: Any Constant, Variable or Function is called a Term.
  • Predicate: A mapping of n-ary terms to a boolean value. e.g. IsTeacher(Socrates, plato)
  • Quantifier: Quantifiers assert a particular property with respect to a bound variable. Q(X, IsTeacher(X, alexander)) Here Q is some quantifier (see below) and X is the bound variable.
  • ForAll (∀): Also known as Universal Quantifier, this asserts that for each value of the bound variable the expression holds True
  • Exists (∃): Also known as Existential Quantifier, this asserts that for atleast one of the values of the bound variable the expression holds True

The above description is far from complete. Below are some more examples. For these examples assume that the program is intelligent enough to know the answer.

# All constants are in lowercase
# All variables are in uppercase
# All functions starts with a lowercase but subsequent words start with an uppercase
# All predicates follow camel casing

# teacher(person) is a function which maps a person to his teacher
>>> teacher(alexander)
aristotle
>>> teacher(teacher(aristotle))
'Socrates'

# f(X) = 2X + 3y + 1
>>> f(X, Y)
Add(Mul(2, X), Mul(3, Y), 1)

# IsTeacher(a, b) is a predicate which returns if a is the teacher of b
>>> IsTeacher(alexander, plato)
False

# Assume that the following properties hold for the relationship teacher
# If A is the teacher of B, and B is the teacher of C, then A is the teacher of C (Transitivity)
# A is always the teacher of A (Identity)
>>> domain = [socrates, plato, aristotle, alexander]
>>> ForAll(X, IsTeacher(X, alexander))
True
>>> Exists(X, IsWarrior(X))
True
>>> ForAll(X, IsWarrior(X))
False
>>> Exists(X, InventedCalculus(X))
False

# f(x, y) = x + 2y
# g(x, y) = 2x + y
>>> ForAll(X, ForAll(Y, Greater(f(X, Y), g(X, Y)))
False
>>> Exists(X, Exists(Y, Greater(f(X, Y), g(X, Y)))
True

I hope the above examples make the definitions clearer. However SymPy is a symbolic library (rather than a numeric one). The implication of the same is that we are much less concerned about the value of the expression than the expression itself. For example, we will be predominantly dealing with functions, variables and predicates without being bothered by their (real life) values. What follows next are examples of conversion of natural language sentences to First Order Logic and their representation in SymPy. One would notice that we have only dealt with the new FoL concepts in the examples above and completely ignored the traditional propositional operators. Using the idea of Atom described at the beginning of this post, if we think of every Predicate (and Quantifier) as an Atom, then we can use every idea of propositional logic here. The examples below describe the complete idea.

# All dolphins all mammals
>>> X = symbols('X')
>>> Dolphin = Predicate('Dolphin')
>>> Mammal = Predicate('Mammal')
>>> ForAll(X, Dolphin(X) >> Mammal(X))
ForAll((X), Implies(Dolphin(X), Mammal(X)))


>>> X, T = symbols('X T')
>>> Person = Predicate('Person')
>>> Time = Predicate('Time')
>>> CanFool = Predicate('CanFool')

# You can fool some of the people all of the time
>>> Exists(X, ForAll(T, (Person(X) & Time(T)) >> CanFool(X, T)))
Exists((X), ForAll(T, Implies(And(Person(X), Time(T)), CanFool(X, T))))

# You can fool all of the people some of the time
>>> ForAll(X, Exists(T, (Person(X) & Time(T)) >> CanFool(X, T)))
ForAll((X), Exists(T, Implies(And(Person(X), Time(T)), CanFool(X, T))))

# You cannot fool all of the people all of the time
>>> ~ForAll(X, ForAll(T, (Person(X) & Time(T)) >> CanFool(X, T)))
Not(ForAll((X,T), Implies(And(Person(X), Time(T)), CanFool(X, T))))

I hope this should provide you with some basic idea about First Order Logic. I will be updating this post with ideas related to Interpretation/Model. Stay tuned.

Aloha!!!

by SD at June 17, 2014 08:54 PM

Akshay Narasimha

Gsoc Week-4

Week 4 of Gsoc has come to an end this week and I have decided the class structure of Plane which I am going to implement in the following week.
Here is the link to the Plane class I have implemented so far.

A plane can be represented with the help of 3 points or with a point and a normal vector to the desired plane. Regarding which implementation to follow , I had a long discussion with Stefan and decided to go with the latter as 3 points will be redundant in the algorithms I am planning to implement.

Implementation:
>>> a = Plane(Point3D(1, 1, 1), Point3D(2, 3, 4), Point3D(2, 2, 2)
>>>a
Plane(Point3D(1,  1, 1), [-1, 2, -1])
>>> b =Plane(Point3D(1, 1, 1), normal_vector=[1, 3, 5])
>>> b
Plane(Point3D(1, 1, 1), [1, 3, 5])
>>> b.equation()
 x + 3*y + 5*z - 9

I plan to continue the work on the Plane class this week an will add more methods to it.

Until then cheers!

by Akshay Narasimha (noreply@blogger.com) at June 17, 2014 01:47 PM

Tarun Gaba

GSoC 14: Fourth Week!

{% include JB/setup %} [ <-Back to posts ](/gsoc14) Fourth week of GSoC'14 has ended. This week has been no less than an adventure for me. A lot of stuff happened in this week, which I will try to elaborate in this post. Meanwhile, something to be proud about: Screenshot ###Accomplishments: - **PyDy Visualizer UI**: I have been working on a refined UI for the visualizer. The critical components of the UI are nearly complete, and are working inside the IPython notebook. It is a great accomplishment, as I was initially skeptical about binding buttons and events _inside_ the notebook's output cell, all thanks to the new features in IPython2.0. Also apart from the UI, the THREE animations are also working inside the output cell. ###Objectives: For the upcoming week, I am supposed to make the UI elements as well as backend functionality for changing the visualization objects via GUI itself. The objects changed should be passed to python side using IPython notebook's Javascript API. This implementation will mean a fully functional visualizer in place, which can be plugged into the PyDy package for usage. ###Issues: There are two main issues which need to be addressed right now: - The animations are buggy, and are not behaving the way they are expected to behave. Probably there is some problem with the object rendering code. - The Javascript Libraries(Prototype.js and jQuery) are conflicting with each other, and with those included in native IPython notebook. This is a major issue which needs to be addressed, as it leads to unexpected behavior when the visualizations are run inside the notebook, leading to notebook freeze. Just a week left before the midterm evaluations start, and I have to get a lot before the deliverables are met. I really hope that the pending issues are addressed before the evaluations, so that I can merge a stable visualizer with the PyDy repository. [ <-Back to posts ](/gsoc14)

June 17, 2014 10:00 AM

June 16, 2014

Sudhanshu Mishra

GSoC'14 progress, week 4: finished medium, working on Fermat's principle and Snell's law

It was a great week! My last to last PR, on enhancing TWave, got merged! I have also finishedMedium and waiting for Sean's comments on it. Here's the link to that PR:
I still have a todo in the tests:
m5 = Medium('m5', permittivity=710*10**(-12)*s**4*A**2/(m**3*kg), n=1.33)
assert simplify(m5.intrinsic_impedance - 6.24845417765552*kg*m**2/(A**2*s**3)) == 0
# XXX: This is supposed to be zero but it turns out to be
# very close to zero and not zero.
# assert simplify(m5.speed - 225407863.157895*m/s) == 0
This m5.speed is equal to 225407863.157895*m/s (or maybe very close to it) but, assertion results in a failure.
Similarly in the constructor of the Medium inconsistency comes due to slight difference in floating point parameters(permittivity and permeability). I've commented out this part of my code. Maybe Sean will suggest a better way to do it.
# XXX: There's issue with precision. Values may be
# different slightly.
if permittivity != u0 and permittivity != e0:
if n != c*sqrt(permittivity*permeability):
raise ValueError("Values are not consistent.")
As soon as it gets merged, I'll make necessary changes in Fermat's principle and Snell's law part of the code and send a PR to review.
Now I realize that the next two weeks are not going to be a joy ride for me. I must finish major part of geometrical optics as written in my proposal.

That's all for now. Cheers!

by Sudhanshu Mishra (noreply@blogger.com) at June 16, 2014 07:32 PM

June 15, 2014

Sushant Hiray

This Week in CSymPy: #4

Week 4 for GSoC just ended and this week I completed implementing the Exponential Module.

Progress

I started the week by implementing ATan2 which was not implemented in the PR of Trignometric Module.

Pull 189 looked into implementing the ATan2 functionality. Since to implement ATan2 we needed to find whether the signs of numerator and denominators. In case of complex symbolic expressions, we are yet to implement numeric-eval in CSymPy, so in those cases we just return the object, in cases where we find the arg in lookup table we simplify and give out the result.
A more detailed version of the assumptions can be found in this comment in PR. It is also documented in the code.

After implementing ATan2, I worked on the Exponential module. Aaron gave some insights regarding the implementation of the exponential module, specifically he mentioned the fact that

It's better to special case E in one place `(pow)` 
than to special case exp everywhere in the code 
that tries to deal with pow objects 

So discussing this with Ondrej we decided to keep exp only as a function which just calls pow(E, x) underneath.

The logarithm class was similar to the current implementation of SymPy barring the part involving complex numbers. The exponential module as a whole was implemented in pull request which has been merged into master.

I have also decided to implement the LambertW function. It wasn’t a part of my original proposal but it seems pretty trivial to implement! Once this is merged into master, the exponential module will cover all major functions implemented by SymPy.

Discussions

As usual most of the actual discussions happened on PR, Aaron also gave some insights on gitter! So its good to hangout there as well.

Week Highlights!

  • This week @isuruf started helping us out in the Number Theory module which was originally added by Thilina.
  • Also I’ve noticed extra traffic in CSymPy gitter. A lot of people have started following the conversations there! That is good to see :)

The Week Ahead

  • Add the LambertW class.
  • Implement the Hyperbolic Module.


Thats all for now :) Will get back next week!

June 15, 2014 06:20 PM

Kundan Kumar

GSoC week 4: Implemented system of three linear equation of first order

I will summarize the work completed until now since beginning of first week. These included

> Implementing classify_sysode to classify all type of equations on the basis of linearity, number of equations, order of equations.

> Implementing system of two linear equations of first order which includes seven type of equation of various forms.

> Implementing system of two linear equations of second order which has eleven type of equation.

> Implementing system of three linear equations of third order which has five type of equation both of constant coefficient and variable coefficients.

> Implementing homogeneous system of n linear equations of first order of constant coefficients.

Those were the work completed till now of which last week work was on implemention of system of three linear equations of 1st order and homogeneous system of n equations. I am currently working on implementing system of two and three nonlinear system of equations of first order and aiming to complete it by next week. Its been a good progress till now and hope to work likely in coming weeks :) .


by Kundan at June 15, 2014 09:26 AM

June 14, 2014

Thilina Rathnayake

[GSoC] Week 4: Fraction Free Algorithms

Hi All,

As the week 4 of the GSoC comes to an end I was able to start implementing Matrix decompositions, starting with the LU decomposition.

Fraction Free Algorithms

Gaussian elimination is the procedure for reducing a given matrix to an echelon form. In Gauss-Jordan elimination, we reduce a given matrix into a reduced row echelon form. I implemented quite a few algorithms for Gaussian elimination and Gauss-Jordan elimination. The conventional algorithm for Guassian elimination is a very straight forward one and can be found in[1]. There are two problems in this algorithm. One is that we have to perform division which is not a good thing to do if the divisor is too small. Also, the above implementation doesn’t allow pivoting, which makes the algorithm not applicable to non-singular matrices.

To avoid the problem with division, I implemented a set of algorithms which results in fraction free matrix entries. These algorithms are described in detail in [2] and [3]. One disadvantage in these fraction free methods is that the resulting matrix entries tends to grow larger than the naive algorithm described in [1]. These can be avoided by observing that certain entries are divisible by other entries in the matrix. This is described in detail in [3].

Pivoting is also a concern in these algorithms. Since we are mostly dealing with symbolic matrices, first row below the current pow which has a non-zero element is chosen as the pivot row.

Also, I implemented a fraction free version of LU decomposition.

Plan for upcoming week

I hope to implement some other Matrix decomposition algorithms like QR and Cholesky in the coming week. Also, I wish to implement a Gaussian elimination algorithm which is more suited to matrices with symbolic entries. For the time being I am studying about it in [4].

References

[1] http://www.cs.rutgers.edu/~venugopa/parallel_summer2012/ge.html#algo
[2] Fraction-free algorithms for linear and polynomial equations (1997) by George Nakos, Peter Turner, Robert Williams
[3] A Simplified Fraction-Free integer Gauss elimination algorithm, Peter R. Turner, Ph.D.
[4] Modified Gauss Algorithm for Matrices with Symbolic Entries, Benno Fuchssteiner.


by Thilina Rathnayake at June 14, 2014 11:51 PM

Sachin Joglekar

GSoC Week 4: Vector framework done

This has been a tricky week. Not hectic, but rather tricky. Most of it was spent fixing the innumerable bugs that I encountered while polishing and ‘packaging’ the basic version of the new vector module – it’s in a PR here.

I cannot stress how helpful the pre-written unit and example tests were, in making things clearer and pointing out the many bugs in the first version of my code. In some cases, I spent quite some (read: a lot) of time fixing the issues, while in some, I just redid the code/API to work around the intricacies of the SymPy core. Its quite easy to get lost in there.

Currently, the following functionalities are supported stably by my code(with appropriate error handling wherever required)-

1. All basic vector operations-

a) Addition/Subtraction

b) Multiplication/Division by scalars

c) Dot/Cross product

2. Including coordinate variables (spatial variables) in vectorial expressions.

3. All basic use cases of the Del operator in vector/scalar expressions, including-

a) Gradient

b) Divergence

c) Curl

d) Directional derivative

I know it seems a little low for something that’s based on code which is already a part of SymPy, but the whole point was to base it directly upon the SymPy core – for users not acquainted with the physics module. Moreover, with some rudimentary timing techniques, I found that on an average, the new module was able to do a set of most-basic vector operations (add, sub, dot, cross) approximately 3-4 times faster than sympy.physics.vector. But I guess it’s too early to judge now, since the new module has no overhead of coordinate systems.

I can say I am successful (hopefully all the tests on the PR should pass soon) in supporting all the functionality mentioned above, in a stable implementation – though I am still waiting for Jason and the SymPy people to review the PR. I hope it gets in soon.

The next step would be to start working on a new branch that would include the classes for coordinate systems and stationary points in 3D space (cartesian system). This is going to be tricky – the code for these classes in sympy.physics.vector is enough to prove that. As usual, the first step would be to write out the ‘expected-to-succeed’ unit and example tests. Hopefully, by that time, I would get sufficient feedback on the current PR too. But since the underlying vector framework is working well, I wouldn’t have to worry about bugs in that area – I can focus purely on the one-level-up code for multiple coordinate systems.

Cheers to a fun week! Hopefully I will be reporting just as much of progress next week too. Have a great week ahead :-).

 


by srjoglekar246 at June 14, 2014 08:15 PM

June 13, 2014

Jim Crist

GSoC Week 4: Linearizing Lagranges Equations

As of last week the Linearizer class implementing the general form discussed in Luke and Gilbert's paper was completed. The methods contained for linearization work for any system that can be expressed by any combination of the following:

\begin{aligned} f_{c}(q, t) &= 0_{l \times 1} \\ f_{v}(q, u, t) &= 0_{m \times 1} \\ f_{a}(q, \dot{q}, u, \dot{u}, t) &= 0_{m \times 1} \\ f_{0}(q, \dot{q}, t) + f_{1}(q, u, t) &= 0_{n \times 1} \\ f_{2}(q, \dot{u}, t) + f_{3}(q, \dot{q}, u, r, t) &= 0_{(o-m) \times 1} \end{aligned}

with

\begin{aligned} q, \dot{q} & \in \mathbb{R}^n \\ u, \dot{u} & \in \mathbb{R}^o \\ r & \in \mathbb{R}^s \end{aligned}

This works for most systems (it was derived with kanes method in mind specifically). However, systems expressed using lagranges method can't be brought into this form.

This week I spent some time rederiving the general form to make it fit lagranges method, as well as kanes method. I plan to write up a formal paper expressing the derivation as a reference and documentation of the class; here I'll just give a brief overview.

In general, Lagrange's Method expresses the system using 3 equations:

\begin{aligned} m_{c}(q, t) \dot{q} + f_{c}(q, t) &= 0_{m \times 1}\\ m_{dc}(\dot{q}, q, t) \ddot{q} + f_{dc}(\dot{q}, q, t) &= 0_{m \times 1}\\ m_{d}(\dot{q}, q, t) \ddot{q} + \Lambda_c(q, t) \lambda + f_{d}(\dot{q}, q, r, t) &= 0_{n \times 1}\\ \end{aligned}

with

\begin{aligned} q, \dot{q}, \ddot{q} & \in \mathbb{R}^n \\ r & \in \mathbb{R}^s \\ \lambda & \in \mathbb{R}^m \end{aligned}

In this case, the first equation encompass the time differentiated holonomic constraints, as well as the nonholonomic constraints. The second equation is then the time derivative of the first equation. The third equation represents the dynamics of the system, as formed by the lagrangian. The lagrange multipliers ($\lambda$) enforce these constraints.

With some rearranging of the above, they can be merged with the previous general form for Kane's Method, forming a set of equations that should be able to contain most equations of motion:

\begin{aligned} f_{c}(q, t) &= 0_{l \times 1} \\ f_{v}(q, u, t) &= 0_{m \times 1} \\ f_{a}(q, \dot{q}, u, \dot{u}, t) &= 0_{m \times 1} \\ f_{0}(q, \dot{q}, t) + f_{1}(q, u, t) &= 0_{n \times 1} \\ f_{2}(q, u, \dot{u}, t) + f_{3}(q, \dot{q}, u, r, t) + f_{4}(q, \lambda, t) &= 0_{(o-m+k) \times 1} \end{aligned}

with

\begin{aligned} q, \dot{q} & \in \mathbb{R}^n \\ u, \dot{u} & \in \mathbb{R}^o \\ r & \in \mathbb{R}^s \\ \lambda & \in \mathbb{R}^k \end{aligned}

Note that the only changes are the addition of a $u$ term in $f_2$, and the $f_{4}$ term holding the lagrange multipliers. For Lagrange's method, $\dot{q} = u$, and $k = m$; for Kanes method $k = 0$, and everything looks the same as it did before.

The returned $M$, $A$, and $B$ linearized form then is:

$$ M \begin{bmatrix} \delta \dot{q} \\ \delta \dot{u} \\ \delta \lambda \end{bmatrix} = A \begin{bmatrix} \delta q_i \\ \delta u_i \end{bmatrix} + B \begin{bmatrix} \delta r \end{bmatrix} $$

where

\begin{aligned} M &\in \mathbb{R}^{(n+o+k) \times (n+o+k)} \\ A &\in \mathbb{R}^{(n+o+k) \times (n-l+o-m)} \\ B &\in \mathbb{R}^{(n+o+k) \times s} \end{aligned}

As before, the $M$ matrix can be inverted, and the square state space matrices $A$ and $B$ calculated.

The functionality described above has been implemented in the LinearizeLagrange branch of sympy on my github. As this is a superset of the functionality I implemented last week, I'm going to hold off on submitting this to master until my current pull request is merged. For now I made a local PR here. Please take a look through, I need all the code review I can get.

Two tests have been implemented linearizing a system generated with Lagrange's Method. I plan on adding more next week, as well as improving the documentation.

Linearizing a Non-minimal Pendulum with Lagrange's Method

A demonstration of the current functionality for a non-minimal realization of a pendulum is below. The pendulum has two generalized coordinates, $q1$ and $q2$. As this is Lagrange's method, the generalized speeds are just the time derivatives of the coordinates (i.e. $u = \dot{q}$).

Pendulum system

# Create the required symbols
q1, q2 = dynamicsymbols('q1:3')
q1d, q2d = dynamicsymbols('q1:3', level=1)
L, m, t = symbols('L, m, t')
g = 9.8

# Compose World Frame
N = ReferenceFrame('N')
pN = Point('N*')
pN.set_vel(N, 0)

# A.x is along the pendulum
theta1 = atan(q2/q1)
A = N.orientnew('A', 'axis', [theta1, N.z])

# Create point P, the pendulum mass
P = pN.locatenew('P1', q1*N.x + q2*N.y)
P.set_vel(N, P.pos_from(pN).dt(N))
pP = Particle('pP', P, m)

# Constraint Equations
f_c = Matrix([q1**2 + q2**2 - L**2])

# Calculate the lagrangian, and form the equations of motion
Lag = Lagrangian(N, pP)
LM = LagrangesMethod(Lag, [q1, q2], hol_coneqs=f_c, forcelist=[(P, m*g*N.x)], frame=N)
LM.form_lagranges_equations()

At this point the equations of motion have been formed, but not linearized. Linearization requires that dependent and independent coordinates be chosen. In this case we'll chose $q2$ as independent, and $q1$ as dependent.

# Choose the independent and dependent coordinates
q_i = Matrix([q2])
q_d = Matrix([q1])
u_i = Matrix([q2d])
u_d = Matrix([q1d])
linearizer = LM.to_linearizer(q_i, u_i, q_d, u_d)

# Compose operating point
q_op = {q1: L, q2: 0}
u_op = {q1d: 0, q2d: 0}
ud_op = {q1d.diff(t): 0, q2d.diff(t): 0}

# Perform the linearization
A, B = linearizer.linearize(q_op=q_op, u_op=u_op, ud_op=ud_op, A_and_B=True)
print(A)
print(B)

Output:

Matrix([
[           0, 1],
[-2*lam1(t)/m, 0]])

Matrix(0, 0, [])

Note that the lagrange multiplier apppears in the linearization. However, for a given operating point, each multiplier has a specific value (i.e. they're not free choices). Using the structure of the system of equations, the values can be solved for, and substituted in:

# Take advantage of problem structure to solve for lams
mass_matrix = LM.mass_matrix.col_join((-LM.lam_coeffs.row_join(zeros(1, 1))))
force_matrix = LM.forcing.col_join(LM._f_cd)
lam_op_vec = Matrix((mass_matrix.inv()*(-force_matrix))[-1:])
lam_op_vec = lam_op_vec.subs(ud_op).subs(u_op).subs(q_op)
lam_op = dict(zip(LM.lam_vec, lam_op_vec))

# Substitute the value for the multipliers at this operating point
A = A.subs(lam_op)
print(A)

Output:

Matrix([
[     0, 1],
[-9.8/L, 0]])

This is the correct linearization for a pendulum linearized about hanging at rest operating point. You can try out the added functionality demonstrated above by cloning my LinearizeLagrange branch of sympy here.

While functional, it still isn't finished. I still need to add documentation, more tests, and finalize the interface. I plan on working on this next week.

<script type="text/javascript"> if (!document.getElementById('mathjaxscript_pelican_#%@#$@#')) { var mathjaxscript = document.createElement('script'); mathjaxscript.id = 'mathjaxscript_pelican_#%@#$@#'; mathjaxscript.type = 'text/javascript'; mathjaxscript.src = 'https:' == document.location.protocol ? 'https://c328740.ssl.cf1.rackcdn.com/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' : 'http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'; mathjaxscript[(window.opera ? "innerHTML" : "text")] = "MathJax.Hub.Config({" + " config: ['MMLorHTML.js']," + " TeX: { extensions: ['AMSmath.js','AMSsymbols.js','noErrors.js','noUndefined.js'], equationNumbers: { autoNumber: 'AMS' } }," + " jax: ['input/TeX','input/MathML','output/HTML-CSS']," + " extensions: ['tex2jax.js','mml2jax.js','MathMenu.js','MathZoom.js']," + " displayAlign: 'center'," + " displayIndent: '0em'," + " showMathMenu: true," + " tex2jax: { " + " inlineMath: [ ['$','$'] ], " + " displayMath: [ ['$$','$$'] ]," + " processEscapes: true," + " preview: 'TeX'," + " }, " + " 'HTML-CSS': { " + " styles: { '.MathJax_Display, .MathJax .mo, .MathJax .mi, .MathJax .mn': {color: 'black ! important'} }" + " } " + "}); "; (document.body || document.getElementsByTagName('head')[0]).appendChild(mathjaxscript); } </script>

by Jim Crist at June 13, 2014 11:00 PM

June 11, 2014

Sushant Hiray

This Week in CSymPy: #3

Week 3 for GSoC just ended and we decided to work on fixing bugs which were high priority.

Progress

I started off the week and decided that I’ll work on the exponential module.

However last week we had discovered a bug in the way CSymPy was handling powers, so we decided to fix it first before moving on to the exponential Module.

The main issue in the bug was the way we were handling the dict. There was no uniformity so we needed up inserting 2 as well as 1/2 into the dict. We discussed with Aaron and the way SymPy handled the issue of powers was to make sure that the power was always positive. If not the reciprocal was added in the dict.

I used the convention SymPy was using and sent a pull request. However, there wasn’t much of help usin this convention as we were ending up with both 2 and 1/2 in some corner cases. To avoid this, I searched the dict first if the reciprocal existed, if it existed, then I added the reciprocal or else added the number to the dict.

However as was expected using a find operation degraded the benchmarks heavily. So after a lot of discussion with Ondrej over the PR as well as on gitter, we figured out a way to solve the problem. The convention was to always enter the number in the form num/den such that num >= den. We updated the PR according to this convention and the issue was fixed. Consuming almost a weeks time!

Discussions

Most of the discussions with Ondrej and Thilina happened over the pull requests and the issues itself. We finalized the convention over gitter.

The Week Ahead

  • Add the ATan2 class.
  • Implement the Exponential Module. After this module we shoule be able to do things like rewrite__as_cos for the exponential class.



Thats all for now :) Will get back next week!

June 11, 2014 07:30 AM

June 10, 2014

Avichal Dayal

In the third week I started with the second phase of my project which is to implement formal power series. I spent some time in designing the class structure. Another important part is to figure out a way to represent infinite power series.
Literature mentions two ways to represent infinite series - one is to have function to compute the kth coefficient while other is lazy evaluation. Many series operations seem to have very simple on-line algorithms if we use lazy evaluation,

To work with lazy evaluation, FPS(formal power series) will be represented as a stream object consisting of head and tail. Head will be the first term of the series while tail is the rest of the series but in an unevaluated form.
So FPS of exp(x) will have its generator as x**n/n!. Head will be 1 and tail will be FPS with generator as 1/(n+1)!. Tail will remain unevaluated. This way any number of terms can be found.

To find generator of any function, Gruntz and Koepf give an algorithm.
It consists of three steps:-
- Finding a simple differential equation for the given function whose coefficients are rational functions.
- Converting the differential equation to a recurrence relation.
- Solving the recurrence relation to get a closed form.
If the original function was a rational fraction, it can split using partial fraction decomposition (already implemented in SymPy as apart) and apply the binomial series to get closed form.
For e.g.:- 1/(1-x) can be expanded as 1 + x + x**2 + x**3 + ... So the generator in this case is x**n.

Some initial results are:-


After completing the algorithm to compute generating function, next step will be to implement series operations like multiplication (cauchy product), inversion, convolution etc.

by Avichal Dayal (noreply@blogger.com) at June 10, 2014 06:03 PM

Sudhanshu Mishra

GSoC'14: Third week

It was another not so productive week but I learnt a lot regarding how to write new classes by extending classes like BasicSymbol and Expr.
Last week I got stuck with the rewrite expression mechanism which turned out to be a bug in the core. Thanks to Sean for the fix! I need to do some more cleanup to get my last PR merged.
I also tried to fix an issue related to assumptions in Function but the complex behaviour of __new__ constructor makes it difficult to follow the inheritance. I also need to understand the working of SymPy's  @cacheit  decorator.

Lately I've been working on the implementation of optical medium. The motivation behind defining aMedium came from the laziness of passing constants like electric permittivity and magnetic permeabhility of the medium everywhere. It will help carry out operations(events) very easily in optics.
With the merger of Point3D by Akshay I planned to implement Fermat's principle and thus felt the need of Medium first.
I've already sent a WIP pull request for Medium. Here's the link:
I also took a quick glance at variational calculus to solve the problem involving Fermat's principle. I'm planning to start this as following:
class FermatsPrinciple(Expr):

def __init__(self, p1, p2, v=Symbol('c')):
if isinstance(p1, type(())) and isinstance(p2, type(())):
if len(p1) == 2 and len(p2) == 2:
self.p1 = Point(*p1)
self.p2 = Point(*p2)
elif len(p1) == 3 and len(p2) == 3:
self.p1 = Point3D(*p1)
self.p2 = Point3D(*p2)
elif isinstance(p1, Point) and isinstance(p2, Point):
self.p1 = Point3D(p1.x, p1.y, 0)
self.p2 = Point3D(p2.x, p2.y, 0)
elif isinstance(p1, Point3D) and isinstance(p2, Point3D):
self.p1 = p1
self.p2 = p2
else:
raise TypeError("p1 and p2 can only be tuple, Point or Point3D")
It's not a very good API and needs more refining which I'll discuss with Sean once I get it working.
My work in this week(fourth week as I am writing this post very late) will be to complete the above discussed things.
On a lighter note, the temperature here is 46 degrees and no sign of rains yet. :(
That's all for now. Cheers!
If you have some feedback or questions regarding this post, please add comments. I would be happy to get some feedback.

by Sudhanshu Mishra (noreply@blogger.com) at June 10, 2014 10:08 AM