The 8^th week has ended and we are now in the middle of phase –III.

Last week was a bit of research-based, understanding matplotlib and how it can be used to plot a beam diagram. I had a conversation with Jason Moore and Jashan where Jason shared a link of a repository, which also was a bit of help as I took some hints from it for the draw() function. After a lot of investigation and analysis, I was finally able to make a draft PR # 17240 which did the work as we intended.

Here is an example of how it would function:

# This example has no prior logic involved. It just tests whether every functionality works or not

>>> E, I = symbols('E, I')

>>> b1 = Beam(50, E, I)

>>> b1.apply_load(-10, 0, -1)

>>> b1.apply_load(R1, 10, -1)

>>> b1.apply_load(R2, 30, -1)

>>> b1.apply_load(9, 5, 0, 23)

>>> b1.apply_load(9, 30, 1, 50)

>>> b1.apply_support(50, "pin")

>>> b1.apply_support(0, "fixed")

>>> b1.apply_support(20, "roller")

>>> b1.draw()

Here we are using matplotlib and numpy by importing them as external modules. Of course, it would be better to have it done via SymPy’s own plot(), but I think that is something we could work on in later stages as SymPy’s plot() is limited to work on equations and stuff (although on can use _backend attribute for further functionalities).  Also to be noted here that SymPy’s plot() is not a replica of matplotib’s plot() but it makes it easier for SymPy equation to be plotted and it uses matplotlib to do so.

Following are the matplotlib modules/classes used:

• matplotlib.patches.Rectangle -to draw the beam
• matplotlib.pyplot.annotate – to draw arrows of load
• matplotlib.markers– to draw supports

Also, considering Jason’s comment in the PR, I will have to work on making SymPy’s plot() to accept a singularity function, so that it would be easier to plot loads which are indeed equations of Singularity function. This is still in consideration, so I will have to look into it and of course will have a discussion on how it is to be done.

Currently, I am not able to determine how to plot parabolic loads. I think this could be added later as we should currently focus on plotting simple parts and certainly work on other complexities later. But we can have a discussion on it.

Other PR’s are still being parallelly worked on.

Next Week:

• Working on the idea of plotting singularity function via SymPy’s plot()
• Plotting parabolic loads
• Writing documentation and tests

Will keep you updated!

Thanks!

Phase-II has been completed and now it’s time to present all the work done during this phase. This week’s blog is little early in comparison to my previous blogs because I’ll not be active for next 2 upcoming days. In the whole phase we worked on Computations with Polycyclic Groups though the tasks mentioned in proposal for this phase were quite different. But let me tell you it worth that much time, Computation with Polycyclic groups(solvable groups) shows the actual development in computational group theory.

Below are the functioalities that were added to the polycyclic group in this phase, Here is the PR sympy/sympy#16991.

Testing Collector

As discussed in week-5 blog at the time of Collector implementation we did not have the implementation of polycyclic presentation so some hand made tests were used. Here is an example of S(4).

>>> from sympy.combinatorics import *
>>> from sympy.combinatorics.free_groups import free_group
>>> F, x0, x1, x2, x3 = free_group("x0, x1, x2, x3")
>>> pc_relators = { x0**2: (), x1**3: (), x2**2: (), x3**2: (),
...                 x0**-1*x1*x0: x1**2, x0**-1*x2*x0: x2*x3,
...                 x0**-1*x3*x0: x3, x1**-1*x2*x1: x3,
...                 x1**-1*x3*x1: x2*x3, x2**-1*x3*x2: x3
...               }
>>> word = x3*x2*x1*x0
>>> relative_order = [2, 3, 2, 2]
>>> group = word.group
>>> collector = Collector(pc_relators, relative_order, group)
>>> collector.collected_word(word)
x0*x1**2*x2*x3

The final word x0*x1**2*x2*x3 is said to be collected. For more information about collected word please look into the docstrings of the method Collector.collected_word().

Computation of Polycyclic Sequence and Series

Polycyclic sequence and series are the building blocks of polycyclic presentation (have a look at week-6 blog) . One thing to note is that, the derived series of a group may change on different course of execution so we may have different pc sequence and series for the same group.

>>> from sympy.combinatorics import *
>>> G = SymmetricGroup(4)
>>> PcGroup = G.polycyclic_group()
>>> PcGroup.pcgs
[Permutation(0, 1, 2, 3), Permutation(3)(0, 2, 1), Permutation(0, 3)(1, 2), Permutation(0, 1)(2, 3)]
>>> 
>>> PcGroup.pc_series[0] == G
True
>>> PcGroup.pc_series[1] == AlternatingGroup(4)
True
>>> PcGroup.relative_order()
[2, 3, 2, 2]

Computation of Polycyclic Presentation

Few approaches were used to compute polycyclic presentation, the current implementation is discussed in week-7 blog. Below is a small example to show the functionality.

>>> from sympy.combinatorics import *
>>> G = SymmetricGroup(4)
>>> PcGroup = G.polycyclic_group()
>>> from sympy.combinatorics.free_groups import free_group
>>> len(PcGroup.pcgs)
4
>>> free_group, x0, x1, x2, x3 = free_group("x0, x1, x2, x3")
>>> PcGroup.pc_presentation(free_group)
{x3**2: (), x2**2: (), x2**-1*x3*x2: x3, x1**3: (), x1**-1*x3*x1: x2*x3, x1**-1*x2*x1: x3, x0**2: x2*x3, x0**-1*x3*x0: x2, x0**-1*x2*x0: x3, x0**-1*x1*x0: x1**2*x3}

As I mentioned above pc_sequence and pc_series may change on different course of execution and hence the pc_presentation changes accordingly.

Testing Presentation

Due to the changing pc_presentation initially, it was difficult to test presentation but later on a method has been developed and a good amount of code is introduced to test the presentation. The details can be found in the module test_pc_groups.py in the above PR.

Additional methods for Polycyclic groups

There were few additional methods added to the polycyclic group.

• exponent_vector() :- It represents the given generator of a polycyclic group with the help of product of pcgs.
• depth() :- Depth of the first non-zero element in exponent_vector.
• leading_exponent() :- It represents the power of polycyclic generator at the above depth.

