GSoC (officially) done!

GSoC is now over, and I have officially passed the Final Evaluations :-).
Looking back, this was an amazingly interesting and exciting process, and helped me learn a lot about SymPy, and also Python. Quite a few things that I learnt in the last three months (some resulting from my own coding, others from just messing around with SymPy code) were so full of Python 'magic' that I wondered whether I would ever use them again in any of my own work!
There were a few problems, like dealing with architectural difficulties, physics-related problems (that I had not though of while writing my proposal) and the radical change in the plan of action mid-way through the GSoC period. This change came because Prasoon's module didn't shape up as we had planned earlier- perhaps we underestimated the difficulties associated with creating such a complex module and integrating it with the rest of the codebase. As a result the code I had written based on our decided API became useless - almost. Thankfully,  while coding all the new stuff (a mechanics core to support fields, electrostatics, a unified class for ReferenceFrame and Particle) I had made a modified version of the mechanics core to test it all on. Hence, Gilbert and I decided around July end that the best thing to do would be to modify the current framework accordingly and base all my work on it.
After the mid-sem evals, I could not work as hard as I had during my vacations due to projects and course-work at college. Inspite of all that, I managed to get two PRs merged...one modifying sympy.physics.mechanics.essential to support scalar and vector fields, and other to add a big function to the functions.py file to calculate motion attributes from time-dependent vectors and boundary conditions.
I have a few more PRs in the pipeline, these would essentially just modify the code I wrote over the summer to work with the current module and add the documentation for the done work. The first one of them is already in the review process :-)
The help from Gilbert, Jason, Stefan and at times even Aaron, has been immense, and I am really thankful to them for it. Gilbert was a great mentor, especially while brainstorming solutions to problems that we faced from time to time. It was almost like solving problems with a college and getting tips to understand how some difficulties could be resolved.
Obviously, I will continue to work for SymPy (mainly sympy.logic) and PyDy - though mostly PyDy for the foreseeable future, since I have to add the E-M module and extend it as I had envisioned earlier. It would be fun to code complex electromagnetic concepts to work with dynamic systems, and having Jason, Gilbert and DL Peterson (and the rest of the PyDy team) to help would be quite the experience on its own.
For me, GSoC has just started my involvement with this community and its codebase, and I aim to be an active developer for them :-D Physics is one of the few things that I miss being in Computer Science, so writing Python code based on it in my free time is something that I am obviously looking forward to!
I will keep updating my blog as and when I get something merged or I work on something that's worth writing about. Once again, thanks a lot to SymPy, Gilbert, Aaron, Ondrej, Stefan and obviously- Google and Carol, for this amazing opportunity!

GSoc: Weeks 12 and 13

Work gained pace in the past two weeks.
First off, I got my first GSoC PR merged into the SymPy base. This PR essentially modified the sympy.physics.mechanics core to incorporate the concept of coordinate variables and manipulations on them.
For example, consider a vector V defined in frame R1, with coordinate variables of frame R2 occurring as symbols in its definition. Now,  while calculating d/dt(V) wrt R1, it would be necessary to substitute values of the base scalars of R2 wrt those of R1(as they will be functions of t). One small point- Coordinate variables, as far as possible should not be used as normal Symbols. This may mislead the 'dt' and 'express' methods to try substitution of variables of frames not linked to each other.

We also added a few optimisations, such as a dcm cache(turned out to be trickier than I thought) to 'remember' the DCMs once calculated, for later use. Also added the 'express' and 'dt' methods to the ReferenceFrame class- this would help re-expressions of scalar as well as vector fields in the said frame. During all these changes, I took care to preserve the old API of SymPy classes and methods- we din't want to make all the code written wrt the previous module version redundant.

Now, I currently have a PR in the pipeline for adding the get_motion functions to the sympy.physics.mechanics.functions file. After this, I will turn to the rather-tricky task of modifying the Point class to incorporate functionality as needed by ParticleCharge from my EM work (all the while, making sure the previous API does not change). Parallely, I will be adding documentation on my previous PR to the Sphinx docs.

However, all this after the 18th- the day my first semester exam ends. After that, I can start working on the project full-time again :-).

Anyways, have a great week :-)

GSoC Week 10 and 11 : Modifying the mechanics core

Unfortunately, as GSoC is drawing to an end, my workload at college seems to be increasing too. Hence, I haven't been able to work for the project as much as I would have liked, though I have a new PR in queue. Since Prasoon's Vector module isn't close to completion yet, Gilbert and I decided to go ahead and modify and hack the current framework to make space for the modifications I require for my work.
Last week, I opened a PR with the said modifications, but as I had expected, there are a *lot* of errors to be resolved. Changing one thing at one place gives rise to a hundred errors elsewhere in the module - my lesson for the week. We decided to go ahead and implement the caching of dcms and variable maps generated wrt pairs of reference frames, and it has led to quite some inconsistencies here and there. I am still struggling with them, I hope I can resolve all of them in a day or two. By this week-end, I should have this PR merged and a new one, with the new classes, opened.
About EM, I dint get to do much, I just implemented ideal dipoles and wrote tests for them. I still haven't pushed them to the PR. Will do that soon, I guess.
That's all there is for this week. Most probably, my GSoC work is going to stretch beyond the timeline, though that's not really a problem. Hope at the end of it, the code is perfectly shippable :-)
Have a good week!

GSoC Week9 : Beginning of electrostatics

Week 8 of GSoC is now over, and just 5 more weeks remain till the hard pencils-down date.
However, my work on the project will most likely exceed that time, since I will be having my first tests around September 17. I am sure I will be able to get all the code in by that time, but to make the entire work shippable, things like Sphinx documentation, doctests etc will also have to be finished, which may take a little more time beyond the actual deadline.

