Showing posts with label python. Show all posts
Showing posts with label python. Show all posts

Sunday, April 02, 2017

Testing Python Unicode Handling the Lazy Way

Background

At work, we've been working on a Python API for IBM's Spectrum Scale (GPFS) filesystem.

The API is well tested, with good coverage. But one thing that's caught us out a couple of times is unicode.

For those not familiar, Python 2.7 has both byte strings and unicode strings. In most cases, the two are interchangeable. But there are some quirks, which can lead to bugs if not handled correctly.

Now it's usually safe to assume that users will work with byte strings - mostly because when you create a string, that is the default format. And internally, we always use byte strings.

But this being an API, it gets used by other people. And what I've come to learn is users do all sorts of weird things. The amount of bugs we've had from frankly bizarre filenames... (why would anyone put a newline character it a file name?!)

And users aside, unicode can also come from interfacing with other code - for example, the json.loads method.

All of which is to say, realistically, it's best not to ignore unicode.


Unittesting

Listen, I understand the value and importance of unittesting. I just find it so dull compared to writing 'proper' code. And frustrating! In my experience, most of the errors and fails raised by unittests come from the tests themselves, rather than the code they're supported to be testing. (Maybe that's on me).

So yeah. The thought of having to write a whole bunch of new unittests - most of which were going to be more or less exact copies of the tests that already existed - didn't appeal to me.

Instead, I wondered if maybe there was a lazier way of doing it.

The easiest way to do that was to take advantage of the tests that we'd already written. All we needed was a way to run those tests with some mechanism for convert any strings passed to API function calls to unicode.


Set Trace

One thing I particularly like about Python is the ability to monkey patch things - playing with existing code on the fly, without having to change the actual source code.

That's why the sys.settrace function is one of my favourites.

I came across set trace when I was looking at ways to do logging in the API. We wanted to be able to log function calls, functions returns, etc. The problem was, we didn't really want to add individual logging statements to each and every function. Aside from the effort, it'd really clutter the code.

After considering a couple other options (decorators, metaclasses), I came across settrace.

settrace isn't well explained in the official docs. You can find a better explanation here.

settrace allows you to set a trace function - a function which will be called for every line of code run, for every function call and return, for every exception - perfect!

A trace function receives 3 arguments - event, frame, and args.

event is a string indicating what the code is doing (line, call, return).

args is some arguments, which vary depending on the event - for example, if event is 'return', then args will hold the function return value.

And frame is a frame object. This is where we get most of our information for reporting.

The frame object holds some interesting stuff - the f_locals attribute holds the frame locals, and in the case of a call event, these locals are any variables that have been passed into the function.

There's also f_code - the code object for the function being called. And from that we can get things like f_code.co_name - the name of the function being called/returned from.

So as a simple example we might have

import sys

def trace(frame, event, args):
    if event == "call":
        print frame.f_code.co_name, "called with", frame.f_locals
    elif event == "return":
        print frame.f_code.co_name, "returned", args
    return trace

sys.settrace(trace)

I ended up using settrace to write an 'strace' style script, which can be used to trace API function calls for a given script or piece of code. Which is pretty cool.


The Solution

So how does this apply to the unicode problem?

As mentioned above, we can get the parameters passed to the function from frame.f_locals. And because f_locals is a dict, it's mutable. That means that we can change it's values, and those changes will persist when the function being traced continues executing.

This is how this solution works - we convert any strings in f_locals to unicode. The code being 'traced' then behaves as if its functions had been passed unicode to begin with.

While we're at it, we have to make sure we also convert any strings in lists, tuples, dicts - in particular because *args and **kwargs are ultimately just a tuple and a dict.

Here's the complete solution

"""Unittest wrapper, which converts strings to unicode.

Check that your code can handle unicode input
without having to write new unittests.

Usage is identical to unittest:

$ python -m unicodetest tests.unit.test_whatever
"""
import atexit
import sys

# mimic the behaviour of unittest/__main__.py
from unittest.main import main, TestProgram, USAGE_AS_MAIN
TestProgram.USAGE = USAGE_AS_MAIN


def unicodify(value):
    """Convert strings to unicode.

    If value is a collection, its members
    will be recursively unicodified.
    """
    if isinstance(value, str):
        return unicode(value)
    if type(value) is dict:
        return {k: unicodify(v) for k, v in value.iteritems()}
    if type(value) in (list, tuple, set):
        return type(value)(unicodify(v) for v in value)
    return value


def unicoder(frame, event, args):
    """For all function calls, convert any string args to unicode."""
    if event == "call":
        for k, v in frame.f_locals.iteritems():
            frame.f_locals[k] = unicodify(v)
    return unicoder


if __name__ == '__main__':
    # make sure unicoder is disabled at exit
    atexit.register(lambda: sys.settrace(None))

    # activate unicoder
    sys.settrace(unicoder)

    # run unittests
    # cli args are passed thru to unittest
    # so the usage is identical
    sys.argv[0] = "python -m unicodetest"
    main(module=None)


I decided mimic the unicode/__main__.py script, so that it works as a drop-in replacement for the Python unittest module - e.g.

$ python -m unicodetest discover -v tests/

This sets the trace to the unicoder function, then calls the usual unittest method to run whatever tests we have pre-written.


Dummy Test

$ cat test_dummy.py

from unittest import TestCase

def dummy(value):
    return value

class Test_Type(TestCase):

    def test_string_type(self):
        self.assertIsInstance(dummy('foo'), unicode)

$ python -m unittest test_dummy

test_string_type (testtype.Test_Type) ... FAIL

======================================================================
FAIL: test_string_type (testtype.Test_Type)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "testtype.py", line 13, in test_string_type
    self.assertIsInstance(dummy('foo'), unicode)
AssertionError: 'foo' is not an instance of <type 'unicode'>

----------------------------------------------------------------------
Ran 1 test in 0.001s

FAILED (failures=1)

$ python -m unicodetest test_dummy

test_string_type (testtype.Test_Type) ... ok

----------------------------------------------------------------------
Ran 1 test in 0.002s

OK


Comments and Caveats

In unicodify I've used  

if type(value) is dict

PEP8 recommends instead using

if isinstance(value, dict)

But this caused issues for us. In the API, we use OrderedDict as one of the base classes for collection objects. But their init functions doesn't have the same function signature as a dict (ordered or otherwise).

So using isinstance causes things to break. But that's fine - in the case of our collection objects, the members don't need unicodifying anyway. This will, however, miss actual OrderedDicts, so you may wish to change the code accordingly. 


When I first tried the code, I kept getting the following error message

Exception TypeError: "'NoneType' object is not callable" in <function _remove at 0x7f4a9aae36e0> ignored

It wasn't fatal - the code worked in spite of it - but it was a little off-putting. With a little Googling, I found someone with the same problem, and a corresponding solution.

Basically, the issue was with the trace function not being disabled properly on exit. That's what the atexit line is for.


If the code you're testing is fairly simple, then you can use the code above as is. If it's a bit more complex, you'll probably find that converting strings to unicode in ALL functions causes problems in code outside of your control - in builtin modules, 3rd party modules.

This was the case for our API - the conversion seemed to upset regex (or something).

In this case, we need to make a slight tweak to make it only affect calls to functions from a particular module (the one we're trying to test).

In this case, we use the inspect.getmodule method with the frame's f_code member. That lets us identify what module the function being called came from, and apply unicoding conversion (or not) accordingly.

module = inspect.getmodule(frame.f_code)
if module and module.__name__.startswith('mymodule'):
    # etc.

Now, I'm not a fan of hard-coding a module name into the script. Ideally, I'd probably add a command line flag to specify the module of interest. But parsing that would make the code more complicated, and would break the 'drop-in replacement' nature of the script. So this is left as an exercise for the reader.


Conclusion

So, rather than writing a whole bunch of new tests, I got away with writing only ~50 lines of code. And I didn't even have to change any of the existing code.

Hurray for laziness.


Oatzy.


[inb4 use Python 3]

Sunday, July 31, 2016

Barber Queue

Time was when I needed a hair cut, I'd go at noon on a weekday, when the barber's is typically empty. One of the perks of being unemployed.

But these days I have to get my haircuts on Saturdays, before noon. Which is bad enough in itself - ideally I'd never see Saturday mornings at all. And to make matters worse, Saturday morning is also when the barbers is at its busiest.

Hence, I found myself sat in the waiter area of a barbershop for the best part of an hour. But this got me thinking about how queuing works at a barbers.



First In Who's Next?

At its core, the barber's queue is just a first-in first-out (FIFO) queue. But it has two interesting features:


1) The queue 'structure' is unordered

