Linear algebra courses at MIT from Prof. Gilbert Strang

Linear algebra is fundamental mathematical knowledge for those who need to perform computational natural sciences. It is a neat formalism to express things in a compact way, and describe precious algorithms to solve computational problems from chemistry, physics, astronomy, and so on.

I found these precious and very clear lectures from MIT professor Gilbert Strang. The lectures are available to the general public under the MIT OpenCourseWare distant learning initiative. I am pretty sure not to say anything excessive when I claim that such initiatives should be declared Intangible Heritage of Humanity due to their important role in scientific advancement and education.

Here is the list of the lectures, linked to the proper material on the MIT OCW website:

  1. The geometry of linear equations
  2. Elimination with matrices
  3. Multiplication and inverse matrices
  4. Factorization into A=LU
  5. Transposes, permutation, spaces R^n
  6. Column space and nullspace
  7. Solving Ax=0: pivot variables, special solutions
  8. Solving Ax=b: row reduced form R
  9. Independence, basis and dimension
  10. The four fundamental subspaces
  11. Matrix spaces; rank 1; small world graphs
  12. Graphs, networks, incidence matrices
  13. Quiz 1 review
  14. Orthogonal vectors and subspaces
  15. Projections onto subspaces
  16. Projection matrices and least squares
  17. Orthogonal matrices and Gram-Schmidt
  18. Properties of determinants
  19. Determinant formulas and cofactors
  20. Cramer’s rule, inverse matrix, and volume
  21. Eigenvalues and Eigenvectors
  22. Diagonalization and powers of A
  23. Differential equations and exp(At)
  24. Markov matrices; Fourier seriesQuiz 2 review
  25. Symmetric matrices and positive definiteness
  26. Complex matrices; fast Fourier transform
  27. Positive definite matrices and minima
  28. Similar matrices and Jordan form
  29. Singular value decomposition
  30. Linear transformations and their matrices
  31. Change of basis; image compression
  32. Quiz 3 review
  33. Left and right inverses; pseudoinverse
  34. Final course review

The third eye of the iguana

Did you know that iguanas (and not only them) have three eyes ?

the feature is called Parietal eye, it is highly advantageous to spot the shadow of a predator coming from above, and it’s not unique to iguanas, but it can be found in frogs, lizards, sharks and other reptiles and fishes. In the Tuatara the eye is way more developed, but it regresses as the animal gets older. All these eyes can influence the pineal gland, apparently responsible for the circadian rhythm regulation.

Offsetting carbon emissions with trees: a little math

I want to do a thought experiment, hoping not to make any mistake.

The world produces 29,888,121,000 metric tons of CO2 per year. According to Wikipedia, one ton of dry wood sequestrates 1.8 tons of CO2. If you put these two together, you quickly realize that you need to produce 16,604,511,666 tons of wood per year in order to offset the worldwide CO2 production.

The largest single stem tree available today is the General Sherman giant sequoia. It is 83.8 meters tall, with a diameter of 7 meters. Its estimated dry weight is 1121 tons. It took thousands of years to reach that size, and it’s a unique tree, but I will use it as a greatly exaggerated upper limit to the size and mass of a tree. I estimate a normal tree as a fifth of that, or 250 tons.  In order to compensate the CO2 emission, we need to plant and grow to full size 66 millions trees. Assuming the tree occupies an area of 10 x 10 square meter this adds up to a required forest size of 6600 square kilometers, which is one sixth of Switzerland. Every year.

Of course, I am not considering the fact that a tree must grow to full size in order to offset the CO2. If we do, and if we assume it takes 20 years for a tree to reach the size given above, we need to plant a little more than three Switzerlands of trees right now, to offset the CO2 production of the next 20 years, for a grand total of 1.3 billion trees. Considering the population of Europe is 852 million people, it means that each of us needs to plant a little more than one tree.

I’m working on that. And you ?

 

A raytracer in python – part 2: rendering multiple objects

A quick addition needed to the raytracer is providing freedom to add more objects to the rendering scene. In Part 1, the design was such that only one object, a sphere, could be drawn. The new code allows much more flexibility. I added a Plane object, introduced assignment of colors to the objects, divided the source into multiple files, and fixed a bug relative to rendering direction. Let’s see more specifically.

The class World is the main interface to user programming of the raytracer. This small program creates and renders a scene containing four spheres, a white one in the origin, the others along the axes, each with different colors.