>>> from sympy.combinatorics import *
>>> from sympy.combinatorics.free_groups import free_group
>>> G = SymmetricGroup(4)
>>> PcGroup = G.polycyclic_group()
>>> pcgs = PcGroup.pcgs
>>> len(pcgs)
4
>>> free_group, x0, x1, x2, x3 = free_group("x0, x1, x2, x3")
>>> PcGroup.exponent_vector(G[1], F)
[1, 1, 1, 1]
>>> G[1] == pcgs[0]*pcgs[1]*pcgs[2]*pcgs[3]
True
>>> PcGroup.depth(G[1], free_group) == 1
>>> PcGroup.leading_exponent(G[1], free_group) == 1

Currently, we are discussing about organizing above methods of polycyclic group. As Kalevi feels that the suitable place for pc_presentation(currently, it’s a method of PolycyclicGroup class) is the Collector class, Perhaps the structure of both the classes should be changed and the same will be reflected in the examples mentioned above.

This was the eigth week meeting with the GSoC mentors which was scheduled on Saturday 20th July, 2019 between 12:30 - 1:30 PM (IST). Me, Yathartha and Amit were the attendees of the meeting. Finally the PR #16890 is merged to Sympy master and the work for Lambert has been completed. It was only possible because of the mentors and especially @smichr for his great suggestions in code.

• Discussing previous week’s progress

In this meeting the goals for the time left for GSoC were decided as follows:-

1. Lambert:- Completed and merged.

2. Solve modular:- What I experienced with lambert’s PR that got stretched so much which I was not expecting. Algorithms made from heuristics takes so much time to be ready for merging so I really don’t know how much time it will take. But It is sure that it will be completed before the final evaluation.

3. ImageSet Union:- This task will be taken after the solve modular. This is a very complex task and will need proper guidance and algorithm. I had searched for some research papers, but what I found out was not that what we want. Before GSoC final evaluation this task will be started to implement, but I am not really sure if it would get merged before final evaluation.

• Next week goals

• Work upon _solve_modular PR

• If time left then find plan for Imageset Union.

Code improvement takes time!!

“Those who plan do better than those who do not plan even though they rarely stick to their plan.” ― Winston Churchill Welcome everyone, this is your host Nikhil Maan aka Sc0rpi0n101 and this week I have some good news for you all. We have got a few plans to get the progress of the project rolling as the second evaluation approaches. The meeting also took place this week, so, after waiting for so long, we have some good insights from the meeting this time.

Week 7 ends.. - Phase 2 of the coding period is smoothly being traversed. I recently had a meeting with Sartaj on 11th of July, Thursday. Here were the minutes of the meeting, along with the deliverables completed over the week. The FiniteFormalPowerSeries class PR needs some changes. Currently, it is taking in an...

The week was successfully completed as planned. The work on Column class has been completed.

The documentation and tests have been written and with some changes in the solve_slope_deflection() and critical_load(), the Column class is now able to handle cases with trivial solutions of the constants ( C1 & C2) which made the deflection equation zero.

Apart from this, another problem that we had with the pinned-fixed end condition, where solve() wasn’t giving the output in the required form, has temporary been handled by making an XFAIL test against it. We can work on it later. Either there has to be some changes in solve() so that we would be able to handle our case or we might have to figure out a way to rewrite it into the desired form.

With the end of this week, PR #17122 and PR #17153 are complete and ready for review. I have made some changes addressing some of the reviews, and we can have further discussions on it.

Now, also moving on to the next phase, I have done a bit of research on it. I will most probably open a discussion to have an initial discussion regarding how work will progress in this stage. This phase is regarding plotting the beam diagrams using matplotlib. I have also considered pyglet plotting module of SymPy, which according to the documentation is capable of plotting geometries, but there has been some problems in this module and it doesn’t seem to be working well. I had earlier made an issue #16537 regarding the same, but there seems to be no improvement here.

So, we will be discussing the rest in an issue-cum-discussion, in the upcoming week.

Next week:

• Working on the Stage-III
• Simultaneously, discussing the leftover PR’s and trying to finish them and make a merge.

Most probably, on successful discussion and planning, I will be opening a draft work-in-progress PR for the draw() function in stage –III.

Will keep you updated!

Thanks!

This week required a lot of thinking before jumping to code the stuff. Interested? Okay move on to next paragraph.

Basically, I worked on three PRs, #17163 for continuous time Markov chains, #17174 for random matrices and #17146 for symbolic Ranges. The first and the last PRs are very much intensive. I developed a new algorithm for the query handler of ContinuousMarkovChain.probability method, because the previous one which I implemented in DiscreteMarkovChain.probability, was not easy to maintain, quite ad-hoc, rigid and difficult to extend. The philosophy behind the algorithm is recursion i.e., boil everything down to Relational query, convert them to sets and then calculate the probability. You can find the complete description here. I am waiting for any critical objections from my mentors and after that I will refactor the code as suggested by oscarbenjamin and jksuom. So, now let’s move on to random matrices. As it was to be implemented from scratch, it required a bit of thinking to reach a decent architecture. Currently, the PR is at a basic level, and some more testing is to be done. Now, coming on to symbolic Range. Let me tell you, it requires a lot of logical thinking to make Range accept symbolic parameters. A lot of tests fail, and a lot of debugging has to be done to make a method work. In fact, we might deprecate xrange support from Range because we are going to drop Python 2 support from SymPy.

This week I learnt to combine the concepts from algorithms and software engineering to develop the stuff I mentioned above. This was the best week of my overall GSoC experience till now.

A lot more lies ahead. Bye!!

The seventh week of coding period has ended and a few methods has been introduced to polycyclic groups, also pc presentation has been modified. Previously for pc presentation we were computing the LHS for both power and conjugate relators via separate methods and then finally their RHS was computed.

Now, the computation of presentation starts from the bottom of the polycyclic generating sequence(pcgs) and polycyclic series. Storing all the previous generators from pcgs and then taking the last generator as the generator which acts as a conjugator and conjugates all the previous generators in the list.

To get a clear picture let’s take an example of S(4) For S(4) we’ll have 4 generators in pcgs say [x0, x1, x2, x3] and the relative_order vector as [2, 3, 2, 2]. Starting from bottom of this sequence the presentation is computed in order as below.