In general the queue 'structure' will be a waiting area with a bunch of seats. When someone new joins the queue, they're free to sit wherever. In fact, odds are, they'll pick a seat in a similar way to how men choose urinals - attempting to maximise personal space.

But the key point is, if someone were to just look at the queue, they wouldn't be able to tell who was next.


2) Each member of the queue knows whether or not they're next

Each member of the queue probably doesn't know who exactly is next, but they do know (with reasonable certainty) whether or not it's them.

The way this works is relatively simple - when you join the queue, you're aware of who was there when you arrived (and of anyone who arrives after you). So when all the people who were there ahead of you have gone, you know that you're next.



O(M G)

Okay, lets break out some Python (2.7)

Just for fun, let's say that the capacity of the queue is fixed - i.e. the waiting area has a fixed number of seats (though in practice, people are free to stand, as I was forced to).

class BarberQueue(object):

    def __init__(self, capacity):
        self._capacity = capacity
        self._queue = [None]*capacity
        self._length = 0

    def __len__(self):
        return self._length

    def __str__(self):
        return ', '.join(str(i) if i is not None else '_'
                         for i in self._queue)

    def push(self, obj):
        pass

    def pop(self):
        pass

So each member of the queue has some awareness of who is ahead of them. But they don't need to know specifically who's who, they just need to keep track of how many are remaining. And in fact, that remaining count is exactly equivalent to the member's position in the queue.

class Member(object):

    def __init__(self, obj, position=0):
        self.value = obj
        self._position = position

    def __str__(self):
        return "%s (%s)" % (self.value, self._position)

    def is_next(self):
        return self._position == 0

    def move_up(self):
        self._position -= 1

So, going back to the push and pop methods

def push(self, obj):

    if self._length == self._capacity:
        raise Exception("Queue is full! Please come back later.")

    for i, m in enumerate(self._queue):
        if m is None:
            self._queue[i] = Member(obj, self._length)
            self._length += 1
            return

Here, we're picking a 'seat' by iterate over the queue looking for the first empty slot (with a value of None). We could implement any seat picking strategy we fancy, this is just the easiest.

Once we find an empty seat, we create a new 'Member' object for the item, with position set to the current length of the queue, then increment the queue length. Also, if the queue has no empty slots, we raise an exception.

def pop(self):
    if self._length == 0:
        raise Exception("The queue is empty")
    for i, m in enumerate(self._queue):
        if m is not None:
            if m.is_next():
                value = m.value
                self._queue[i] = None
                self._length -= 1
            else:
                m.move_up()
    return value

Here we iterate over the queue looking for the member who is 'next' (has position 0). While we're looking for the next person, we also de-increment the positions of the other members.

Of course, this isn't how things work in practice. The barber doesn't go to each person and say "are you next?", "how about you?". They simply say "who's next?", and the person who believes they are next steps forward. Though arguably, that's just equivalent to asking every member concurrently. But let's not complicate things.

>>> b = BarberQueue(3)
>>> for i in xrange(3):
 b.push(i)
 
>>> print b
0 (0), 1 (1), 2 (2)
>>> b.push(4)

Traceback (most recent call last):
  File "<pyshell>", line 1, in <module>
    b.push(4)
  File "<pyshell>", line 36, in push
    raise Exception("Queue is full! Please come back later")
Exception: Queue is full! Please come back later
>>> b.pop()
0
>>> print b
_, 1 (0), 2 (1)
>>> b.push(4)
>>> print b
4 (2), 1 (0), 2 (1)
>>> for _ in xrange(4):
 b.pop()
 
1
2
4

Traceback (most recent call last):
  File "<pyshell>", line 2, in <module>
    b.pop()
  File "<pyshell>", line 46, in pop
    raise Exception("Empty queue!")
Exception: Empty queue!
>>> print b
_, _, _

So there we have it. Of course, this type of queue isn't really useful from a programming perspective.

All of insertion, deletion, and lookup are worst-case \(O(N)\), where N is the queue's capacity (not the number of people in the queue). Which is pretty much worse than all other types of queue.



Why do items keep disappearing from my queue?

Okay, let's move away from computer-sciencey queues. There are certain behaviours in real-world queues that don't apply or wouldn't make sense to programmatic queues.

For one, as I alluded to earlier, the capacity of the queue is not enforced - there's room for overflow, even if it means people have to stand.

But on the other hand, when a place does get that full, people are less likely to stick around.

In particular, we have two situations:


1) There's a non-zero probability that a person will not stick around if the queue is full or close to full. This probability will tend to be related to the length of the queue when that person arrives

\[p(not join) \sim f(capacity, length)\]

For example,

\[p(not join) = A\left(\frac{length}{capacity}\right) - C\]

where A and C are some constants relating to how likely the person is to stick around if the queue is 'full', and at what point they consider the place to be 'too full'.


2) There's a non-zero probability that a person already in the queue will leave before they're served.  This probability will typically depend on how long the person has been waiting, and how many people are still ahead of them

\[p(leave) \sim g(wait, position)\]

It's interesting because as time passes, wait increases, but position decreases. So how the probability evolves depends on how those factors balance against one another.

In particular, the probability evolution will likely depend on how each particular person responds to the sunk-cost fallacy - i.e. are they the sort to think "well, I've waited this long, I might as well see it through to the end", or do they think "this is taking too long, I've got better things to do with my time"?

For the sake of arguing, lets go with an exponential function for the general form.

For a sunk-cost person we might have

\[p(leave) = B \cdot \exp\left(-d\cdot \frac{wait}{position}\right)\]

This is a function where p goes to zero as position goes to zero or wait goes to infinity (B and d are some arbitrary constants).

Whereas for a non-sunk-cost person we might have

\[p(leave) = B \cdot \exp\left(-\frac{d}{wait \cdot position}\right)\]

This is a function where p goes to zero as position goes to zero, but goes to one as wait goes to infinity.

This gives us an interesting graph

Because 'position' is discrete you get this nice step function, with intervals of the probability steadily rising, then suddenly dropping. We can also see that there's a point at which probability of leaving is maximum, around the time you're in the middle of the queue. Which seems plausible.


We also have situations where a person will leave and come back later. But since, when they come back, they have to join the back of the queue, they're mathematically indistinguishable from a someone arriving for the first time.


One other complicating situation is people in groups. For example, if there's a parent and child ahead of you in the queue, the child is getting their hair cut by one member of staff, the parent is waiting; another member of staff asks 'who's next?' - is it you or the parent?

Situations like these add uncertainty into a member's queue position, and by extension their knowledge of whether they're next.

In that situation, we might wait to see if anyone else steps up, and if not we can assume it's our turn.

So we have \(p(next)\), which is a Bayesian probability that updates over time to reflect whether anyone else has stepped up yet. The longer we wait with no-one stepping up, the closer our probability gets to one.
Of course, if you wait too long, someone behind you might assume you're not in the queue after all and try to go ahead of you. But that's a topic for another blog.



Nothing's so simple that it can't be made complex

I've written about queuing before, in the context of a mathematical model of a cafe.

The long and short of it is this - people arrive at random (following a Poisson distribution) and join the queue with some probability (see above). Each iteration, some people are served, some join the queue, some get tired of waiting and leave, etc.

I'm going to iterate in 5 min time-step and say that a haircut takes 10-25 minutes (i.e. 2-5 steps). The exact duration is randomly generated for each customer.

