analytics

Monday, December 29, 2014

Shortnet definition and source code

I'm defining the Shortnet of a group of points to be the subset S of the collection of every possible line connecting pairs of points C, such that every line crossing a line shorter than itself has been removed.
Example: You have your points...
 Then you create all the lines connecting every point pairwise...
Then you remove any line that crosses a line shorter than itself...
Here is one with 50 points:

There will never be any lines that cross left as one or the other of each pair of lines that cross is removed.
Source Code using PIL drawing library and Python:
import random
import Image, ImageDraw
def cross(a1, b1, c1, d1, a2, b2, c2, d2):
    if (-a2*b1+a2*d1+c2*b1-c2*d1+b2*a1-b2*c1-d2*a1+d2*c1)== 0:
        return False
    if (-a2*b1+a2*d1+c2*b1-c2*d1+b2*a1-b2*c1-d2*a1+d2*c1) == 0:
        return False

    s = 1.0*(a1*d1-d2*a1+d2*c1-c1*b1+c2*b1-c2*d1)/(-a2*b1+a2*d1+c2*b1-c2*d1+b2*a1-b2*c1-d2*a1+d2*c1)
    t = 1.0*(-b2*c1+a2*d1+b2*c2-a2*d2+d2*c1-c2*d1)/(-a2*b1+a2*d1+c2*b1-c2*d1+b2*a1-b2*c1-d2*a1+d2*c1)
    if s > 0 and s < 1.0:
        if t > 0 and t < 1.0:
            return True
    return False
def gen(num):
    points = []
    for i in range(0, num):
        point = [random.randint(100, 700), random.randint(100, 700)]
        points.append(point)
    return points
def shorter(a1, b1, c1, d1, a2, b2, c2, d2):
    d1 = ((c1 - a1)**2.0 + (d1 - b1)**2.0)**.5
    d2 = ((c2 - a2)**2.0 + (d2 - b2)**2.0)**.5
    if d1 <= d2:
        return True
    else:
        return False
def makelines(points):
    lines = []
    for i in range(0, len(points)-1):
        for j in range(i+1, len(points)):
            lines.append([points[i][0], points[i][1], points[j][0], points[j][1]])
    return lines
def picklines(lines):
    marked = []
    for i in range(0, len(lines)):
        marked.append(True)
    for i in range(0, len(lines)-1):
        for j in range(i+1, len(lines)):
            if cross(lines[i][0], lines[i][1], lines[i][2], lines[i][3], lines[j][0], lines[j][1], lines[j][2], lines[j][3]):
                if shorter(lines[i][0], lines[i][1], lines[i][2], lines[i][3], lines[j][0], lines[j][1], lines[j][2], lines[j][3]):
                    marked[j] = False
                else:
                    marked[i] = False
    for i in range(0, len(lines)):
        if marked[i] == False:
            lines[i] = False
    return lines
def drawlines(d, lines):
    for line in lines:
        if line != False:
            d.line(line)
def drawpoints(d, points):
    for point in points:
        d.ellipse([point[0]-10, point[1]-10, point[0]+10, point[1]+10], fill=(255,0,0))
def main():
    points = gen(10)
    plot = Image.new("RGB", [800, 800])
    d = ImageDraw.Draw(plot)
    drawpoints(d, points)
    plot.save("dots.png")
    lines = makelines(points)
    drawlines(d, lines)
    plot.save("all.png")
    plot = Image.new("RGB", [800, 800])
    d = ImageDraw.Draw(plot)
    drawpoints(d, points)
    lines = picklines(lines)
    drawlines(d, lines)
    plot.save("shortnet.png")
main()

Friday, December 19, 2014

Parallelizable convex hull algorithm

Suppose you have a number of points in any number of dimensions, here I will just show 2.
First find the centroid, and the point farthest from that centroid Pc and construct the circle or n-dimensional sphere to it from the centroid... (Centroid shown in blue)
Now consider all the straight lines C to P(i) to B(i) from the Centroid to each Point and continuing on to a point on the Bounding region. Make a list of Pc and all points P(i) such that P(i) to B(i) is shorter than any other line from a different point P(j) to B(i).
This last step is the part that can be done in parallel, passing the entire set to as many processors as there are points and each one finding the solution for its given point. The points that are closest to their respective bounding point than any other point is to that bounding point are the hull. 

