I recently migrated my blog from Wordpress to Jekyll, a fantastic website generator that's designed for building minimal, static blog posts to be hosted on Github pages. Firstly, I wanted to use Ghost but then simplicity of Jekyll's theming layer and writing workflow had me.

Lets get the introduction out of our way first. I am Sahil Shekhawat, a sophomore majoring in Computer Science. I am Google Summer of Code intern at PyDy and have worked with Backpack for about one year. Thats pretty much it!

Once again I got accepted into Google Summer of Code! I'll be working on assumptions system of SymPy. This time, SymPy is participating under Python Software Foundation.

SymPy is a Python library for symbolic mathematics. It aims to become a full-featured Computer Algebra System while keeping the code as simple as possible in order to be comprehensible and easily extensible.

Here's what ideas page says about the project:

My mentors are Aaron Meurer and Tim Lahey.

Currently SymPy has two versions of mathematical assumptions. One is called "old assumptions" because a new implementation has been carried out recently. Since "old assumptions" were developed a long back, they are more mature and faster. However, because of its design, it is not capable of doing some interesting things like assuming something over an expression e.g. x**2 + 2 > 0.

Old assumptions store assumptions in the object itself. For example, the code x = Symbol('x', finite=True) will store the assumption that the x is finite in this object itself.

Both systems expose different APIs to query the facts:

Old:

In [1]: from sympy import *

In [2]: x = Symbol('x', imaginary=True)

In [3]: x.is_real
Out[3]: False

New:

In [4]: y = Symbol('y')

In [5]: ask(Q.real(y), Q.positive(y))
Out[5]: True

My work includes but is not limited to:

• Identifying inconsistencies between old and new assumptions and eliminate them.
• Improving performance of the new assumptions.
• Making new assumptions read old assumptions.
• Removing assumptions from the core as much as possible.
• Making API of old assumptions call new assumptions internally.

That's all for now. Looking forward to a great summer!

TL;DR I created a gallery for IPython/Jupyter notebooks. Check it out :-)

A couple of months ago I put online a website that displays a collection of IPython/Jupyter notebooks. The is a website that collects user-submitted and publicly available notebooks and displays them with a nice screenshot. The great thing about this website compared to other similar efforts is that this one updates and categorizes the notebooks (by date and views) automatically. You can even search for anything in this database that contains already more than 400 notebooks!

Vision

I would like to make a website where it is possible to

1. Find the notebook you are looking for. There is precious information (examples, documentation, tutorials, etc.) contained within IPython/Jupyter notebooks. It should be easy to find this information.

2. Discover new notebooks. I would like to see new notebooks as they are submitted in order to discover new notebooks, possibly about new technologies that I did not know about.

How it works

Under the hood there's a django app, for which the source code lives here (don't hesitate to use the issues feature in github to suggest features or to report bugs). To propose new notebooks there's a tab where anyone can leave the URL of a notebook. If all goes right, the django app will take a screenshot will be taken and incorporated into the database. The only non-trivial part of this process is to take the screenshot, for which I use a bit of javascript around phantomjs.

Future plans

What you see is just the beginning!. As time permits I would like to implement user authentication so that registered users can bookmark their favorite notebooks, up and down-vote notebooks etc. Cathegorization of notebooks (e.g. in bins such as Math, Phisics, Machine Learning, R, Python, Julia, etc.) is also high on my list. Leave me a comment if you would like to see some specific feature!

Last Friday was PyData Paris, in words of the organizers, ''a gathering of users and developers of data analysis tools in Python''.

The organizers did a great job in putting together and the event started already with a full room for Gael's keynote

My take-away message from the talks is that Python has grown in 5 years from a language marginally used in some research environments into one of the main languages for data science used both in research labs and industrial environment.

My personal highlights were (note that there were two parallel tracks)

• Ian Ozsvald's talk on Cleaning Confused Collections of Characters. Ian gave a very practical talk, full of real world examples. The slides have already been uploaded on his website. Many tips and many pointers to libraries. In particular, I discovered fixes text for you.

