A Python Riddle

From the book Fluent Python (which you can get from the Humble Bundle right now):

What would the following piece of code do?

t = (1, 2, [30, 40])
t[2] += [50, 60]

Four choices:

    1. t becomes (1,2,[30,40,50,60])
    2. A TypeError is raised because tuples does not support item assignments?
    3. Neither
    4. Both 1. and 2.

The answer is the link, which is quite surprising.

It turns out that t[2].extend([50, 60])doesn’t break Python, and this riddle is really a super esoteric corner case…

Vectorization over C

The title is probably misleading, but this is a lesson I needed to talk about.

I wrote out some simple code for the quadrature over the reference triangle last time, which involves a double loop. To my chagrin, my immediate reaction to speeding up the code was to put it into Cython, and give it some type declaration.

This did speed up my integrals, but not as much as vectorization. By simply condensing one of the loops into a dot product, and using vector-function evaluation, I sped up my code a substantial amount, especially with higher order integration of “hard” functions.

def quadtriangle_vector(f, w00, w10, x00, x10):
    total = 0
    for i in range(len(w00)):
        total += w00[i] * np.dot(w10 / 2, f([(1 + x00[i]) * (1 - x10) / 2 - 1, x10]))
    return total

To see what I mean, consider the following function

from scipy.special import eval_jacobi as jac
def f(x):
    return jac(2, 1, 1, np.sin(x[0] - x[1]))
p = 20
x00, w00 = rj(p + 1, 0, 0)
x10, w10 = rj(p + 1, 1, 0)

The speedup I get is staggering.

%timeit quadtriangle(f, w00, w10, x00, x10)
3.11 ms ± 24.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

348 µs ± 1.67 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Also, I tried to fully vectorize by removing the outer-loop. This actually slowed down the code a bit. Maybe I did it wrong? But for now, I’m decently happy with the speed.

The (lack of a) Matrix

I think I finally understand why software packages like PETSc has an option for an operator when doing something like conjugate gradient. Why isn’t having a matrix good enough for everyone?

Well turns out that while all linear operators can be translated to a matrix, it may not be the best way to represent the operator. As an example, consider a basis transformation from Bernstein polynomials to Jacobi (or vice versa). It’s certainly possible to find and construct a matrix which does the operation, but it’s ugly.

On the other hand, it’s not that bad to write a code which utilizes the properties of the polynomials and convert it within in O(n^2) time. The key is that a Jacobi polynomial is a sum of Bernsteins, and Bernsteins can be degree raised or dropped at will.

This function will outperform the matrix in many sense. For one, there’s no need to construct a matrix, which will take n^2 operations in the first place. Next, matrix multiplication will take a n^3 operation, so if we optimize enough, we will always beat it. Finally, it’s really less painful to code, because each line of the function serves a visible purpose.

Anyways, I’m sold.

(I’ll eventually publish the code in the summer)

Notes: SSD edition

Some notes from the past week:

  1. It is incredibly easy to be an impostor in a more academic party. First of all, most of the people will be already intoxicated to the point where bullshit science can’t be discerned from actual science. This is good as I can just say random facts I remember from Popular Science.Another acceptable thing to do is to just ask questions upon questions. “What’s your research? … Oh that’s so cool! Tell me more about it! … So does this connect to insert scientific news here? Wow.” That’ll burn around 5 minutes minimum.The main problem comes when you run out of questions in the initial barrage. It also fails when the person is laconic or can’t speak English.

  2. Installing a SSD is extremely easy, but installing operating systems are hard. Right now, I have around 8 entries on my GRUB menu before I migrate everything over to my new distro.I followed the mount guide provided here, which seems intuitive enough on where to put mount points. I’ve also learned that
    mount

    and

    df -h

    are my friends. There’s also that good GParted software.

  3. The Lloyd Trefethren numerical linear algebra book is quite good for a quick overview of the subject. It doesn’t get bogged down with the analysis, and generally refers to other books (mainly the Van Loan) throughout.
  4. Holy shit URF mode.
  5. I need to be more brave in a certain subject….

Ray Casting with JOGL

I won’t post the entire code here, because it’s pretty damn ugly. But here’s what I ended up doing:

  1. I used the
    gluunproject

    statement to find the beginning and end points to extrapolate a line from.

  2. Now that I have a line, I use the formula provided by Wikipedia using the vector formulation of the line.
  3. Simply do a loop over the vertices and find the minimum.

Sorry for note posting recently… I got caught up in things… 🙁

VTune Profiler Error: “The Data Cannot be displayed, there is no viewpoint available for data ”

The solution to this if you’re using the GUI on Linux can quite possibly be that ptrace_scope was set to 1.

From the Intel forums:

Note: In Ubuntu 11.4, you may need to disable ptrace_scope.

cat /proc/sys/kernel/yama/ptrace_scope; “0” is expected, if it is “1”,

Then do

$echo 0 | sudo tee /proc/sys/kernel/yama/ptrace_scope

Took me a while to find… hope this saves someone some time.

Hash

Doing another SPH implementation for parallel computing. It’s amazing how fast adding a hash table, instead of looking at all particles does for one’s speed. 1429.22 seconds to only 327 seconds. That’s 5 times as fast!

Wine IE

  • Installs Wine to play Hearthstone.
  • Made a gif for a project.
  • Double clicked gif to see how it turned out.
  • gif opened in IE under Wine.

Life’s Update

I’ve been working on a few things:

  1. Position-based fluids for my CS5643 class, which is a pain to debug. Normally, I have tons of intuition on what numbers are suppose to be but my intuition is nil here. I’m unsure as we’re guessing the correct parameters or even the right formulas and procedure.It’s a great class, just the project are incredibly tedious. For example, my particles right now just seem to float of into the x direction and stay there. Why would it do that?

    Hopefully, what ends up happening is something like the following, except less pretty.

  2. I finally installed Hearthstone for my Elementary OS Linux box. It was incredibly easy following this tutorial.
  3. In terms of research, I’ve been reading up on different ways to generate the Krylov basis for GMRES as described in this paper. I’ll probably write up and refined the old  “idiot’s guide” again for this paper information (and fix the bad latex in there).
  4. This soundtrack (not allowed to embed).
  5. My card handeling has suffered :(.