I'm going to say that on average one person shows up every 10 minutes (0.5 per step). Here's an example of how that might look over 12 steps (one hour), using the Poisson distribution: [2, 0, 2, 0, 0, 0, 1, 0, 0, 1, 0, 0]

I set up the simulation so that you specify some number of steps for the barbers to be considered 'open'. After that, no more people are added to the queue, but the simulation keeps running until the queue is empty.

To begin with, I made it so that everyone who arrived stayed.

With 2 workers, capacity 10, and an arrival rate of 0.5, the queue length typically stayed below 3. The highest I saw it go was 7, which is still comfortably within capacity.
Above is an example of how the queue size varied over time in a particular simulation.

Increasing the arrival rate to 0.7, the average maximum queue length goes up into the low teens. And when the rate goes up to 1, the maximum queue length goes all the way up into the 30s.

Homework question - how does maximum queue length vary as a function of number of workers and arrival rate?

As I mentioned, those simulations assumed that everyone stuck around. Once you turn on probabilistic leaving, things get a bit more interesting.

I tried both versions of the leaving probability. The result was largely the same, except that sunk-cost people tend to leave sooner - average wait before leaving 2.6 for sunk vs 7.4 for non-sunk. This is what we'd expect - people adhering to sunk cost will tend to leave before they get too invested.

In the above example, orange is a time step when someone left the queue, and red is when a newcomer decided not to join the queue. I tuned the probability constants so that people don't start leaving until we're close to or at capacity, as we'd expect in real life.



You and I have different ideas of what constitutes 'interesting'

Going back to the original description of the barber's queue, here's an example of a full queue (bracketed numbers are queue positions)

542b5af6 (3), 33e323cf (5), 69e60241 (6), b3f12010 (0), fc89732e (7), 991f8709 (1), a0cb93cc (2), 17186c75 (4), 57269a3e (8), 8f1b61ca (9)

Notice how, even with the basic seat picking strategy (take the first available), the members aren't in a predictable ordered.


When we look at the waiting times, we can see some interesting things. For example
...
1c7081d1 waited 7.0
34 1
35 3
50be68cf waited 9.0
36 3
37 4
26ca11b7 waited 3.0
38 3
39 3
a571a7da waited 5.0

...

Here we have a person who had to wait 9 steps (45mins) to be served, followed by someone who only had to wait 3 steps (15mins). Which just goes to show, how long you wait in a queue is very much a matter of timing and luck.


It's also interesting that you can run the simulation multiple times with the exact same settings, and one time the queue will never go higher than 4, while in the the next it'll go as high as 13. This is complexity at work - various small random factors in the model interacting to produce wildly different outcomes.


So yeah. If you're interested, you can see the full code here. I may have gotten carried away with the object-orienting.


Oatzy.