• Chloe-Agathe gave a short talk on DREAM challenges. In her talk she mentioned GPy. One year ago, I visited Neil Lawrence at his lab in Sheffield and at that point they were in the process of migrating their Matlab codebase into Python (the GPy project). I'm very glad to see that the project is succeeding and being used by other research institutions.

• Serge Guelton and Pierrick Brunet presented “Pythran: Static Compilation of Parallel Scientific Kernels”. From their own documentation: “Pythran is a python to c++ compiler for a subset of the python language. It takes a python module annotated with a few interface description and turns it into a native python module with the same interface, but (hopefully) faster”. The project seems promising although I do not have had experience as to judge the quality of their implementation.

• Antoine Pitrou presented: “Numba, a JIT compiler for fast numerical code”. I must say that I'm an avid user of Numba so of course I was looking forward to this talk. One thing I didn't know is that support for CUDA is being implemented into Numba via the @cuda.jit decorator. From their website it looks like this is only available in the Numba Pro version (not free).

• Kirill Smelkov presented wendelin.core, an approach to perform out-of-core computations with numpy. Slides can be found here.

• Finally, Frances Alted gave the final keynote on “New Trends In Storing And Analyzing Large Data Silos With Python”. Among the projects he mentioned, I found particularly interesting bcolz, his current main project and DyND, a Python wrapper around a multi-dimensional array library.

My latest research paper[¹] deals with the estimation of the hemodynamic response function (HRF) from fMRI data.

This is an important topic since the knowledge of a hemodynamic response function is what makes it possible to extract the brain activation maps that are used in most of the impressive applications of machine learning to fMRI, such as (but not limited to) the reconstruction of visual images from brain activity [²] [³] or the decoding of numbers [⁴].

Besides the more traditional paper that describes the method, I've put online the code I used for the experiments. The code at this stage is far from perfect but it should help in reproducing the results or improving the method. I've also put online an ipython notebook with the analysis of a small piece of data. I'm obviously glad to receive feedback/bug reports/patches for this code.

One of the lesser known features of the memory_profiler package is its ability to plot memory consumption as a function of time. This was implemented by my friend Philippe Gervais, previously a colleague at INRIA and now working for Google.

With this feature it is possible to generate very easily a plot of the memory consumption as a function of time. The result will be something like this:

where you can see the memory used (in the y-axis) as a function of time (x-axis). Furthermore, we have used two functions, test1 and test2, and it is possible to see with the colored brackets at what time do these functions start and finish.

This plot was generated with the following simple script:

import time

@profile
def test1():
    n = 10000
    a = [1] * n
    time.sleep(1)
    return a

@profile
def test2():
    n = 100000
    b = [1] * n
    time.sleep(1)
    return b

if __name__ == "__main__":
    test1()
    test2()

what happens here is that we have two functions, test1() and test2() in which we create two lists of different sizes (the one in test2 is bigger). We call time.sleep() for one second so that the function does not return too soon and so we have time to get reliable memory measurements.

The decorator @profile is optional and is useful so that memory_profiler knows when the function has been called so he can plot the brackets indicating that. If you don't put the decorator, the example will work just fine except that the brackets will not appear in your plot.

Suppose we have saved the script as test1.py. We run the script as

$ mprof run test1.py

$ mprof run test1.py
mprof: Sampling memory every 0.1s
running as a Python program...

$ mprof plot

Phew. Sorry about the late blogpost (and a missed one last week). Google Summer of Code is now officially over, and I am happy to report that I was successful in building a basic vector-manipulation/vector-calculus module for SymPy. The code is now in SymPy’s master, which means it will be released with SymPy 0.7.6. Its been a great 3 months, with intensive software designing, coding and bug-fixing going hand-in-hand to being out ‘sympy.vector’.

I would like to thank Jason Moore, my mentor during this year’s GSoC. Initially, I had begun on my usual path of the ‘lets start coding immediately’ philosophy, but he made me wait it out, figure out extensive use cases, and give stress on the software designing process. This was the first software written by me, with the test-driven development approach. And must say, it does benefit a lot. There were quite a few tricky places, but I am glad I got through it all :-). I couldn’t implement all I wanted to, but atleast I am happy that all I wrote is stable, and works perfectly :-D.

