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.]

Sharing a Burger Between Three

There are several ways to divide a burger evenly between three people.

Probably the best is to cut it radially (like a pizza). It can be tricky working out exactly where to make the cuts, but if you can pull it off, then all the pieces will be roughly identical (topping distribution notwithstanding).

But for the sake of arguing, lets say you want to divide the burger by making two parallel cuts: Where, then, do you make the cuts so that all three people get the same amount of burger?

NB/ This gets quite maths-heavy, so if you're not interested in that sort of thing, feel free to skip right to the end for the solution.



Geometry

For simplicity, we're going to consider the burger as a circle, and make the cuts so that each chunk has the same area. The two cuts are going to be the same distance from, and parallel to, the central axis of the burger, so we only need to consider the position of one of the cuts.

Here's the set-up

We work out the area of the cut-off as the area of the circular segment, minus the area of the triangle.


- Aside: Radians

Radians are basically an alternative way of measuring angles. For maths and physics they're generally more useful than degrees.

They're relatively easy - there are 2pi radians in a full circle, so 2pi radians = 360 degrees

1 radian = 180/pi = 57.3 degrees
1 degree = pi/180 = 0.017 radians, etc.

*    *    *

Back to the circle; with angle x in radians, the area of the circular segment is

The area of a triangle is half base times height..


- Aside: Area of the Triangle

We start by splitting the triangle down the middle, so that we have two identical right angle triangles

The height, l = r*cos(x/2)

The base, b = 2*(r*sin(x/2))

So the area of the triangle is (l*b)/2 = r^2 sin(x/2)cos(x/2)

Finally, use the identity sin(2x) = 2sin(x)cos(x)
to get

*    *    *

So the area of the cut-off is

and it needs to equal a third the area of the circle = 1/3 pi r^2.

So first, we need to find x satisfying

or

Once we have a value for x, we find where to make the cut from

Intermission

The thing about this equation is it doesn't have an exact, analytical solution - to find the solution you have to use numerical methods. Well, I say you have to use numerical methods; these days you can just type the equation into WolframAlpha, and you'll get a solution like *snaps fingers*

Which is nice. I even have the WolframAlpha app on my phone. But when I thought up this question I was on holiday in Sherwood forest, where there was literally no mobile singal.

So that was out of the question. And since I'm not in the habit of carrying a scientific calculator around with me, I was stuck with the basic calculator on my phone. It looks like this:

No trig functions, no square roots, no pi button. It doesn't even do brackets, or have a memory function. Luckily, I am in the habit of carrying around a notepad and pen.

Anyway, there are two ways of working this out with only a basic calculator. The first is 'easier', but only if you know some stuff, and the numbers happen to be nice (in this case, they kind of are). The second is harder, in that it requires more number crunching, but it'll work with any numbers, and can be more precise.

Again, feel free to skip to the solution if you're not interested in the gritty details.



Method One

First of all, here's a graph of the two sides of the equation

hand-drawn with Skitch

We want to find the point at which the two graphs cross. We can see that that happens somewhere between 2pi/3 and pi (120 and 180 degrees). So, lets make a guess that it's exactly halfway between these two values: 5pi/6 (150 degrees).

For the right hand side of the equation: 5pi/6 - 2pi/3 = 0.52

NB/ I'm using pi=3.1416 (rounded to 4 decimal places). If you prefer, you could use the approximation 22/7. The result should be roughly the same.

For the left hand side of the equation, we need to work out sin(5pi/6)

At A-level, we were expected to memorise sin() and cos() of angles 0, 30, 45, 60, 90, and 180 (degrees). We were also expected to know the formulas for sin() and cos() of sums of angles. For sin(), it works like

Why is this important? Well 150 degrees = 180 - 30 (5pi/6 rads = pi - pi/6)

So

And since 0.5 is pretty close to 0.52 - less than 5% error - we can accept the convenience of that answer and say it's close enough.

So our approximate value of x is 5pi/6 = 2.618

