The code is available on github in case anyone is interested. It requires sage to run though. https://github.com/syamajala/KnotContac … /braids.py

]]>Another thing we realized is that there is a way to linearize, which we are interested in. Anytime we have a product like a_ij*a_kl we can replace that by -2*a_ij*-2*a_kl - 4. I haven't used this fact yet in my code.

Although it reduces the problem to something that's no longer interesting, it seems apparent to me that you can use this fact to rewrite phi_s*k* so that every invocation, no matter how deeply nested, always returns a data structure of fixed size. If that's possible, you'd certainly save a lot of time by doing so in the original calculations rather than by calculating a huge tree of lists and then running it through an algorithm to simplify it.

It's also possible I've misunderstood something, so if you have any follow up questions feel free to ask.

]]>`-1a_3_1*a_1_2*a_2_1*a_1_4*a_4_1*a_1_2*a_2_1*a_1_2*a_2_1*a_1_4*a_4_1*a_1_2*a_2_1*a_1_3*a_3_1*a_1_2*a_2_1*a_1_4*a_4_1*a_1_2*a_2_1*a_1_4*a_4_1*a_1_2*a_2_1*a_1_3*a_3_1*a_1_2*a_2_1*a_1_0`

So it's definitely big, but not unmanageable. The main source of slowdown is probably Python, and if you rewrote it in a more appropriate language you might have been able to generate the sequence in considerably less time (but still less than instantaneous -- say, 30 minutes instead of 3 hours). On the other hand, if you're ok with leaving it to crunch numbers overnight, maybe Python's speed isn't a problem here. Still, if you're planning to deal with much bigger expressions, you will soon find yourself leaving it over the weekend, or longer.

Like I said before, it does generate sets of like terms that need to be combined, so you'll need to account for that in the code that makes the matrix Phi_L_B. I haven't read the paper and you didn't explain exactly how that matrix is calculated; I'm assuming the matrix you're going for is much smaller than the entire expression, which is why generators work with you (since you don't need the entire expression at once).

@jakob: I was thinking the same thing when I started, but having finished I'm not sure tail-call optimization can help here. The recursion depth isn't that great -- only around 25 -- and I'm not even sure it's possible to do TC optimization in this case; the logic is just too complex. However, a Lisp would definitely be a more natural choice for this kind of problem.

]]>Taviantor, I need the function to build the Phi_L_B matrix I described in my third post. If you mean more generally what do I need this function for, basically, I'm trying to reimplement the algorithms described in this paper, http://arxiv.org/pdf/math/0302099.pdf The guy who wrote that paper wrote a mathematica program to do the computations but only 1 person in our group has mathematica and he said that the program was written for version 5 but the latest version is 9 so he had some trouble getting it to work.

Anyway, it turns out that there is another method described in the paper for generating the Phi_L_B matrix using a "so called" chain rule. I've also been playing with implementing that.

Another thing we realized is that there is a way to linearize, which we are interested in. Anytime we have a product like a_ij*a_kl we can replace that by -2*a_ij*-2*a_kl - 4. I haven't used this fact yet in my code.

Trent, I'll take a look at your code tomorrow when I am less tired. Also, I gave cython a try on my "optimized" expand function in my second post, but all that happened was that it became faster at expanding small terms, according to the profiling I did.

]]>http://paste.pound-python.org/show/34425/

The classes are just containers for variables and terms ("term" being a sequence of variables multiplied together in order, plus a numeric coefficient). I wrote them to ease the complexity of the nesting and to make it easy to print out an expression for debugging. If you need the result in something other than a list-of-terms, it shouldn't be too difficult to write a conversion function. **Edit**: Custom classes removed; interpret() now returns a sequence of dicts, which represent terms. Closer to the original problem and less extraneous fat.

The interpret() function works on the recursive data structure you defined (e.g. [(1, 0), (2, 1), [(0, 1), (1, 2), (2, 0)]] ) and returns a generator. The generator yields terms, one after the other. Like terms are not combined, because to do so you would have to keep track of all previous terms, which is not possible for obvious reasons. I found that it was surprisingly hard to write this in a way that didn't bleed a *lot* of memory, but I think it works now; I can write a loop to churn through the generator output and leave it running for some time without the process being killed or locking up my desktop.

However, that brings me to the other caveat: the example input at the end of the code you originally posted? It's really big. Really, really, really big. The expansion contains millions of terms. In the last 10 minutes my computer has generated 32 million terms and I have no way of knowing how close it is to "finished". I imagine the time required to do so could be much, much longer than you could conceivably wait. If my gut instincts are right, rewriting it in C wouldn't help you much; you need to find another approach that doesn't involve enumerating every term in the expression. Perhaps there is some shortcut you can take when generating the expression lists.

To give you an idea of the scale of the problem, the following relatively short input has an expansion 33 terms long:

`[(1, 0), [(1, 4), (1, 2), [[(2, 4), (2, 1), (1, 4)], (2, 1), (1, 4)]], [[(4, 0), (4, 1), (1, 0)], [[(4, 2), (4, 1), (1, 2)], (4, 1), (1, 2)], [(2, 0), (2, 1), (1, 0)]]]`