The docs are expected to go in soon, after which you (if interested) may read through them and understand how to use the new module. All in all, its been a great summer!

Cheers!

The GSoC program officially ended this Monday, and so my work for SymPy has concluded. I got a lot done in these last two weeks though. Here's what's new. In order of completion:

Complete overhaul of the codeprinting system

I wasn't happy with the way the codeprinters were done previously. There was a lot of redundant code throughout ccode, fcode and jscode (the main three printers). They also had a lot of special case code in the doprint method for handling multiline statements, which I felt could be better accomplished using the visitor pattern that is used by all the other printers. The issue is that some nodes need to know if they are part of a larger expression, or part of an assignment. For example, in C Piecewise are printed as if statements if they contain an assignment, or inline using the ternary operator if they don't.

After some thought, this was solved by adding an Assignment node that contains this information, and then dispatching to it in the printer just like any other node. Less special case code, and allowed the base CodePrinter class to contain a lot of the redundancies. For those implementing a new code printer (perhaps for Octave?) all you'd need to do is add how to print certain operators, and a dictionary of function translations. Everything else should just work. I may add little cleanups here and there, but I'm pretty happy with the refactor.

Code printers now support matrices

This was the original goal, but got put aside to do the previously described refactor. The codeprinters now support matrices - both as inputs and outputs. For example, the following now works:

# Expressions inside a matrix
x, y, z = symbols('x, y, z')
mat = Matrix([x*y, Piecewise((2 + x, y>0), (y, True)), sin(z)])
A = MatrixSymbol('A', 3, 1)
print(ccode(mat, A))

A[0] = x*y;
if (y > 0) {
   A[1] = x + 2;
}
else {
   A[1] = y;
}
A[2] = sin(z);

# Matrix elements inside expressions
expr = Piecewise((2*A[2, 0], x > 0), (A[2, 0], True)) + sin(A[1, 0]) + A[0, 0]
print(ccode(expr))

((x > 0) ? (
    2*A[2]
)
: (
    A[2]
)) + sin(A[1]) + A[0]

# Matrix elemnts in expressions inside a matrix
q = MatrixSymbol('q', 5, 1)
M = MatrixSymbol('M', 3, 3)
m = Matrix([[sin(q[1,0]), 0, cos(q[2,0])],
    [q[1,0] + q[2,0], q[3, 0], 5],
    [2*q[4, 0]/q[1,0], sqrt(q[0,0]) + 4, 0]])
print(ccode(m, M))

M[0] = sin(q[1]);
M[1] = 0;
M[2] = cos(q[2]);
M[3] = q[1] + q[2];
M[4] = q[3];
M[5] = 5;
M[6] = 2*q[4]*1.0/q[1];
M[7] = 4 + sqrt(q[0]);
M[8] = 0;

There even was a Piecewise inside a Matrix in there. As long as there is an assignment between two compatible types (matrix -> matrix, scalar -> scalar), the new codeprinters should print out valid expressions.

codegen now supports matrices

This is more of a continuation of the above. The code generators have been modified to recognize instances of MatrixSymbol as array variables and act accordingly. There actually wasn't that much to change here to make this work. The biggest change that happened is that all C functions that have a return value (non void functions) allocate a local variable of the same type. This is to cover a larger set of expressions, while still generating valid code. So now, when performing codegen on "sin(x)" you won't get "return sin(x)", you'll get:

result = codegen(('sin_c', sin(x)), "C", "file", header=False)
print(result)

double sin_c(double x) {

   double sin_c_result;
   sin_c_result = sin(x);
   return sin_c_result;

}

This isn't as pretty, but handling return inside expressions is a tricky problem, and this solves it without much work. Modern compilers should remove the variable assignment if it's unnecessary, so there shouldn't be a resulting speed loss in the code.

Cython wrapper for autowrap now works

There was a code wrapper for Cython in the codebase, but I don't think it has ever worked. It does now:) It can do everything f2py can do, and I plan on adding more useful features. In it's current state it can:

• Handle both scalar and matrix input, input-output and output arguments
• Internally allocate output arguments
• Pull inferred variables (such as matrix dimensions) out of the function signature
• Create a multiple return value tuple