[Incidentally, the identity for sin(2x) is just a special case of the above, with a=b=x; i.e. sin(2x) = sin(x+x) = 2sin(x)cos(x)]


Now we just need to work out l/r = cos(x/2) = cos(5pi/12)

For this one, 5pi/12 rads = 75 degrees = 45 + 30, so we can use

So

And we just have to evaluate that. But we don't have a square root button. Now, I just happen to know that sqrt(3) ~ 1.73 and sqrt(2) ~ 1.41.

But I'm just weird like that. Let's say you don't. How do you work it out?


- Aside: Square Roots

There are several ways of working out square roots with just basic operators. For two easy examples:

The first is 'Trial and Improvement' - pick a number, square it, does that give the right answer? If not, pick another number based on whether the last guess was too big or too small.

For example: sqrt(3)

1.5 -> 2.25 -> too small
1.7 -> 2.89 -> too small
1.8 -> 3.24 -> too big
1.75 -> 3.0625 -> too big
1.73 -> 2.9929 -> too small
1.74 -> 3.0276 -> too big
1.735 -> 3.010225 -> too big
1.7325 -> 3.00155625 -> too big
1.732 -> 2.999824 -> too small
etc.

The second method is the "Babylonian Method". It's more systematic, and can converge to the correct answer quicker than guessing. But it can be irritating if your calculator doesn't have a memory function.

It uses the recurrence relation

Basically, you make a guess xn. Divide the number you want to square root (S) by xn. If xn is lower than the actual square root, then S/xn will be greater than it. That means the actual root will be between xn and S/xn, so we make the next guess xn+1 the average of these two values. Repeat until x is sufficiently accurate.

For example: sqrt(2)

x0 = 1.5 -> 2/1.5 = 1.33
x1 = (1.5 + 1.33)/2 = 1.4166.. -> 2/1.4167 = 1.41176..
x2 = (1.4167 + 1.41176..)/2 = 1.41421.. -> 2/1.41421 = 1.41421..

*    *    *

Whatever way you do it, you repeat the process until you get the degree of accuracy you're happy with.

You should get the answer around l/r = 0.259



Method Two

We go back to the equation sin(x) = x - 2pi/3

We still have to find the solution numerically, we still don't have a calculator with a sin() function, and this time the numbers don't work out nicely.

So, the question is, how do we calculate sin(x)?


- Aside: Taylor Expansion

The Taylor Expansion of a function is a way of fitting a polynomial (sums of powers) to a more complicated function. It works like this

Basically, it gives a way of converting a function we can't calculate into an infinite sum of powers of x, which we can calculate.

It's usually expanded around the origin (x0=0), since the equations work out neater. But you can do it around any point, x0=a. This is useful if the value you are trying to calculate is far from x=0. The closer x is to x0=a, the quicker the sum converges.

Even though the expansion is an infinite sum, it's usually sufficient to just take the first few terms, since each additional term makes a smaller and smaller contribution to the sum.

So the trick is working out how many terms you need to include to get some desired level of accuracy.

*    *    *

In this case, I'm going to use the Taylor Expansion of sin(x) around x0=pi, since the approximate value (2.6) is nearer to pi than 0.

Here's what the expansion looks like

So, how many terms do we need to include?

Here's what the graph looks like for different numbers of terms

plotted with WolframAlpha

For a value around 2.6, it can be shown that including the first two terms is correct to ~3 decimal places; the first three terms is correct to ~5 decimal places; the first four terms to ~7 decimal places, etc.

So I would probably go to the third term (for 5dp), but only take the result to 3dp.

NB/ We shouldn't get too hung up on getting an extremely accurate value for x, since we're already getting rounding errors from the factors of pi in the expansion. Also, since we're calculating x to 5dp, we should use pi=3.14159

That means we want to solve

which can't be solved exactly.

So, for finding the correct value (without WolframAlpha), we can use any root-finding method. For what it's worth, I used Trial and Improvement; the other methods are easier with a computer.

