Go deh!

Mainly Tech projects on Python and Electronic Design Automation.

Sunday, July 16, 2017

Python investigation of the Shoelace Formula for polygonal area

Shoelace formula for polygonal area

Sawtooth generator

In [33]:
def sawtooth_wave(teeth, base=2, height=3):
    xy = [(0, 0)] # start
    for n in range(teeth):
        xy += [(base * (n + 1), height), (base * (n + 1), 0)]
    xy.append((0, 0))  # return to start
    return zip(*xy)   # x then y
In [34]:
from pprint import pprint as pp
from bokeh.plotting import figure, output_notebook, show
from bokeh.models import ColumnDataSource
from bokeh.charts import Scatter
output_notebook()
 BokehJS 0.12.5 successfully loaded.
In [35]:
teeth = 3
sawx, sawy = sawtooth_wave(teeth)
sawx, sawy = zip(*saw)
pl = figure(plot_width=400, plot_height=400, title=f'Sawtooth polygon with {teeth} right-angled triangle teeth')
#help(pl.line)
pl.line(sawx, sawy, line_width = 1, line_color="navy")
show(pl)

Calculate area by shoelace formula

In [36]:
def area_by_shoelace(x, y):
    "Assumes x,y points go around the polygon in one direction"
    return abs( sum(i * j for i, j in zip(x, y[1:])) + x[-1] * y[0]
               -sum(i * j for i, j in zip(x[1:], y)) - x[0] * y[-1]) / 2 
In [37]:
print(" Area is:", area_by_shoelace(sawx, sawy))
 Area is: 9.0
We know that the area of a right angled tringle is half base times height.
So each tooth of the saw has area 1/2 * 2 * 3 = 3 units and we have three teeth for a total area of 9.

Why Shoelace?

It's called the shoelace formula because of a visualization of the calculation.

Get the coordinates, in order.

list the x,y coordinates of each point of the polygon, going around the polygon in one direction, and going back to the starting point. (For more complex polygons, no line segments should cross).
Split the ordered points into two identically ordered lists of the x coordinates and the y coordinates.
Visualize the x coords as numbers beside one side of a row of eyelets of a shoe, and the y coords adjacent to the eyelets on the other side.

First "lacing":


lacing1 = sum(x[0]*y[1] + ... x[n]*y[n+1]) + x[N]*y[0]

Second "lacing"

lacing2 = sum(x[1]*y[0] + ... x[n+1]*y[n) + x[0]*y[N]

Complete formula

area = abs(lacing1 - lacing2) / 2
All tied up!

Wednesday, September 21, 2016

Non-transitive dice in Python



All dies are fair and six-sided but with differing numbers on each face.
One die is said to be stronger, > than another if that die has a higher probability of rolling a higher number than the other and so winning.
In [29]:
from itertools import product
from collections import namedtuple
In [35]:
Die = namedtuple('Die', 'name, faces')
dice = [Die('A', [3,3,3,3,3,3]),
        Die('B', [4,4,4,4,0,0]),
        Die('C', [5,5,5,1,1,1]),
        Die('D', [2,2,2,2,6,6])]
In [44]:
def cmp_die(die1, die2):
    'compares two die returning <, > or ='
    # Numbers of times one die wins against the other for all combinations
    win1 = sum(d1 > d2 for d1, d2 in product(die1, die2))
    win2 = sum(d2 > d1 for d1, d2 in product(die1, die2))
    return '>' if win1 > win2 else ('<' if win1 < win2 else '=')
In [45]:
dice
Out[45]:
[Die(name='A', faces=[3, 3, 3, 3, 3, 3]),
 Die(name='B', faces=[4, 4, 4, 4, 0, 0]),
 Die(name='C', faces=[5, 5, 5, 1, 1, 1]),
 Die(name='D', faces=[2, 2, 2, 2, 6, 6])]
In [46]:
cmp_die(dice[0].faces, dice[1].faces)
Out[46]:
'<'
In [47]:
shortform = []
for (n1, f1), (n2, f2) in zip(dice, dice[1:] + dice[:1]):
    comparison = cmp_die(f1, f2)
    print('%s: %r %s %s: %r' % (n1, f1, comparison, n2, f2))
    shortform.append('%s %s %s' % (n1, comparison, n2))
print('\nIn short:', ', '.join(shortform))
    
A: [3, 3, 3, 3, 3, 3] < B: [4, 4, 4, 4, 0, 0]
B: [4, 4, 4, 4, 0, 0] < C: [5, 5, 5, 1, 1, 1]
C: [5, 5, 5, 1, 1, 1] < D: [2, 2, 2, 2, 6, 6]
D: [2, 2, 2, 2, 6, 6] < A: [3, 3, 3, 3, 3, 3]

In short: A < B, B < C, C < D, D < A
Notice that it forms a circle of winning die!

Monday, March 14, 2016

Mathematicians Discover Prime Conspiracy: Python check


In the article Mathematicians Discover Prime Conspiracy essentially is described the discovery that the last digit of the next prime number is less likely to equal the last digit of its predecessor.
To check this I decided to download a list of the first million prime numbers from here and check.
In [27]:
from collections import defaultdict

this2previous = defaultdict(int)
previous = None
with open('primes1.txt') as f:
    # Skip text file header
    print(f.readline())
    for line in f:
        line = line.strip()
        if not line:
            continue
        for primelastdigit in (prime[-1] for prime in line.split()):
            if previous is not None:
                this2previous[(primelastdigit, previous)] += 1
            previous = primelastdigit
                 The First 1,000,000 Primes (from primes.utm.edu)

In [28]:
def pretty_print(counts):
    lastdigit = None
    for (this, previous), count in sorted(counts.items()):
        if this in '25' or previous in '25':
            continue
        if lastdigit != this:
            print()
        lastdigit = this
        print('%2s %s after %s frequency: %6i' 
              % (('#' if this == previous else ''), this, previous, count))
In [29]:
pretty_print(this2previous)
 # 1 after 1 frequency:  42853
   1 after 3 frequency:  58255
   1 after 7 frequency:  64230
   1 after 9 frequency:  84596

   3 after 1 frequency:  77475
 # 3 after 3 frequency:  39668
   3 after 7 frequency:  68595
   3 after 9 frequency:  64371

   7 after 1 frequency:  79453
   7 after 3 frequency:  72827
 # 7 after 7 frequency:  39603
   7 after 9 frequency:  58130

   9 after 1 frequency:  50153
   9 after 3 frequency:  79358
   9 after 7 frequency:  77586
 # 9 after 9 frequency:  42843
As you can see, the frequency of the last digit being followed by the same last digit in consecutive prime numbers (for this sample of the first million primes), is markedly lower.

Followers

Subscribe Now: google

Add to Google Reader or Homepage

Go deh too!

whos.amung.us

Blog Archive