The last thing I want to do to make this really nice is to add support for informative docstrings. Even so, this is already usable:

x, y, z = symbols('x, y, z')
mat = Matrix([x*y, Piecewise((2 + x, y>0), (y, True)), sin(z)])
func = autowrap(mat, 'c', 'cython')
func(1, 2, 3)

array([[ 2.        ],
       [ 3.        ],
       [ 0.14112001]])

For some reason the Fortran/f2py has around a 2 microseconds faster than the C/Cython code. I think this has something to do with array allocations, but I'm not sure. For larger expressions they're pretty equal, so this shouldn't be that big of a deal. I still plan to look into code optimizations I could make in the Cython wrapper.

Project Status

Overall, I accomplished most of what I set out to do this summer. Some things (pre-solution linearization) were nixed from the project due to changing goals. Here's a short list of what was done:

1. General linearization methods added for both KanesMethod and LagrangesMethod.

2. Code cleanup and speedup for KanesMethod and LagrangesMethod.

3. Creation of msubs - a specialized subs function for mechanics expressions. This runs significantly faster than the default subs, while adding some niceities (selective simplification).

4. Complete overhaul of codeprinters. Fixed a lot of bugs.

5. Addition of support for matrices in code printers, code generators, and autowrap.

6. Overhaul of Cython codewrapper. It works now, and does some nice things to make the wrapped functions more pythonic.

7. Documentation for the above.

The Future

I had an excellent summer working for SymPy, and I plan on continuing to contribute. I have some code for discretization that I've been using for my research that may be of interest to the mechanics group. I also want to get common sub-expression elimination added to the code generators, as this kind of optimization may result in speedups for the large expressions we see in mechanics. My contributions will unfortunately be less frequent, as I need to really focus on research and finishing my degree, but I still hope to help out.

I plan on writing another post in the next few days about the GSoC experience as a whole, so I won't touch on that here. Let me just say thank you to Jason, Luke, Oliver, Sachin, Tarun, Aaron, and all the other wonderful people that have offered me guidance and support throughout the summer. You guys are awesome.

Hi All,

During this week I completed more functionalities regarding to CSR Matrices. I mainly implemented a generalized binary operation between two matrices. It’s similar to the SciPy function csr_binop_csr_canonical. This method applies a given binary function, to each corresponding elements of the two matrices which are given as input arguments, i.e. entry of the output matrix is equal to where and are the entries of input matrices respectively.

Also, I implemented methods for CSRMatrix class. It took me while to get this completely right as I was discarding some of the possible cases.

I originally planned to start writing python wrappers but I couldn’t do it as I was a bit busy with a project in my university. But I hope to get started on this in upcoming week. Although the GSoC is officially over, I will work on the project for two additional weeks as I took a two weeks break during the programme.

Hello everyone,

This GSoC post has been long pending, but to be honest the last two weeks have been simply been so hectic that I haven’t managed to find the time to write a lot of code, or blog about it. One of the major things I have worked on is greatly improving the resolution code I wrote earlier. This revision fixes a small bug and greatly improves the speed of the resolution process. Initially if a clause was not entailed by the KB, then it took almost forever for the method to return False. In the scenario where it was entailed, it was still performing lots of useless computations. This has been completely removed from the new resolution, and the method quite efficiently gives the correct value. The other major task has been to fully integrate the Constants class. Now, any non-boolean object that will internally get converted to a Constant and the user can be simply put in strings and numbers and the system will still play nice with it. Except for this, I have finished allSAT (all models for a given formula). I will be sending out a PR for this as soon as I finish testing it properly. I think that is all for now.

Bye!!!

Hi All,

I am writing a blog post for GSoC after three weeks. I couldn’t work on the project for two weeks as I participated in a competition. However I am back again and started working again. Now, I am mainly concerned about finishing my work with Sparse Matrices. Then we can write some wrappers and do some benchmarks.

We decided to use SciPy’s sparse matrices code and I have ported most of the important functionalities. I created a CSRMatrix class and created a method to create a CSRMatrix from a normal COO matrix. I ported the methods for testing the canonical form and sparse matrix multiplication.  Also, I implemented three methods mainly to be used in extracting diagonal elements and scaling.