x3**2	       	---| ---> using only [x3] from pcgs and pc_series[1]
x2**2		   |
x2**-1*x3*x2	---| from bottom up because pc_series[0] is an identity.

x1**3		---| ---> using [x3, x2] from pcgs and pc_series[2]	
x1**-1*x3*x1	   | 
x1**-1*x2*x1	---|      from bottom up(which have both the gens).

x0**2		---| ---> using [x3, x2, x1] from pcgs and pc_series[3]
x0**-1*x3*x0	   |
x0**-1*x2*x0	   |      from bottom up(which have all three gens).
x0**-1*x1*x0	---|

There were 3-methods which were added namely:

• Exponent vector
• Depth
• Leading Exponent

>>> from sympy.combinatorics import *
>>> from sympy.combinatorics.free_groups import free_group
>>> G = SymmetricGroup(4)
>>> PcGroup = G.polycyclic_group()
>>> pcgs = PcGroup.pcgs
>>> group, x0, x1, x2, x3 = free_group("x0, x1, x2, x3")
>>> PcGroup.exponent_vector(G[0], group)
[1, 0, 0, 0]
>>> exp = PcGroup.exponent_vector(G[1], group)
>>> g = Permutation()
>>> for i in range(len(exp)):
...     g = g*pcgs[i] if exp[i] else g
... 
>>> g == G[1]
True
>>> 

For the details of these methods one can look into the docstrings and doctests of these methods in the PR sympy/sympy#16991.

Tasks I hope to complete next week:

• Get the polycyclic group pr ready to be merged.
• Get started with quotient groups.

Till then good byee..

This was the seventh week meeting with the GSoC mentors which was scheduled on Monday 15th July, 2019 between 6:00 - 7:00 PM (IST). Me and Yathartha were the attendees of the meeting. In this blog I will be describing the lambert equation solver for Sympy and what problems it faced before and how PR #16890 will solve the problems. It is preassumed that you know what lambert type equations are, so I will not be explaining that.

Explaining the function _solve_lambert (main function to solve lambert equations)

Input - f, symbol, gens
OutPut - Solution of f = 0 if its lambert type expression else NotImplementedError

This function separates out cases as below based on the main function present in the main equation.

For the first ones:
1a1) B**B = R != 0 (when 0, there is only a solution if the base is 0,
                   but if it is, the exp is 0 and 0**0=1
                   comes back as B*log(B) = log(R)
1a2) B*(a + b*log(B))**p = R or with monomial expanded or with whole
                            thing expanded comes back unchanged
   log(B) + p*log(a + b*log(B)) = log(R)
   lhs is Mul:
       expand log of both sides to give:
       log(B) + log(log(B)) = log(log(R))
1b) d*log(a*B + b) + c*B = R
   lhs is Add:
       isolate c*B and expand log of both sides:
       log(c) + log(B) = log(R - d*log(a*B + b))

If the equation are of type 1a1, 1a2 and 1b then the mainlog of the equation is taken into concern as the deciding factor lies in the main logarithmic term of equation.

For the next two,
   collect on main exp
   2a) (b*B + c)*exp(d*B + g) = R
       lhs is mul:
           log to give
           log(b*B + c) + d*B = log(R) - g
   2b) -b*B + g*exp(d*B + h) = R
       lhs is add:
           add b*B
           log and rearrange
           log(R + b*B) - d*B = log(g) + h

If the equation are of type 2a and 2b then the mainexp of the equation is taken into concern as the deciding factor lies in the main exponential term of equation.

3) d*p**(a*B + b) + c*B = R
   collect on main pow
   log(R - c*B) - a*B*log(p) = log(d) + b*log(p)

If the equation are of type 3 then the mainpow of the equation is taken into concern as the deciding factor lies in the main power term of equation.

Eventually from all of the three cases the equation is meant to be converted to this form:-

f(x, a..f) = a*log(b*X + c) + d*X - f = 0 which has the
solution,  X = -c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(f/a)).

And the solution calculation process is done by _lambert function.

Everything seems flawless?? You might be thinking no modification is required. Lets see what loopholes are there in it.

What does PR #16890 do?

There are basically two flaws present with the this approach.

1. Not considering all branches of equation while taking log both sides.
2. Calculation of roots should consider all roots in case having rational power.

1. Not considering all branches of equation while taking log both sides.

Let us consider this equation to be solved by _solve_lambert function.

-1/x**2 + exp(x/2)/2 = 0

So what the old _solve_lambert do is to convert this equation to following.

2*log(x) + x/2 = 0

and calculates its roots from _lambert. But it missed this branch of equation while taking log on main equation.

2*log(-x) + x/2 = 0

Yeah you can reproduce the original equation from this equation.So basically the problem was that it missed the branches of equation while taking log. And when does the main equation have more than one branch?? The terms having even powers of variable x leads to two different branches of equation.

So how it is solved? What I has done is that before actually gets into solving I preprocess the main equation and if it has more than one branches of equation while converting taking log then I consider all the equations generated from them.(with the help of _solve_even_degree_expr)

How I preprocess the equation? So what I do is I replace all the even powers of x present with even powers of t(dummy variable).

Code for targeted replacement
lhs = lhs.replace(
            lambda i:  # find symbol**even
                i.is_Pow and i.base == symbol and i.exp.is_even,
            lambda i:  # replace t**even
                t**i.exp)
Example:-
Main equation -> -1/x**2 + exp(x/2)/2 = 0
After replacement -> -1/t**2 + exp(x/2)/2 = 0

Now I take logarithms on both sides and simplify it.

After simplifying -> 2*log(t) + x/2 = 0

Now I call function _solve_even_degree_expr to replace the t with +/-x to generate two equations.

Replacing t with +/-x
1. 2*log(x) + x/2 = 0
2. 2*log(-x) + x/2 = 0

And consider the solutions of both of the equations to return all lambert real solutions of -1/x**2 + exp(x/2)/2 = 0.

Hope you could understand the logic behind this work.

2. Calculation of roots should consider all roots in case having rational power.

This flaw is in the calculation of roots in function _lambert. Earlier the function_lambert has the working like :-

1. Find all the values of a, b, c, d, e in the required loagrithmic equation
2. Then it defines a solution of the form

   -c/b + (a/d)*l where l = LambertW(d/(a*b)*exp(c*d/a/b)*exp(-f/a), k)
   

   and then it included that solution. I agree everything seems flawless here. but try to see the step where we are defining l.

