Showing posts with label Maths. Show all posts
Showing posts with label Maths. Show all posts

Wednesday, October 2, 2019

Generate an STL 3D model with a surface described by an equation

I recently needed to generate a 3D model for printing based on an equation and couldn't find software to do the job.  In the end the easiest thing to do was to write a quick Python script called "stl-surface.py" to do the job.

At first the task seemed daunting, but generating the STL file is rather easy with the use of the numpy-stl library. The equation for the surface is calculated over a structured square/rectangular grid and each little cell of the grid is split in two to form triangles. You use many many of these small triangular faces to construct the model. All you have to do is to generate the coordinates for the 3 vertices that make up each triangle and numpy-stl does the rest. The bottom and sides are flat which makes things rather easy, and the top is defined by the equation.

When constructing the model it's important to list the coordinates of each face in a counter clockwise direction when looking at the model from the outside. This allows numpy-stl to later calculate things like volume of the model and centre of gravity and also test if the surface is closed (although they use a slightly buggy way to test this).

I've shown a few examples of the generated models below.

10 - 2*(1 - math.cos(2 * x * math.pi / 20)) * (1 - math.cos(2 * y * math.pi / 20))
Dimpled 3D model
10 + 2 * math.cos(math.sqrt((x - 50)**2 + (y - 50)**2)/2)
Rippled 3D model

Here you can see how the triangles are assembled to construct the model.
The triangular faces that make up the model

The triangles that form the edges of the model
Edge of the model

From the top you can see the grid points where the surface function is evaluated.
Top of the model

I hope this code can help someone else. I may not be the easiest thing to use and if you need help getting started get in touch.
Get the code!


Sunday, September 15, 2019

Self Normalising Price Markdowns For Short Life Products

When selling products that have a relatively short shelf life, it becomes necessary to reduce the price on items as they approach their use by date so they sell and aren't thrown out.  This is done for several reasons.

  • Sending plastic packaging and food that's safe to eat to landfill is a waste of resources, and in 2019 we can do a lot better than that
  • If a business can reduce its waste, they can save money on garbage collection
  • It's an acknowledgement to customers that the product isn't as flexible as one with a longer use by date
  • It adds a little excitement for customers if they can buy a product that they normally wouldn't buy because it's normally too expensive

In my experience, supermarkets perform markdowns in an incremental way to make sure that stock isn't wasted, and to get as much money for the product that is possible.  This was usually accomplished by reducing the price by 20% on the 2nd last day of sale, 30% on the morning of the last day of sale, 40% at lunch time on the last day of sale and then going 50%, 60%, 70%, 80%, and 90% until the close of the store.  In essence, it's a Dutch auction. This works well but has problems.  People buy the most popular items first leaving an assortment of unappealing products towards the end.  As fixed percentages are used across the board, it means you might be reducing some items more than you need to and some less than you need to.  It's also reasonably labour intensive.

If you track markdowns and plot a graph of the most common ones you end up with a graph like the one below, which I'll use to demonstrate a point.
markdown percentage graph
Markdowns shouldn't waste time or money

In the graph above you can see that average markdown is around 50%.  Traditionally you'd start with a 20% markdown, but the graph shows that in this case almost no one will buy the product at that discount.  This is the "Waste of Time" section of the graph.  Most likely you'll have to come back and do another markdown on the product.  Conversely if you go straight to an 80% markdown you've wasted money because you mostly likely could have sold the product for more.  The "Waste of Money" section.  Ideally most of the markdowns should be in the 30% to 70% range, the "Sweet Spot". But how do you calculate the markdown percentages?

I'm going to demonstrate a method that doesn't necessarily generate optimum markdown percentages, but does generate percentages so that popular and unpopular products sell at roughly the same rate with a reasonable chance of selling on the first markdown.

Let's start with some Lamb Hearts, and to keep things simple, from most to least recent the markdown percentages are as follows, 98%, 66%, and 6%.  This process uses a method called Kernel Density Estimation to generate a smooth curve of markdown probabilities.  This means that at each of those percentages above you place a predefined shape, in this case a Gaussian distribution.
markdown percentage graph
Initial placement of kernels

You can see the distributions listed above.  You may notice that because the grey and blue curves are truncated at the sides of the graph the area under them is less.  This means that they'll have less of an influence on the result.  This can be corrected with a simple scaling process though.
markdown percentage graph
Compensation for truncated kernels

The areas under the graphs above are now all the same, but is this what we want?  Most recent markdowns should have a greater impact on the process than ones that were performed long ago.  To correct this, each peak is weighted by 3 for the most recent, 2 for the next most recent, and 1 for the oldest markdown.
markdown percentage graph
Weighting the kernels for recency

Adding each of these graphs gives the final markdown distribution in yellow.
markdown percentage graph
Summation of the kernels

Now we perform a process called integration to get the light blue line.  For those of you unfamiliar integration it basically records the area under a graph. For example at the 20 point on the bottom axis the blue line records the area under the yellow line up until that point.
markdown percentage graph
Integration of the distribution to produce the final curve