Tuesday, December 9, 2014

recursive formula for points on the unit circle of regular 2^n gon's

I found that formulas of this type:


describe the points on a circle corresponding to regular 2^(n+1)-gons where n is the number of square root symbols in the formulas and a,b,c... are each 1 or 2. For example a 16-gon
I haven't quite figured out where each point will be on the unit circle by the signs of the terms though. 

Monday, December 8, 2014

Fast parallelizable way to generate approximate points for a circle without using trig functions

This method starts with the four points of the bounding square...
Then replaces each point with two points (1- root 2/2) away from the sides... This is basically like truncating every corner to the right proportion...
Then iterates this process but using .25 instead of (1-root 2/2) that I found by experiment works very well. I suppose for perfect accuracy there is a recursive series involving root 2 but it approaches .25 so quickly that I didn't bother with it.


I found 4 iterations is enough for fairly large circles...
And in the end you have an list with all the points in counter clockwise order around the circle. 
Python source code (Using PIL for drawing)
import Image, ImageDraw
def plot(points, draw):
    for i in range(0, len(points)-1):
        draw.line([points[i], points[i+1]])
    draw.line([points[len(points)-1], points[0]])
def truncate(points):
    newpoints = []
    param = .29
    if len(points)>4:
        param = .25
    for i in range(0, len(points)-1):
        x1 = points[i+1][0]*(param)+points[i][0]*(1.0-param)
        y1 = points[i+1][1]*(param)+points[i][1]*(1.0-param)
        x2 = points[i+1][0]*(1.0-param)+points[i][0]*(param)
        y2 = points[i+1][1]*(1.0-param)+points[i][1]*(param)
        newpoints.append((x1, y1))
        newpoints.append((x2, y2))
    x1 = points[0][0]*(param)+points[len(points)-1][0]*(1.0-param)
    y1 = points[0][1]*(param)+points[len(points)-1][1]*(1.0-param)
    x2 = points[0][0]*(1.0-param)+points[len(points)-1][0]*(param)
    y2 = points[0][1]*(1.0-param)+points[len(points)-1][1]*(param)
    newpoints.append((x1, y1))
    newpoints.append((x2, y2))
    return newpoints
def main():

    points = [(100,100),(100,700),(700,700),(700,100)]
    
    for i in range(1, 5):
        p = Image.new("RGB", [800,800])
        draw = ImageDraw.Draw(p)
        points = truncate(points)  
        plot(points, draw)
        p.save("plots"+str(i)+".png")
main()

Thursday, December 4, 2014

Using physics to try and find a Hamilton Circuit of a graph with odd number of vertices

Supposing you have a simply connected planar graph of even number of vertices with no loops edges labeled like so:
I found one can set one edge leading to one vertex, say A arbitrarily to 1, then write equations from Kirtchoff's laws for the graph...
http://en.wikipedia.org/wiki/Kirchhoff%27s_circuit_laws

Kirchoff's laws say that the sum of the edges connected to each vertex is 0 and the sum around each face is 0.

It's kind of interesting that after you set one edge to 1, this will always be the right number of equations because of the graph equation V+F=E+2, and we have V+F equations to solve for E unknowns, The plus two in this formula is canceled by the fact that usually the formula counts the region outside of the graph as a face, but we don't use write an equation for that, and we've set one edge value to 1 so there is one less unknown edge.

Note there will be no solution if c,g,l, or m in this example is the edge chosen to be 1. But for any the other edges it will either be this solution or -1*all the values.

Then we can plot these solutions...
Now we can follow the flow of the graph around through every vertex exactly once by starting with A which is 1, then following the path from that vertex that is -1, and from there to the vertex following the path of 1, and so on alternating values. Since there are an odd number of vertices this works out evenly.
If you had an odd number of vertices to start out with, you could break the edge that you set as 1 into two edges one -1 and 1...

**Reasoning**
I don't know how to prove it, it just seems to work in every case. It makes sense to me electrically that if there is a nonzero flow around each face's edges, that some faces will combine and form one flow around the outside. I don't know how to prove that that will always go through every vertex though.

Friday, November 28, 2014

Branching poisson disk algorithm