Let us suppose a hypothetical algorithm just like algorithm used in _lambert in which equation to be solved is

x**3 - 1 = 0

and in which we define solution of the form

x = exp(I*2*pi/n) where n is the power of x in equation

so the algorithm will give solution

x = exp(I*2*pi/3) # but expected was [1, exp(I*2*pi/3), exp(-I*2*pi/3)]

which can be found by finding all solutions of

x**n - exp(2*I*pi) = 0

by a different correct algorithm. Thats y it was wrong. The above algorithm would have given correct values for x - 1 = 0.

And the question in your mind may arise that why only exp() because the possiblity of having more than one roots is in exp(), because if the algorithm would have been like x = a, where a is some real constant then there is not any possiblity of further roots rather than solution like x = a**(1/n). And its been done in code like this:

code
num, den = ((c*d-b*f)/a/b).as_numer_denom()
p, den = den.as_coeff_Mul()
e = exp(num/den)
t = Dummy('t')
args = [d/(a*b)*t for t in roots(t**p - e, t).keys()]

Thank you! Thats all it was!!.

With this ends the seventh week of the official coding period. During the end of 6th week and the beginning of 7th week, I was mostly travelling, so I was not able to write a blog for the sixth week. Instead, I will try to summarize my work during the last two weeks here.

For the last few weeks, I have been focused on optimizing the code of new assumptions to enhance its performance. Most of my work has been exploratory, as Aaron says 😅. Indeed I have dry-run, backtracked, profiled, and ran the same code with a debugger too many times to understand the slow parts and the improvements I can make here and there. Mostly the code is optimized given the class structure of SymPy. But it is also the class structure that is adding up to the performance issues. Already noted in my last blog, classes like And and Or sorts their args, hence take a great amount of time. But other SymPy class constructors also take significant time.

With the success of fifth week’s attempt, I have been desperate to bring down the execution time 😅. Some of the attempts I have made, which are included in #17144, are:

• I have modified CNF class, which essentially is a low-level implementation for the cnf of any boolean expression. CNF object holds a set of clauses. These clauses are themselves frozenset of Literal objects. Literal class is being implemented just to reduce the unnecessary creation of Not objects (It takes significant execution time and is called many times).
• I have also modified the code of sympify(), it appeared to take more time than expected when the argument is a SymPy object already. Consider this, almost one-third execution time.

In [2]: x = Symbol('x')

In [3]: timeit sympify(x) # before change
601 ns ± 14.6 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

In [3]: timeit sympify(x) # after change
239 ns ± 11.8 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

• Finally, I rewrote to_cnf() for CNF objects. By using mostly Python’s built-ins and removing any SymPy object construction during its execution. The performance gain is quite subtle 😎.

In [1]: from sympy import *

In [2]: from sympy.abc import a,b,c,x,y,z

In [3]: k = (x & y & z) | ( a & b & c)

In [4]: from sympy.logic.boolalg import CNF

# Before
In [5]: timeit CNF.from_prop(k) # It is using to_cnf()
1.41 ms ± 18 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

# after
In [5]: timeit CNF.from_prop(k) # It is using the new to_CNF()
31.5 µs ± 1.48 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

There is definitely a limit to performance we can get with Python. But implementing most of the things in Python built-ins we can definitely make things much faster.

For the upcoming week, I will try to complete the following:

• The CNF objects are still not simplified. I have to implement simplification to reduce the number of clauses. These have to be fed into satisfiable at the end, which can benefit from lesser number of clauses.
• rcall used with sathandlers is another major portion having high execution time. I will try to work it out.
• Clean up the code and make it ready for reviewing.

Apart from that, I will also try to shift my focus towards the other part of my project and write some theory solvers. After all, I need to enhance the assumptions not just make it faster 😎.

Peace

I spent most of this week on extending wildcard support for matrix expressions, along with some more explorations in printing array contractions.

Matrices and Wildcards

As I’ve probably mentioned in the last two blog posts, SymPy’s support for matching matrix expressions through the Wild class is currently severely limited (when it works). While it is possible to construct a non-commutative Wild, it isn’t able to match expressions in a matrix multiplication:

>>> W, X = symbols('W, X', cls=Wild, commutative=False)
>>> from sympy.abc import  N
>>> A, B = MatrixSymbol('A', N, N), MatrixSymbol('B', N, N)
>>> type((A * B).match(W * X))
<class 'NoneType'>

It’s also currently not possible to combine matrices and wildcards in expressions, since wildcards don’t have a defined shape and so may only function as scalars:

>>> W + A
TypeError: Mix of Matrix and Scalar symbols
>>> W * A
NotImplementedError: noncommutative scalars in MatMul are not supported.

MatrixWild

I spent most of this week working on #17177, which implements a MatrixWild class that functions as both a wildcard and a matrix expression. In order to construct the wildcard, we need to give it a shape:

>>> from sympy.abc import N
>>> from sympy.matrices.expressions.matexpr import MatrixWild
>>> W, X = MatrixWild('W', N, N), MatrixWild('X', N, N)

Unlike in the example above using Wild, compound expressions are able to match against a matrix multiplication:

>>> A, B = MatrixSymbol('A', N, N), MatrixSymbol('B', N, N)
>>> (A * B).match(W * X)
{W_: A, X_: B}

Note that in order for matrix wildcards to match, their shape must match with the target expression:

>>> x = MatrixSymbol('x', N, 1)
>>> e = A * x
>>> e.shape
(N, 1)
>>> type(e.match(W))
<class 'NoneType'>

However, if we don’t care about dimension, we can include another wildcard in the matrix wildcard’s shape:

>>> M = MatrixSymbol('M', 3, 3)
>>> w = Wild('w')
>>> Y = MatrixWild('Y', w, w)
>>> M.match(Y)
{w_: 3, Y_: M}

While this is a good first step to the matching functionality I was looking for with unify for rewriting matrix expressions, there is still quite a bit of functionality (and tests) to be implemented, along with an unknown number of bugs to fix.

Printing Indexed Bases

