Author: Marshall Jiang

Posted on

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…

Posted on

Submitted!

Good news; my advisor and I submitted our first paper this past Monday. I’ve uploaded a draft here for those curious.

 

Posted on

Privilege

NYTimes Article

I think I always forget how privileged I am. I benefited from many points laid out in the article by Mr. Levy: my parents explicitly moved to one of the best school zones in Tallahassee just for me to get a better shot at college applications, I also had the benefit of them expediting the citizenship process so I could apply as an American citizen, and they provided enough money for all the extracurricular I needed. All this after living in graduate student housing for 2 years.

I like to think that I can do what I am doing without my parents’ support, but the truth is I can’t. So much of success is dependent on having a good family. I don’t think any of the people currently in my graduate program came from an under-privileged home… and to be honest, all of my friends are quite affluent.

Do I agree with the article? Yes. Do I think it’s enough? No. College is the end of the road which starts from pre-K. Certainly universities can offer more aide, or allocate more slots for under-represented students, but will this really pull them up? I don’t think so; the system right now is too fundamentally broken for a top-down approach.

Posted on

Some Updates

  1. Went to RPI applied math days this weekend, and gave a talk about preconditioning. It went pretty darn well, but I really do have to feel bad about the international grad students who don’t have a mastery over the English language… because they can’t bullshit out of anything like I can.
  2. On the way there, I saw for the first time a cop pulling over someone going too slow in the left lane. I rejoiced.
  3. On the way back from the gym, I saw a cancer patient in a car by the traffic lights. Chemo really does take the life out of people’s eyes.
  4. Being ghosted sucks ass, especially when you thought you clicked well with the other person.
  5. Taxes suck ass.
Posted on

Subsequences of Heights

I really liked this problem that my friend Shamile gave during dinner:

Consider 50 people lined up in a random order. Show that there exists a subsequence of length 8 such that the height is non-increasing or non-decreasing.

The solution is quite elegant. Assign an ordered pair (a_i, b_i) to person i in the line such that the first coordinate is the length of the longest subsequence which is increasing from person i, and the second coordinate is the length of the longest subsequence which is decreasing from the same person.

Let’s assume by contradiction that the longest such length is less than 8; this implies that the set of acceptable coordinates is from (0, 0) to (7,7). There is only 49 total combinations, so by pigeonhole principle, there exists a pair of the same number. But this cannot be! The rest is left to the reader.

Posted on

Gradient of the Barycentric Coordinate in 2D

I always forget this fact, so hopefully typing this out will help.

Let T be a triangle with vertices v_1, v_2, v_3, then there exists unique linear functions \lambda_i such that \lambda_i(v_j) = \delta_{ij}. This is the so called barycentric coordinates.

The fact here is that \nabla \lambda_k = -\frac{|\gamma_k|}{2|T|}n_k where |\gamma_k| is the length of the edge opposite vertex v_k, and |T| is the area of the triangle. The vector n_k is the unit outward normal.

The derivation is simple. The direction of the gradient is towards the highest growth, and as the barycentric coordinates are linear functions, it’s clear that -n_k is the correct direction. The scaling is to simply note the area as d|\gamma_k| = 2 |T| where d is the shortest distance from the edge to the vertex. This gives us the inverse of the slop as the barycentric coordinate has a value of 1 at the vertex.

Posted on

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.



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

Quadrature Rules on Triangles

Quadrature rules on the interval [-1, 1] is pretty well-documented in many different textbooks (see for example Cheney and Kincaid). By transformations, we can easily integrate across the whole real line spectrum. Going into two dimensions, we can take tensor products to integrate rectangles easily. On the other hand, integration on the 2D simplex (i.e. triangles) is a bit harder. There really exists two types of quadrature rules for the triangle:

  1. Tensor products + transformations
  2. “Symmetric quadrature rules”

The latter offers better convergence with smaller points, and has better spread-out points than the tensor product version, but offers less technical support through software packages; we refer the reader to the following paper. We will discuss the first, and exposit a little more; I will be borrowing heavily from Karniadakis and Sherwin’s fluid dynamics textbook.

Suppose our triangle at (-1, -1), (1, -1), and (-1, 1) (so called reference triangle) uses x, y as it’s coordinates. If we use the change of coordinates u = 2\frac{1 + x}{1 - y} - 1, v = y, it’ll map the triangle to the square defined by (-1, -1), (1, -1), (-1, 1) and (1,1). Hence, we can express the integral using a change of variables (details omitted)

\int_T f(x, y) \, dx dy = \int_{-1}^1 \int_{-1}^1 f(u, v) \frac{1 - v}{2} \, du dv.

We note that the weird factor is simply the change of variables Jacobian. Now it’s easy to use Gaussian quadrature on this! But there are better version, as we have the Jacobian term; looking up the quadrature rules, we see that Gauss-Jacobi fits our bill. We refer the reader to page 144 of the Sherwin textbook to see the full-details.

A simple implementation is below (this is far from the best). I guess I’ll talk about this at the end of the week.

from scipy.special import roots_jacobi as rj
def quadtriangle(f, w00, w10, x00, x10):
    total = 0
    for i in range(len(w00)):
        subsum = 0
        for j in range(len(w10)):
            subsum += w10[j] / 2 * f([(1 + x00[i]) * (1 - x10[j]) / 2 - 1, x10[j]])
        total += w00[i] * subsum
    return total
def f(x):
    return x[0]**2 + x[1]**2
p = 4
x00, w00 = rj(p + 1, 0, 0)
x10, w10 = rj(p + 1, 1, 0)
print(quadtriangle(f, w00, w10, x00, x10))
Quadrature rules