My program outputs a distribution like this:
Every dot is at least a distance of 10 pixels from any other.
The way it works is to start at the middle, and branch out in 3 evenly spaced directions at a random angle from the center dot, then from those 3 points branch in a similar way generating more points and continue in that fashion for a certain number of iterations, above is 25 generations. I think this picture helps explain it a lot.
The speed is fairly quick, generating these graphs in python took only about 1.5 seconds for around 1300 points, I think it could be speeded up considerally using parallel processing for example on a gpu with a faster language. 
**Python source code**
import math
import random
import Image, ImageDraw
def check(d, pointlist, point, gridsize):
    tooclose = False
    if point[0] < 10 or point[0] > gridsize[0]-10 or point[1] < 10 or point[1] > gridsize[1]-10:
        return True
    for p in range(0, len(pointlist)):
        if ((point[0]-pointlist[p][0])**2.0 + (point[1]-pointlist[p][1])**2.0)**.5 < d-1:
            return True
    return tooclose
def findspot(pointlist, added, d, gridsize, draw):
    newlist = []
    newadded =[]
    for point in pointlist:
        newlist.append(point)
    for point in added:
        tx = point[0]
        ty = point[1]
        theta = random.random()*2*3.14159
        for i in range(0, 3):
            px = tx + d*math.cos(theta+i*(2*3.14159)/3)
            py = ty + d*math.sin(theta+i*(2*3.14159)/3)
            if check(d, newlist, [px,py], gridsize) == False:
                newlist.append([px,py])
                newadded.append([px,py])
    return newlist, newadded
class graph():
    def __init__(this, d, gridsize, draw):
        this.pointlist = [[400,400]]
        this.added = [[400,400]]
        this.d = d
        this.gridsize = gridsize
        this.draw = draw
    def tree(this, limit):
        newlist, newadded = findspot(this.pointlist,this.added,this.d, this.gridsize, this.draw)
        if limit < 25:
            this.pointlist = newlist
            this.added = newadded
            this.tree(limit+1)
def main():
    gridsize = [790,790]
    d = 10.0
    point = [400,400]
    plot = Image.new("RGB", [800,800])
    draw = ImageDraw.Draw(plot)
    G = graph(d, gridsize, draw)
    G.tree(0)
    for point in G.pointlist:
        for i in range(-1, 2):
            for j in range(-1, 2):
                plot.putpixel((int(point[0]+i), int(point[1])+j), (255,255,255))
    plot.save("plotdot.png")
main()

Sunday, November 23, 2014

recursively moving around the unit circle with constant angle to the origin

Starting with solving this general formula for two vectors [x,y] and [a,b]:
The above says the dot product between them equals to c and the magnitude of the cross product equals to d (counterclockwise). And then it is solved for a,b.
So we begin by thinking about if [x,y] starts as the vector [1,0]... then we can find possible values for c and d like so, picking any value as close or far away from 1 as we'd like for x and solving...
So plug those in for d and c respectively...Or alternatively use a rational value from the rational parameterization of the circle. These are the values for the dot and cross product from the vector to the generated point to the vector [1, 0] because the dot product will just be the x coordinate of the point we generated and the cross product will just be the y.   Let's call this Formula 1.
Now c and d is also our first vector after [1,0] so we can plug that in for x and y as well and get...
This vector [a,b] is the same angle to [.9801,.141] as [.9801,.141] was to [1,0], because we made the dot and cross products stay the same.