Anyways, coming to this week's work, I wrote some code for electrostatics. I added some functions and the ParticleCharge class (inheriting from Particle of mechanics module).
However, the E-M module's handling of charged particles required me to modify my earlier implementation of Particle- most importantly, adding a 'set_motion' method. Like I mentioned earlier, Particle earlier had a frame attached to it. Gilbert and I discussed at length on this, and we agreed that making the user initialize a new frame everytime, just to attach it to a Particle, was quite cumbersome.
Hence, I modified the class in the mechanics core- It still has a frame as an attribute, but its now 'hidden' from the user (nothing in Python in hidden from the user, but what I mean here is that the user can do everything with a Particle without caring about the frame). The set_motion method I mentioned earlier just initializes a new frame according to the motion parameters set by the user (similar to MovingRefFrame).

The main reason I want a frame attached to a Particle is the ease of inner working. Since a Particle will not have an 'angular acceleration', those parts of the MovingRefFrame code will never really get executed in Particle. However, Gilbert has a point in pointing out that its unnecessarily adding a complex attribute to the class. We havent finalized anything yet, but re-writing the motion-setting methods for Particle would be tedious, and more importantly, will lead to quite a lot of code duplication. Lets see what we end up doing.

Anyways, I also finished writing tests for the code I wrote this week, and I am waiting for a review from Gilbert. Jason and he also encouraged me to add example-tests to ensure everything works smoothly.
Here are the example tests I wrote (they are really basic, but its better that way, since strenuous tests only make things worse)-

 def test_example_1():  
   """  
   This is a sample example meant to test the basic functionality of  
   electrostatics.py  
   
   Consider a frame R, and a particle charge P(mass m, charge q) situated  
   at its origin initially.  
   Now, at time t=0, a field R.x+R.y+R.z is 'switched on' and kept that way  
   for time 't0'.  
   The user needs the x-coordinate of the particle charge as a function of  
   time, after the field is switched off at time=t0.  
   """  
     
   #Basic initialization  
   m, q = symbols('m q')  
   P = ParticleCharge('P', m, q)  
   P.set_motion(R, pos_vector = 0)  
   field = R.x + R.y + R.z  
   time = dynamicsymbols._t  
   t0 = Symbol('t0')  
   #The acceleration is equal to the electrostatic force experience by the  
   #particle charge, divided by its mass  
   acceleration = P.electrostatic_force(field)/P.mass  
   #Use get_motion_acc from the mechanics core to find velocity and position  
   #parameters using acceleration function and boundary conditions  
   translation = get_motion_acc(acceleration, \  
                  0, P.pos_vector_wrt(R), frame = R)  
   #Calculate the motion parameters of the particle charge at time t0  
   acceleration = translation[0]  
   velocity_at_t0 = translation[1].subs({time:t0})  
   position_at_t0 = translation[2].subs({time:t0})  
   #Set motion of P accordingly  
   P.set_motion(R, trans_acc = 0, trans_vel_b = velocity_at_t0,  
          pos_vector_b = position_at_t0)  
   #assert that the dot product of P's pos_vector wrt R, with R.x  
   #matches the actual result  
   assert P.pos_vector_wrt(R).dot(R.x) == \  
       q*t*t0/m + q*t0**2/(2*m)  
   
   
 def test_example_2():  
   """  
   This is a sample example meant to test the basic functionality of  
   electrostatics.py  
   
   Consider a set of 3 particle charges of equal mass and charge placed  
   symmetrically on a circle of unit radius around the origin, with one  
   point lying at position vector of R.y  
   We calculate the energy required to assemble this system of charges from  
   infinity(with no kinetic energy in any particle) and find out values of  
   the electrostatic fields and potentials at the origin  
   """  
   
   #Basic initialization  
   from sympy import sin, cos, pi  
   m, q = symbols('m q')  
   P = ParticleCharge('P', m, q)  
   Q = ParticleCharge('Q', m, q)  
   S = ParticleCharge('S', m, q)  
   #Place particle charges as per required configuration  
   P.set_motion(R, pos_vector = R.y)  
   Q.set_motion(R, pos_vector = -cos(pi/6) * R.x + -sin(pi/6) * R.y)  
   S.set_motion(R, pos_vector = cos(pi/6) * R.x + -sin(pi/6) * R.y)  
   #Check outputs of required values  
   assert charge_assembly_energy(P, Q, S) == 3*3**(-0.5)*k*q**2  
   assert S.electrostatic_field(R, 0) + \  
       Q.electrostatic_field(R, 0) + \  
       P.electrostatic_field(R, 0) == 0  
   assert P.electrostatic_potential(R, 0) == \  
       Q.electrostatic_potential(R, 0) == \  
       S.electrostatic_potential(R, 0)  

If you think the API is a little complex, I think so too, but those modifications are for later. It would be good, according to me, to have something like an 'environment' infrastructure. Basically, the user will keep adding components (fields, particles etc etc) to it, and the constituents will keep getting modified accordingly. However, it will be good to have the module at a certain level of working before thinking on those lines.

That's all there is to report for this week. The coming week, I will be reading up on magnetostatics from Griffith's and start coding for the same. However, I am afraid things are only going to get more and more complex hereon. Lets hope for the best.

Have a great week :-)

GSoC: Week 8

Week 8's over.
This week was more of 'studying' and 'preparing' rather than coding. Since the entire infrastructure behind my initial proposal has now shifted considerably, I am drawing out the plans anew and getting the basic things done before I start working on sympy.physics.em (tentative name).

First off, I finally got the logic docs up at the SymPy development docs. You can see them here.
Second, I wrote the basic field functions for electrostatic theory (including the tests).
Apart from that, I was mostly busy with college stuff (first week of classes are usually a little difficult for adjusting) and going through Griffith's book. So hopefully, I will get a lot done by the end of the coming week (code + further plans).
Prasoon's code is shaping up well, and I may help him with the testing and debugging part to speed up the development of sympy.vector. I feel quite optimistic about the projects (Prasoon's and mine) now.

Well, thats all (short post, I know) for now.
Will come back with more next time.
Have a great week :-)

