Monte Carlo solver

Just wanted to post a quick video of a coupled Monte Carlo transport-Poisson equation solver for a 2D MESFET. The colored background is a heat map of electric field strength, while the white dots represent individual electrons being accelerated and scattered in the device. Sorry for the low quality. YouTube decided that was acceptable, despite having a higher-quality source available. Each frame is 2e-15 seconds, and the total length corresponds to 8e-12 seconds.

IPython dependencies

IPython has some nifty features, but since IPython is still in development, they are awfully hard to get running, at least on OS X. Compounding this is the dependency on Qt for the cool new features to work. Qt isn't the most seamless package to get installed and running. I'm not complaining at all. I just want to list all the steps that I went through to get the Qt console working.

1. Try starting the IPython qtconsole:
   python qtconsole
   and learn that I'm missing all kinds of dependencies.
2. Try to install pyqt using pip:
   pip install pyqt
   and learn that the egg for pyqt is broken.
3. Try to build pyqt using the files that pip installed
   python configure.py
   and learn that pyqt depends on sip.
4. Try to install sip using pip(!):
   pip install sip
   and learn that the sip egg is broken in the same way as pyqt
5. Install sip using the files that pip downloaded:
   python configure.py
   make
   make install
   This actually works.
6. Try to build pyqt again:
   python configure.py
   and learn that I need the Qt SDK. This one is my fault.
7. Install Qt from the Nokia website
8. Repeat step 6, find that for some reason qmake was not installed.
9. Uninstall Qt
10. Uninstall sip
11. Install Homebrew
12. Try to install pyqt:
    brew install pyqt
    which successfully installs Qt, but fails because it can't link to the right version of Python, which I have somehow accumulated 4 different installations of
13. Install python through homebrew:
    brew install python
    This fails sip has hard coded the location of the old, broken framework.
14. Uninstall sip.
    brew remove sip
15. Try to install pyqt:
    brew install pyqt
    This still fails. No idea why.

MATLAB may not be completely terrible

I use MATLAB a lot, but I wish I could use it less. I have a few reasons for not liking it:

• I prefer to use VIM and the command line for everything. The MATLAB IDE isn't my favorite.
• It's closed source (not optimal, but not necessarily a deal-breaker) and expensive (unless you're a student).
• Lots of useful code is sold as addons. The parallel computing toolbox, for example, is not included in the standard version of MATLAB.
• The language itself isn't as nice as Python is, at least for me. Language choice is a highly personal matter, however.

For these reasons, I am trying to install Numeric Python as a replacement for MATLAB. So far, it's a huge struggle. The value in having a working product out of the box is not insignificant. Python support on OS X isn't as strong as could be hoped, and the recent upgrade to Lion has also broken some dependencies. I am eschewing pre-built binaries for locally compiled versions to get around these issues. I'll update this post with resources.

Rookie mistakes

Working on some Poisson solver code today. Made a rookie mistake: used 'x' and 'X' as variables within the same segment of code. They are even related: X is the extent of the device I am simulating, while x is the current position along the device within a loop. I rewrote ~300 lines of code trying to figure out what had gone wrong. Nothing worked.

INCORREC

Then I walked away from the computer and thought about the output. I was seeing this, which isn't totally wrong, but is wrong enough that I was losing my mind trying to figure out.

The key to finding my error was that the error seemed to increase on the right hand side. You can see what look like some 'tearing' up near the top of the image on the right hand side.  The error seems to increase left to right.

CORRECT

The answer isn't that shocking: instead of looping, I was using a fixed value (x versus X), which was causing wrong values to be substituted. Using the right variable yields the right answer. Much better. The Poisson equation is solved!

JavaScript is the scourge of the universe

I want to win my fantasy baseball league, but I think ESPN's interface leaves something to be desired. I want to take the daily stats from my team and put it into a more useful format. Sadly, the hardest part is getting the data.

ESPN doesn't want any of the raw statistics to get out, I suppose for good reason. However, I really just want a daily feed of how each player performed, points-wise. ESPN has a simple web page to handle this, which it terms the 'quick box score'. Even better, these box scores aren't protected by the authentication interface that the rest of the fantasy baseball functions are behind, so you can just query a pretty straight-forward URL and get the quick box score for any day of scoring. The first day of play (yesterday) is day 1, and ESPN simply increments that paramater in the URL each day. These box scores contain nicely formatted tables with all the data I want, as shown.