The next step is to introduce the concept of a markdown level which ranges from 0 to 10.  0 corresponds to a 0% markdown and 10 corresponds to a 100% markdown, but in between, things are very different.  As a demonstration we'll calculate the markdown percentages for markdown levels 2, 5, and 8.  The process is quite simple.  Find the levels on the left hand side of the graph and project them across to the blue curve and then down to the bottom of the graph.  This gives percentages of 53, 80 and then 94.
markdown percentage graph
Calculating markdowns percentages from the graph

For another demonstration let's try something popular like chicken breast that has recent markdowns of 30%, 10%, and 20%.  This leads to percentages of 10, 22, and 33.
markdown percentage graph
Calculating markdowns percentages from the graph

Like I said above, this isn't designed to generate an optimum markdown, it's designed to create a markdown specific to each product so that they'll sell at similar rates.  In the demonstration above the situation may warrant a level 2 markdown.  This means that a product that's hard to sell like lamb hearts gets a 53% markdown, while a popular item like chicken breast only gets a 10% markdown.

So a strategy may be to use a level 2 markdown on the second last day of sale, and then taper throughout the trading day from level 3 to 10 on the last day.

The demonstrations above have been quite simple and don't give a full picture of how powerful this method can be.  Let's put together a complicated situation with fake data to test it.

A new flavour of sausages is introduced and we are unsure how they'll perform.  They initially use an assumed prior that isn't mentioned above, but it quickly becomes obvious that they are a good seller and don't require large markdowns, with the first level 2 markdown being in the region of 15%.  After 500 recorded markdowns a new similar cheaper line is introduced and the sales of the original line suffer.  Before long the system calculates that the first markdown at level 2 needs to be about 55%
markdown percentage graph
Animation of the markdown calculation process

The animation above is generated from 500 random markdowns with an average of 20% and then 500 random markdowns with an average of 80%.  The calculations take into account the last 100 markdown and are weighted from 1 to 100.

There are a couple important points to consider.  Firstly the process uses a lot of statistical methods and equations but isn't really grounded in theory.  It uses the properties of the methods to create a stable system that fits the requirements that are specified. The next thing to remember is that this process will follow the markdowns.  If you start the markdown process at level 8 without even trying something less than a level 5.  The system will generate larger and larger markdowns.  A stabilising term can be added to prevent this.

It's not a perfect method but it's a lot better than some of the optimised markdown systems I've seen.  I actually trialed this process a while ago(I have no deep access to the computers and had to do it manually) in my store and it's fascinating to see it converge on higher and lower markdowns based on the saleability of the product.

Saturday, September 14, 2019

Reducing Evaporation Rates With Rippled Surfaces

A theory that I want to test out is that adding a ripple to a surface makes it easier to clean in certain situations.  I work in retail, and sometimes the cleaning of non food safety related issues are delayed due to more important tasks.  In particular I'm talking about drips.  You might see these on metal trays under bottles of milk in the fridge, or on plastic sheets under the shelves in the meat department.  In both cases drips are meant to be caught on a surface that is easy to remove and clean.  The problem is that if you leave a spill too long it dries and become hard to remove and requires vigorous cleaning which may damage the surface and waste time.

To avoid this, you want to slow the rate of evaporation.  This can be done by decreasing the surface area by increasing the depth.  A rippled surface is perfect for this.  You essentially make little cups to hold the spill.  Another requirement that would be nice to include is no tight internal corners.  Anything that is designed to be cleaned shouldn't contain a concave surface that you can't get a finger into.  Anyone that knows me well understands that I think in equations, and the one below matches our criteria perfectly.

$$$z=10 - 8 \left(\dfrac{1-\cos(2\pi x / 20)}{2}\right)\left(\dfrac{1-\cos(2\pi y / 20)}{2}\right)$$$

Let's talk about what this equation means.  The two large sections in the brackets are periodic terms that create ripples in the x and y directions.  The 20 means that this ripple will repeat every 20mm.  By subtracting the cos term from one and then dividing that result by 2, the term in the bracket will range from 0 to 1.  Multiplying the two brackets together will also give a result between 0 and 1.  By multiplying this by 8 we now have a function that ranges from 0 to 8.  The 10 describes the maximum of the equation.  By subtracting the rest of the equation from 10 we now have a surface that ranges from 2 to 10 above zero.  The important things to remember are, 20 specifies how wide the ripples are, and 8 describes how deep they are.

Sometimes though it helps to have the real thing in your hand to test so I created a 3d model to send away for 3d printing.  There doesn't seem to be anything out there to create a 3d model from an equation so I had to write some software to do that.  I'll post that when I tidy it up and comment it properly.

I'm not made of money so the model is only 100mm x 100mm x 10mm (a volume of 100mL) with 20mm wide ripples that are 8mm deep.  By a stroke of luck, I happened to write the software is a way that easily calculates the volume of the model.  In this case the model is 80 mL.  As the bounding box of the model is 100 mL this means the volume of the 25 little cups is 20mL.  Each one holding 0.8 mL.