import raytrace
from raytrace import objects

w=raytrace.World()
w.add_object(objects.Sphere(center=(0.0,0.0,0.0),
                            radius=10.0,
                            color=(1.0,1.0,1.0)
                           )
            )
w.add_object(objects.Sphere(center=(50.0,0.0,0.0),
                            radius=10.0,
                            color=(1.0,0.0,0.0)
                           )
            )
w.add_object(objects.Sphere(center=(0.0,50.0,0.0),
                            radius=10.0,
                            color=(0.0,1.0,0.0)
                           )
            )
w.add_object(objects.Sphere(center=(0.0,0.0,50.0),
                            radius=10.0,
                            color=(0.0,0.0,1.0)
                            )
            )
w.render()

The resulting image is the following

Clearly, the white sphere is not visible, as it’s hidden by the blue sphere. This test allowed me to discover a problem with orientation: the green sphere was on the wrong side. I therefore had to analyze a bit the orientation and the different coordinate systems into play here

  1. The measure of the screen is given in pixels. We are used to this when it comes to screen size, for example. An image which is 320×200 means that it’s 320 pixels wide (horizontal resolution) and 200 pixel high (vertical resolution). Don’t fall into “thinking matrix”, where NxM means N rows x M columns. It is the exact opposite.
  2. The geometry uses the cartesian system, which has the origin in the center of the picture. The x axis is oriented towards the right, the y axis towards the top, and the z axis towards the observer. This is not different from a traditional cartesian layout: for the z=0 plane, the top left pixel correspond to a (-,+) coordinate, the bottom right to a (+,-) coordinate. Points closer to the observer have positive z, honoring the right hand system. The camera in the above picture is at z=+100.0
  3. The raytracer uses pixel coordinates for the viewplane. Pixel 0,0 is at the bottom left. Changing the first index moves horizontally (along the row) from left to right. Changing the second index moves vertically (along the column) from bottom to top. As a consequence, the point at the bottom right is (hres-1,0), and at the top left is (o, vres-1). Note that this is equivalent to a cartesian (x,y) system with origin on the bottom left corner. Not a surprise, since there is a direct mapping between the pixels and the rays’ origins.
  4. Finally, the pixel image indexing of pygame and PIL. For them, pixel 0,0 is top left. Like the case above, incrementing the first index also moves along the row from left to right. However, incrementing the second index moves vertically from top to bottom, which is the opposite of the raytracing index. Bottom left is (0, vres-1) and the bottom right is (hres-1, vres-1).

A remapping is therefore needed from the pixel coordinate of the rendering and the pixel coordinate of the display (e.g. pygame). The transformation is trivial, of course, but it must be kept into account, otherwise pics will be flipped horizontally.

Another interesting fact is that, according to “Ray Tracing from the ground up” it’s useful to perform the raytracing operation starting from the bottom and working our way up, rendering pixels from left to right. According to them this is for coding symmetry and future convenience, so we stick to it.

Finding the foremost object

In the book, finding the foremost object is made more complex by the language used, C++. In python, you can use functional programming style to obtain the same in a very concise statement. The idea is to cast a ray, then go through all objects in the world to find the intersection point (if any) between the ray and the object. If more than one object is hit, the one with the intersection point closer to the observer commands the pixel color, since it’s in front.

I achieve this with the following code

def hit_bare_bones_object(self,ray):
   def f(o):
     shadeRec = o.hit(ray)
     if shadeRec:
       return (shadeRec.parameter, o)
     else:
       return None

   try:
     foremost=sorted( \
                filter(lambda x: x is not None, \
                  map(f, self.objects)
                ), key=lambda x: x[0]
              )[0][1]
   except IndexError:
     return None

   return foremost

What does this code do ? I defined a simple internal function f which accepts an object, performs the hit and returns a tuple containing the hit point position (as a parameter of the ray, so it’s a single number, not a xyz coordinate) and the object.