But, note that the function is decreasing

So if the guess gives a value greater than zero, you need to increase the value of x (and vice versa).

If you run through all that (I won't go into detail), it gives a value around x=2.605


Alternatively, you could expand around x0=5pi/6 (if you know/can work-out sin and cos of 150 deg without a calculator).

In this case you'd only need up to the term in x^2 (correct to ~4dp). Using this expansion would mean solving a quadratic equation, which is easy. But using this expansion can introduce more rounding errors from the factors of sqrt(3). It's a matter of preference, I guess. The answer should be about the same.


Finally, we need to calculate cos(x/2)

Again, we use the Taylor Expansion to calculate cos(). In this case, we're doing the expansion around x0=0; the expansion is

In this case, you just keep adding terms until the result remains approximately constant to some desired degree of accuracy (3pd).

This gives a value around l/r = 0.265



So What is the Real Answer?

Once I got to somewhere where I could get at WolframAlpha, I checked the real numbers; here are the results:

The approximation of x from Method One (2.618) is an over estimate by ~0.5%, which is relatively acceptable. The approximation from Method Two (2.605) is correct to 3 decimal places, which is definitely acceptable.

And for the value of l/r

From Method One (0.259), the approximation is an under estimate by 2%, and correct to 2 decimal places, so is probably acceptable. The approximation from Method Two (0.265) is, again, correct to 3 decimal places. So that is also acceptable.

So, if the numbers happen to be convenient and you know some trigonometry, you're probably as well using Method One. If not, or if you just want more accuracy, then go for Method Two.



Applying the Results

The results are actually quite nice, in terms of practical application (dividing up a burger). The ratio of the radius (0.265) being close to one quarter, you find the cuts like this

That is, find the central axis, then find the (imaginary) line halfway between the centre and the edge - make the cut halfway between the centre and this imaginary line (maybe cut an extra hair's breadth towards the edge). Repeat on the other side.

Easy.


So, now you know. Obviously, all this applies to dividing any circular thing evenly between three people. You could probably even adapt the methods for sharing between even more people.

And in theory, You could do all this with just pen and paper (no calculator). Though you probably wouldn't want to. I know I wouldn't..


Oatzy.


[Wow, I really managed to stretch that one out.]

Where are the Carriages?

<rant>

So it's around 9am, and you're waiting for a train. The trains in your area are a bit scummy, but whatever, you can't drive, and you've got to get to work/uni somehow. It's pretty busy, what with it being 9am, and when the train finally pulls up.. it's a single carriage.

Now a single car has around 50 seats, and this train is packed tight, with people standing in the aisles and the doorways, being forced to get intimate. And with everyone on board the conductor can barely get through the door. There are clearly enough people on this train to fill two cars, with people still having to stand.

So what gives? You're pretty pissed about having to stand for half an hour, and you've been inadvertently touched by strangers in ways you're not comfortable with. So you get your complaining hat on, and you turn to the internet to unleash the fury.

The train company's website directs you to their Twitter feed, where a poor public relations person is taking a barrage of vitriol from other angry commuters, while doing their best to remain polite and professional.

Several other people have already made your complaint, but all they're getting in return is "Apologies we try to avoid this where possible", or that the short trains were "due to operational reasons". And that's not really a satisfying response. You're not even sure it means anything.

</rant>


So Where Are The Carriages?

Disclaimer: I have no idea how the train company actually operates. I just like to write blogs about applied maths.

Here's the setup - you're in charge of logistics for the train company. You have a fixed number of trains/cars, and, for each day, a list of services and expected numbers of passengers.

How do you apportion the carriages between the services?

Obviously, every service needs at least one car, and services with more passengers need more cars.

So if you had, for example, three services with 50, 100, and 150 passengers respectively, and 6 cars, then it's easy to divvy them up - one for the first, two for the second, and three for the third.

But in general, it won't be possible to divide up the cars exactly like that; so what do you do about remainders?

This is actually similar to the problem of apportioning parliamentary seats between states in the US.

Currently there are 435 seats in the House of Representatives, which need to be shared out between the 50 states according to each state's population - the idea being that each seat should represent roughly the same number of voters.

There are various methods for doing this, but I'm only going to go over two of them.


The Hamilton method

Hamilton's Method is the easiest and most intuitive.

Say you've got only two services - service A has 85 passengers, service B has 115 passengers (200 passengers total). And you happen to have 4 carriages - a total 200 seats. So seats for everyone! Except service A needs 1.7 cars, and service B needs 2.3 cars. And you can't divide a car into two chunks.

But for starters, you can give one car to service A (50 seats), and two to service B (100 seats). So what about the fourth? Regardless of which service you give it to, some people are going to end up having to stand. So what you want to do is minimise that number.

For service A, 35 people need a seat. For service B, 15 people need a seat. So it makes sense to give the last car to service A. 15 people still have to stand, but - short of building more cars - there's really nothing you can do about that.

This is also known as the largest remainder method.

What if the passenger numbers were 75 and 125? Well in that case, I guess it would have to be a judgement call.

Anyway, that's the basic idea of the Hamilton method - divide up the cars as far as you can, then give the remaining car(s) to the service(s) which 'need them most'.

Things get trickier when you're dealing with larger numbers of passengers and cars, but ultimately it's pretty straight-forward.


The Huntington-Hill Method

Where trains differ from apportioning of parliamentary seats, is that a carriage has a fixed capacity, whereas the number of people a seat can represent is free to change.

And this means that you won't encounter some of the 'quirks' that arise from the Hamilton method - such as the Alabama Paradox (where increasing the total number of seats can mean a state losing seats), or the Population Paradox (where increasing a state's population can result in that state losing seats).

[You won't encounter them, but they're worth mentioning 'cause they're pretty cool.]

Huntington-Hill's Method - the one currently used by the House of Representatives - is, generally, a more powerful method than Hamilton's, and isn't susceptible to the various paradoxes. Another benefit is that it can be set up to guarantee that each state will get at least one seat.

It works by assigning seats, one by one, based on each state's 'priority quotient' - itself based on the state's population, and a 'modified divisor' based on the number of seats already allocated to that state.

But one advantage Hamilton's Method has over HHM, is that it will always give each service its ideal number of cars, either rounded up of down to the nearest whole number. So if a service's ideal number is 3.67, then Hamilton's method will assign either 3 or 4 cars to that service.

HHM, on the other hand, can 'violate quota' - that is, it can result in a given service being assigned more or fewer cars than it's ideal number (e.g. the 3.67 service might end up with only 2 cars). But this problem only occurs when the number of cars is fixed prior to apportioning. And sadly, it usually is.


As an aside, it turns out it's impossible to find a 'perfect' method of apportioning - one with neither paradoxes nor quota violating. Incidentally, the mathematical theory surrounding voting is quite fascinating. You know, if you're into that sort of thing.


Unfortunately..

Either method would work perfectly well. The problem is, while rounding up or down to the nearest car may seem trivial, the 50 seats that that one car represents can be a significant gain/loss for a service.

The lack of precise information can cause problems as well. You can only estimate how many passengers a given service will have in advance, and if you under-estimate, people are gonna be pretty pissed.

And on top of that, the numbers have to be recalculated from time to time as passenger numbers change..

Basically, getting the number of cars right can be tricky.


Of course, this all assumes that Northern Rail (or whatever company) has enough carriages to adequately satisfy its needs in the first place. Some how I doubt that's the case.

So ideally, NR needs to work out how many cars it actually needs, and build them. But I can't see that happening. And even if it did happen, it'd probably mean higher fares. And fares are pretty bad as they are.

Really, though, NR could do with a complete makeover just in general, given how shitty it is. Seriously, you should see the difference between them and the London Midland service. Ridiculous.

Incidentally, I also have strong feelings about the buses; don't even get me started. Trams are alright though - one every 10mins, seldom short of seats, sufficient leg-room.. bliss.


tl;dr - Where are the carriages, then? As it turns out, it's probably just that they don't have enough to go around.

This is what we get for privatising the trains...


Oatzy.


[My complaining hat is a trilby. I like to look bad-ass while I'm complaining.]

Gift Wrapping

In these times of austerity, one has to be economical with one's wrapping paper.

Yeah, I know, this would have been much more useful about a week ago. But I just never got around to it. And I use the word 'useful' loosely; any saving one might get from these 'techniques' will most likely be insignificant.

Nonetheless...


Rectangles

All gifts - no matter shape or size - can be wrapped with a rectangular piece of paper.

This is easily proved. However doing so may not be the neatest or most efficient approach. For example, it's easy to efficiently wrap a cuboid (book, DVD, box) with a rectangular piece of paper, it's less easy to wrap a bike (efficiently) with a single, and in this case, very large piece. However, in the case of a bike, one could use several, smaller pieces of paper for more efficient wrapping.

For a ball, one might be tempted to say that a more unusual shape - maybe a truncated icosahedral shell - would be more efficient. However, cutting out such a shape would result in lots of little slivers of waste, and would take a lot of time and effort. Instead, it's usually easier, and not greatly inefficient to just use a rectangular piece, and some creative folding and scrunching.

So here's the point - for a given gift, we can define a rectangle (or set of rectangles) which most efficiently wraps that gift.

For a cuboid of sides L(ength), W(idth), H(eight), the minimum paper would be: (H+L) by 2(H+W), with a bit extra for overlap.

For a cylinder of dimensions H(eight) and R(adius), the paper would need to be: 2PI*R by (H+2R), and a bit.

[Proof that these cuts are optimal is left as an exercise for the reader.]

For other shapes, this may be more tricky to determine, but for the sake of arguing, we will assume we know all the required rectangle's dimensions.


Bin Packing

So for a given collection of gifts to be wrapped, we have defined a set of rectangles. We now have to determine the most efficient way to cut these out from rolls of wrapping paper - i.e. a large rectangle, ~ 1m x 4m.

This is similar to solving a two-dimensional version of the Bin Packing problem.

In the Bin Packing problem, we have several bins (cuboids) of fixed dimensions, and a collection of smaller cuboidal objects. The problem is to pack the objects into as few bins a possible - i.e. find an efficient way of fitting cuboids into boxes.

The problem is NP-Hard, meaning that finding the optimal solution can take a very long time. However, there are algorithms which are fast and give adequately optimal solutions.

One, in particular, is called the First Fit algorithm, and it works something like this:

1) Sort the rectangles into size order
2) Cut out the largest pieces still in the set
2 i) If there is not enough paper left on the current roll, move on to the next/start a new roll
3) Repeat (in descending size order) until there are no more pieces left to cut.

There are also instructions for how to place the rectangles, for example first fit from bottom and right, to left and top.

This article gives a more precise description of a similar algorithm, though that one is defined for a 'roll' of non-restricted dimensions.


Realistically

That's all well and good mathematically. But in the real world, it's likely more trouble to solve the problem than the saving in paper is worth.

However, the general principle of the algorithm can be loosely applied (if indeed it isn't the approach you already use) to improve on efficiency.

In general, you can easily sight-guess which gift in going to require the largest amount of paper. And from that you can apply the First Fit algorithm.

One change I make is, cut out and wrap the largest one left, first. Then wrap the gift(s) which best fits the off-cut.

This is particularly good for rolls, where you want to unroll and use it a 'slice' at a time (rather than unrolling the whole thing).


Anyway, something to think about next time you happen to be wrapping a large number of gifts.

Even better if you can get some of that Tesco wrapping paper with the grid on the back.


Oatzy.


[No wasted paper from my wrapping this year. Just saying.]