3D model of a rippled surface
Rippled Drip Tray

I don't have the print yet, but it has been done and photos sent to me.  In theory this surface should hold 20 mL of liquid and covers an area of 100 square centimeters, so as a test I poured 20 mL of water on a flat surface and it spread to cover 180 square centimeters.  So already the ripple pattern has reduced the surface area by 45%

3D printed model of a rippled surface
3D Printed Drip Tray Top

That may not sound like much, but the ripples can be made deeper.  If they were 3x times deeper (24mm) the surface would hold 60mL.  Once they get too deep though, you would need to make the ripples wider to make cleaning easier.  Changing the width of the ripples doesn't effect the volume that the surface would hold though.

3D printed model of a rippled surface
3D Printed Drip Tray Top

In this demonstration I've shown the surface as a solid block with depressions.  In reality you'd use something like a polypropylene sheet moulded to this shape.  It would give an object shaped similar to an egg carton.

3D printed model of a flat surface
3D Printed Drip Tray Base

Anyway, this is just a thought that I wanted to explore.  Maybe it'll work out, maybe it won't.  Either way the process was enjoyable.

Sunday, October 15, 2017

ABC Logo Lissajous Curve

Over the last couple of days I've seen a few people mention how similar the ABC Australian logo is to the new Disney Movies Anywhere service. I don't know much about the legal side of things, but I thought an explanation of why the ABC logo looks the way it does might be interesting.