The next input I tried, which is about 4 times as long, has an expansion 737 terms long, and the next input, which is about twice as long as that, has an expansion 64,864 terms long and takes a few seconds to calculate. The entire expression is over 3 times longer than the 65k-term expression. There's no rigid correlation between input size and number of terms, which just makes it hard to reason about performance. I tried to figure a way to calculate how big the expansion would be, but then I realized I'm not getting paid for this.

tl;dr -- here's some code that works better for shorter inputs and doesn't lock up for longer ones, but I can't guarantee it will finish churning the test data this millenium. 70 million terms and no end in sight. Hope you didn't need the whole expression at once because that just ain't going to happen.

]]>From my experiences you need recursion very very rarely in python.]]>

Also, I simplified what expand returns. Its now [(-1, (0, 2)), (-1, (0, 1, 1, 2))] instead of a list of tuples for the second entry in the tuple. Now, when iterating through we just take 2 elements a time instead of 1. I might possibly be able to simplify it more and just have it be a list of tuples. [(-1, 0, 2), (-1, 0, 1, 1, 2)]

]]>[(1, 0), (2, 1), [(0, 1), (1, 2), (2, 0)]] are the indices into the matrix Ga and represents -a1s - a21*(-as1 - a12*a2s)) we want to expand it into this -a1s + a21*as1 + a21*a12*a2s, but we need to keep track of certain things like whether we are subtracting or adding an element and which elements we are multiplying together, so, we use a list like [(-1, [(1, 0)]), (1, [(2, 1), (0, 1)]), (1, [(2, 1), (1, 2), (2, 0)])] where in a tuple, the first element is the sign, and the second element is a list of elements we multiple together.

That's just too much nesting, both for Python and for mortals. Plus the meaning of each data structure depends heavily on both its contents and context. Only a mathematician would come up with something like this...

Okay, okay, here's what I'm thinking: you have two options. Preferably, you would rewrite your phi_s*N* functions so that they simplify as they go, and always return a simple (non-recursive) data structure like a list-of-terms; that would probably cost much less than trying to simplify the entire massive expression after the fact.

The other option is to store the whole mess into a tree-like data structure, and simplify it in-place using postorder traversal. Which is not entirely different from what you're doing now, but you need to store it in such a manner that the context doesn't determine what you do with a subexpression (i.e. is this tuple part of a list, or is this list the 2nd member of a tuple, etc.) -- make the data structure more symmetric, and you'll be a good step closer to taming that complexity. It has to be postorder and in-place to keep your memory usage under control.

Actually, if you can modify your existing algorithm to do the simplification in-place (rather than appending to a new data structure) that might help enough to make it work. I'm not really sure that's possible though. In any case, I doubt I'll be much more help.

]]>```
Ga = [as0, as1, as2]
[a1s, a11, a12]
[a2s, a21, a22]
```

We then want to compose a function phi_sk (where k goes between 1 and n-1) with either itself or some other function phi_sk, n times and apply that function to each as* and a*s element. Doing that we get lists like

`a = [(1, 0), (2, 1), [(0, 1), (1, 2), (2, 0)]]`

The tuples in the list represent the indices into the matrix Ga. By definition of phi_sk, we'll either get a single element (in this case for example phi_s1(a2s) = a1s), or we'll get 3 elements (phi_s1(as1) = -as2 - as1*a12) that will always look like -1*(some element) - (the product of two others). So, I represented the equation -as2 - as1*a12 by the list [(0, 2), (0, 1), (1, 2)].

So, now say we've applied all the phi_sk's to all the as* and a*s elements n times. We'll get back lists like [(1, 0), (2, 1), [(0, 1), (1, 2), (2, 0)]] and even much more complicated ones like this [[(2, 0), (2, 1), (1, 0)], (2, 1), [(1, 0), (1, 2), [(2, 0), (2, 1), (1, 0)]]] What we want do now is build another matrix Phi_L_B where the ij entry are the coffeicients of the ajs terms of the composition of the phi_sk's on a_is. In order to get those coefficients, we need to take a list like [(0, 2), (0, 1), (1, 2)] and associate a sign with each element and figure out the products and sign of the products.

For example looking at this [(1, 0), (2, 1), [(0, 1), (1, 2), (2, 0)]] this represents -a1s - a21*(-as1 - a12*a2s))

if we expand that out we get -a1s + a21*as1 + a21*a12*a2s. Its important to note that the structure we are working over DOES NOT have a commutative multiplication so we cannot change the order. So, I choose to represent the expanded out version in the following way, [(-1, [(1, 0)]), (1, [(2, 1), (0, 1)]), (1, [(2, 1), (1, 2), (2, 0)])] where in a tuple, the first element is the sign, and the second element is a list of elements we multiple together.

**TLDR**

[(1, 0), (2, 1), [(0, 1), (1, 2), (2, 0)]] are the indices into the matrix Ga and represents -a1s - a21*(-as1 - a12*a2s)) we want to expand it into this -a1s + a21*as1 + a21*a12*a2s, but we need to keep track of certain things like whether we are subtracting or adding an element and which elements we are multiplying together, so, we use a list like [(-1, [(1, 0)]), (1, [(2, 1), (0, 1)]), (1, [(2, 1), (1, 2), (2, 0)])] where in a tuple, the first element is the sign, and the second element is a list of elements we multiple together.

]]>