Sadly, however, these pages are not generated by the server. Instead, some sort of AJAX query is used to generate them in the browser. Unfortunately, the code that underpins ESPN's site is as ugly as the site itself. Looking at the source reveals piles of indecipherable JavaScript to generate the tables. This is a problem because any automated method of fetching a box score won't be able to parse the JavaScript. Thus, none of the important data will ever load.

And that's pretty much where I gave up. I know I can have my script open the page in a GUI web browser, then save the resulting static page back out and have the script continue working, but that amount of interactivity is too much. I really wanted the script to handle this for me, running the stats once a day and generating the report I want. I will work on this more eventually, but for now I just don't want to have to deal with all of these terrible web issues. Unless someone else has a suggestion for parsing the JavaScript, I'll leave this alone for a few days.

Saguaro Cluster

So this Monte Carlo simulation is pretty cool, but extremely computationally expensive (compared to a normal analytical drift-diffusion simulation, that is). I'm going to try to get some time on the ASU Saguaro Cluster to run the simulator because the code is extremely easy to parallelize. If I have an ensemble of N particles, each particle is independent of the others, so for large values of N, throwing more cores at it results in a roughly linear speedup. In my case, N=5000 or 10000, so that's good. A simulation that takes 50 minutes on my computer should be done much, much faster on the cluster. I'll let you know how it goes.

Optimization

I've been working on a Monte Carlo method particle simulator for a class I'm taking. I've been wondering about the most efficient way to write it. I know premature optimization is the end of all things, but I can still see a few problems ahead.

Right now, I am using MATLAB. I have decided to go the object-oriented route, because each particle has a number of properties that I'd like to keep together for clarity's sake, such as energy and momentum. So, I defined a class to contain all of these variables. The program needs to track a fairly large number of particles (~10000 or so) in order for the statistics of the Monte Carlo method to work out. The question I have regards the efficiency of the numerical methods applied.

As I understand it, most linear algebra libraries understand that if a bunch of independent operations can be applied in parallel to a large set of data, those operations will be applied as efficiently as possible, using multiple cores as needed. However, there is a chance that the library may not recognize my code as being a case where parallel execution is possible. In other words, if I had an array of numbers that I wanted to perform an operation on, the library would recognize that it could operate on the array in a parallel manner. However, if I have an array of objects, and want to operate on all the members of those classes, the library may not optimize this in the same way.

Basically, I wonder if MATLAB will know to treat my array of particles in an efficient method. I would like to rewrite this simulator in Python for code clarity anyway. Wonder how the linear algebra methods in NumPy would handle this.

Euler 007

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.
What is the 10001st prime number?

My solution below.

Euler 006

Next!

The sum of the squares of the first ten natural numbers is,
1² + 2² + ... + 10² = 385
The square of the sum of the first ten natural numbers is,
(1 + 2 + ... + 10)² = 55² = 3025
Hence the difference between the sum of the squares of the first ten natural numbers and the square of the sum is 3025 - 385 = 2640.
Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

Didn't brute force this one. Read about math instead: 

Summation

Square Pyramidal Numbers

Going to move the solutions inside the fold.

Euler 005

Sorry for three in one day, but this one was awesome in Haskell.

2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.
What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?

Neil, you might punch me for this one.

Euler 004

Problem 004: 

A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 X 99.
Find the largest palindrome made from the product of two 3-digit numbers.

BOOM:

Euler 003

The prime factors of 13195 are 5, 7, 13 and 29.
What is the largest prime factor of the number 600851475143 ?

My response:

Euler 002

Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.

My solution, in Haskell:

Euler 001

The first Project Euler problem asks

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.
Find the sum of all the multiples of 3 or 5 below 1000.

Here's my solution, implemented in Haskell:

My own Watson?

So, IBM has told us all the secrets to building Watson. How hard could it be?

OK, pretty hard. But the whole system was built out of some fundamental building blocks. From a high level perspective, the most interesting building block to me is the Apache UIMA. This is the engine that parses natural language documents, looking for facts. IBM used some pretty heavy duty annotators to find the things they were looking for, but the default annotators are pretty neat to begin with. A short example page shows off the general idea.