The shape of logo is called a Lissajous figure or curve. The shape is generated by a parametric equation where the x and y coordinates are sinusoidal.  The frequency and phase relationship between the two equations for x and y determine the shape.  In the case of the ABC logo, the frequency of the y coordinate is 3 times that of the x coordinate and there is a 90 degree phase shift (I'll clarify this later) applied the y equation.

$$$x=cos(t)$$$ and $$$y=cos(3t+\pi/2)$$$

To make things clearer I've put together an animation. As the vertical bar sweeps across the screen it will intersect the x and y equations. The y coordinate is projected across and the x coordinate is projected across and up. The Lissajous figure is drawn where the project lines intersect.

Tracing a 3:1 Lissajous Curve x=cos(t)  y=cos(3t + $$$\pi$$$/2)
The phase shift is also very important to the shape. In the animation below you can see how the it changes as the phase shift is cycled across all possible values from 0 to $$$2\pi$$$.

Changing the phase relationship of a Lissajous Curve x=cos(t) y =cos(3t + $$$\delta$$$)
These curves aren't just a mathematical curiosity, they have a real world application. They used to be a very important tool for broadcast engineers. If two signals are feed into an oscilloscope (a tool that plots electrical waveforms) while it's in x-y mode, the Lissajous curve on the screen will reveal things about the signals.

The first thing to notice is that the number of horizontal and vertical lobes indicates the ratio of the frequencies. If the ratio of frequencies is rational (can be expressed as the ratio of two integers) the curve will be stationary. If not, it will slowly rotate like the second animation. If the two signal are meant to be locked together so that one is exactly 3 times the other like in the example above but the curve rotates, you know there's a problem. The rotation rate of the curve tells the engineer the deviation from the desired frequency. There are simple lookup tables like the one below that show what the curves should look like for a given phase shift and frequency ratio.


Frequency Ratios and Phase Differences

Earlier I said I'd elaborate on phase shift. The main thing to note is that you are working with two different frequencies.  1 degree of phase shift on one signal will take a different amount of time to 1 degree of phase shift on the other.  So it's important to not just note that there is a phase shift of 30 degrees, you have to specify what waveform you are referring to.  That's why a lot of the table results may be different from what you measure.  In a mathematical sense, it also makes a difference if you are talking about sine or cosine signal as one is a phase shifted version of the other.

Just to clear up another thing as well. The second animation is generated by plotting a waveform frame and then changing the phase an repeating this process. This is what makes it rotate. By chance though, this is exactly what you would see if the frequencies weren't locked together. A time varying phase is no different to making a small deviation to the frequency.
 
The code to generate the animations above is located here.  I did it as a quick project so didn't comment it.
 

Friday, May 19, 2017

Testing All The PINs On A Lock Box With A Forgotten Code

A long time ago I did a blog post on PIN coded lock boxes that you put on your house to hold a spare set of keys. Their main function is to give easy access for emergency services in case they need keys to get in if an elderly relative has a fall and can't get up. The point of my post was to demonstrate that although they have 10 buttons, look secure, and you can choose how many digits are in the PIN, they only really have 1024 PIN combinations. This is because the numbers can only be used once in each PIN.

Lock
PIN coded key lock box
The math is explained in my initial blog post but the results are below. I opined it would take about 4 hours to try all the combinations, and in this case it would require pressing the open button 1024 times, the clear button 1024 times, and pressing digits (0*1 + 1*10 + 2*45 + 3*120 + 4*210 + 5*252 + 6*210 + 7*120 + 8*45 + 9*10 + 10*1) = 5120 times.  In total there are 7168 button presses. That works out at about one button press every 2 seconds.
Table
Number of PINs vs PIN length
At some point in the 5 years since I first wrote about these locks (the time flies doesn't it) it occurred to me that I was being inefficient. Testing if the box opens doesn't clear the current code. So I  can test multiple PINs at once by just chaining them together. It's similar to the De Bruijn sequence, a mathematical tool that can be used to brute force another type of PIN. In this case however you have to reset the lock after a few buttons are pressed. To illustrate the process I'm trying to explain, imagine that the clear button has been pressed and I then enter the numbers 0 through 9, testing if the lock will open in between numbers. I've actually just tested the codes (Null) (0) (01) (012) (0123) (01234) (012345) (0123456) (01234567) (012345678) (0123456789). I've tested 11 codes with one press of the clear button, 11 presses of the open button, and just 10 digit presses. The whole process won't be this efficient but it'll be better than nothing.

But where do we start? In a perfect world, you may see by looking at the chart below and table above, that the best we can do is 1 sequence testing PINs of length 0-10,  followed by 9 sequences testing PINs of length 1-9, 35 sequences testing PINs of length 2-8, 75 sequences testing PINs of length 3-7, 90 sequences testing PINs of length 4-6, and 42 sequences testing PINs of length 5.

This will still result in 1024 presses of the test button, only 252 presses of the clear button, and (1*10)+(9*9)+(35*8)+(75*7)+(90*6)+(42*5) = 1646 presses of the number buttons. For a new total of 2922 button presses, or about 41% of the original estimate.
Graph
Distribution of PIN lengths
There's no guarantee that these sequences exist though and I had no deep understanding of how to create them, so I tried the first thing that came to mind and got lucky.

First, generate all possible combinations for each PIN length and then sort the numbers in each PIN. Then sort the PINs for each length comparing them element by element.  Start in the middle with the PINs of length 5 as there are more of these than the others. Then take the PINs of length 4 and go through them one by one from the start and place them to the left of first 5 digit PIN that could follow it. For instance (1, 5, 6, 7) might go to the left of (1, 2, 5, 6, 7) as only a 2 would have to be pressed to get to the 5 digit PIN. Do the same for the 6 digit PINs and place them on the right of the 5 digits. In this case (1, 2, 3, 5, 6, 7) might go to the right of (1, 2, 5, 6, 7) as only a 3 needs to be pressed. Repeat this left right procedure until all PINs are processed. It will look like a mess, but if you sort the sequences by length you will get a spreadsheet that looks like the one below (zoomed and rotated to fit). Look familiar? It's reflects the distribution graph above, and it shows that the sequences can be generated.

Spreadsheet
Rotated spreadsheet of the sequences
I've placed my code in a Github Gist for you to generate your own sequences in case there are a different number of buttons on your lock. If you have a lock with 10 buttons I've already generated the file for you. It looks a little different from the output of the Python file because I've done some some find and replace operations in Notepad++. I've used the word test instead of unlock as well.

Update1 There is a MP3 file of the instructions in this blog post, Linux Text To Speech With Saved Audio.

Update2 I've sorted the sequences by how often the 4 digit PINs appear in the 2009 Rock You data breach. There are two new files that cover the all the PINs. One of them has the sequences sorted by the 4 digit PIN popularity but still groups the sequences lengths, while the other sorts by PIN popularity only. An Excel file with all this data is also included so you can sort the data however you want. The associated files are in this Google drive folder. It could be argued that Rock You users and the users of lock boxes are a completely different demographic, and it's true. However, their users will exhibit similar behaviours like using birthdays and years that make the effort worthwhile. Although the instructions for the lock suggest selecting a PIN between 4 and 7 digits long, 4 digit PINs have been targeted as they are what people tend to think of when you mention PINs, even though 5 digit PINS are more secure.

So what about those locks you see in banks and other offices that have 14 buttons.  I've done the maths so you don't have to, but they can be broken with 50316 button presses. So 4 extra buttons buys you an increase in security by a factor of about 17. Not that it matters, whenever I've seen these used, people aren't too discreet about entering the PIN.

Door Lock
PIN Door Lock

Friday, April 21, 2017

Efficient Centroid Calculation for Discrete Areas

A project I'm working on requires the repeated calculation of weighted centroids for arbitrary regions in an image.  Calculation of centroids is fairly straight forward but computationally expensive.  To calculate the x and y coordinates you first need to calculate three other values, the denominator, x-numerator and y-numerator.  The denominator is simply equal to the sum of all the values in the region, the x-numerator is the sum of all the values in the region multiplied by their x coordinate, and the y-numerator is the sum of all the values in the region multiplied by their y coordinate.  From these values the x coordinate of the centroid is equal to the the x-numerator divided by the denominator likewise for the y coordinate of the centroid.  

You can see from this description that there are a lot of time consuming multiplication operations to perform for each area that you want to calculate the centroid for.  As it turns out though, you can perform two cumulative summations along the rows of the input data as a first step, store this and then perform the centroid calculation by only accessing values at the boundary of the region. For reasons I'll explain later a column of zeros with x coordinate -1 is added to the start of the data.  A cumulative summation is performed along each row to obtain $$$P$$$. A cumulative summation along each row of $$$P$$$ is performed to generate $$$Q$$$ The following equations describe the process.  Before we get started, the derivation of these  equations can be found in this small article I wrote.
Equations
Centroid Calculations
It may not be immediately obvious what these mean so I'll give a quick explanation.  First of all each of the above calculations are done row by row and added together. This is what the summation symbol over $$$y$$$ represents.  For the denominator the value to be summed is the difference between two values on a row of $$$P$$$. For the y-numerator the calculation is the same as the denominator but each row is multiplied by the y coordinate.  The calculation of the x-numerator is a little different. It's similar in that values of $$$P$$$ at the boundary are subtracted and multiplied by $$$x+1$$$ but the additional step of subtracting the difference of the values of $$$Q$$$ at the boundaries is now added.  Maybe an example will help.

Spreadsheet
Example Calculation

In the first data set on the top left an area is marked green.  The centroid of this region is to be found. A column of zeros is added to the start of the data.  This is seen in the data set on the top right.  These are added because when performing the operation in software you end up accessing array outside of their bounds.

Typically when calculating centroids each green cell would have to be accessed for each calculation. However when using the new method described above only the cells in blue are accessed.  You may also notice that the width of these doesn't change the number of operations.  The complexity of the calculation is only dependant on the height of the region.  I've included this spreadsheet for you to play around with and get a better understanding of the process.

One last this to note is that the addition of the column of zeros while maintaining the original $$$x$$$ indices is silly as it creates a column with an index of negative one.  This is where the adjusted x index comes in.  Using this and the following adjusted equations allow the calculation to be performed easily on computers.
Equations

Let's for example calculate the centroid of the green section on only row 5.  Its $$$x$$$ bounds are $$$x_1=4$$$ and $$$x_2=7$$$, but in the above equations $$$P$$$ and $$$Q$$$ are accessed via the adjusted $$$x$$$ index at 4 and 8 (8 because of +1 on the upper index).  This means the denominator is equal to $$$(356-170) = 186$$$, the y-numerator is equal to $$$5(356-170)=930$$$, and the x-numerator is equal to $$$(8 \times 356-4 \times 170)-(1560-420)=1028$$$.  This leads to centroid coordinates of (5.526, 5). This is what you would expect as it's a single row the y value is obviously 5 and the x value is a little to the right of centre due to the increasing trend in the numbers. The centroid coordinates calculated are given in the original x index to allow accessing the original data.

Saturday, March 4, 2017

Converting Binary to Gray Code with XOR

Over the last week I've been researching Gray Codes and how to generate them. I understand the basic concept but I wanted a deeper understanding.  Most of the explanations out there are a little dense so I'm going to try to make the available information a little bit more accessible to everyone. So let's start with a basic explanation of what a Gray code is.

A Gray Code is a special way of counting, where adjacent entries differ by only one digit. In practice this almost always refers to a binary counting sequence and more specifically a special code called the Binary Reflected Gray Code (BRGC). From here on, this is the Gray code that I'll be talking about.

These code have practical uses in many areas.  The one that I'm most familiar with is in rotary encoders where it's advantageous that only one digit changes at a time.  In a normal binary counting sequence you may have multiple bits change at once.  For example from $$$0111$$$ to $$$1000$$$. The problem with this is that if all the bits don't change at exactly the same moment, the value $$$1010$$$ may be read from the sensor.  If Gray codes are used, as you progress through the series only one bit changes at a time.  This leaves no ambiguity as to what the value is.

Let's start with a formal definition of a Gray code and try to explain things in easier to understand terms.

$$G_0=\epsilon \\ G_{k+1}=0G_k,1G_k^R$$

The series for zero bits is shown as $$$G_0=\epsilon$$$ and indicates an empty series.  The next equation, $$$G_{k+1}=0G_k,1G_k^R$$$ indicates that the series for a Gray code series of $$$k+1$$$ bits is created by taking the the Gray code of $$$k$$$ bits and joining a reversed copy to its end, with zeros added before entries in the first half of the series, and ones added before entries in the second half of the series.  This is how the series ensures that consecutive entries only differ by one digit.  This may be clearer in the diagram below.

Equation
Constructing $$$G_{k+1}$$$ from $$$G_{k}$$$

If $$$G_k$$$ is a Gray code where consecutive entries differ by one digit, then appending a mirror of itself to its end won't change that, except for where the two series join, as those two entries will be the same. Adding zero to the first half of the new series and one to the second half of the new series won't change it either but it does mean that the entries where the series join now differ by one digit. This means that the new series $$$G_{k+1}$$$ is now also a Gray code.  It also proves that the entries in the output are unique.  Starting with $$$G_0$$$, we know that all its entries are unique.  If we then assume that all the elements of $$$G_k$$$ are unique, then the elements of $$$G_{k+1}$$$ must also be unique as $$$G_k$$$ with zeros added before it and a reverse copy of $$$G_k$$$ with ones added before it create a new unique series. It's then provable by induction that the each entry in a Gray code is unique.

This gives us the instructions needed to construct a Gray code of any length.  Start with the minimal Gray Code that contains one empty element and keep doubling its length when adding a new digit.  

Gray Code
Constructing $$$G_3$$$ from $$$G_2$$$
Getting a picture in your mind of what a Gray code looks like will be helpful for the next step. For some situations the above representation can help, but in others, the vertical table representation with the index $$$n$$$ of each element is better.
Gray Code
Gray code indices in decimal

Now that I've explained the basics of what a Gray code is, I want to detail a proof of a well known formula to generate them with an xor operation.

$$H[n] = n \oplus (n>>1)$$

Where $$$H[n]$$$ is the $$$n$$$th element of the Gray code. $$$H$$$ is used at the moment as it's a hypothesis that this formula is equivalent to $$$G[n]$$$ described above.

I found a well thought out explanation of this proof in a Quora post but it could do with some further explanation.  So I'll repeat the proof here with some commentary to help those new to the subject. To be clear all of these operations are performed on the binary representation of the numbers involved.

First, let's define a helper function.

$$F[n,k]=(2^k-1) \oplus n$$

$$$F[n,k]$$$ is a function that inverts all the bits of a binary string $$$n$$$ of length $$$k$$$. $$$2^k$$$ is a binary string of length $$$k+1$$$ that is all zeros except for the most significant bit. Subtracting one from this will create a binary string of length $k$ that is all ones.  When the string $n$ is xored with this all of the bits are inverted. This is a function that essentially reverses the direction of  a binary count. For example, when $$$k=3$$$, and for a value of $$$n=001$$$, $$$F[001,3] = 110$$$.  Where $$$001$$$ is the second entry when counting in binary, $$$110$$$ is the second last.  It's a mathematical representation of reflecting/mirroring a binary count.

Now we have that out of the way we will prove the following

$$H[x] \oplus H[F[x,k]]=2^{k-1}$$

What does that mean?  It's saying that the output of the hypothesis function xored with the hypothesis function with the input ordering reversed will give a string of length $$$k$$$ that is all zeros except for the most significant bit.  This is a fancy way of saying that the output of $$$H$$$ is symmetrically except for the most significant bit.  This is a property of the Gray code that we are trying to prove. From the beginning of the series each element of a Gray code is equal to the element the same distance from the end except for the most significant bit.

As $$$F[x,k]$$$ is an inverted version of $$$x$$$ the most significant bit of one or the other, but not both, will be equal to one. Let's assume that it's $$$x$$$. It doesn't matter which one is chosen, the proof works either way.  We can now write

$$x=2^{k-1} \oplus y \quad where \quad y<2^{k-1}$$

It then follows that $$$F[x,k] = F[y,k-1]$$$.  This says that inverting $$$x$$$ which we have defined to start with one, will now start with zero and is equal to the inverted value of the new shorter bit string $$$y$$$.  For example if $$$x=1100101$$$ then $$$y = 100101$$$.  So $$$F[1100101,7]=0011010$$$ and $$$F[100101,6]=011010$$$.

Now onto the proof.  It's all maths at this point, you can ignore conceptual ideas here and just work through the process.

\begin{align}H[x] \oplus H[F[x,k]]&=(2^{k-1} \oplus y) \oplus (2^{k-2} \oplus (y>>1)) \oplus H[F[y,k-1]] \\ &=2^{k-1} \oplus 2^{k-2} \oplus H[y] \oplus H[F[y,k-1]] \\ &= 2^{k-1} \oplus 2^{k-2} \oplus 2^{k-2} \\ &= 2^{k-1}\end{align}

Now we've proven this, for strings of length $$$k$$$ using strings of length $$$k-1$$$ it just needs to be proved for strings of length 1 ie $$$k=1$$$ for the inductive proof to be complete.

\begin{align} H[0] \oplus H[1] &= 0 \oplus 0 \oplus 1 \oplus 0 \\&=1\end{align}

All of that above just proves that $$$H[n] = n \oplus (n>>1)$$$ results in an output series that is mirrored with the most significant bit inverted.  Now we need to prove the following

$$H[2^k \oplus x] = 2^k \oplus H[F[x,k]] \quad where \quad x<2^k$$

This states that if we add a one to the input to the hypothesis function, it's equivalent to inverting the input and then adding a one to the output.  It's hard to see, but this describes the process of going from a Gray code of $$$k-1$$$ bits to a Gray code of $$$k$$$ bits.

\begin{align} H[2^k \oplus x] &= 2^k \oplus x \oplus 2^{k-1} \oplus (x>>1) \\ &= 2^k \oplus 2^{k-1} \oplus H[x] \\ &=2^k \oplus H[F[x,k]]\end{align}

As we did before, this needs to be proved for $$$k=1$$$ for the induction proof to work.

\begin{align}H[10 \oplus 1] &= 10 \oplus H[0] \\H[11] &= 10 \oplus 0 \oplus 0  \\ 11 \oplus 01 &= 10 \\ 10=10 \end{align}

OK, that was a confusing rabbit hole we just fell down, but it proves that $$$H[n]=G[n]$$$ for all $$$n$$$.  I'm not entirely sure we need both parts of that proof, but it can't hurt to prove it twice. 😅 

We can now be sure that $$$G[n] = n \oplus (n>>1)$$$  How simple is that!  A Gray code can be generated by xoring its index with itself but right shifted by one.  This means that each digit of the Gray code is the xor of two consecutive digits in the input.

Generating a Gray code from an index - input in green output in purple
In the image above it's kind of like differentiating the bit stream.  If two consecutive bits are different it generates a one in the output.  So if we go from index to Gray code with something like a differentiation you might think we can go the other way with integration.  You sure can.  I won't prove it but you can look it up.

Now all this seems a little simple and useless apart from the use case I gave at the start, but this goes into some deep concepts.  We're looking at mathematics over a finite field, in this case $$$GF(2)$$$. Now here's the really cool part.  If you use each digit of a Gray code as a coordinate and plot it, you visit each vertex of a square for a two bit Gray code.  For a three bit Gray code you visit each vertex of a cube travelling along the edges. More generally for a k bit code you travel along the edges of a k-dimensional hypercube visiting all vertices.  This is called a Hamilton cycle.

Sorry if I'm ranting but I find all this fascinating, but it hasn't solidified in my mind yet.

If you've made it this far you deserve a joke.

Q. What's the best way to visit all the vertices of a hypercube?
A. On a Hamilton bicycle.

Hey give me a break.  I didn't say it'd be funny.

In case the equations don't render properly, here is a PDF of this page.
A PDF of the original solution on Quora.

Wednesday, February 22, 2017

Hilbert Curve Generation With Lookup Tables

I've been researching something that requires a Hilbert curve and I thought I'd share how to generate the path and also move between the index and coordinates of points.

For those of you unfamiliar with Hilbert curves let me quickly explain what one is and why you would want to use one.  The curve covers every point in a square with side length 2^n by moving up, down, left and right, starting at one corner and ending at an adjacent corner.  You could accomplish the same requirements by just scanning back and forth across the rows, but the advantage the Hilbert curve has is that nearby points on the 2D grid are generally also near each other on the curve.  If you were to scan back and forth you get a lot of points that are directly beside each other in 2D space but very far apart on the curve.

So why would you want to use one?  They're great for turning 2D areas into 1D streams of data while maintaining locality.  This comes in handy in image processing.  Depending on what you want to do, this may make the processing a lot easier.

Generation of the curve can be done recursively by first selecting an initial shape type and then using tables of sub-types and locations to generate the next stage.  In the animation below the initial shape is type 1.  When you go to the sub-type table you see that type 1 is replaced with types 2, 1, 1, and 4 in that order.  The position of each sub-type is defined by the parent type.  As the parent is type 1 the locations of the sub-types from the location table are:

2 at (0, 0)
1 at (0, 1)
1 at (1, 1)
4 at (1, 0)

The location table can then be used to find the coordinates of the sub-types.  As binary coordinates are used, the x and y coordinates of the sub-type can be appended to its coordinates to find the new sub-coordinates.   For example, the coordinates of the type 4 sub-type at position (1, 0) above are (1, 1) (0, 1) (0, 0) (1, 0).  Appending these to the coordinates of (1, 0) you find the new sub coordinates (11, 01) (10, 01) (10, 00) (11, 00).  Trust me it's confusing at first but after a while it all makes complete sense.  I'm still trying to come up with appropriate terminology for the process.  These tables are from a paper


It's a little confusing at first but I recommend reading that to help understand the process.


Recursive Hilbert Curve Generation

What if you don't want to generate the full curve and just want to know where along the curve a point lays or where a point on the curve is in 2D space?  Let's use the below example to demonstrate this.

44th Point at location (111, 101)

The above image shows the 44th point of the third stage at (111, 101).  It's actually the 45th point but we start counting at  0.  We add an index row to the tables provided in the above paper to make the process easier.

Generation Tables

We start by converting the index to a binary number with 2 bits for every stage.  In this case it's the third stage, so we get 6 bits.  These bits are placed two at a time on rows under the index column.  We also define the initial shape as type 1.  From this, the sub-type and location for the first row can be found.  The sub-type becomes the type for the next row.  This process is repeated two more times until we have three locations.  To get the final location coordinates, just append the x and y coordinates downwards to get (111, 101).

Converting an Index to a location

To go from coordinate to index, split the location up by removing the highest bit from the x and y coordinate for each line.  As the start type is known, the index can be found.  This index is then used to find the sub-type for the next row.  This is repeated two more times.  The indices are then concatenated to get 101100 in binary which is 44.

Converting a location to an index

I hope this helped to explain the process.  Trust me, even if it's still confusing, having diagrams and a couple of hopefully easy to follow examples will help.

Monday, November 21, 2016

How Can Polls Be Wrong and Right At The Same Time?

The results of the 2016 Presidential election have surprised many and caused people to call into question how accurate polls are.  There are endless think pieces about how the polls were wrong and how people were lulled into a false sense of certainty, but what seems to get overlooked in all of these conversations is the probabilistic nature of polling.  That's what I want to delve into a little deeper and understand better.

Before going too far, it should be said that conducting an unbiased poll is an incredibly hard thing to do.  You need to sample a representative cross section of the electorate, and that's becoming harder as the available ways to contact people are fragmenting due to technology.  There are many other causes for polls giving unexpected results, but what I'm really interested in is what the underlying mathematics tells us about the results of a poll when all these other effects are stripped away.

From this point on we're going to look at a simplistic poll where we ask people if they prefer Apples or Bananas, and we'll consider a response for an Apple as positive response.  We also need to ask the question "what are we trying to discover with this poll?".  The answer to that is we are trying to determine the proportion of the public that prefer Apples to Bananas.  We don't know what that number is.  It's a hidden parameter of the population that we are trying to discover through random sampling.

Let's define some variables, do some math and take a poll.  It won't be anything too formal.  I'm not labelling my graphs, this is just a casual chat between friends.

p is the probability of a positive response for an Apple
n is the number of people surveyed for the poll
k is the number of positive responses for an apple
P(k;p) is the probability of k responses given the probability of a single response is p
Binomial pmf
Binomial pmf range
The probability mass function for a series of Bernoulli trials is the the binomial function.  In this case it is telling us the probability of a certain result given that we already know the underlying probability parameter, p.  However, we don't know p, so it's more useful to look at the same equation but this time with the results of a series of trials.  Given that we know the number of positive results, what is the likelihood of the parameter p being equal to a value?  The new equation L(p;k) is referred to as the likelihood function, and it doesn't follow the same rules as a probability density function, but it does allow comparison of relative likelihoods.
Likelihood Function
For example: 100 people are polled and 35 people say that they prefer apples to bananas.  This gives a likelihood function of:
A poll of 100 people with 35 positive responses
When the above likelihood function is graphed we get the following.

A poll of 100 people with 35 positive responses
From this we can see that the relative likelihood that p is 0.35 compared to 0.3 is around twice as likely but the possibility isn't remote.  We can visually see that p could be anywhere between about 0.27 and 0.43 and you wouldn't be surprised.  What if we conduct more polling?  Let's say we poll 300 people.

A poll of 300 people with 105 positive responses
You can now see that the peak of likely values has now tightened, and we can more confidently say that the value of p likely lies between 0.3 and 0.4.  This shows how additional polling allows the margin of error to be reduced.  If you're astute you may have noticed that the results of the polls don't rely on the size of the population being polled.  As long as you select a representative cross section of society the size of the population is irrelevant.  You may also be interested in what happens if a relatively small number of people are polled.  This can be seen when you poll 20 people.

A poll of 20 people with 7 positive responses
It wouldn't really surprise you if p was between 0.18 and 0.56 a pretty wide range. and doesn't reveal too much about the actual value of p.  It's also interesting to note that all of the example above have 35% positive response rates and therefore p peaks at 0.35 for all surveys.  Polling more people only  increases confidence in our results.  In reality the peak would move around a little bit with each new bit of data.  However, the more people that you poll, the more confident you would be that about one in three people prefer apples to bananas.

There is a rule of thumb for margin of error in polling that states that 100 divided by the square root of the number of respondents is equal to the percentage margin of error.  So for our 3 examples above the margin of error for 20, 100, and 300 people respectively are 22%, 10% and 6%.  It's not uncommon to see a political poll of 2000 respondents, and using our rule of thumb above, it would still have a margin of error of 2.2%.

I suppose what I'm trying to get at is that even if you could remove all of the external human factors affecting the poll and you chose a perfectly random cross section of the population, your sample size will limit how accurately your poll determines the underlying parameters the population.

For all the talk about how every polling model got the result of the Presidential election wrong, remember that the 538 model gave this outcome a one in three chance of occurring.  That isn't an insignificant probability.  It's true that polling methods need improvement, but in this case I think the main problem is the way probabilities are communicated to the general public are too confusing.  People don't intuitively understand what percentage probabilities mean.  The IPCC have this problem and even have defined terms in everyday language to communicate probabilities to the public.

IPCC likelihood scale
I actually have some modelling problems in mind that I would like to work on one day, and until now I had always envisaged communicating them as percentages, but given the intended audience I think I might change my approach.