[Post-Script

My boss recently pointed out that it'd been over a year since my last blog post. That was another perk of being unemployed - more time to come up with dumb blog posts. Anyway, here's a quick update on some stuff.


Pirate Game

The last blog post was about the making of an Android game - The Pirate Game.

The game is now finished-ish and has been released in 'beta' on the Play Store.

In the previous post, I mentioned that the game would eventually get a less utilitarian design. I ended up making that design myself (because I'm a control freak). I'm pretty pleased with how it turned out.
Also, following a... less than positive review, I added some new game play modes.

I never did figure out multiplayer, though. If I ever get the time or inclination to go back to the game, that'll be on the todo list. But for the time being, I don't anticipate any updates to the game. Certainly not any time soon.


New Job

So yeah. I finally got a job. I'm now a Software Developer at a company called Pixit Media.

The company sells large scale 'storage solutions' to companies primarily in the VFX industry, as well as universities and other such people that do high performance computing.

What I personally work on is primarily a Python API for the IBM Spectrum Scale (GPFS) filesystem. You can see the API docs online. I wrote a decent amount of the documentation (and the code that's being documented).

In particular, the 'Getting Started With List Processing' guide. Admittedly the topic is a bit niche - I doubt many readers of this blog even know of GPFS, let along have a cluster with it installed. But you might still find it interesting; you can learn some stuff about MapReduce - a technique for taking advantage of parallelism when processing large data-sets.

There's also the 'example scripts' repository - scripts written to use the API, some of while I wrote. But, again, they're a bit niche.

]

Monday, February 23, 2015

Designing a General Relativity Based Casual Game

Let's suppose, for the sake of arguing, that I'm unemployed. It'd pretty great if there were a way I could make a quick buck (or a quick quid, in my case). Obviously, the solution is to make a super addictive mobile game. Then I can just sit back and watch the cash come pouring in.

Okay, so that's a pretty terrible plan - making a mobile app of any kind is no trivial thing, and there's no guarantee of making any significant amount of money (or any money at all). Still, I thought about it, and wondered - if I were to make a game, what would I do. This is one of the ideas I came up with.


Game Concept

A common type of casual game is the physics-based game - that is, a game where objects obey a form of classical Newtonian mechanics: your Angry Birds, your Flappy Birds, your games that don't involve birds. So, for example, when you catapult an angry bird, it follows a parabolic path like a real bird would (drag not-withstanding).

And that's all fine, but I wondered if you could make a game based on non-classical physics - quantum mechanics, relativity, that sort of thing.

Now, I don't want to call this 'Interstellar: The Game' (not least because of copyright). That is, however, a good short hand for the idea behind the game.

The idea is this: you've been off exploring space, and now it's time to go home. So the main objective is to get back to Earth as quickly as possible, whilst dodging obstacles like planets, stars, black holes. But these massive objects have a gravitational pull, which can make maneuvering a little trickier, especially if you get too close to a black hole.

On the other hand, you can use massive objects to get a speed boost, by doing for example a powered flyby. But being close to massive objects for too long (black holes in particular) will cause time dilation, which might make you late home.


Game Construction

The way I imagine the game is as a 'side-scroller', where you move right to left (assuming the phone/tablet is held in landscape). I figure it would be a top-down view, moving in the x-direction with the spaceship - i.e. the sprite would be fixed in the x-direction, with freedom of movement up and down. All the obstacles would then move towards the sprite in the negative x-direction.

Here's what a (non-relativistic) binary orbit looks like from a frame of reference moving in the x-direction with the particle (blue).


In the code, this is done by calculating the particle velocity as normal, but instead of adding the velocity in the x-direction to the particle, you subtract it from all the obstacles.

For the sprite's up-down movement, you could either constrain the player to within the screen, or else if they drift off the edge of the screen, it would be a 'lost in space' game over.

Level construction can either be done by hand, or levels can be generated randomly / procedurally by adding massive objects of varying type/size/mass etc. at various points along the game map. I'd say go for random/procedural, since that makes creating levels easier, and would mean there are effectively unlimited levels. Though it might be worth storing generated levels for replayability.

For the representation of black holes, you could show an accretion disc (like they did in Interstellar), have a background star field that's gravitationally lensed, or show a dotted line for where the event horizon is. Alternatively, you could do nightmare mode - give no visual indication of a black hole, and let the player infer its presence from its gravitational pull.

I tend to think that controls should be kept simple, and should be appropriate to the medium - in other words, no on-screen d-pads on touchscreens (if you can help it). Instead, I'd say use four directional swipes for a speed boost in whichever direction. In terms of code, this would be done by adding some constant to the velocity in the direction of the swipe.

One thing to remember is that in space there is no drag - nothing to slow you down (gravity notwithstanding). Steering like this can be surprisingly tricky. If you accelerate left, you will keep moving leftwards until you accelerate right to counter balance.

On the other hand, no drag means the player could keep accelerating forward until they're going arbitrarily fast. To stop this, you could put a limit on how many times the player can accelerate - for example, say that there's a limited amount of fuel. This also means the player would have to use their fuel wisely - baring in mind that breaking counts as accelerating. This makes tricks like gravitational assists even more important.

Anyway, having an idea is all well and good, but really an idea isn't worth much if you can't prove that it's viable. So...


Game Mechanics

Preamble

I studied theoretical physics at university, and that included a module on General Relativity. But that module only covered the basics. As the lecturer pointed out, General relativity is a graduate level topic.

I read a bunch of articles for this blog, and have tried to get as scientifically accurate as possible - or at least I've tried to avoid making massive errors. But even so, there are numerous assumptions, approximations, and inaccuracies in this formulation. So just bear that in mind - I don't necessarily know what I'm talking about. Do feel free to offer corrections.

Prototyping - For coding the prototype/simulations I used Python with PyGame.

In PyGame, coordinates are defined with the origin (0,0) in the top right corner. This isn't a problem mathematically, it just means the graphical layout is 'upside-down'.

To make the videos in this post, I added the line
pygame.image.save(screen, 'output/frame%05d.png' % framenum)
to the code's main loop to save individual frames as a series of images (you'll need to iterate 'framenum'). I then used 'ImageJ' to convert those images to (avi) video.

Terminology - I'm going to refer to moving objects as 'particles', 'planets', or as a 'spaceship' in the context of the game. I'll refer to static, gravitating objects as 'obstacles', '(massive) bodies', or 'stars'. In images and videos, planets are blue, and stars are yellow.

Units - I chose to use 'pixels' as the unit of length, and 'frames' as the unit of time. In this case, velocity is measured in 'pixels per frame' (px/fr).

This is useful, because it means the velocity of our particle is updated as \(v(t) = v(t-1)+a(t)\) and the position is updated as \(x(t) = x(t-1)+v(t)\). For my simulation, I used a frame rate of 50fps.

Physical Constants - We could use geometrised units, where the Gravitational Constant (G) and the Speed of Light (c) are both set equal to one. In that case we wouldn't need to include them in the maths/code. However, my inner-mathematician likes generality, so I'm going to keep them around.

The constant 'G' only ever appears alongside obstacle masses (as GM), so we can set it to one and say that it's value is absorbed into the (otherwise arbitrary) mass values. It's worth keeping G around though, so we can easily tweak the strength of gravity, if need be.

The speed of light controls the strength of relativistic effects - larger c, weaker relativity. In my code I set the speed of light to c = 12 px/fr. This seems to make things behave in a way that looks 'right'.

Object Properties - We could try to go for a certain level of realism - try, for example, to work out a conversion between real world distances and pixels, try to get everything to scale. But that's too much faffing. Besides, if everything were to scale, a 1px Earth would orbit at a distance of 23000px around a 100px sun.

For obstacle masses, I set M = 500 (arbitrary units). The particle mass isn't really important since it can be cancelled out in all the equations. From a theoretical perspective, all that matters is that it's much smaller than the obstacle masses. I set both the obstacles and particle radii to 12.5 pixels (25x25px sprites). Object radii aren't important outside of collision handling.


Newtonian Gravity

For particle acceleration, we start with Classical Newtonian gravity.

\[ a = \frac{G M}{r^{2}} \]

where G is the gravitational constant, M is the mass of the gravitating body (the obstacle) and r is the distance between the centres of our particle and the obstacle. We're assuming that the particle is moving in a 2D plane, so we have \(r = \sqrt{dx^{2} + dy^{2}}\), where \(dx = x_p - x_{ob}\) and \(dy = y_p - y_{ob}\) are the x and y distances between the particle and obstacle.
This gives us the magnitude of the acceleration, pointing from the centre of the particle to the centre of the obstacle. For our purposes, we need to resolve this acceleration into x and y components.

The components are given by

\[ \ddot{x} = -\frac{G M}{r^{2}} \cos(\theta) \equiv -\frac{G M}{r^{3}} dx\\ \\
\ddot{y} = -\frac{G M}{r^2} \sin(\theta) \equiv -\frac{G M}{r^{3}} dy\]

Where \(\theta = \arctan\left(\frac{dy}{dx}\right) \) is the angle between the x-axis and the position/acceleration vector. The double dots mean 'second derivative with respect to time', i.e. acceleration.

If there are multiple gravitating objects, we have to calculate the contributions from each, and add them all together to get the total acceleration.

\[ \ddot{x} = -\sum_i \frac{G M_{i}}{r_{i}^{3}} dx_{i} \\ \\
\ddot{y} = -\sum_i \frac{G M_{i}}{r_{i}^{3}} dy_{i} \]

As mentioned above, once we have the acceleration we add that to the velocity, then add the velocity to the position to work out where the particle moves to. Below is an example of a Newtonian orbit using the above.



General Relativity

Since we want the game to include General Relativistic effects like time dilation, we also have to take into account the other effects of relativity - in particular, how it modifies particle motion/acceleration.

For simplicity, we're going to assume our obstacles (planets, stars, black holes) are uncharged, and non-rotating. In that case, we use the Schwarzchild metric, which describes the gravitational field in the vicinity of an (uncharged, non-rotating) massive object.

\[ c^2 d\tau^2 = \left(1 - \frac{r_s}{r}\right) c^2 dt^2 - \frac{dr^2}{\left(1- \frac{r_s}{r}\right)} - r^2 d\theta^2 \]

We're assuming the particle is moving in the 2D plane around the equator of the massive object. \(r_s\) is the object's Schwarzchild radius (also known as the 'event horizon' in the context of black holes), defined as

\[ r_s = \frac{2GM}{c^2} \]

Without going into too much detail, we can derive from the Schwarzchild metric the particle's acceleration in Cartesian-like coordinates as

\[ \ddot{x} = - \frac{G M}{r^3} dx - \frac{3 G M L^2}{c^2 r^5} dx \\ \\
\ddot{y} = - \frac{G M}{r^3} dy - \frac{3 G M L^2}{c^2 r^5} dy \]

Where the first term is the classical Newtonian gravitation, as seen above.

The second term is purely relativistic, and will usually only have a significant effect when a particle is sufficiently close to a massive body's Schwarzchild radius.

The main effect of this term is to increase gravitational acceleration in the vicinity of the massive object, and to cause close orbits to precess, as can be seen below


Notice that the point at which the particle is closest to the 'star' (the perihelion) moves over time. This doesn't happen in the classical limit of Newtonian gravity. In fact, explaining the anomalous precession of Mercury was one of the first pieces of evidence supporting the theory of General relativity.

In the acceleration equation, L is the angular momentum (per unit mass) of our particle, and c is the speed of light. Angular momentum (per unit mass) is calculated as

\[ L = \left|r\right| \left|v\right| \sin(\phi) \equiv (dx\ v_y - dy\ v_x) \]

where \(\phi\) is the angle between the radial vector (r) and the velocity vector (v)

The angle can be calculated as

\[ \phi = \alpha - \theta = \arctan\left(\frac{v_y}{v_x}\right) - \arctan\left(\frac{dy}{dx}\right) \]

where \(\alpha\) is the angle between the velocity vector and the x-axis, and \(\theta\) is the angle between the position vector and the x-axis.

In the two body case, angular momentum is a constant of motion, so only needs to be calculated once - say, once the particle's initial position and velocity has been set.

Caveat - Strictly speaking, the radial length 'r' in the Schwarzchild metric is not the same as the Euclidian distance \(\sqrt{dx^2 + dy^2}\). At least, not when relativistic effects are significant. Rather, r is defined as the circumference of a sphere surrounding the massive body, divided by \(2\pi\). Defining r this way means the length of an interval dr isn't affected by the curving of spacetime near the massive body.

By comparison, an observer falling towards a black hole would see lengths stretching longer and longer the closer they got to the black hole event horizon - an observer at a distance r would measure the interval dr to have a length of \( \left( 1 - \frac{r_s}{r} \right)^{-\frac{1}{2}} dr \).

All this to say, the particle distances displayed in the game (and the simulation videos) are distances in Schwarzchild coordinates, projected onto a Euclidean plane, not distances as viewed by an observer such as our particle/spaceship.


Relativity for Multiple Bodies

This is where things get dicey. In the General relativistic limit, combining the fields of multiple massive objects is not so straightforward.

The lazy way of doing this is to just vector sum the relativistic accelerations, like we did for the Newtonian case. This is problematic, though, because it assumes the particle's angular momentum around each massive body is constant. This is not true.

Instead, we can make a simplification - we can sum the Newtonian terms as before, but we'll only consider the relativistic term when we're sufficiently close to any given body. This cut-off is going to be fairly arbitrary. Ideally, we want the 'radius of influence' to be as big as possible, but small enough that the particle can never be within the relativistic limits of more than one body at a time. The main problem here is that the way the game is set up means all the massive bodies are unrealistically close together, so it's hard to make the radius of influence big enough to be ideal. In my code, I chose the cut-off to be \(10 r_s \simeq 70px\).

Now, the angular momentum around this close body is still not constant (although, arguably, it may be sufficient to assume it is). The forces from other local bodies can cause torque, which will change the particle's angular momentum.
Torque is calculated as

\[ \Gamma = \left|r\right| \left|F\right| \sin(\gamma)\ \equiv r_x\ F_y - r_y\ F_x\]

where \(\gamma = \theta_1 - \theta_0 = \arctan(\frac{dy_1}{dx_1}) - \arctan(\frac{dy_0}{dx_0}) \) is the angle between the position vector (r) and the force vector (F), and \(\theta_0\) and \(\theta_1\) are the angles between r and the x-axis, and F and the x-axis, respectively.

For the force from a gravitating body, the torque can be written out as

\[\Gamma = \frac{G M_1 }{r_1^3} \left(dy_0\ dx_1 - dx_0\ dy_1\right) \]

We only need to take into account Newtonain gravitation in calculating torque, since we're already assuming the particle is outside the relativistic limit of these other bodies.

So once the particle enters a massive body's radius of influence, we calculate it's initial angular momentum, as above. Then, in each time step (for as long as the particle is in the radius of influence), we calculate the total torque from all local bodies, and add that torque to the particle's angular momentum (before calculating the new acceleration).

This is a big simplification. But it's good enough, and correct in the limiting case of a single/well isolated body. Especially if the radius of influence can be made sufficiently big.

While we're on torque - if our spaceship accelerates (fires its thrusters) while it is close to a massive body, we will need to take into account any torque from that as well

\[ \Gamma = dx_0\ a_y - dy_0\ a_x \]

Where \(a_x\) and \(a_y\) are the x and y accelerations caused by the thrusters. This is ignoring the details of how a spaceship actually maneuvers. If you're interested, you'll have to look into that for yourself.

So putting it all together we have

\[ \ddot{x} = -\sum_i \frac{G M_i}{r_i^3} dx_i \ - \frac{3 G M_0 L_0^{2}}{c^2 r_0^{5}} dx_0\\ \\
\ddot{y} = -\sum_i \frac{G M_i}{r_i^3} dy_i \ - \frac{3 G M_0 L_0^{2}}{c^2 r_0^{5}} dy_0 \]

where \(M_0\) is the body whose radius of influence the particle is within (if any), and

\[L_0 = L_{0}(t-1) +\sum_{i\ne 0} \frac{G M_i}{r_i^3}(dx_i\ dy_0 - dx_0\ dy_i) + (dx_0\ a_y - dy_0\ a_x)\]


Time Dilation

For a single massive object, the time dilation can be derived from the Schwarzchild metric (see above). Dividing through by \(c^2 d\tau \) and rearranging we get

\[  \left(1- \frac{r_s}{r}\right)^2 \left(\frac{dt}{d\tau}\right)^2 = \left(1- \frac{r_s}{r}\right)\left(1 + \frac{r^2 \dot{\theta}^2}{c^2}\right) + \frac{\dot{r}^2}{c^2} \]

where the dots mean 'derivative with respect to \(\tau\)'. And since \( \dot{r}^2 + r^2 \dot{\theta}^2 \equiv v_x^2 + v_y^2 = v^2 \), we can re-write and rearrange further to get

\[ dt = \frac{d\tau}{\left(1- \frac{r_s}{r}\right)}\sqrt{1- \frac{r_s}{r}\left(1 + \frac{r^2 \dot{\theta}^2}{c^2}\right) + \frac{v^2}{c^2}} \]

where dt is the 'coordinate time' - time as measured by an observer at rest, far away from any gravitational fields. For the sake of the game we can say that this represents time as measured on Earth. In reality, the Earths gravitational field does cause it's own time dilation effect. \(d\tau\) is the 'proper time' - time as measured by a clock on our spaceship.

Notice, in the limit of a stationary particle, that is \(\dot{r}=r\dot{\theta}=0\), we get the purely gravitational time dilation

\[ dt = \frac{d\tau}{\sqrt{1 - \frac{r_s}{r}}} \]

Similarly, in the limit of \( r \gt\gt r_s\) (i.e. when our particle is far from any gravitating bodies), the time dilation equation reduces to the Special Relativistic case

\[ dt = d\tau \sqrt{1 + \frac{v^2}{c^2}} \]

Note - in this, the velocity v is the 'proper velocity' - velocity with respect to proper time. This is different from coordinate velocity - velocity with respect to coordinate time.

While a particle can't have a coordinate velocity greater than the speed of light 'c', because of time dilation it can have a proper velocity greater than 'c'. This doesn't, however, mean that a particle can travel faster than light, since light has a proper velocity of infinity. Coordinate velocity and proper velocity are related by

\[ v_c = v \left(\frac{d\tau}{dt}\right) \equiv \frac{v}{ \sqrt{1+\frac{v^2}{c^2}}}\]

where the equivalence is true in the Special relativistic limit. Notice that when proper velocity equals the speed of light (c), the coordinate velocity equals \(\frac{c}{\sqrt{2}}\) - less than the speed of light.


Time Dilation for Multiple Bodies

Here, we once again have the problem of combining the effects of multiple gravitating masses. Some would argue, at this point, that trying to be scientifically accurate is more trouble than it's worth. Still, I'm going to at least try.

To start, we can replace the Schwarzchild radius terms with the local (Newtonian) gravitational potentials as

\[ \frac{r_s}{r} = \frac{2GM}{c^2 r} \rightarrow  \sum_i \frac{2 G M_i}{c^2 r_i} =  \sum_i \frac{2U_i}{c^2} = \frac{2U}{c^2}\]

Where \(U_i\) are the individual gravitational potentials of the various local gravitating bodies.

For the term in \(r^2 \dot{\theta}^2\), we note that \(r^2 \dot{\theta}^2\ = \frac{L^2}{r^2} \), where L is angular momentum (per unit mass). So we can use the same trick we did for combining accelerations - that is, only take this term into account when the particle is within a body's radius of influence.

So putting it all together, we have

\[ dt = \frac{d\tau}{\left(1 - \frac{2U}{c^2}\right)}\sqrt{1 - \frac{2U}{c^2} - \frac{2U_0}{c^2}\frac{L_0^2}{r_0^2 c^2} + \frac{v^2}{c^2}} \]

With \(L_0\) defined as above. Alternatively, we could just omit the angular momentum term since it's of order \(c^{-4}\). In my simulation, it amounted to a ~0.17% decrease in the dilation ratio.

The typical dilation ratio in my simulation was \(\frac{dt}{d\tau} \approx 1.07\), which admittedly isn't significant.

In the movie Interstellar, a lot of the time dilation came from the fact that their black hole (Gargantua) was spinning very fast. If you want to do that, you'll have to work it out for yourself - for that you would use the Kerr metric. With rotating black holes, you'd also need to take into account other effects like frame dragging.

Of course, we've chosen arbitrary masses, an arbitrary speed of light, distances not to any sort of realistic scale. We can always tweak the specific time dilation to our liking. In particular, if the realistic dilation isn't dramatic enough, we could artificially inflate it. So long as we keep key behaviours, like dilation going to infinity when the particle reaches a black hole event horizon, etc.


Collision Handling

Collision handling is important for figuring out when the game is over - for figuring out if the player crashed into a planet, or fell into a black hole. Now, PyGame has in-build collision handling. But where's the fun in that.

The easiest way to check if two circular objects have collided is to check if the distance between their centres is less than the sum of their radii

In other words, we have the condition - a collision occurred if \( dx^2 + dy^2 \le (r_p + r_{ob})^2 \).

Easy. There is a special case though - because the particle moves in discrete steps from frame to frame, if the particle is moving fast enough, it could move from one side of an obstacle to the other without ever overlapping.

The basic collision handling would miss this. So this is a little trickier to catch. In the above diagram, I've traced the path going from the particle's position in the previous step to its current position.

We can say that the particle collided with the obstacle if there's a point on that path that's closer to the centre of the obstacle than \( R = r_p + r_{ob} \). In other words, if the particle had moved continuously along the path from \(x_0\) to \(x_1\), would there have been any point(s) where the particle and the obstacle overlapped.

To figure this out, first we find the equation for the line passing through the particles' two positions as

\[ y = m x + c = \left(\frac{v_y}{v_x}\right) x + \left(y_0 - \frac{v_y}{v_x} x_0 \right) \]

(in this case, m and c are gradient and constant, respectively, not mass and speed of light).

Second, we're going to define a function for the distance between a point on the particle's path and the obstacle's centre

\[ f(x) = (x_{ob} - x)^2 + (y_{ob} - y)^2 = (x_{ob} - x)^2 + (y_{ob} - m x - c)^2 \]

Now we're going to look for the point along the path of closest approach to the obstacle. To do that we look for x such that the separation f(x) is minimised. In other words \( \frac{d}{dx} f(x) = 0\)

If you work through the maths, you get

\[ x_{min} = \frac{x_{ob} + m y_{ob} - m c}{1 + m^2} \]

Which gives us the condition - a collision occurred if \(f(x_{min}) \le R^2\) and \(x_0 \le x_{min} \le x_1\) (assuming \(x_0 \le x_1\)).
Both these collision handlers can also be used to check if the particle crosses a black hole's event horizon - just replace \(r_{ob}\) with \(r_s\), the Schwarzchild radius.

Of course, in the game, the spaceship sprite wouldn't be circular. But, you know. You could create an 'imaginary' circle around the sprite, or give the sprite a radius of zero (so it collides when its centre overlaps the obstacle). Or you could just use a third party library. Either way.


Bonus: Elastic Recoil

For my own amusement, I had the particle elastically recoil off of the obstacles, as well as the sides of the screen, so that I could just watch it drifting and bouncing about endlessly. It's weirdly hypnotic.

For screen boundaries, you have, for example
In this case, we flip the particle velocity in the x-direction as \( v_{x}' = -v_{x} \), and move the particle's x position to \( x_{p}' = W - (x_{p} - W) = 2 W - x_p \), where W is the pixel width of the game window. You do the same sort of thing for the other 3 boundaries.

For recoil off obstacles, things are not so straight forward. The procedure works like this: First, check for a collision. Then, find the point along the particle's path where it first makes contact with the obstacle - that is where \(f(x) = R^2\).
Finding this point is just a matter of solving a quadratic equation to get

\[ \begin{align*}
x' &= \frac{x_{ob} + m y_{ob}-mc}{1 + m^2} \pm \frac{\sqrt{(x_{ob} + m y_{ob}-mc)^2  - (1+m^2)(x^2_{ob} + y^2_{ob}+c^2 - 2 c y_{ob} - R^2)}}{1 + m^2} \\ \\ &\equiv a \pm b
\end{align*} \]

If \(v_x \gt 0\) then x' = a - b, if \(v_x \lt 0\) then x' = a + b

So move the particle back to the point of first contact. Then we want to flip the particle velocity so that the angle of incidence \(\phi\) between the velocity vector (v) and the position vector (r) equals the angle of reflection
The new velocity is given by

\[v'_x = -\left|v\right| \cos(\beta)  \\ \\ v'_y = -\left|v\right| \sin(\beta)\]

with \(|v| = \sqrt{v_{x}^{2} + v_{y}^{2}} \), and

\[ \beta = \alpha - 2 \phi = 2\arctan\left(\frac{-dy}{-dx}\right) - \arctan\left(\frac{v_y}{v_x}\right) \]

where \( \alpha \) is the angle between the x-axis and the old velocity vector (v), and \(\beta\) if the angle between the x-axis and the new velocity vector (v').

Note - in the code I used 'arctan2(y,x)' from the Numpy library, for which the signs of 'x' and 'y' are important. dy and dx are negative in the above because I have the position vector pointing from the particle to the obstacle (opposite to how it's defined).

Finally, we want to move the particle to where it should be from recoiling.

In other words, we want to move the particle to the point along it's new velocity vector, such that it's the same distance from the contact point as it would be if the particle hadn't collided. To do this, we can just resolve the distance \( r' = \sqrt{(x_1 - x')^2 + (y_1 - y')^2} \) in the direction of the new velocity vector, as

\[ x_1 ' = x' - r' \cos(\beta) \\ \\ y_1 ' = y' - r' \sin(\beta) \]

Admittedly, this doesn't work perfectly - sometimes it behaves a little screwy. Especially when the particle spirals in on an obstacle. Of course, none of this recoil stuff is necessary for the game, since you want collisions to end the game. Like I said, this was more for my own amusement.


Conclusion

If you were feeling really ambitious, you could do the game from the first-person perspective of someone on the spaceship - see for example 'Falling into a Black Hole', or the Interstellar papers.

There are almost certainly things wrong with this 'formulation' of General Relativity. But then, this is meant to be for a game, so it doesn't need to be 100% scientifically accurate. I did feel like I should make an effort, though. And I'm willing to call this good enough. If you're someone who knows what they're talking about, feel free to offer your thoughts/corrections.

Anyway, I've done my job as a theoretician - now it's up to someone else to actually make the game. If you do, please send me link to the completed game. And maybe give me a name-check? A little kickback would be nice too... ;)

You can have a look at my prototype/simulation code here.


Oatzy.


[Gravity, don't mean too much to me.]

Thursday, January 01, 2015

Lights Out (Game)

I was watching this year's RI Christmas Lectures, and it reminded me of this game I had when I was much younger. It was, I think, a Christmas present from my grandparents, back in '96. It was called Lights Out, and it looked something like this.

Basically, the console had a 5x5 grid of these rubbery, transluscent buttons with little red lights underneath them. When you press a button, that button and the four adjacent buttons switch state - light on to light off, or off to on. So, in the example below (where blue square are lights), if you press the centre square, the centre square turns off and the four adjacent squares turn on.


In the game, you're given a pattern of lights, such as the one above, and the task is to press buttons until all the lights are out.

Interestingly, in the pattern above (on the left), you can clear the board by just pressing the squares that are on at the start (in any order).

Anyway, I was reminded of that. Meanwhile, I'd been learning some Pygame for this other thing (that I'll talk about in a later blog). Now, Pygame is a Python library for making games (amongst other things). So you can probably see where this is going...

Coding a 'Lights Out' clone was straightforward enough - the code is short and not too complicated. The above images are actually screenshots. You can look at my code/download it here. You'll need Python and Pygame to run it.

This is a very threadbare version of the game. You can load a level as a text file of five lines of ones and zeros, like here. Or else you can just play from a blank board - make patterns then try to get rid of them again (not necessarily as easy as it sounds).

If you want something more fancy, or if you don't want to install Python and all that, there are probably loads of playable clones online.

Now, at this point you're probably thinking "Yeah, that's great. But what about the maths of Lights Out?". And you'd be damn right to ask. If you're interested, there's a good summary, as well as links to articles/papers, on the wikipedia page. I don't want to just repeat what's written there. Basically, it all comes down to linear algebra.


Anyway, have fun with that.


Oatzy.


[Call this a late Christmas present.]

Tuesday, October 14, 2014

Dark Matter and Machine Learning

If you're a long time follower of this blog you'll probably have noticed I haven't posted much in the last two/three years. This is because I've been busy with university stuff. This blog is sort of related to that. The following section is an edited version of the background I wrote for my final year project report (hence why it's so formal). Enjoy.


Background - Dark Matter and DM-ICE

According to current theories, dark matter makes up 23% of the mass-energy density of the universe. However, to date, it hasn't been conclusively observed directly. Instead, it's existence is inferred from large-scale gravitational effects, like the rotation curves of galaxies, and gravitational lensing. Based on rotation curves, we find that the outer edges of galaxies move as if there is more matter present than what we can directly observe. Consequently, it has been theorised that there is a halos of dark matter around the outer edges of galaxies, including our Milky Way.


As our Solar System orbits around the galactic centre, it moves though this halo, experiencing an effective dark matter wind. This, combined with the tilt of the Earth's orbit around the sun, means we expect to see an annual modulation in the dark matter flux hitting the Earth. This modulation should have it's maximum around June when the Earth is moving into the wind, and it's minimum around December when it's moving away from the wind, regardless of location on Earth.

The current best candidate for dark matter are so-called Weakly Interacting Massive Particles (WIMPs). In direct detection experiments, we look at interactions between WIMPs and the nuclei in some target, such as Sodium Iodide (NaI) crystal. The WIMPs scatter elastically off nuclei in the target, and the recoil of the nuclei cause scintillation photons to be emitted, with summed energies of ~10keV range. Direct detection is hard, however, since WIMP interactions are rare events, where we expect << 1 count/day/kg. So identifying WIMP events, especially at low energies, requires significant background and noise suppression.

Various direct detection experiments (DAMA, CoGeNT, CRESST, CDMS, etc.) have already been run, but none yet have conclusively detected dark matter. Of particular interest from these experiments are the results from the DAMA/NaI and DAMA/LIBRA experiments. The DAMA experiments (Gran Sasso, Italy) ran direct detection using NaI(Tl) (Thallium doped Sodium Iodide) crystals, and observed an annual modulation in the 2-6keV energy region, which they claim is the result of the dark matter wind. However this interpretation is disputed. It is argued that the modulation could have come from some other unaccounted for signal – for example muon flux also has an annual modulation with peak around June (in the Northern Hemisphere).
Source - ArXiv:0804.2741
DM-ICE is a joint venture between University of Sheffield (UK), and University of Wisconsin (Madison, USA). The aim, along with other independent direct detection experiments being run around the globe, is to test the DAMA result by looking for this same modulation, and by ruling out other possible sources of modulation. By running at the South Pole, DM-ICE is in the opposite hemisphere to DAMA, so seasonal effects – such as temperature and muon flux – are reversed, while the dark matter modulation should stays the same.



Context

For my final year project I was looking in particular at the performance stability of the detectors (loss of light yield, etc.) and how to correct the data for those effects. I also did some preliminary modulation analysis - looking for evidence of dark matter in the data.

This blog post isn't really to do with any of the work I did on my project, though. Rather, this is looking at how to remove noise from the data. For my project, this work had already been done, so I didn't have to worry about it. So why am I looking at noise removal now?

When the academic year was over, my project supervisor asked if I wanted to carry on working on the project (funded) for the summer. It made sense, since there was still work to do, and since I already knew the project well. So mostly I was doing more of what I did for the project. But the supervisor also wanted to adapt the DM-ICE project into a simplified version for groups of 3rd year students to do.

Specifically, the students would have to figure out the best cuts to remove noise, then do a modulation analysis to decide if there was evidence of dark matter.

So I set up a simplified version of the data and wrote up a description of the project for the students. The supervisor then asked if I'd work through the project and write up a model answer so that he'd have something to mark against. Hence, it was my turn to try the noise removal stuff.



Pulse Classification

In a DM-ICE detector, we have an 8.5kg NaI(Tl) crystal between paired photo-multipliers (PMTs). These two PMTs - DM0 and DM1 - record scintillation events as voltage pulses (in ADC units). When the raw data is processed, we look for pulses from both of the PMTs which occur within 100 ns of each other. These pairs of pulses are assumed to have come from the same scintillation event.

At low energies, pulses look something like this


where the various spikes represent single photo-electron events (SPE). Unfortunately, at low energies there is also a significant amount of noise events - specifically, electromagnetic interference (EMI) that look something like this


and 'thin peaks', that look something like this


So the task is to identify and remove these EMI and thin peak events.

From looking at plots of the pulses, it's fairly straightforward to tell the different event types apart. But computers can 'look' at the pulses. Instead, we have to define some parameter - values that we can calculate that will tell the computer something about the shapes of the pulses.

Take, for example, the EMI pulses. From the plot we can see that EMI pulses oscillate rapidly between positive and negative ADC values, where the SPE and thin peak events don't. So we can define some parameter, which we'll call 'emi', which basically counts how many times the pulse goes from positive to negative (or negative to positive) in the first 40 bins of the pulse.

So since the emi value is going to be much higher for EMI-type pulses than for SPE or thin peaks, identifying and removing EMI events is relatively easy.

Separating the SPEs from the thin peaks is less straightforward. For this, we calculate various other parameters - for example, since thin peaks will typically only have one peak, while an SPE event will have several, we can look at the peak number. Other things we can look at are how quickly the pulse decays to zero (log-meantime), or how the pulse energy is distributed. The details of how these parameters are calculated aren't important to this blog.



Human Learning

What you can do, then, is look at a bunch of pulses - plot it, what type of event is it? What are it's parameter values? You collect together the results and you ask, what sort of emi values do EMI-type events have? What sort of peak numbers do SPEs and thin peaks have. How can you use that information to tell the peaks apart?

And that's fine. But it's also boring. You have to plot, and calculate, and record. Repeat. Analyse. Boring.

I subscribe to the philosophy of "why do yourself what you can make a computer do?". In practice this usually means making the work more 'complicated' in the short term, but once it's done it makes life a whole lot easier. And while the work is more complicated, it's at least more interesting as well.

So what I did was I wrote a program. It looks something like this
(I used Python with Numpy, Matplotlib for plotting, and PyQt4 for the GUI.)

Basically it plots the pulses and presents you with three buttons - signal (SPE), noise (thin peak), and EMI. On the front-end, all the user has to do is look at the plot and decide what type of event the pulse is. Meanwhile, on the back-end, the program calculates the parameters, and sorts them according to the user's classifications.

Once you've classified a set of pulses, you're presented with histograms of the different parameters, like this (for the emi value parameter)


where the different event types are plotted in different colours. This makes life a lot easier. For example, in the above, it's immediately apparent that EMI-type events can be removed from the data by cutting any event with an emi value above 25.



Machine Learning - Naive Bayes Classifiers

Lets make this more interesting. Sure, we can click through a bunch of peaks telling the program what's what. But wouldn't it be cooler if we could teach the program to tell the events apart itself?

Yes. Yes it would.

This is where machine learning comes in. Now there are various machine learning techniques, but the one I used in the program is called a 'Naive Bayes Classifier'. This is the sort of thing that's used in spam filters - it's a way of calculating, for example, what is the probability that a particular email is spam given that it mentions, say, penis enlargement or Nigerian princes.


Quick Introduction - Bayes' Theorem

Imagine there is a person. Sight unseen, what is the probability they are female? Well, given that there are roughly equal numbers of males and females in the world, the probability is about 50%. This is called the prior probability.

[Aside: in this example, we're assuming for simplicity that gender is a binary state.]

Now, we're given a new piece of information - this mystery person's height is 5'4 (five feet and four inches). Now the question is, what is the probability the mystery person is female, given they are 5'4? Well, females tend to be shorter than males (on average) so we might say that the mystery person is more likely to be female that male. But how do we quantify this new probability?

This is where Bayes' Theorem comes in. It looks like this

Basically, this calculates the probability of a hypothesis (H) given some observation/evidence (O). p(H|O) is called the posterior probability.

For our example, we have the hypothesis "female" (F) and the observation "h=5'4". So the terms for this example are:
  • p(F)  -  The probability of a person being female. This is our prior.  [ 50% ]
  • p(h=5'4|F)  -  The probability of a person being 5'4 given they are female. Or to put it more plainly, the probability of a female being 5'4.   [ ~15.1% ]*
  • p(h=5'4)  -  The probability of any person being 5'4. This can be calculated as the (weighted) average of the probabilities of a female being 5'4, and a male being 5'4. In other words, the probability can be expanded as

    p(h=5'4)  =  p(h=5'4|F) p(F) + p(h=5'4|M) p(M)

    where p(h=5'4|M) is the probability of a male being 5'4   [ ~1.8% ]*

Putting it all together, we calculate the (posterior) probability of the mystery person being female given their height is 5'4, as


As we expected, this probability is much greater than 50% - the mystery person is more likely to be female than male, given their height.

*Aside: These figures come from average heights for Americans aged 20-80 for 2007-2008.

Obviously the probabilities would vary depending on where the mystery person is from, how old they are, etc. Since all we know about this mystery person is their height, we shouldn't really make assumptions about their background. But for the sake of demonstration (and convenience), these values are sufficient.


Gaussian Naive Bayes - When the Training Set is Small

So how does this apply to noise removal?

Let's go back to EMI - we can ask, what is the probability that an event is EMI-type given it has an emi value of 30?

Using Bayes' Theorem, we have


Now we need to calculate the probabilities. We're going to assume we have a training set of (31) events that have already been classified. From this we can calculate/estimate these probabilities.

The prior, p(E), is pretty straightforward - what fraction of the events in the training set are EMI? This turns out to be around 22.6%

The probability of an EMI-type event having an emi value of 30 - p(emi=30|E) - is a little trickier to find. If our training set is small we may have no events with an emi value of exactly 30. But there could, at the same time, be several events with emi values slightly above or below 30.

Here then we have to make an assumption - we're going to assume that the distribution of emi values (for each event type) can be approximated by a normal distribution. This approach is called the 'Gaussian Naive Bayes', and is usually a reasonable approximation - if you look at the plot of emi values higher up, you'll see that they are roughly normally distributed.

So to find p(emi|E) we need to calculate the mean and standard deviation of the emi values for the EMI-type events in our training set. Obviously, when the training set is relatively small, the mean and standard deviation are going to have a large margin of error. So it's important to make sure we have a sufficiently large training set.

For p(emi=30|!E) - that is, the probability that an event that isn't EMI-type (signal or thin peak) has an emi value of 30 - we do the same procedure of calculating the mean and standard deviation of emi values for non-EMI events, again assuming a normal distribution.

For my training set of 31 events, the probability works out at p(E|emi=30)  =  97.96%

In other words, an event with an emi value of 30 is almost certainly an EMI-type event.


Machine Peak Classification

For the full Naive Bayes Classifier, you just repeat the process above for all the parameters, then combine the probabilities. For example, the probability that an event is noise given its log-meantime value and peak number, would be given by


It's 'naive' because it assumes all the parameters are independent of each other. This isn't strictly true, but it's an acceptable approximation.

In fact, in practice, rather than calculating the probability, we calculate the 'log-likelihood ratio'

where Pi are the various parameters. From that it's pretty straightforward to extract the probability that an event is noise given its parameters. Alternatively, we can just say that if the log-likelihood is greater than 0 - if the probability of being noise is greater than the probability of being signal - then the pulse is classified as noise.

The point here is that, rather than investigating the differences between the different pulse types ourselves, then hard coding how to tell pulse types apart - e.g. explicitly writing in the code that an event is EMI-type if it has an emi value greater than 25, etc. - we instead show the code some examples of different pulses, and it 'learns' to tell the difference itself. This saves us a lot of work. In fact, we never even have to know what the differences between pulse parameters are.

When you're using the program, you can have it display what the algorithm thinks pulses are (and the associated probabilities). So you can click through a bunch of pulses, telling the algorithm what they are, and the algorithm will tell you what it thinks they are. As you 'teach' it, the algorithm will get better at telling the difference between pulses. And when you're happy with how well its classifications match up with your own, you can press a button and it will classify the rest of the pulses for you. How cool is that?



Finding the Best Cuts

Once the algorithms is adequately trained, we could just set it to work, going through the data, classifying and removing noise events. But for their project, the students aren't looking at machine learning, or statistics, or anything like that. Instead, they're asked to look for cutoff conditions for each of the parameters - for example, a pulse is noise if it has emi > 25.

To get a better idea of the problem we can look at our histograms of the parameters, with the different pulse type plotted in different colours. This is where the pulse classification comes in useful - the more pulses we classify, the clearer the distributions of parameters.

When you're looking at the emi values, as in the plot higher up, the cutoff is pretty easy - there's a clear gap between emi values for EMI-type pulses, and those for signal and thin peaks. For other parameters, the distinction is less clear, and often we have to make some trade off between leaving behind noise events and removing signal events. For example


Here, if you want to keep all the signal, you leave behind 23 (thin peak) noise events (34%). If you want to remove all the noise you end up also removing 6 signal events (54%).

Aside: Since EMI-type events are easily distinguished and removed with the emi value cut, they are not taken into account when we look at cuts for the other parameters, which are meant to separate signal from thin peaks.

So the question is, how can we determine the 'best' cuts?

This is basically an optimisation problem - we want to find the cut for which the maximum amount of noise AND the minimum amount of signal is removed. To do this, we need to come up with some metric/way of scoring how effective a given cut is, then we look for the cut that scores best.

We'll define Ni as the number of noise events and Si the number of signal events removed by some cut Ci. Because there are generally more noise events than signal, we're going to score based on the fractions of noise and signal removed - (Ni/N) and (Si/S) - where N and S are the total number of noise and signal events respectively. So the optimisation problem is finding the cut that gives (Ni/N) closest to one and (Si/S) closest to zero, simultaneously.

One approach to this optimisation is a 'nearest neighbour search'. To make this clearer, we can look at a scatter plot with points representing the fractions of signal and noise removed for different cuts. In this set-up, optimising means looking for the point that is nearest the bottom-right corner of the graph (1,0).


We can find the distance of each point from (1,0) using Pythagoras, making the scoring function
which we want to minimise.

There is one caveat though - a cut that, for example, removes 10% of the signal and 60% of the noise will have the same distance score (0.41) as a cut that removes 40% of the signal and 90% of the noise. The latter removes more noise (almost all of it) but also removes more signal. In this case we have a choice - do we want to remove more noise, or save more signal?

To differentiate the two cases we can use the angle between the point and the noise-axis
So if, for example, we want to save as much signal as possible then we want to minimise the angle. Or perhaps for a better balance between removing signal and noise, we should prefer the cut that gives an angle closest to 45deg. For myself, I prefer to save as much signal as possible. In practice, however, having two or more cuts with the same distance score is rare.

When it comes to reporting the (distance) scores, they are re-formatted as
This gives a score of 1 for a cut that removes all noise and no signal, a score of -1 for a cut that removes all signal and no noise, a score of 0 for a cut that removes no signal or noise at all, etc. This score could be used in the optimisation algorithm (where it would need to be maximised), but the results would be the same, and the Pythagorean distance is easier to calculate.

Below is an example cut from this algorithm.
Notice that it sacrifices one signal event for 15 noise events, but prefers to keep one more signal event, rather than removing an extra 5 noise events.

Ultimately, this technique of removing noise with parameter cutoffs is less effective than the Naive Bayes, largely because it considers all the parameters individually/independently. However, there are things called 'support vector machines' (SVM) which would probably be ideal for this problem. They look for a dividing 'line' between signal and noise, as above, but they consider all of the parameters together. OpenCV has an implementation of SVM, so I might give that a try at some point.



Last Word

As it started out, I just wanted to throw together a quick GUI to make my life easier, without worrying about error handling or bugs or any of that other fiddly stuff. Then I decided I wanted to try out some machine learning. And the more I worked on the program, the more features I thought of, and added.

So now, I have a program designed specifically for removing noise from DM-ICE data. But since DM-ICE has already had its noise removed (as far as possible), the program is basically of no use to anyone. Well, except maybe to any students working on this project. But that would be cheating - hence why I'm wary of posting the code.

Still, it was a fun little project.


Oatzy.


[Keep an eye out for news on DM-ICE.]