GSoC Midterm Report (Week 7)

Week 7 of GSoC is over, which means I am now in the 'second half' of the coding period. The midterms were conducted by Google from 29th to 2nd, and I am happy to report I passed them :-).

As I had predicted in my last blog post, this week wasn't very exciting, atleast as far as the project is concerned. I am now at Goa, and classes have started. However, I did get a few things done - The Dyadic framework is now almost complete, with all the classes functioning as expected (from point of view of dyadic operations). I also added all the special operations like cross, dot, express and time-differentiation. For now, just like last week, they use SymPy Symbols for temporary working.

Moreover, I finished writing all the tests for my first PR. Gilbert reviewed them, and I made a few changes (mostly API-wise).
I also tested my MovingRefFrame implementation of vector time-derivatives on my hacked mechanics module, and I am happy that it worked fine. Thus, I now have the proof-of-concept for all the work I have done so far.

I have also submitted a PR to get the logic module docs into the online development documentation, but I am quite sure I am going wrong somewhere. I have pinged Aaron, and I hope it gets solved by the end of this week.

To sum up my work till now, this is what I got done-

1) MovingRefFrame class - The basic class to represent moving reference frames in 3-D space. Inherits from the CoordSysRect class of the new vector module.

2) Particle class - Class to represent particles (mass but no volume) in space. Each particle is associated with a frame of its own, and all it's parameters are calculated wrt it.

3) Miscellaneous mechanics functions like the ones to obtain position vector given a time-dependent vector and boundary conditions.

4) Dyadic framework - Classes Dyadic, BaseDyadic, DyadicAdd, DyadicMul. All these classes are based on relevant SymPy core classes and build dyadic-based functionality on top of them.

Currently, I am reading Griffith's book (very helpful suggestion from Aditya - my collegemate) and taking notes for the EM module. Since the development of the new mechanics core is now nearing completion, I can focus time on working on the new addition to sympy.physics.

Hence, the coming week will be spent on polishing off the mechanics work, making sure all required parts are in place, maybe helping Prasoon with some functions of vector fields (so that he can concentrate on vector-space integration).
Plus, I have to design the vector-dyadic interfacing. Basically, how dyadics will be instantiated from vector 'outer' functions. Obviously, having a dyadic framework is useless without fitting it to the vector module.
And last but not the least, I plan to code the electrostatic functions for the EM sub-module by this weekend. Hope I succeed.

Have a great week :-)

GSoC Week 6: First half(almost) done

First, about the work done this week. It took me some time to set up my programming environment and get used to Ubuntu's interface. I fumbled around a bit messing up some things on the way, but now I can say I am quite comfortable...at least with the basics.
For one, I got the docs for my work on the logic module merged into the master docs. Aaron also pointed me to the method for merging those changes with the development version  of the SymPy docs. I will try and get it done in a day or two, before I leave for my campus.
Second, Gilbert and I submitted a PR to Prasoon's branch with some changes that we made to his code. Since the number of merge conflicts are high, he would most probably be merging the relevant parts manually.
Third, I started working on the code for dyadics. Initially, I submitted a slightly modified version of the earlier Dyadic class, but Gilbert soon made me realise that compatibility with SymPy's core would mean a lot more than just a few changes here and there. Sigh.
So, I am now busy working on classes for dyadics, compatible with SymPy's global architecture-
1) Dyadic(Expr), the super class - all operations and initializations will be handled from here
2) BaseDyadic(Dyadic) - the class to represent basic components of dyadics, things like (R1.x|R2.y)
3) DyadicAdd(Add, Dyadic) and DyadicMul(Mul, Dyadic) - additions and multiplications of dyadics

The last two classes are a little shaky for now, and the code I have submitted at the PR at the time I write this makes it quite evident. Inspite of that, I have got the basic operations- add, sub, mul, div and some others like 'and' to work as expected(conceptually). To see how it looks as of now, you can see this (real SymPy session with mock arguments). However, it's quite obvious that getting these new classes to behave exactly according to SymPy's way of doing things and return accurate results is going to take quite some time, maybe a week more. Not to forget docs, doc examples and tests.

Anyways, I can't believe I am almost at the end of the first half of my GSoC period. I have learnt *quite* a lot, well, that's obvious from my blog posts- not just vector math, but a bunch of Python as well. I am happy with the progress till now, though nothing has been merged yet. But well, Prasoon and my projects are such that when things will get merged, a huge chunk will go in together.
I plan to get the mechanics core done perfectly before I start with the EM module..and by 'done perfectly', I mean getting the main code merged, along with a lot of commentary-style documentation on the new module. I also have to enquire about putting deprecation signs on the old module- will have to consult the mailing list for this.

As mentioned earlier, I will be moving to Goa on Tuesday, so that period is going to be quite busy with getting things in order @the insti. Then I can start working again.
The mid-term GSoC evaluations will be done from 29th to 2nd of August, after which I will get my first big payment in 4-5 days :-)

Thats all for now!
Have a great week, and 'stay tuned' :-D.

GSoC: Week5

A late post this week.

In terms of progress, last week was a bit slow, but things are going to move pretty fast now.

The past seven days, I wrote the tests for the code of MovingRefFrame, Particle and all the helper functions in the functions.py file. Trust me, if writing tests are the only thing you are doing, it can get pretty boring. Testing your code and experimenting with it is a lot of fun, but sadly, writing unit tests is not. Anyways, I finished writing most of them for the classes built till now. Some more may have to be added, but that will be after Prasoon's code gets pushed into SymPy's master repo.

Apart from that, I did some groundwork on Dyadics, which would be the focus of my GSoC work for the upcoming week. I will be taking hints from the old module's Dyadic class, but I will have to modify and hack it considerably to fit SymPy's global architecture. Plus, some methods like the __eq__ one need improvement, which I plan to do. I hope I can finish the class (operation methods as well as those for printing, etc) by the end of this week. Tests later.