This raises the question 'how can I use this?' I was trying to think of a novel, game-changing application for this, but then I realized I wasn't a computer scientist, so I should probably think a little bit smaller for now. I'd like to know how many unique names come up in every paper ever published in an IEEE journal. This is probably still unreasonable, because I can't just download all of the IEEE corpus, as cool as that would be. ASU probably wouldn't enjoy me abusing their library access by screen scraping all the PDFs, either. New plan needed. Thinking even smaller, I attended an IEEE conference last summer that just happened to provide all of the proceedings of the conference for the last 30 years on a CD. This might be a reasonable amount of data to start with.

This is another one of those projects that will probably be slow to happen, but it doesn't seem insurmountable. The hard work is all done. I just have to put other people's work to use. I'll let you know if I make any progress.

functional programming

So I've decided to think smaller. Reading up on MapReduce has made me interested in learning a functional programming language, at least a little. What is a functional programming language? Well, if I knew, I'd tell you. I have a very vague idea in my head, but it's so vague I can't even delude myself into thinking I know what it is.

I'm going to need some projects to learn this new method, however. Luckily, Project Euler has a large selection of simple problems to work on. They increase in difficulty monotonically. Eventually I'd like to solve them all, but for now I'm going to start at the beginning.

I'm still broke, so I'll have to use online resources for the time (which might be the best kind anyway), but if anyone has any suggestions for learning Haskell, let me know. I plan on trying to solve each problem using both a functional and an imperative language. This should really show the differences between the two approaches. In the end, maybe this exercise will result in a better understanding of the functional concepts that underlie something like the MapReduce juggernaut.

Data!

The next thing I plan on doing, when I have time, is learning something about http://en.wikipedia.org/wiki/MapReduce. To do this, I'm going to need a lot of data. Luckily, someone has collected a bunch of publicly-available datasets. I'm not sure which one sounds coolest. I realize I may be getting ahead of the game here, since I have absolutely no experience with this sort of programming, but I can at least ask for suggestions here.

So, anyone seen any fun data sets in there? I feel BLS data might be sort of interesting. I could use it to prove all sorts of irresponsible political points. It's hard to argue with a few hundred megs of raw numbers. I like the idea of the Usenet corpus from 2005-2009, but it's 28 gigs. Woof.

I really like the idea of this kind of massively parallel data manipulation. Had I known this kind of stuff existed back in high school, maybe I would have done more with computer science in college. Never too late to teach myself, though.

Debugging is hard sometimes

So my last post finished with me asserting, boldly, that there was a minor issues with projecting coordinates that was causing my fractals to come out warped. I went to bed thinking that there'd be no problem. I spent about 6 hours today writing and re-writing projection code because I thought for sure that that was the part of the code that was causing me problems.

I was wrong. It turns out that unlike a <div>, <canvas> tags do not like to have their dimensions set through CSS attributes. Rather, they are sized directly using HTML attributes. It's the difference between

<canvas id="fracMap" style="width:512px; height:512px;"></canvas>

and

<canvas id="fracMap" width="512" height="512"></canvas>

Switching from the first form to the second cleared up my problem. Sadly, this was the very last thing I tried (Obviously. You always find your keys in the last place you look, right?). In the meantime, though, I half-heartedly cleaned up all of the rest of the code. I haven't written this much code in a while, and it feels good. I will never love JavaScript, however. Google has done some pretty awesome things with such an eccentric language. Whenever I move to independent hosting, I'll upload what I have so far.

 

Fractal coding

The fractal coding is actually going quite well so far. I've implemented the code to use the canvas tag to draw the Mandelbrot set. Blogger doesn't really offer an easy way for me to upload the JavaScript, so I'll just attach a screenshot.

I've also started working on the Google maps interface for zooming and panning. I've replaced the tile getting code with my fractal drawing code, but it still needs more work. I think I'm having issues mapping between two slightly different coordinate spaces. Luckily, there's a way to mess with the projection function, which should allow for an easy mapping between the Google coordinate system and the Mandelbrot set coordinate system.

Hopefully I can finish this in the next few days. I have some other engineering work to do, but this is much more fun right now.

Well, crap

I should have worked harder, because Google beat me to it.