Now, I use this function to map all the objects defined in the world. I will obtain a list with one entry per each object, either a None (not hit) or a 2-tuple containing the parameter and the hit object. I filter out the None entries, leaving only the 2-tuples and then sort according to their first element. The 2-tuple with the lowest parameter is now at index 0, and the [1] element of this tuple is the foremost object. At any time, the list may be empty (such as if you don’t have any object,  or no object is hit. In that case, a IndexError will be raised and that will indicate that the ray hit nothing. I may rework on this function later on, but for this second round, it suits my needs.

It’s now time to move on to samplers. Given that the code is growing in size, I created a git repository you can clone from. The release of this post is available here. The code is under BSD license.

Robotic overlord domination!

Remember when I posted about the lovely but fake Titan robot, and I told you the future was close? Well, some time passed, and the future is now even closer. New robots are among us, and they are learning. Some of them prefer to stay at home, folding our clothes, washing our dishes or bringing drinks and staying with us, acting as tables. If you are at home wasting time in front of the computer, you could just ask one to bring you a beer.

Some of them decided to find a job outside, so they pick strawberries, administer weddings, or wash our hair with the power of sixteen fingers. If we are sick, they may go to school for us, and relay the lesson. Sometimes during breaks they can lose themselves in the big city, so you could meet one of them asking for directions. Maybe they will take a bike ride to get to the pub, and play pool.

Robots can also act as our trusted pet friends, and although they will never replace them, they can become quite interesting marvels of technology. Sony Aibo was the first robot dog commercially available. I played with one, and it’s an amazing little thingie. These two, however, bring the doggy robot to a whole new level: small dog or big dog, the choice is yours. And speaking about animal robots, you can also have a tree-climbing snake, a crawling snake, a bug-like thingie that both walks and rolls, a fish, a spider.

They are also into performing arts. They can dance the Bolero, and they are damn good at it (really, take your time with that video, it’s pure mechanical beauty), or some other dance maybe? With violin they are not yet Uto Ughi, but they are improving. With trumpet they are already pretty good, apparently. As for sport and outdoor, they can learn archery skills, or and some of them very soon will ask us to become a real boy.

Why I won’t talk about neutrinos

Because I don’t know enough on the subject to say anything about it. What I do know however is that the finding is exciting, regardless of the final result. Science thrives in being challenged, and shows humility in front of facts.

Switching to a once-a-month schedule

During the course of the last years, I was able to keep a twice-a-month post schedule. In the following months however, I will relocate, my work will become more time consuming, posts are becoming more demanding to be produced (as I also have to write code for them) and my interests are spreading to other fields. Hence, I decided to reduce the throughput that devours my scheduled posts queue, which, at the moment and to be fair, is rather long.

Note that this does not mean that I will post once a month. It means that I will have a scheduled post once a month, instead of two. Shorter, “spur of the moment” posts may still happen.

As you noticed, I started a series of posts on raytracing. This will carry on for a while, depending on the depth level I will be willing to reach. I still have to post the last “molecule that changed the rule of the games” and I say, don’t worry, it will come.

A raytracer in python – part 1: basic functionality

Some time ago I visited Sydney, and I made a tragic mistake: I entered the University bookshop. Why a mistake, you say? I am book maniac. As soon as I enter a book shop (live or on web) I end up spending up to a thousands euro every time. This time, it was not the credit card, but the limit of 20 kg on my luggage to put a limit on what I could buy. Needless to say I got really good stuff, in particular this 761 pages of pure awesomeness: “Ray tracing from the ground up”. This book teaches you how to write a raytracer, step by step from a first basic skeleton to incredibly complex optic effects. It also provides code, available for download under GPL. I haven’t downloaded it, but according to the book it’s in C++. Since I want to keep myself fresh with python, I will write my own version in this language.

The first objection I may hear is performance. Native python is much slower than C++. I agree, and I do it on purpose. My objective is in fact not only to write a raytracer for fun (and the world is full of already awesome raytracers), but also to perform some infrequent python exercises: interfacing with C, parallelization and, hopefully, some OpenCL programming. This is my hope, at least. I will try to spend some time on it and see how far I can get.

What is raytracing and how does it work ?

Raytracing is a technique to produce a photorealistic image. It works by projecting rays from the observer to the scene, and coloring pixels on a viewplane for every ray that intersects an object.

This mechanism resembles how vision works, although in the opposite direction. Light rays from a lamp hit objects and their reflection happens to scatter around. Some of these reflections will enter our eyes and allow us to see the world around us. Doing the opposite, tracing from the observer, is clearly more efficient as we don’t care about the rays not hitting the observer (at least in its basic implementation), only those who do.

Performing raytracing (or to be more accurate, for now just its basic form raycasting)  involves the following steps:

  1. define a geometric object in space, like for example a sphere
  2. define a view panel made of pixels
  3. shoot one straight line (ray) from the center of each pixel
  4. if the ray intersects the object, mark the pixel colored, otherwise mark it with the background color

That’s it. Basically, it’s an exercise in geometry: finding intersections between lines and 3D objects in 3D space. This first program, ray.py, does exactly that. You will need to install the Python Imaging Library, pygame and numpy. The result is intriguing:

A rendered sphere

Ok, I have a very loose definition of “intriguing”, but it’s a start.

Pac-mecium and other games

This is a rather peculiar use of paramecia: by directing their preferential swimming movement by means of an applied electric field (a property called galvanotaxis) we can use them as player characters in videogames. I don’t expect this to become common gaming, though.

Computational chemistry development in research

Imagine you are a professor in organic chemistry. You received financial support for a project, and you are ready to hire a Ph.D. student to make it happen. The project requires the synthesis of a new compound.

Imagine you interview your best candidate. At the whiteboard, you present him with various problems of how to synthesize different products, and you find he is very qualified in planning these strategies. He also knows how to estimate the yields of the products, and how much initial reagents are needed. He knows temperatures and catalyzers required, as well as conditions that may break down the product. He can also write excellent articles. Seems a very promising person, so you hire him.

After a few weeks, he joins the group and enters your lab for the first time. Here you find out the following about him:

  • It is the first, maybe second time he enters into a lab.
  • He does not know the names of lab glassware. He never heard terms like “flask” or “test tube”. He actually looks at you with a strange face when presented with them, and say that the only time he did an experiment during his course, he used a pan and it was enough.
  • He is unable to keep his bench and glassware clean, and most often than not, he throws away used glassware instead of washing it.
  • He constantly breaks glassware, scattering glass and reagents around, damaging the work of his colleagues who have to clean his mess or just deal with it in the most improbable ways.
  • He is unable to assemble glassware (such as he is unable to build a distillation setup). He puts clamps in the wrong way, does not put silicone grease in the joints, nor any clamp to keep these joins sealed to prevent “popping” due to pressure buildup.
  • He never performed a chromatography, nor an Infrared Spectrum. He just vaguely heard about them, but his course professor said that they are useless and you can identify the substance with a sniff, which is enough most of the time.
  • He pretends to prepare and use his own very impure reagents, even when high-purity, standard reagents are available from a reagent reseller. He claims that it takes too much time to call them and make the order. Also, he is unable to understand the codes and numbers written on the bottle’s label.
  • He labels his flasks with cryptic names, and pass them to colleagues who constantly have to ask him what they contain. Most often, he does not remember and have to check on his notebook. When he remembers, it turns out that the product inside is either polluted or damaged.
  • He disposes of his byproducts in the drain, instead of using the proper disposal units. When addressed about it, he says “who cares? it takes less trouble”
  • When asked about his low quality work, he claims that it’s not his task nor its core expertise to produce top quality laboratory work, nor to know how to use laboratory equipment, as long as he can manage to get the required product at the end. He also points you at a very efficient synthetic strategy he just devised at the whiteboard, claiming that he is doing an excellent work.
  • At the end of his employment, he leaves a notebook containing the process to create the product. The notebook is completely disorganized, the pages are mixed and not numbered. The handwriting is close to impossible to understand, it is written in four different languages, and the quantities are specified as “a pinch of”. However, the glassware setup he made and left, when fed with the contents of an unlabeled reagent bottle you find in a bench downstairs, gives you the product, but only if the humidity in the laboratory is at exactly 75 %.

My question for you, the professor of the group is the following: would you keep this person in your group, or would you dismiss him?

Same story, different chemist

The scenario presented above is considered the norm in a different branch of chemistry, theoretical and computational chemistry. While “white coat” chemists use whiteboard synthesis strategies, a laboratory, glassware tools, spectroscopic instruments, chromatography and distillations, theoretical and computational chemists use mathematics, computers, editors and software development. Both these sets are tools for the job, and in order to practice the discipline, proficiency with them is expected. Yet, in the discipline of theoretical and computational chemistry, the general level of competence and proficiency with its tools can be rightfully compared one-to-one with the scenario given above for the case of organic chemistry.

Organic chemists are supposed to be able to use a laboratory environment with high quality standards and protocols. By analogy, computational chemists are supposed to be able to use a software development environment to high quality standards and protocols. In both cases, to perform their research they should be required to be masters of the basic tools, and proficient in a broad set of specialized and modern tools. It is inexcusable not to be.