One great thing that happened this week- I finally decided to give up Windows as my coding environment. Its not like you cannot code in Windows; if you are good enough for what you are trying to do, you can pretty much do whatever you want there...but in Linux, its 10x easier. Especially since any library you want is just one small command away.

My laptop's LAN card got replaced, and it prompted me to set up a dual boot system on my machine.

Anyways, I am very happy with the linux(Ubuntu 12.04 LTS) experience!

Yeah...one more thing! SymPy 0.7.3 finally got released, and this will be the first official release to include code from me (the work on sympy.logic)! However, the documentation for that code hasn't been added to the website docs, so I probably need to learn the basics of Sphinx documentation for that. I hope I can get it done this week. I tried working with Sphinx on Windows earlier, but it quite frankly was a little tricky for an impatient guy like me.

That's all the 'progress' that happened this week.

The next week's post will be the last GSoC update I'll be posting from Mumbai. I'll be going to my campus at Goa on next Tuesday, that week's work may be a little bit delayed around the Tue-Wed period.

But thats for later.

Have a great week! 

GSoC Week4: Improvements

This has been a slow week, with a lot of ups and downs, work-wise.
Firstly, my laptop's LAN card decided to stop working properly, keeping me busy talking to the Dell customer service and convincing them that something is wrong with the thing. After a lot of troubleshooting and failed tests on their part (they even tried remote desktop connection to see what was going wrong), they finally agreed that something needs to be done. I hope it gets fixed before I go back to campus.

As far as the GSoC work is concerned, firstly, I got the Particle code finalised on Tuesday after a review by Gilbert.
Next, I spent the whole of Wednesday hacking the old mechanics framework to act as the base to test my own code on. I felt the need to do this, as I was anxious to see whether my own work was functioning as I thought it would, and it would take some time before Prasoon could get his framework in a completely consistent state (thats natural, considering the complexity of the work involved). It took a lot of time, as I had to tweak the old framework to support coordinate variables, use substitution on them during re-expression and time-differentiation, and so on...
Finally, I got it to work and tested my own code. It worked perfectly! There were a few errors that had to be rectified for the tuple-args format to function, but the rest of it went smoothly.

This is a real SymPy session using the hacked module and my own work-

1:  >>> O = MovingRefFrame('O')  
2:  >>> F1 = MovingRefFrame('F1', trans_vel = 3 * O.x + 4 * O.y,  
3:                    pos_vector_b = O.x + O.y + O.z, parentframe = O)  
4:  >>> F1.pos_vector_in(O)  
5:  (3*t + 1)*O.x + (4*t + 1)*O.y + O.z  
6:  >>> Q = Symbol('Q', positive = True)  
7:  >>> F2 = MovingRefFrame('F2', ang_vel = Q * O.z, parentframe = O)  
8:  >>> F3 = MovingRefFrame('F3', pos_vector = F2.x, parentframe = F2)  
9:  >>> F3.pos_vector_in(O)  
10:  F2.x  
11:  >>> O.express(F3.trans_acc_in(O))  
12:  - Q**2*cos(Q*t)*O.x - Q**2*sin(Q*t)*O.y  
13:  >>> P = Symbol('P')  
14:  >>> F4 = MovingRefFrame('F4', orient_type = 'Axis', orient_amount = [P * t, F1.x], \  
15:                    parentframe = F1)  
16:  >>> F4.ang_vel_in(F3)  
17:  P*F1.x - Q*O.z  
18:  >>> F5 = MovingRefFrame('F5', orient_type = 'Axis', orient_amount = [Q, F1.x], \  
19:                    trans_vel = ((0, 0, 1), 0), parentframe = F4)  
20:  >>> F3.express(F5.pos_vector_in(O))  
21:  (-t*sin(Q*t)*sin(P*t + Q) + 4*t*sin(Q*t) + 3*t*cos(Q*t) + sqrt(2)*sin(Q*t + pi/4))*F3.x + (-3*t*sin(Q*t) - t*sin(P*t + Q)*cos(Q*t) + 4*t*cos(Q*t) + sqrt(2)*cos(Q*t + pi/4))*F3.y + (t*cos(P*t + Q) + 1)*F3.z  
22:  >>> F5.pos_vector_in(F4)  
23:  - t*sin(Q)*F4.y + t*cos(Q)*F4.z  
24:  >>> F5.trans_vel_in(O)  
25:  (-P*t*cos(Q) - sin(Q))*F4.y + (-P*t*sin(Q) + cos(Q))*F4.z + 3*O.x + 4*O.y  

Line 3 means that frame F1 should have an initial position of O.x + O.y + O.z wrt O. Substituting t = 0 in the output at line 5 will confirm this. To put in the boundary conditions at some other value of time, user just has to enter the keyword arg 't'.
Line 21 adequately shows how complicated things can get after a while.
Line 19 is an example how a tuple argument can be used in the new framework to position/orient new frames wrt their parents. Line 19 essentially means that F5 will have a translational velocity of 1 * F5.z after being oriented as per the user's conditions.
(All these things are adequately explained in the docs)

This work with the old module convinced me that dynamicsymbols should also be a part of my own work, since it makes handling of time variables considerably simpler. Hence, I added the function in my last commit to the PR.

Last but not the least, I planned out some basic vector/scalar field functions for my E-M module. But more on that later.

Now that I have my entire workspace up and running on my old laptop (temporarily), I will continue with the following work in the coming time-
1) Enable users to enter initial orientation in the same way they enter the time-dependent orientation. As of now, they get to do so only by entering the initial rotation in terms of a vector. Awkward, I agree.
2) Write tests for the MovingRefFrame and Particle classes as well as all the helper functions. Now that I have a system to test these things on, it will be much easier to do so.
3) Start working on the last part of the new framework - dyadics.

Thats all for now.
Have a great week ahead :-)

Singletons in Python (and GSoC update)