I’ve also been working on a small pull request to improve the functionality the printing IndexedBases so that it instead uses intermediate values (represented through the new code generation classes) to accumulate the values of contractions. Currently, this does nothing but break existing compatibility (Fortran versions older than Fortran 95 don’t support variable declarations in arbitrary locations, and the variable currently defaults to a 32-bit floating point number), though I think this is a good first step for supporting the printing of more complex contractions.

Next Steps

For this week, I plan to finish with the implementation ofMatrixWild (and hopefully get started with using it for rewriting matrix expressions), along with making some more progress on the indexed bases pull request.

“It’s not what we do once in a while that shapes our lives. It’s what we do consistently.” ― Anthony Robbins Welcome everyone, this is your host Nikhil Maan aka Sc0rpi0n101 and this week I’ve tried to be more consistent with my schedule, work and communications with the Organization. I’ve taken a few small steps and plan to improve as I keep working. Consistency!? Merging the Pull Requests?

The sixth week has ended with a lot of work to be done ahead.

Last week the work was majorly focused on the work in progress PR #17122. I have included the critical load function which makes the Column class capable of determining the critical load. Some problems still came up in solving some equations. I have made an issue related to those.

An equation similar to tan(x) – x comes up while determining the critical load for the pinned-fixed end-condition. SymPy’s solve() won’t be able to solve such an equation, and as per the solution given in the issue, I think that nsolve() would surely help in this case. So I will be going ahead to solve it using the approximation returned by nsolve() to handle this condition.

Another problem that I faced was determining deflection and critical load for the pinned-pinned end-condition. Here, the deflection comes out to be:

C1*sin(sqrt(P)*x/(sqrt(E)*sqrt(I))) + C2*cos(sqrt(P)*x/(sqrt(E)*sqrt(I)))

Now on solving it for constants C1 and C2, using initial boundary conditions, both come out to be 0, making the deflection zero. This implies that no buckling occurs, which is not the case.

Even when solving it manually, this situation occurs, we deal with it by putting C2 = 0 and instead of putting C1 = 0, we consider the sin term equal to zero and then solve for P (critical load). So, I will be adding a few more lines of code to deal with this situation.

Apart from working on this module, I have also opened another PR #17153 which implement methods to determine section modulus and polar modulus of any polygon (more precisely a cross-section). Initially it was a draft PR, but now the work has been completed on it. Once I get the approval, I will also be adding the same for the Ellipses module. Also, if cut_section() gets successfully implemented I will be adding another method to determine the first moment.

I am pretty sure the work on Column class will be successfully completed before the end of the next week. Also, we will be heading towards the next stage which intends to plot beam diagrams using matplotlib. Till then we can have an initial discussion regarding the same.

Next Week:

• Improving the critical_load() to handle the above problems
• Completing the Column class (documentation and tests)
• Starting with the next phase

Will keep you updated!

Thanks!

This week was a mix of discussion on design and extending previous work. I also got to know about some new cool features of SymPy.

According to the plan proposed in Week 4, I have completed my work on DiscreteMarkovChain via PR #17083. I used the as_set and as_relational methods which helped me to cover many miscellaneous cases and probably, now DiscreteMarkovChain is quite dynamic and can handle various generic probability and expectation queries. I have also started the PR #17163 for adding ContinuousMarkovChain and I am observing that it’s a bit tricky to maintain both the performance and result quality while working on it. Now, moving on to symbolic Range,well, the work has been started in the PR #17146 and I have figured out one disparity between Range and python’s range(details available at this thread). I will try to fix it by making minimal changes to the code. The tensorflow related PR #17103 which I started in the previous week is also almost complete and is waiting for Tensorflow 2.0 release. I am also studying a bit about the architecture of the above framework to make changes to lambdify. Regarding random matrices, I believe that discussion has reached its final stages and I am waiting for the comments from Francesco for improvements at the issue #17039.

Let me share with you about my discoveries and learnings in this week. Well, thanks to Francesco for telling me about, sympy.multipledispatch. It helps in implementing operator overloading like in C/C++. I liked it very much. I also read about continuous Markov chain and discovered about generator matrix, forward and backward equations. Adding one interesting fact, that Poisson process and continuous Markov chain are very closely related via generator matrices it will make the implementation of the former much easier :D.

Leaving you for now, Bye!!

The sixth week of coding period has ended and a good amount of work has been done on polycyclic groups. Polycyclic presentation, Polycyclic generating sequence(pcgs) and it’s series is implemented which for sure need some improvement sympy/sympy#16991.

The polycyclic series is computed starting from the bottom of the derived series of a group by adding the missing generators in the subgroups, and collecting these missing generators provide us the polycyclic generating sequence.

As we discussed last week here that to compute conjugate relators of a polycyclic group we were missing the RHS term, which was of the form x[i]**-1*x[i+1]*x[i] == RHS. So, starting from the bottom of the polycyclic generating sequence forming the subgroup and finding all the generators of the RHS using generator_product, mapping these generators with the free group elements and forming a word, finally collect the above formed word which will give us the collected RHS.

Below is an example to compute polycyclic presentation for S(9).sylow_subgroup(3)

>>> from sympy.combinatorics import *
>>> from sympy.combinatorics.free_groups import free_group
>>> F, x0, x1, x2, x3 = free_group("x0, x1, x2, x3")
>>> S = SymmetricGroup(9)
>>> G = S.sylow_subgroup(3)
>>> pc_group = G.polycyclic_group()
>>> group = F
>>> pc_group.pc_presentation(group)
{x3**3: (), x2**3: (), x1**3: (), x0**3: (), x2**-1*x3*x2: x3, x1**-1*x3*x1: x3, x1**-1*x2*x1: x2, x0**-1*x3*x0: x2**2*x3**2, x0**-1*x2*x0: x3, x0**-1*x1*x0: x1*x3}

One problem that we’ve encountered is that the generators in pcgs may change for the same group on executing it several times which makes it difficult to test pc_presentation but, Kalevi advised me to initalize random.seed with some chosen value and then it will result in the same repeatable result, will try it by today!

The tasks that I’m hopping to accomplish next week are

• Add tests for plycyclic presentation and be sure that it works properly.
• Include more functionalities to pc groups like exponent_vector, element_depth, leading_coefficient.
• Add documentation for all the functions.

Till then, good byee..