I am going to implement binary operations with CSR matrices during the coming week and also want to get started with writing wrappers so we can use the Matrix classes from python and SymPy.

The GSoC is about to end. We have the suggested deadline tomorrow and the firm deadline next week. We are wrapping up the project now, cleaning up the code, writing the documentation and turning the project into a usable state. Since the new solvers cannot be called compelete I have to make sure that the next person working on it can pick it up easily. I plan come up with a post with the summary of the work that was done, that challenges we faced and thing that needs to done and how can they be tackled.

Phew, that was close I just had full week of written tests by recruiting companies. Well, not good for me. I suck at electrical and that too when we have to study ancient electrical subjects which has no longer any use in this modernized world. My college, maybe all Indian colleges, lag behind in electrical subjects. Though I always find online courses of electrical in U.S colleges fascinating, Oh but I am in India, can’t help with it. Let’s see how good I will do in GRE and TOEFL, maybe after getting in better colleges I will get to discover more.

Oh, I got driven by emotions, this is my GSoC blog so I will start with concluding my work I have done till now. Well before I conclude I will admit I wasn’t doing much work from last few days. In my 10th and 11th week I just had my checksysodesol PR merged. And two PRs for doc and nonlinear one are still behind. In the beginning of next week I will complete that and get it merged. This week I started the doc PR to complete all documentation and will add few tests to make coverage_report of ode to increase which is 85% now. I will discuss on nonlinear one with Sean and Tim about what modification is required.

My work in GSoC is on system of ODEs. In this I included linear system of ODEs for first order of two, three and `n` equations and second order of two and three equations and non-linear system of ODEs for first order of two and three equations. In order to solve these I modified `dsolve` function to incorporate solving these equations, made a new function `classify_sysode` to classify these system of ODEs and also made a function `checksysodesol` to check for solutions by substituting those solutions in equations and simplifying to get zero.

As this will be the last official blog for GSoC, as Aaron mailed, I would like to say few things here, offcourse we will be on gitter. This is the best experience of work I had in life. Well working continuously without having any care of anything else, it just feels awesome. I learned a lots and lots of stuffs from Sean, Tim, Aaron, Smith, Ondrej and some from Harsh and Sudhanshu too. I just started working on Sympy from 27 Feb, that was too late :), and then I didn’t knew git, github, opensource and many more things, just knew continuously working hard just pave the way to success but not without you guys. Open source is a great place to learn new stuffs and even better due to interaction with lots of people from other places. I will always be contributing to open source. It has lot of free software which otherwise would have to hard to buy in college life. This is going to my most memorable summer.

Thanks you guys!

Jim brought my attention to a very crucial bug in my work this week. Well, SymPy now automatically detects folders to test by looking at the __init__.py files. That is, a given folder in the SymPy code directory would _only_ be tested if it has __init__.py. Since I hadn’t put any such file in _sympy/vector/tests_, they were being skipped till now! Thank God this got fixed.

I started out this week with the seemingly ‘easy’ task of implementing different kinds of printing functionalities for sympy.vector. Well, the basic methodology of doing things is pretty simple, if you just follow the documentation in the _sympy.printing.printer_ file. It will tell you all you need to know about – either implementing a new printer, or adding custom code for printing of your class (to an existing printer).

Some things I had to improve/word around-

1. BaseScalars had to be handled differently from BaseVectors, since BaseScalars are a part of other Exprs (that may be measure numbers in vectorial expressions), and they have their own dedicated printing methods. Hence, I had to implement the _latex and _pretty methods in BaseScalar itself.

2. The pretty _as well as_ latex printing for _both_ Vectors and Dyadics with a single _print_BasisDependent method in the respective printers.

3. I used the printing code in _sympy.physics.vector_ as reference, though quite a few things had to be changed keeping in mind the inheritance from Expr, the structure and nature of the args being different, etc.
I did improve on a few things. _physics.vector_’s pretty printing messes up when dealing with vectors whose measure numbers have varying ‘heights’.