This week, there isn't any interesting stuff to write about as far as the project is concerned, as I pretty much cleaned up what I had done in the previous two weeks. By 'clean up', I mean PEP-8 modifications to the existing code and rectification of some tiny(but serious) conceptual errors in the MovingRefFrame methods. I also spent quite some time reviewing, or rather, trying to understand Prasoon's code. I pulled his vector branch on my machine and started executing some tests on it. I have left some comments on his PR, and hopefully they will get sorted out soon. One good thing is, I also wrote the first draft of my code for Particle class-including methods for angular momentum, translational motion calculation, etc. This wasnt much of a problem, as I 'attached' a dedicated reference frame to each Particle that would be initialized. Hence, all of Particle's methods would call the relevant methods in MovingRefFrame.

I got to learn about a pretty neat method of implementing Singleton classes in Python. How? I saw the __new__ method being called in Prasoon's code, and after an hour of going from website to website, I somehow ended up learning this.

Have a look at the following code-

1:  class MutableSingletonClass(object):  
2:    instance = None  
3:      
4:    def __new__(cls, *args, **kwargs):  
5:      print ("In __new__")  
6:      instance = cls.instance  
7:      if instance is None:  
8:        print ("No instance found")  
9:        instance = cls.instance = \  
10:              object.__new__(cls)  
11:      else:  
12:        print ("Instance found")  
13:      return instance  
14:      
15:    def __init__(self, *args, **kwargs):  
16:      print ("In __init__")  
17:      pass  
18:    
19:  class Class1(MutableSingletonClass):  
20:    def __init__(self, x):  
21:      print ("In Class1 __init__")  
22:      self.x = x  

MutableSingletonClass's subclasses (and that class itself) will have only one instance being generated per session (the first time, that is). The next time you try to initialize a new instance of Class1, that same instance that was created earlier will be modified. Mind you, that one common instance WILL be modified everytime you call the constructor with new arguments. The reason this happens is, every time __new__ returns an object (which __new__ ensures is only one per session in this case), the __init__method is called. You cannot avoid that.
To get a better idea, look at the output-

1:  >>> c1 = Class1(3)  
2:  In __new__  
3:  No instance found  
4:  In Class1 __init__  
5:  >>> c2 = Class1(4)  
6:  In __new__  
7:  Instance found  
8:  In Class1 __init__  
9:  >>> c1.x  
10:  4  
11:  >>>c1 == c2  
12:  True  

Now, if we want to implement a real Singleton type, the Python docs show a nice method-

1:  class FixedSingletonClass(object):  
2:    instance = None  
3:      
4:    def __new__(cls, *args, **kwargs):  
5:      print ("In __new__")  
6:      instance = cls.instance  
7:      if instance is not None:  
8:        print ("Instance found")  
9:      else:  
10:        print ("No instance found")  
11:        instance = cls.instance = \  
12:              object.__new__(cls)  
13:        instance.init(*args, **kwargs)  
14:      return instance  
15:    
16:    def init(self, *args, **kwargs):  
17:      #This is where the subclasses should implement their  
18:      #initialization  
19:      pass  
20:      
21:    def __init__(self, *args, **kwargs):  
22:      #Subclasses should only define this as a dummy so users no which  
23:      #arguments to pass for initialization  
24:      print ("In __init__")  
25:      pass  
26:    
27:  class Class2(FixedSingletonClass):  
28:    def init(self, x):  
29:      print ("In Class2 init")  
30:      self.x = x  
31:    
32:    def __init__(self, x):  
33:      print("In Class2 __init__")  
34:      pass  

As you can see here, this implementation reduces the __init__ method to a mere formality/dummy (just to tell the users which values to enter for initialization, in the shell). Actually speaking, all the work is done by the 'init' method, which can now be controlled by the program. Hence, this real initializer (not __init__) is called only once..the first time an instance of Class2 is created. The rest of the times, __init__ is called, but it makes no changes to the existing instance. This ensures all references of Class2 point to a common instance.
Hence, the output-

1:  >>> d1 = Class2(5)  
2:  In __new__  
3:  No instance found  
4:  In Class2 init  
5:  In Class2 init  
6:  >>> d2 = Class2(10)  
7:  In __new__  
8:  Instance found  
9:  In Class2 init  
10:  >>> d1.x  
11:  5  
12:  >>> assert d1 is d2  

I hope this was useful.
That's all for now.
Have a great week ahead :-)

GSoC: Week 2 : Reference Frames

This is an update for the second week of my GSoC coding period. Its been long and quite tedious (though very interesting at times), as far as the code was concerned. Got a lot of stuff done in the commits this week, and I guess the major components of the reference frame class are now in place. Well, I can't be sure whether its perfect till Gilbert reviews it, but all that we discussed is more or less there. The things that got done this week-

• The initializer for the MovingRefFrame class. Its huge and dense as I had said earlier, but it gets the job done well, and prevents any 'special cases' from wreaking havoc while at the same time giving all the needed functionality. As mentioned in the last week, the functions for boundary conditions have now been implemented in a file of their own, and they get called in this class' 'init' function. There were some other issues as well, like allowing velocities to be defined in frames that don't 'exist' yet, etc. But thankfully, they got solved after some thinking(read: thinking hard).
• I rewrote the tree algorithm for finding relative motion parameters, since the one I had implemented earlier went via the global frame everytime- something that made the whole process very inefficient(and computationally impractical) at times. To get an idea of what I am talking about, imagine a 100 reference frames, each defined at some orientation/position wrt its predecessor. Now I want to know the motion of the 98th frame wrt the 96th. Earlier, I would find the motion of the 98th wrt the first(global) frame, do the same for 96th, and then 'process' the info to get the required motion params. Some tests with the current framework made me realise how stupid that was. So now, it just goes this way - 98th->97th->96th, the way its supposed to be.

A small mock SymPy session with the use of this class is shown below-

The conditions are-

1) A frame R1