Week 6 ends.. - Phase 2 of the coding period is smoothly being traversed. I recently had a meeting with Sartaj on 4th of July, Thursday. Here were the minutes of the meeting, along with the deliverables completed over the week. The meeting commenced with the discussion that we need to wrap up the...

See the previous post for Week 5.

This week I’ve made some progress on matching and tensors, though I haven’t filed any pull requests.

Unification

I have a working implementation of rewriting non-commutative expressions using SymPy’s unify. It works by generating a ReplaceOptim object that applies the rewriting rules to any term it’s called with. Here’s how we specify the rewriting rules:

>>> from sympy import Symbol, MatrixSymbol
>>> n = Symbol('N_matcher', integer=True)
>>> X = MatrixSymbol('X_matcher', n, n)
>>> Y = MatrixSymbol('Y_matcher', n, 1)
>>> variables = [n, X, Y]
>>> matcher = X**(-1) * Y
>>> goal = MatrixSolve(X, Y)

Here, the combination of matcher and variables specifies that we’re looking for any expression of the form X^− 1Y, where both X and Y can be any compound matrix expression. The inclusion of n in variables imposes the additional restriction that the matrix expression matched by X must be square (i.e. n × n) while the expression matched by Y must be a vector (i.e. n × 1). goal specifies what the matched expression should be replaced with, where X and Y serve as stand-ins for the matched terms.

After specifying our goals, we can construct the object and apply the replacement to some expressions:

>>> replacer = gen_replacement_operator(matcher, goal, variables)
>>> A, B, x = MatrixSymbol('A', 3, 3), MatrixSymbol('B', 3, 3), MatrixSymbol('x', 3, 1)
>>> replacer(A ** (-1) * x)
(MatrixSolve(A, vector=x))
>>> replacer(A ** (-1) * B)
A ** (-1) * B

The first term was replaced since the dimensions of A and x agreed with what was specified in matcher, while the second expression was left untouched since B is not a vector.

While the matcher does work, I haven’t filed a pull request because of some problems which don’t seem like they could be easily addressed:

• I had to give add the suffix _matcher to the variable names to avoid variable capture, since SymPy symbols are considered equal if they have the same name. unify does not support Dummy symbols as variables.
• Some compound expressions are not matched. I’ve narrowed this down to the way the variables are being passed to unify, since they need to be converted to symbols. It seems like this conversion sometimes causes expressions to no longer be unifiable.
• Unification doesn’t seem to work for a mixture of commutative and non-commutative expressions. I’m not sure if this is a problem with unify itself or the way that I’m using it, since the only test of unify in the SymPy codebase involving matrix expressions is on matrix multiplication.

As I mentioned in my last blog post, SymPy already supports this sort of pattern matching through Wild, though it currently does not support expressions involving matrices. Before trying to address these issues, I think it would be worthwhile to look into extending the functionality of Wild as an alternative.

Tensors

I’ve made some progress in low-level code generation of matrix expressions. I tried seeing if instances of classes in the array_utils module could be converted to SymPy’s AST representation before being passed off to the code generators. This doesn’t seem possible at the moment, since the AST has a number of limitations (such as not supporting variables in for loop ranges). The IndexedBase printer already has some of the functionality that I’m trying to implement, so I’ve settled on extending the printer to support arbitrary contractions. This same functionality can probably be reused for the array_utils printers. The implementation will hopefully be straightforward.

Next steps

My goal for this week is to have a pull request for the tensor code generation ready, along with a plan for what to do with matching.

This was the sixth week meeting with the GSoC mentors which was scheduled on Sunday 7th July, 2019 between 1:45 - 2:45 PM (IST). Me, Yathartha and Amit were the attendees of the meeting. This meeting was short.

• Discussing previous week’s progress

In this meeting both the mentors were convinced by the code for Lambert’s. And few modifications in documentation and code clean up were suggested by them. In this week the whole idea of power_list was droppped because @smichr suggested code for replacing the symbol more targetted as we wanted by which the whole code was improved. And it was decided to work upon on _solve_modular mainly now onwards.

• Next week goals

• Getting merge existing PR for Lambert

• Work upon _solve_modular PR

• If time left then find plan for Imageset Union.

Code improvement takes time!!

“A language that doesn’t affect the way you think about programming is not worth knowing.” ― Alan J. Perlis Welcome everyone, this is your host Nikhil Maan aka Sc0rpi0n101 and this week was one of the most chaotic weeks till now if not the most chaotic. I had to move places as the Uni opens next week, had the summer classes exams and had to convert the C parser to use SymPy’s Codegen AST from Python AST all in one week.

With this the fifth week and the first phase of the official coding period has ended. I will try to give a brief summary of my work during this week.

I spent most of this week learning the inner working of satask module. It took a lot of debugging to understand the on-going processes efficiently 😅 . My major job was to reduce the unwanted slowing portions in the code. I had #11789 for reference. Some of such performance reducing portions of code were:

• A major segment of code was creating And objects unnecessarily. As suggested in #17087 the sorting in the And constructor takes up a significant amount of time. These have to be reduced.
• Using SymPy objects is itself a bottleneck for performance. Having a system built over SymPy objects slows things down. Python’s built-in types should be used as much as possible.
• A specific segment (used many times in the code) calls rcall over propositions. The rcall which is a recursive process also takes a significant of time.
• Also, I have tried to pre-compile results as much as possible.

With the above and some other changes, the overall performance of satask has improved much. I have made a PR over this #17144. For an instance,

There is still scope for improvement in performance. In the coming week, I will try to work these out. I will also try to improve the performance of ask module.

Also, the first evaluations are over now and I feel happy to announce that I passed it. During the first phase I learnt a lot. In last few weeks I got to explore profiling and got to understand how small segments can influence performance. Before this I felt that I already know the codebase, but in reality I had much to explore. My mentors always gave me a good starting point and a direction over the course of this phase. With the hope to work much better in the coming phases, I take your leave now 😄 .

The evaluation results for phase 1 are out, and I am very glad to share with you that I have passed with flying colors. I received, “Well done so far.” as the feedback for my work till now.