Overall, it was quite a learning experience and tricky coding – to implement printing, compared to my initial impression of it being a boring and tedious job.

Well, I have implemented all the requried methods and tests (all of which pass). I have pushed them to my third GSoC PR, so thats in line. Hopefully, the dyadics PR will soon show all tests passing, after which Jason or I will get it merged. It will add the OOP structure, and the dyadic code, as well as fix some crucial bugs.

A demo of the new printing implementations for a tedious vector expression-

Assume N is a CoordSysCartesian instance. Then,
(a**2 + b)*N.i + (Integral(f(b)))*N.k pretty-prints as


u'\u239b 2 \u239e N_i + \u239b\u2320 \u239e N_k\n\u239da + b\u23a0 \u239c\u23ae f(b) db\u239f \n \u239d\u2321 \u23a0 '


and LaTeX-prints as-


(a^{2} + b)\mathbf{\hat{i}_{N}} + (\int f{\left (b \right )}\, db)\mathbf{\hat{k}_{N}}


I missed my post last week due to my research for grad school all going to hell at the worst possible time :(. There wasn't much to report on last week, so I'm not too perturbed. Turns out even with PRs done and tested, it still takes a long time to get them merged and have everyone agree on them. Fortunately I made up for it this week; as of now I am officially working on code generation! This is behind schedule on my original timeline, but should have been expected, as per Hofstadter's law:

Hofstadter's Law: "It always takes longer than you expect, even when you take into account Hofstadter's Law."

Things are moving ahead now, and I have some hope left that I can accomplish (most of) everything I set out to do.

Project Status:

I set out to accomplish a number of things this summer. Here's the status of each project goal:

• General Linearization Form:

  Done!

• Linearization methods for Kane's and Lagrange's Methods:

  Done! Lagrange stuff got merged last week.

• Documentation for the above:

  Done? This PR is waiting for review, but I think it's good.

• Pre-linearized solution methods:

  Nixed from the plan due to lack of interest.

• Code Generation:

  In Progress...

I also accomplished a few other things, that I found necessary for my project:

• Refactored KanesMethod and LagrangesMethod:

  This resulted in more readable, pythonic code, and also a speed/memory improvement. PR for this is still awaiting review.

• Faster subs function for mechanics:

  My msubs implementation is several orders of magnitude faster than subs for the kinds of expressions seen in mechanics (even without the benefit of cacheing that subs has). This is in a PR awaiting review.

• Soooo much work on the nan and oo issue:

  Still not solved for all cases... :(

TODO:

There are only 3 weeks left of GSoC! In my last remaining weeks, here's what I plan to get done:

• Get the Bicycle Example working:

  After my work everything runs faster, results in smaller, more tractable expressions, and uses less memory. Except for the bicycle example. For some unknown reason I can not get this thing to result in anything except nan. This is a regress in performance (even though everything else runs better), and needs to be solved.

• Code generation:

  I've already got some stuff working, and it's really exciting. More on this below.

• Get all my current stuff merged:

  All that works needs to get into Sympy. As not everyone else is being paid to do this, it can take some time and effort to get things through the review process and into master, but I have hope that my remaining contributions will eventually make it in.

I think I can do it, but it'll be a stretch.

Code Generation

Sympy currently contains some facilities for code generation, but they lack support for the matrices that are necessary for working with dynamics problems. I hope to remedy that, as well as to make general improvements to the entire codegen module.

Code generation in sympy has three levels:

1. Code Printers ccode, fcode, and the like

   These are printers that know how to print simple sympy expressions using functionality and syntax found in that language. For example, ccode will print exponents using pow, which is found in the math library in C. These printers don't have any knowledge of functions, multiple statements, or header files. They simply print a single expression out on one line.

2. The codegen submodule

   This submodule contains facilities for representing generalized routines, and generating functions in various languages (currently C and FORTRAN) that can be compiled as a library without any changes. They know about function and variable declarations, header files, library imports, and multi-line statements. However, they have no idea how to make this generated code work with python.

3. Code wrapping, usually accessed through autowrap

   This is where the functionality for wrapping the generated code lives. Using the functionality provided here, one can compile and wrap generated code, and then call it from python. The autowrap function is the main entry point, allowing for all 3 steps to be done in one call.

The first thing I wanted to fix was getting code generation to work with matrices and matrix expressions. This turned out to be harder (and more confusing) than I expected. There is currently support for a "matrix like" object named sympy.tensor.IndexedBase, but I really don't understand the purpose behind it. Reading through the code and examples though it seems to be for representing indexed loop operations in a concise form. This unfortunately has nothing to do with the indexed types (matrices) that I plan on implementing.

I spent a long time reading through the code and playing around with it using pdb trying to figure out the control flow in the codegen function, and am still a little lost. Most of what's there seems to be for supporting the Indexed operations. After some time trying to bend them to work for matrices, I changed plans and now am supporting Matrix and MatrixExpr types for matrix operations only. Indexed types can be used elsewhere, but they shouldn't be used for representing matrices with expressions inside them.

I currently have this "working", but am not happy with it yet. The current layout of the module made for some hacky work adding in matrix support. I plan on doing some refactoring to make this implementation cleaner. Currently, on my codegen branch the following is supported:

• Generating C code for a matrix with expressions in each element:

  Matrices are set as input-output type arguments, and are modified in place before being returned.

• Passing in a MatrixSymbol as an argument:

  Here the plan is to use matrices to pass in a large number of arguments. You can think of this kind of like a vector. There's another symbolic vector type as well in Sympy (DeferredVector). I may end up supporting it, but I'm not really sure what it's for. In its current implementation, the following works:

  q = [q1, q2, q3, q4, q5]
  q_vec = MatrixSymbol('q', 5, 1)
  sub_dict = dict(zip(q, q_vec))
  # Replaces each q with elements from q_vec
  expr = msubs(expr, q_vec)
  # Generate a function that takes a numpy array (q) and returns expr
  # This works if expr is an expression, or a matrix
  func = autowrap(expr)
  

After I clean this up, I plan to add support for:

• Common Subexpression Elimination (cse):

  Even though modern compilers do this already, experimentation shows that the large expressions generated in mechanics benefit from generated code having cse performed. This will be implemented as a boolean kwarg (default False). When True, sympy's cse function will be run on the expression, and the code for each subexpression will be generated, followed by the final expression. I actually don't think this will be too difficult to implement, and should give some speed improvements on the compiled code (at the cost of slower generation).

• A ctypes code-wrapper:

  Currently the only code wrappers supported are f2py and cython, neither of which is in the standard library. While the wrappers generated with those functions may be more robust, a ctypes wrapper is also possible, with the added benefit that ctypes is in the standard lib.

• Support for matrix expressions:

  In an ideal world, I'd implement the excellent work done by Matthew Rocklin, discussed in this video from SciPy 2013. The idea here is that we have some knowledge about each of the matrices involved in an expression (for example $A^{-1} B$). We may know that A is positive definite, or symmetric, or upper triangular, etc... For each case, there may be a faster inversion routine that we could take advantage of rather than using a one-size-fits-all inverse function. As I don't have time to implement support for all possible operations, and the many BLAS/LAPACK routines that support them, I'll focus just on the inverse, as it's commonly found in expressions in mechanics. The thought is, we should be able to run:

  func = autowrap(Inverse(M)*F)
  

  And have code generated that solves the expression in a fast manner, without having to symbolically find the inverse of M and combine it with F into one matrix beforehand.

Of course this is a wishlist, and it's unlikely all of this will be accomplished in the next 3 weeks. Still, I plan to keep supporting sympy after my GSoC ends, so if it's not done by then it will eventually get there.

Other exciting news of the week:

I got accepted to the GSoC reunion at the end of October! As this is the 10th annual GSoC, Google is throwing a big reunion shindig for past and present students. As there are lots of us, only a few were chosen based on a lottery, and I made it through! I'm very excited to meet other students that completed the program, listen to some interesting talks, and see the GooglePlex. I also bought my tickets to get there a day early so I have some time to explore the bay area. Last time I was out there I was 14, and I didn't get to see much of the area. If you also got accepted/live out and would be interested in meeting up, let me know! I'll be in San Jose/San Francisco October 23-26.

This week I worked on the solvers for the equations with radicals. Suppose you have to solve

$$ \sqrt{x} + x - 2 = 0 $$.

then you have to move the radical to the right hand side of the equation.

$$ x - 2 = - \sqrt{x} $$

and then square at both sides

$$ x^2 - 4x + 4 = x $$

Now the equation is a polynomial in \( x \) can be solved with usual polynomial solving methods. Note that squaring both sides produce some extra solutions and we will have to check all the solutions obtained against the original equation. If there are more than one radicals involved we may have to apply the method recursively. For example in solving \( \sqrt{2x + 9} - \sqrt{x + 1} - \sqrt{x + 4} = 0 \) the method will recurse twice.

To implement the method I tried a pattern matching approach. The squaring part is easy the tricky part is identifying which part to move to the right hand side. First I tried to match the expression with the form sqrt(p) + q but it failed even for case like 4*sqrt(x) + x - 2 because no pattern matched to it. I had to use a*sqrt(p) + q with the condition that the expression matched to a shouldn't be zero. Now I can simply move the expression matched with p and terms multiplicated with it to the RHS and square both the sides.

Notice that this method for solving sqrt equation can work with any radical equation, if it were cube root instead of sqrt I just had to cube both the sides. OK so how do I mathch that expression? I tried to pattern matching with assumptions on the wild symbols but it doesn't work. I tried to match with somthing like a*p**Rational(1, m) + q but this also didn't work out because Rational(1, m) raises TypeError no matter what the assumption on the variable are. There is a proposal for a new pattern matcher, I have not closely checked the details but it will be able to work with assumption. You can see the proposal on the wiki here and if it is implemented then things will be good but I can't wait for it. I had no other option to check term by term for rational power. Here's the implementation

def _has_rational_power(expr):
    a, p, q = Wild('a'), Wild('p'), Wild('q')
    pattern_match = expr.match(a*p**q)
    if pattern_match is None or pattern_match[a] is S.Zero:
        return (False, None)
    elif isinstance(pattern_match[q], Rational):
        if not pattern_match[q].q == S.One:
            return (True, pattern_match[q].q)

    if not isinstance(pattern_match[a], Pow) or isinstance(pattern_match[a], Mul):
        return (False, None)
    else:
        return _has_rational_power(pattern_match[a])

This was a tiring week, with me being sick with a bad case of sinus for 3 days. Phew. Even then, I did manage to get quite a bit of work done in the time I could manage to sit in front of the computer without getting a headache.
Some main points of progress –

1. I finished off the Orienter classes, with a nice OO-structure, along with docs and tests. I basically changed the API of the orient_new

I am finally in college, going through a long journey from home. Its final year for me now, got a very tight schedule.

I have updated the checksysodesol with test cases. The test cases includes the equations from first and second order linear and non-linear system of ODEs. I debugged few methods of system of ODEs, find out through checksysodesol.

The updated func method PR is merged which closes issues #7736 and #7723. The PR of non-linear of system of ODEs and PR of checksysodesol is open.

Today I will add all the test cases in checksysodesol and start a PR for non-linear system of ODEs of 2nd order.

Hello Everyone,

I finished up the changes remaining in FOL and will send in the almost final pull request this week. This was only a minor part of the week however. Before jumping on to the implementation of faster SAT, I wanted to revisit faster CNF conversion. I have tweaked things further, but it still seems like very close to CNF turns out to be quite slow. If I can’t imporove the situation any further I will post the changes I made to Tseitin. Except fo that, I worked on tweaking the parameters of SAT, basically EVSIDS, without making any actual code changes. I will be posting the results of the same later.

Thats all for now. Hopefully this weekend will be quite productive.

See ya all later!!!

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.

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 where is a column matrix. I had to enhance the algorithm so that it can be used to solve several systems at once, i.e. can be a collection of column matrix.

This was not a problem in fraction free LU method because and factors can be used to solve column vectors after calculating and 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.

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.

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!

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 worst¹. 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 documentation². 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.

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

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.

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.

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!

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.

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 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 . 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

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!!!

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

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 \).

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.

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 :-).