2) A frame R2 defined wrt R1, such that it initially coincides with R1, but has angular velocity of k̂ wrt it

3) A frame R3 defined wrt R2, such that it has a constant position vector of  î  wrt it.

Since R2 is rotating wrt R1, R3's position vector wrt R1 will be a function of time, as expected. Heres the working-

 >>> from sympy.physics.mechanics import MovingRefFrame  
 >>> R1 = MovingRefFrame('R1', parentframe=None)  
 >>> R2 = MovingRefFrame('R2', parentframe = R1, ang_vel = R1.basis(2))  
 >>> R3 = MovingRefFrame('R3', parentframe = R2, pos_vector = R2.basis(0))
 >>> R1.express(R3.pos_vector_in(R1))  
 - cos(t)*R1.basis(0) - sin(t)*R1.basis(1)  

Thats all for now. You can have a look at the PR if you wish. The coming week will majorly be spent fine-tuning my current code and making it perfectly compatible with Prasoon's API, adding a few other helper functions if needed, and writing unit tests/example docs for all that's implemented so far.
Have a great week :-)

GSoC: Week 1

Hmm. Week 1 of GSoC ends today. Its been a bit of a confused week, with most of my time spent thinking and  trying to determine how I am going to code the initialization of reference frames in the new mechanics core. Its going to be a little dense, with all the motion parameters being checked, analysed and stored, and the constraint of immutability making things a little difficult, implementation-wise. But I guess its better to make the classes immutable to gel well with the rest of SymPy, as the whole point of this part of my GSoC work would be to make sympy.physics.mechanics seem 'friendlier' to the rest of the codebase.

What was gained? Well, I did finish (today) determining how exactly I am going to go the whole initialization procedure, after a lot of brain-racking (its quite true when they say that programmers do most of their thinking not at their desk, but while doing things like eating and showering :-) ). A lot of math was involved, along with considerable thinking given to the exact implementation details. What did I learn? - Always plan out what you are going to implement, on paper, before you start typing code(especially if math is involved). Taking time and understanding what you are going to do, and then doing it, is much better than constantly changing your code as you keep finding faults. Anyways, I did finish coding the rough API and some elementary functions like those of time derivatives, relative motion etc., but all this will be useful only once I finish the __init__ method.

One good thing we decided to do was to code helper functions in the a separate file in the core to calculate motion parameters given differential equations and boundary conditions- using SymPy's integration and equation-solving functions, ofcourse. I proposed this way of doing things since having this in the reference frame code would be confusing and chaotic, not to mention unnecessary. At the same time, a user would want to create a new frame instance based on a velocity vector and initial position boundary values-or something on those lines. So this we felt was a good middle ground to enable both. This is going to take time, but I guess having a good foundation is always good.

In terms of code and API, I will have something concrete in my PR by the end of this week. Will try and show a mock SymPy session in my next post :-). Anyways, thats all for now. I wont be able to do much tomorrow, since I have to travel to a relative's place for lunch..though I guess I could work at night for an hour or two. A long week ahead, wish me luck :-D

Smart caching in Python

Prasoon and I are currently looking for ways to reduce the computational load of calculating the DCM(Direction Cosine Matrice) of one frame with respect to another. This DCM is useful to mathematically understand how one reference frame is oriented with respect to another. The reason why we need to find a way is, if we are talking about something like a 100 frames each oriented in some random way with respect to one of its predecessors, it gets pretty messy- for your computer AND you.

We thought of various ways, but apart from some very inefficient methods like 'World' frames (to store these DCMs for a particular run of the program) and using a global frame as an intermediate, we really couldn't find anything much. Then Gilbert suggested looking at the way SymPy does caching of some of its methods.

We may or may not use this technique eventually, but I learnt a cool way to implement caching in Python-using dictionaries and method decorators. If you don't know what method decorators are, have a look at this - decorators are essentially ways to mutilate(or rather, add to) the functioning of a function/class.

I am going to try and show a simplified version of what I found in the SymPy source, their original method is close but much more detailed- for purposes of immutability and debugging and stuff.

The components involved in this technique are these - 

1) CACHE = {}

This is a Python dictionary (empty intially). This is what we will use as our cache.

2) The caching decorator -

def globalcached(function):

    @wraps(function)

    def wrapper(*args):

        tup = tuple([function] + [args])

        if tup in CACHE:

            print('In cache')

            return CACHE[tup]

        else:

            print('Not in cache')

            temp = function(*args)

            CACHE[tup] = temp

            return temp

    return wrapper

Understand what the above will do. wraps is a decorator from the functools Python library- its there by default. It enables the 'decorated' function to retain the docstring of the original function- we would probably want to retain the docs.

As you might have understood, 'wrapper' is the modified function that will get executed in place of the original one.

Now, tup is a tuple (tuple is immutable in python and hence hashable. You can only use hashable objects as dictionary keys) containing info about the following-

i. The function being called

ii. The arguments of the function

Thus, I know exactly what function it is, and what it is going to work on. Hence, the code will not directly let the original function go computing again- that may be a waste of resources. It first checks in CACHE, whether it has that particular tup as one its keys. If its not there, then this is the first time that function is being called with those arguments. So, I will calculate the needed quantity, but before returning it, I will store the calculated stuff with tup as the key. Hence, the next time the same function gets called with the very same arguments, the decorator will have a 'deja vu' and return the value- from CACHE. Ingenious, isnt it?

I don't know how many other projects use this method, or something similar. But its quite cool, it you want speed and you know the computation is going to be intensive.