So now let us move to the work done in the gap between phase 1 and phase 2. Firstly, both of my open PRs of the previous phase, i.e., #16962 and #16934 have been merged. Though for symbolic dimensions some more work has to be done to make sympy.stats.frv more efficient and maintainable. I have also started my work, PR #17083, to extend the scope of queries for DiscreteMarkovChain and the system has become a bit smarter. In fact, during this week, while working on the PR, #17103, I came across the news that Tensorflow has changed a lot of APIs while migrating from 1.x to 2.x. AFAIK, they are moving towards Function approach from the previous Session approach, and due to that, SymPy’s lambdify faced some issues which I will be fixing soon with the help of other members. The Tensorflow details can be seen here.

Now, let’s move to the learning part. During the transition period I learnt about the dependencies of SymPy. Moreover, I came across, how, some bugs can be unnoticed when left untested. Thanks again to oscarbenjamin for letting me know about the bugs related to variance of finite random variables. I also got to know that, how bare except can even catch keyboard interrupt and that’s what makes it quite vulnerable. Thanks to sidhantnagpal for helping me with this.

So, that’s all for this, see you next week. Bye!!

This week was mostly about testing the collection of a word and fixing small bugs in the implementation pointed out by Kalevi. The major challenge was to construct the polycyclic presentation of a group to test the Collector since we don’t have the implementation of polycyclic presentation and it’s generating sequence yet. So, we decided to form some hand made tests and we started with SymmetricGroup(4) and further we also tried with S(3) the details can be found in the test file of the PR(here).

Now, the next step is to implement polycyclic presentation and polycyclic sequence. In the presentation we’ll need generators which we can easily get and the relators. There are two types of relators needed for the presentation:

• Power relations (ex. x^re = x')
• Conjugate relations (ex. x[i]**-1*x[i+1]*x[i] = RHS and x[i]*x[i+1]*x[i]**-1 = RHS)

For every pair of generators we’ll form above conjugate relations but the tough part is computing that RHS which should be collected and for now we don’t have that much idea about how to get that RHS.

But, let’s hope that in upcoming days we’ll be able to figure it out, till then Good byeee…

A lot of things happened this week and I am happy to inform you that PR #17055 has been successfully merged. The beam module now supports the cross-sectional shape of the beam as an alternative parameter to the second moment. With this, the aim of the stage-I to integrate the geometry module with beam module has been accomplished.

Although we need to add some examples in the docs, to make it easier for the user to understand how to use this new feature.

Coming on to stage-II, I had already, initiated a discussion to finalize the API of the new Column class that is to be implemented as a part of the continuum mechanics module in this stage.

We concluded that it would be much better if the Column class remains non-mutable i.e. unlike the beam class where a beam is formed in a piecewise form, the new Column class would take all its required input data during the declaration and then one can call different methods to calculate different things.

I have made a work-in-progress PR #17122 implementing the Column class which performs the required buckling calculations. Currently, I have not included a method to calculate the critical load as there was a bit of problem with the form of the equation which the dsolve() returns after solving the differential equation of buckling. dsolve() is SymPy’s differential equation solver.

In general, if we solve the general equation of buckling manually, we might apply the method of undetermined coefficients, which of course even dsolve()  is capable to apply, but it gives the answer in an exponent form, while we need it in a trigonometric form (for ease of further calculations). So after seeking different methods trying to convert this equation in terms of sin(x) and cos(x), I finally had to put that problem in the discussion, where Oscar Benjamin, gave an idea to declare the variables as positive in order to get it in terms of sin and cos. I tried that it works well for our case. I will have to figure out the further calculation of the critical load.

Hopefully will be updating the code with a new method to calculate critical load, soon.

Also, I have planned to have a method to solve the unknown reactions and reaction moments, which would use the boundary conditions to get their values.

With all these things going on, this week we also had our first evaluations, and I am very happy to say that I have passed it. Thanks to the mentors!

 Next Week:

• Completing Column class with all its methods
• Adding tests and documentation.
• Starting discussions for the next stage.

I will try to finish working on the Column class this weekend.

Will keep you updated!

Thanks!

Week 5 ends.. - Phase 2 of the coding period has started. This week has gone in wrapping up the left-over work of FormalPowerSeries. I had a meeting with Sartaj on Tuesday 25th of June, about the work left to be done on FormalPowerSeries module. We agreed that some minor changes need to be...

Hello, the first phase is ended and I am happy to pass the first evaluation. I was struggling with my academic projects and final exams during the last two weeks. After talking about my difficulty of spending time contributing on my project with my mentors, Francesco allowed me to have a one-week break in condition that I should make up one week in the next phases. The goal is to have 40 hours work per week on average by the end of this program.

Thanks to the comprehension of my mentor, I could successfully pass the exams. I am going to work more over the second phase in order to have more contributions to the community.

This was the fifth week meeting with the GSoC mentors which was scheduled on Saturday 29th June, 2019 between 11:30 - 12:30 PM (IST). Me, Yathartha and Amit were the attendees of the meeting. I passed my first evaluation, Amit gave his feedback and told me some very important points to take notes on. I do personally believe that his suggestions are the best a mentor could gave to his student after practicing his suggestions in my real life.

• Discussing previous week’s progress

In this meeting both mentors suggested me to work upon the code improvements and documentation improvement. And make it more readable to user. ALthough somehow @smichr had some doubts on the logic that we were implementing. Although a lot of progress has been there. So I decided to create and discussion for thinking new logic for implementing Lambert all solutions and work on the current PR as goes on.

• Next week goals

• Improving existing PR for Lambert

• Improving _solve_modular PR also

• If time left then find plan for Imageset Union.

Code improvement takes time!!

“On two occasions I have been asked [by members of Parliament!]: ‘Pray, Mr. Babbage, if you put into the machine wrong figures, will the right answers come out ?’ I am not able rightly to apprehend the kind of confusion of ideas that could provoke such a question.” — Charles Babbage Welcome everyone, this is your host Nikhil Maan aka Sc0rpi0n101 and this week was the last week of Phase-1 and the evaluations for the first phase.

See the previous post for Weeks 3 and 4.

This week I’ve been mostly doing background reading. This post is mostly a summary of what I learned.

Unification

In short, unification is the process of finding substitutions of variables within two terms two terms to make them identical. For example, if we have the expressions x + 2y and a + 3b, the substitution {x ↦ a, y ↦ 3, b ↦ 2} is a unifier, since applying the substitution to both expressions makes gives us the identical expression of a + 3 ⋅ 2. While this particular substitution includes variables from both expressions, we’re mostly interested in rules involving substitutions of variables from just one expression (a case of unification known as matching). Several well-known algorithms for unification already exist.

Unification in SymPy

SymPy also has an implementation of a unification algorithm that is able to take the commutativity of operations into account. Suppose we wanted to unify the matrix expressions A^TB²C^− 1 and XY^− 1. This is essentially the problem of finding a substitution that makes these two expressions equal. Using the sympy.unify.usympy module, we can discover what this substitution is:

>>> from sympy.unify.usympy import *
>>> from sympy.abc import N
>>> m = lambda x: MatrixSymbol(x, N, N)
>>> A, B, C, X, Y = map(m, ['A', 'B', 'X', 'Y'])
>>> e1 = A.T * B**2 * C.I
>>> e2 = X * Y **(-1)
>>> next(unify(e1, e2, variables=[X, Y]))
{X: A.T*B**2, Y: C}

We’ve reduced this to a matching problem in which the variables are specified only in e2. What’s important to note here is that the matching rule within e2 we specified (XY^− 1) was a compound expression. This is something that is currently not possible for non-commutative expressions (such as matrix multiplication) using SymPy’s Wild interface. unify allows use to express substitution rules that are able to match across sub-expressions in matrix multiplication.

Through unification, we can express substitution rules for optimization as a simple term-rewriting rule. In my previous blog post, I mentioned rewriting the matrix multiplication Ax as a solving operation of MatSolve(A, x) under certain assumptions. The actual implementation is restricted to cases where both the A and the x are matrix symbols, and the optimization can’t identify cases where either the A or the x is a compound expression. With unification, we can identify the same pattern in more complex subexpressions. If we’re given the matrix expression A^T(AB)^− 1xy, a unification based transformer can produce MatSolve(AB, x), provided that the shapes of the matrices match the given rule.

Codegen Tensors

I also looked into generating C and Fortran code from SymPy matrix expressions. For the purposes of code generation, SymPy has a relatively new array_utils module. The AST nodes in this module express generalizations of operations on matrices, which require a bit of background in tensors.

Many array operations (including matrix multiplication) involve contraction along an axis. Contractions are a combination of multiplication and summation along certain axis of a tensor¹. In assigning the matrix multiplication AB to the n × n matrix C, we can explicitly write the summations (using subscripts for indexing matrix elements) as

$$C_{ik} = \sum_{j = 1}^{n} A_{ij} B_{jk}$$

The index j is contracted, as it is shared between both A and B, and describing this summation operation as a whole boils down to which indices are shared between the matrices. This is essentially what the array_utils classes do. This is what happens when we use array_utils to convert the matrix multiplication to an equivalent contraction operation:

>>> from sympy.codegen.array_utils import CodegenArrayContraction
>>> from sympy.abc import N
>>> A = MatrixSymbol('A', N, N)
>>> B = MatrixSymbol('B', N, N)
>>> CodegenArrayContraction.from_MatMul(A * B)
CodegenArrayContraction(CodegenArrayTensorProduct(A, B), (1, 2))

We’re given a newCodegenArrayContraction object that stores, along with the variables A and B, tuples of integers representing contractions along certain indices. Here, the (1, 2) means that the variable at index 1 and index 2 (indices start at 0) are shared. We can confirm this by looking at the above summation, since both the second and third indices out of all indices that appear in the expression are j.

Next Steps

For next week, I’ll try to re-implement the rewriting optimization in terms of unify. This will both make it easier to express rules and extend to sub-expressions as well. I’ll also start with implementing additional printers for the C and Fortran printers. The printer will probably just print naive for loops to keep things simple (and it would probaly be better to use something like Theano for highly optimized code).

The fourth week of coding period has ended and now it’s time for phase-I evaluations. Below, is a brief progress report of the project.

The tasks that were proposed in the proposal for phase-I consists of:

• Implementation of Abelian Invariants
• Implementation of Composition Series
• Computation with Polycyclic groups

Abelian Invariants

Implemented a function to compute the abelian invariants for a given permutation or free group. These are given as a list of prime-powers and describe the stucture of G/G' as a direct product of cyclic groups of prime power order.

• PR link Added method to calculate Abelian Invariants

Composition Series

Implemented a function to compute the composition series. It provides a way to break up a group into simple pieces. Composition series of a group G is defined as the maximal subnormal series G = H_0 > H_1 > H_2 ... > H_k = 1 where every factor group H(i+1)/H(i) is simple.

• PR link Added method for Composition Series computation

Polycyclic Groups

The work on polycyclic group is in progress. For now, collection algorithm has been implemented which needs to be tested and a lot of discussions were made on the polycyclic generating sequence and its presentation and may be in a week we’ll be ready with the stucture of polycyclic groups and collection of words.

• PR link Added Polycyclic Group Class

Documentation

Some documentation is done to increase the sphinx coverage of SymPy.

• PR link Increase Accessibility of docstrings from Sphinx

To follow the discussion on above topics and the further progress of the project one can check Gitter room sympy/GroupTheory

With this the fourth week and the first phase of GSoC 2019 is over. Here I will give you a brief summary of my progress this week.

I started this week by setting up my workspace for profiling the code related to new assumptions. I am using pyinstrument for that. The results of profiler suggests that a significant amount of time is spent in the to_cnf() function which converts the logical expression into their CNF counterparts, to be used by the SAT solver. Also, since this system is built over the SymPy core, a large amount of this time is spent in the core itself (See the graph here). A possible solution to this is to use constructs at a level lower than the SymPy objects, hence removing the overheads.

Also, as suggested in the last blog, there are various ideas proposed for improving the new assumptions mechanism. Some of them have been implemented to some extent in some PRs. Before proceeding for any new strategies, I need to look into these ideas first. I have started an issue-tree to gather them.

Over the end of the week, I also pushed my work on First Order Logic module at #17069. This work is based on #7608 and extends it by adding Equality to it. Currently, there are test failures and some points to decide. I will try to get it done within this week.

I spent most of this week exploring the profiling and benchmarking of code, and I learnt a lot during this. For the coming week, I will focus on speeding up the code in to_cnf. As suggested by Aaron, this seems a good point to start with. Also, I will be working on the FOL module.