Now we can plug this [a,b] in for x and y in formula 1 and get...
And in this way work our way around the unit circle by that small angle (which we actually don't know) without ever having to figure a cosine or sine! An approximation of the angle can be found by seeing how long it takes for this process to get close to a known value like [0,1]

Friday, November 21, 2014

wave interpolation (with natural number harmonics of different amplitudes)

I noticed that it's possible to find solutiond for waves that are 0 at the endpoints of the interval and that pass through a certain number of points.
Like for instance, say you want to know the waves that pass through these points:
[[.5, 1], [1, 0], [1.5, -1]]
Since that is three points, I start by considering waves of the form:
Maple can solve this by entering...
You can see that there are many waves satisfying these conditions, b and c can have any possible value and a just has to be 1+ whatever c is chosen to be... for example a=5.5, b=3, c=4.5
The thing I noticed is that there are two degrees of freedom in the solutions for this wave because the t values interpolated through .5, 1, 1.5 are all multiples of the first value of .5. We can remove a degree of freedom in the solutions by interpolating through .5, 1.2, 1.5, here see that 1.2 is not a multiple of .5 but 1.5 still is...
The solutions only have one degree of freedom now so that c can be chosen as anything. Interestingly the solutions involve the golden ratio! 1.618033... is usually known as phi or the golden ratio. I'm not sure exactly why that is.

So the next step of course is to look to see that a case where the t values are all not multiples of each other yields only one solution...Let's try [.7, 1.2, 1.6]
There is only one solution for those heights through those t values.

** Finding unique solutions for an unknown wave though a number of points**
The trick for this is as you can see above to sample your unknown wave so that none of the t values are multiples of each other... There might be many ways to do this. It's kind of an interesting problem in itself.
One way is to use prime numbers...if you have n points you would use fractions of the (n)th prime like for example if you have six points you could use for t values:
[1/13, 2/13, 3/13, 5/13, 7/13, 11/13]
For example, if our 6 heights at those t values are [1, 1.5, .5, 0, -.65, -.8]


Sunday, November 9, 2014

fanless turbo compressor

An idea I had for an alternative to a turbo compressor with fan blades like they usually have. I thought because the center of gravity is so close to the spindle on this design you could spin it easily really quickly and get better performance. 

yellow on an RGB computer screen vs from the spectrum



The bottom graph shows an overlay of the yellow you get from mixing pure red and green light and what you would get from a light source of the right frequency for yellow. As you can probably see it's off not just by a scaling factor vertically, but the second half of the wave is closer than the first half.  This phenomenon would persist with different colors even if a yellow led was added and after that an orange, etc.Though it would improve with each addition. Maybe one day they'll make an led that can actually be any different color from the spectrum by altering the voltage source to it, so they don't have to add different colors of light to try and make all the colors.  


Symbol gas cryptography

First two computers are initialized to the same configuration of the symbol gas...
Every symbol is given an initial position and an initial motion vector, above I've shown the motion vector for just the k.
The algorithm for running the simulation is made so that no matter how long it is run of simulated time, the final position of the symbols is deterministic based on the starting values. Though it is very psuedo random where every symbol will be after a fair amount of time, because they are all bouncing off of each other and off the walls.  One computer has a message it wants to send to the other computer such as...
"this is a secret message..."
It first sends the position of the symbol "t" rounded to some accuracy, then the simulation is run for some predetermined amount of time such as 1000 seconds of simulation time, then it sends the new position of the symbol "h"  and repeats for the whole message.
To decode the other computer looks for which symbol is at the first position sent from the other computer, then runs the simulation for the predetermined amount of time then which symbol is at the second position sent, etc..
I think because of the simulated entropy that a third party listening in on the message, but not knowing what the initial parameters of the simulation were, will never be able to decode a message. Even if they have an encoded message and a plaintext decoding of the same message, because only positions are given and not vectors or the predetermined time scale it would be impossible to generalize to decode other messages.
I think it's better than other methods because all have a sort of formula using a  key that might be very complex but then there's always the possibility of finding the inverse. Here there is no formula for where one of the symbols will be after a certain amount of time and the simulation simply has to be run to find where a symbol will be so there is no inverse either.


Thursday, November 6, 2014

4r's genetic problem solving

The problem I developed this technique for is to fill a square with 10 points so that the minimum distance between any selection of two points is maximized with wrap around. It produced the following solution after 10000 generations.
The 4 r's pneumonic is for:
1. Randomly produce a population of solutions
2. Reduce the population to the best performing subset of the population.
3. Reintroduce some new random solutions
4. Reproduce the solutions with one another to make offspring solutions.

For the example above the first step was to produce 120 random solutions each of 10 random points in the square. Then the top 10 best solutions were kept and the rest discarded. Then 10 more random solutions are added. This step turned out to be very important as without it the population stagnates and doesn't continue improving. Then the reproduction step was to take every pairwise combination of solutions, and take the average of each of the 10 points location in each to produce an offspring and add that to the population. Then steps 2 through 4 are repeated indefinitely.
**Note
It actually converges quite quickly to nearly the best solution it will be able to find, 100 generations takes only a couple of seconds while 10000 takes around 10 minutes, but the difference between the top solution they found was very small. The great improvement is in only the first handful of generations and then gradual improvement after that.Also as for the performance, after 100 generations a total of around 1000 random solutions are added in, a different approach that simply took the best of 1000 random solutions was no where near as good, so it is in the process not simply the number of solutions considered that the improvements are found. It even did better than looking at 100000 random solutions which took considerably longer to generate.

**Python source code
 import random
def populate(n):
    p = []
    for i in range(0, n):
        p.append((random.random(), random.random()))
    return p
def distance(p1, p2):
    return ((p1[0]-p2[0])**2.0 + (p1[1]-p2[1])**2.0)**.5
def fitness(p):
    mindistance = 2.0
    for i in range(0, len(p)-1):
        for j in range(i+1, len(p)):
            m = distance(p[i], p[j])
            l = distance(p[i], (p[j][0]-1, p[j][1]))
            r = distance(p[i], (p[j][0]+1, p[j][1]))
            u = distance(p[i], (p[j][0], p[j][1]+1))
            d = distance(p[i], (p[j][0], p[j][1]-1))
            if m < mindistance:
                mindistance = m
            if l < mindistance:
                mindistance = l
            if r < mindistance:
                mindistance = r
            if u < mindistance:
                mindistance = u
            if d < mindistance:
                mindistance = d
    return mindistance
def breed(pops):
    pops2 = []
    
    for p in pops:
        pops2.append(p)
    l = len(pops2)
    for i in range(0, l):
        p = populate(len(pops[0][1]))
        pops2.append((fitness(p), p))
    for i in range(0, len(pops)-1):
        for j in range(i, len(pops)):
            p = []
            for k in range(0, len(pops[1][1])):
                p.append(((pops[i][1][k][0]+pops[j][1][k][0])/2.0,(pops[i][1][k][1]+pops[j][1][k][1])/2.0)) 
            pops2.append((fitness(p), p))
    
    return pops2
def main():
    pops = []
    n = 10
    for i in range(0, n*n+1):
        p = populate(n)
        pops.append((fitness(p), p))
    pops = sorted(pops)
    for g in range(0, 10000):
        pops2 = []
        for i in range(len(pops)-10, len(pops)):
            pops2.append(pops[i])
        pops2 = breed(pops2)
        pops2 = sorted(pops2)
        pops = []
        for i in range(len(pops)-10, len(pops)):
            pops.append(pops2[i])
        pops = breed(pops)
        pops = sorted(pops)
    print pops2[len(pops2)-1]
main()
  

Wednesday, October 22, 2014

Regular Polyhedra formula

I found this formula that you put in S for how many sided polygons and P for how many meet at a vertex and gives the number of Edges, Faces, and Vertices it must have...




**Proof
There are 3 formulas describing such polyhedra...
Euler's formula for polyhedra:
F = E-V+2
Another one saying that the number of Faces times the number of sides on each face equals twice the number of edges, because every edge connects two faces.
F*S = 2*E
And one more saying that the number of vertices times the number of edges connecting to it equals twice the number of edges, because every edge connects two vertices.
V*P = 2*E

These solve in terms of S and P as the formula for E, F, and V given.

For example if we want triangles so S=3 meeting 5 at each vertex:

It says there must be 30 edges, 20 faces and 12 vertices, which is the regular icosahedron. 
I find it interesting that the assumptions are quite general, Euler's formula for polyhedra can be proven even over topologies with holes in them like a donut, and the other two would certainly still be true. But it says you could never have say a 1000 vertex donut polyhedra where the whole surface is divided into triangles and 5 meet at every vertex, because the only solution is for there to be 12 vertices for that configuration of S and P.

**Note**
I was thinking about how a toroid can be broken into almost any number of 4 sided figures with 4 edges meeting at each vertex.. there is a caveat that these formulas end up dividing by 0 for S=4 and P=4...



Tuesday, October 21, 2014

easy way to find the volume of a polyhedron

An easy formula for finding the volume of a polyhedron when you know the coordinates of the vertices...


where the bars are absolute value.



The sums in this formula are over all of the faces, Maxheight is the maximum height for any of the vertices, MinHeight the minimum value of the height of all vertices. F(i)Aveheight is the average of the height coordinates of the ith face and F(i)Area is the area of a face polygon after all the height coordinates are replaced with 0.

To see why it works first think of the 2 dimensional version of this formula:

We start finding the area B underneath each face down to the minimum value the height coordinate takes on. The red area would only be counted once as it is below the faces along the top and above any of the faces along the bottom, but we count the area underneath the faces along the bottom as well so the purple area ends up being counted twice. The area of each face is calculated as the average value of the height of the face, times the length when the height coordinates are replaced by 0 for the two vertices in the 2d case this is just the width.  

Now we do almost the same thing but find the area T above the faces of the polygons up to the maximum value of the height coordinate (Note that this area will be negative). 
Now we can solve that the total area of the polygon will be abs(T-B)/2 because it is the symmetric difference of the two overlapping regions.  It's hard to do a direct proof because of the absolute values but one way to reason about it is if you have a cube so that T and B are each the whole volume of the cube, but B negative, then it would be abs(T-B)/2, but that always has to be the solution because it is always the same system of linear equations relating T,B, the red area and the total area. 

So this same sort of thing is presented above for 3 dimensions, but it is average values of faces times the area of face polygons when the heights are replaced by 0 to make shafts of volume, but otherwise the same concept.  

  

Wednesday, October 15, 2014

Finding a path through every point of a graph

For a simple connected graph with coordinates I believe this algorithm finds a path visiting every point once if it is possible.
The first step is to go to a point on the graph X, and find two points out of those that X is connected to, P and N, that maximizes this formula:
This will ensure that the vectors from P to X and from X to N are going around clockwise, and that the angle between them taken counterclockwise is maximized.
For the example above, this will make W=P and N=Y.
So we do this for every point on the graph:
And then eliminate any edges such that every point after removal is still connected to two other points. If this isn't possible than, you either already have the solution or it isn't possible to visit every point once and return to the original point.
**Idea for proof**
The formula being maximized is individual terms from the shoelace formula for polygons, so when each term is maximized the formula as a whole is maximized, and that must mean you have a closed polygon ordered clockwise which will visit every point once and return to the start.
 http://en.wikipedia.org/wiki/Shoelace_formula


Tuesday, October 7, 2014

Generating pythagorean triples from the positive rationals [0..1]

consider this formula:
For rationals [0..1] these (x,y) pairs are on the unit circle in the first quadrant, as proven here: http://mathnow.wordpress.com/2009/11/06/a-rational-parameterization-of-the-unit-circle/
Another way to write the above formula since m is a rational between [0..1] is to substitute in p/q as that rational to get...

For any rational number [0..1] a pythagorean triple can be generated as follows...
Put the rational number in for m (or as p/q in the two other formulas), here I'm using as an example 7/12...
A, B of the triple appear in the numerators of the simplified fractions and C is the denominator common to both...


*Simple direct proof

The above shows that x,y lie on the unit circle as the sum of the squares equals 1. In the first link it shows that p/q [0..1] covers the top right quarter of the circle.
And the above shows that the numerator of x squared plus the numerator of y squared equals the common denominator squared.
This proof really admits a great deal more generality than just for pythagorean triples, considering it doesn't even require p,q be rational numbers. One could have p=pi and q = -sqr(2)/2.

*Alternate more explanatory Proof
Consider some right triangle with whole number hypotenuse and legs such that it is contained in the first quadrant. This has an a,b,c in the usual way for right triangles. We can divide a and b by the distance to that point (a,b) to normalize the vector to that point to somewhere on the circle. Because it is a pythagorean triple triangle, we will have to divide by c to normalize it (because (a^2+b^2)^1/2=c. The formula I gave after substituting in p/q  once simplified can be thought of as such a  normalization because it is two rational numbers with the same denominator and is on the unit circle, hence we clear the common denominator from x,y to get a,b and that denominator becomes c. 
* Note
I just find that interesting that it puts an ordering on primitive pythagorean triples.
I guess it's also interesting that because rational points are dense on the unit circle, if you draw a line through somewhere on the circle from the origin you will pass arbitrarily close to an infnite number of pythagorean triple points, I don't know enough topology to know what to conclude from that.