Anyways, here is a snippet using the above decorator-

 CACHE = {}  
 def globalcached(function):  
   @wraps(function)  
   def wrapper(*args):  
     tup = tuple([function] + [args])  
     if tup in CACHE:  
       print('In cache')  
       return CACHE[tup]  
     else:  
       print('Not in cache')  
       temp = function(*args)  
       CACHE[tup] = temp  
       return temp  
   return wrapper  
 def printcache():  
   for x in CACHE:  
     print 'Function: ' + str(x[0])  
     print 'Arguments: ' + str(x[1:])  
     print 'Value: ' + str(CACHE[x])  
     print()  
 @globalcached  
 def myfunction(x, y):  
   """ myfunction docs"""  
   return x ** 2 - y ** 2  
 class SomeClass(object):  
   def __init__(self, x):  
     self.x = x  
   @globalcached  
   def function1(self, otherinstance):  
     """function1 docs"""  
     return otherinstance.x ** 2 - self.x ** 2  

Using it in a Python IDLE session-

 >>> myfunction(8, 3)  
 Not in cache  
 55  
 >>> myfunction(8, 3)  
 In cache  
 55  
 >>> i1 = SomeClass(4)  
 >>> i2 = SomeClass(7)  
 >>> i1.function1(i2)  
 Not in cache  
 33  
 >>> i1.function1(i2)  
 In cache  
 33  
 >>> myfunction(6, 7)  
 Not in cache  
 -13  
 >>> printcache()  
 Function: <function myfunction at 0x026B38B0>  
 Arguments: ((8, 3),)  
 Value: 55  
 ()  
 Function: <function function1 at 0x026CD330>  
 Arguments: ((<__main__.SomeClass object at 0x026C9B10>, <__main__.SomeClass object at 0x026C9B50>),)  
 Value: 33  
 ()  
 Function: <function myfunction at 0x026B38B0>  
 Arguments: ((6, 7),)  
 Value: -13  
 ()  
 >>>   

I hope its clear now. Its not rocket science, but its a clever way of doing things.

Anyways, that's all for now. And btw, thanks to Google for giving us a free membership to ACM for one year. The GSoC welcome package will also be arriving soon, and I am excited :-D

Some basic math

As suggested by Gilbert, I am uploading some hand-written notes I made about the entry-level math involved in coding framework Prasoon and I are talking of building. This DID help to some extent as a form of revision and cleared up some basics in my own brain.
I guess I should do this more often. I have practically lost the skill of writing things down since my tenth grade.
(By the way, this is a good article to refer to, if you are interested in knowing the conceptual differences between Reference Frames and Coordinate Systems)

sympy.logic - Part 2

(contd from last post)
Ok. A little late on this one. As promised, the part-2 about the features of sympy.logic, one of the lesser-used modules of sympy.logic. Today's post contains information about some stuff that was there before I worked on the module, and some stuff that we worked in. Anyways, here goes-

IV) Equality of boolean expressions

This is a function that became quite a headache for us while were trying to implement it. We thought of many algorithms and kept refuting them, until we settled on a good one. I coded it, and Christopher improved it (based on the SymPy architecture) considerably. It may still have a few bugs (one is right here), but it gets the job done nonetheless.

Lets take an example to see how it works
Take the function that I discussed in part I) of last post. Its SOP form was -
Or(And(Not(a), d), And(c, d))
and its POS form was-
And(Or(Not(a), c), d)
Lets take the POS form to be And(Or(Not(x), y), z), just to get a change of variables.

So, function1 = Or(And(Not(a), d), And(c, d)) and function2 = And(Or(Not(x), y), z)
According to the docstring of the bool_equal function, it "Returns True if the two expressions represent the same logical behavior for some correspondence between the variables of each (which may be different)"
In the above functions, the mapping a : x, c : y, d : z will give the SAME outputs for every input for the variables (a,c,d) or (x,y,z). How do we check this equality?

>>> from sympy.logic import *
>>> from sympy import symbols
>>> a, c, d, x, y, z = symbols('a c d x y z')      <- Remember all SymPy variables are called Symbols?
>>> function1 = Or(And(Not(a), d), And(c, d))
>>> function2 = And(Or(Not(x), y), z)
>>> bool_equal(function1, function2)

Output - True

Want to check the mapping?

>>> bool_equal(function1, function2, info = True)
(And(Or(Not(a), c), d), {c: y, a: x, d: z})

If there is more than one mapping which may make the expressions equal, any ONE is returned. The function that is returned as the first return is the simplified version of the functions (in terms of one set of variables). As the function is already the simplest here, it gets returned as it is.

V) Inferences from knowledge bases

This is a very handy feature of sympy.logic. It basically lets you 'tell it' a set of True boolean expressions, and then 'ask' it if a given input is True based on the 'knowledge' it has. First, let me introduce to two other logic classes in sympy.logic-

Implies

Implies(some expression, some expression)

Equivalent

Equivalent(some expression, some expression)

They behave exactly like the corresponding operators in logic theory. VERY useful when converting natural language statements to logic ones. They also let you convert them to And, Or and Not instances using simplify_logic or to_dnf or to_cnf functions.Ok, now about KBs-
The API is as follows-

>>> from sympy.logic.inference import PropKB
>>> from sympy.abc import x, y
>>> l = PropKB()
>>> l.tell(x & ~y)
>>> l.ask(x)
 True

Did you see what it did there? It returns True for any statement thats satisfiable given the statements it is 'told'.

One more thing.

One inefficient way to check if a statement y is valid given a set of statements x1, x2, ... xn is using the Implies class and simplify_logic.
Just use the following function -

def check_validity(statements, check):
      statements = simplify_logic(And(*statements))
output = simplify_logic(Implies(statements, check))
if output != True:
output = False
return statements, output

>>> check_validity([x|y], x)
(Or(x, y), False)

It will return the simplified version of the ANDing of all the statements in the first input, and the validity of the 'check' statement.

How did I use the above function in Database Systems? Put all the known functional dependencies in the 'statements' input and the dependency to be checked in the 'check' input (in the form of Implies objects). If the second output was not True, that would mean the functional dependency does not hold. Ofcourse there are better ways of doing this, but well, I had little time and I wanted to check some stuff pronto.

Anyways, thats all I have about sympy.logic. If you want to know more, the source is here. Do have a look.

GSoC update - 

The whole of last week was spent deciding the API of our work on frames, coordinate systems and vectors with Stefan, Prasoon and Gilbert. We have a consensus on some things now, and I am likely to start working this week on a new class to depict ScalarFields in n-dimensional space. I will post the progress by the end of next week, and I hope I will have substantial stuff to write about. Excited. Definitely. Cheers!

sympy.logic - Part 1

As I have mentioned before, my first work for SymPy was for its logic module. Most of it took inspiration from the digital design and CS-Logic classes I took in my 3rd semester at college. In this post and the next I will try to explain some of the cool features of this module that I have used applied in useful ways.

Note - Quite a bit of this work was done just a few months back by Christopher Smith(smichr) and me, and hence is NOT included in the lastest version of SymPy(0.7.2), but will find its way into SymPy 0.7.3. Therefore, there is no documentation on the official website also. However, interested people can fork the SymPy source from its github repo - There is enough of in-code docs to help you.

I) Using the truth table to find SOP/POS form 

This is quite useful in digital design, and I myself used this quite a few number of times during my course.
Suppose you have to find the SOP (sum of products) form of a boolean function that has the following specs-
1. The variables (literals) are a,b,c and d.
2. The function MUST return True for the following combinations (in order-a,b,c,d)-
0,0,0,1
0,0,1,1
0,1,1,1
1,0,1,1
1,1,1,1
 where 0-false and 1-true
3. The function may give True or False (a dont-care) for the following combinations-
0,0,0,0
0,0,1,0
0,1,0,1

For those of you who may ask WHY we won't care. Simple- These combinations will never occur in our application if used rightly, OR we dont' care about the output even if they do. For example, a function that deals with binary representations of digits from 0-9 will not care about the inputs from 1001, as they will not occur in operation.
4. All the rest of the combinations must return False

Python session for this? Here it is -

>>> from sympy.logic import SOPform
>>> minterms = [[0, 0, 0, 1], [0, 0, 1, 1], [0, 1, 1, 1], [1, 0, 1, 1], [1, 1, 1, 1]]
>>> dontcares = [[0, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 1]]
>>> SOPform(['a','b','c','d'], minterms, dontcares)

Output - Or(And(Not(a), d), And(c, d))

POS form?

>>> from sympy.logic import POSform
>>> POSform(['a','b','c','d'], minterms, dontcares)

Output - And(Or(Not(a), c), d)

II) Simplification of boolean expressions

Suppose you want to simplify a cumbersome boolean expression to something simple. How?

Use simplify_logic.

Example expression --> "AB(A + B)" [Symbols take their usual meaning]

On paper, you will have to take the following steps-

AB(A + B)
(A + B)(A + B)	DeMorgan's Law
A + BB	Distributive law. This step uses the fact that or distributes over and. It can look a bit strange since addition does not distribute over multiplication.
A	Complement, Identity.

SymPy console session?

>>> from sympy.logic import simplify_logic
>>> to_be_simplified = '~(A & B) & (~A | B)'
>>> simplify_logic(to_be_simplified)
Output - Not(A)

Period.

III) Conversion between cnf and dnf forms

Conversion of boolean expressions to CNF or DNF forms is also possible.

The code to do this is quite simple -

>>> from sympy.logic.boolalg import to_dnf

>>> from sympy.abc import A, B, C, D

>>> to_dnf(B & (A | C))

    

Output - Or(And(A, B), And(B, C))

sympy.abc generates 'symbols' - Variables that can be used to denote any value. The entire framework of SymPy is based on the concept of Symbols (hence, Symbolic-Py)

Thats all for now. In the next post, I will continue with some other features like boolean equality mapping and checking for true-ness of a statement given a set of statements (and using it to find redundancies in databases).

Hello World!

Aha. My first ever blog post. The inspiration for this blog came from the Google Summer of Code program, which will be my primary focus till September.

I came to know that I was selected only yesterday night, and what a day that was! Phew... I was going round in circles thinking about my chances, what I had done right, what I had done wrong, blah blah. Never has been a wait of 24 hours been so torturous. I have to give credit to Gilbert Gede, Sean Vig, Jason Moore, DL Peterson, Angadh and all the other SymPy physics community mentors who helped me shape my proposal to a presentable, well-thought out condition. Earlier it was just one haphazard list of things that I intended to do over the summer(mainly monsoon, due to the late arrival of GSoC), and my chances wouldn't have been half as good without these people giving their input time and again.

Gilbert Gede, from UC-Davis, will be my mentor throughout the whole GSoC period, and I hope to learn a lot from him....not just physics, but also coding and coding styles in general. He is one of the authors of sympy.physics.mechanics (PyDy).

My co-mentor will be Stefan Krastanov, who has done awesome work for SymPy during his GSoC project last year. Considering both my mentors are former GSoCers, their experiences with the program will go a long way in making sure I don't stumble around too much.

I will be working alongside Prasoon Shukla, a student of IIT-R, as his project is closely related to mine in terms of dealing with vector fields.

Some others I need to thank are Christopher Smith, who helped me immensely during my first PRs for SymPy, Aaron, the SymPy head, who took the pains to comment on all my stupid coding mistakes in my initial work, Manoj Kumar, Priyans Murarka, Jay Rambhia and Divyansh Khanna who were great supports during my application period. Thanks to all of you, I am penning my first ever blog post :-)

I hope I have a great 3 months coding for GSoC. I do realise it's not gonna be very easy, and I am bound to get stuck at times, but I hope to live up to my proposal(https://github.com/sympy/sympy/wiki/GSOC-2013-Application:-Sachin-Joglekar:-Electromagnetism-module) and finish well what I intend to do, to begin with.
The past few months have been very eventful and anxious at the same time, I hope I can work in peace now. 

A lot more blog posts are going to come, mostly about my GSoC work, but some also about the other interesting coding things I do here and there. Hope I live up to the expectations of the SymPy community :-).