gpu grid computing

#### GPU Grid Computing for Gene Sequencing

It's always stimulating to have conversations with people outside your field of expertise and recently during a very educational conversation with some friends in biology, I noticed some of the bottlenecks they face might be solved in an interesting way. So as a side project, I'm creating an open source application that uses spare GPU cycles for grid computing. There are obviously a lot of applications for this but given the magnificent work their labs do, it is most practical to start with a gene sequencing package. You can see the stub/concept landing page here (gpyou.org).4D rotation

#### Axis vs Plane of rotation

The idea of an axis of rotation maybe intuitive in 3D but it is only a coincidence that axis of rotation just happens to work. And indeed it works for no other number of dimensions. So we need to think of rotation in another way. As it turns out we can; planes of rotation. In 3D the are three planes of rotation, the XY, XZ, YZ planes. If we call the four coordinates of 4D w,x,y,z then in 4D there are three additional planes of rotation: the WX, WY, WZ planes.#### Matrices

Since there are rotation matrix for point in 3D we simply need to find rotation matrix for those other planes. And here they are (I present the derivation elsewhere)
XY rotation matrix

cos(θ) | sin(θ) | 0 | 0 |

-sin(θ) | cos(θ) | 0 | 0 |

0 | 0 | 1 | 0 |

0 | 0 | 0 | 1 |

XZ rotation matrix

cos(θ) | 0 | -sin(θ) | 0 |

0 | 1 | 0 | 0 |

sin(θ) | 0 | cos(θ) | 0 |

0 | 0 | 0 | 1 |

YZ rotation matrix

1 | 0 | 0 | 0 |

0 | cos(θ) | sin(θ) | 0 |

0 | -sin(θ) | cos(θ) | 0 |

0 | 0 | 0 | 1 |

#### 4D rotation Matrices

XW rotation matrix

cos(θ) | 0 | 0 | sin(θ) |

0 | 1 | 0 | 0 |

0 | 0 | 1 | 0 |

-sin(θ) | 0 | 0 | cos(θ) |

YW rotation matrix

1 | 0 | 0 | 0 |

0 | cos(θ) | 0 | -sin(θ) |

0 | 0 | 1 | 0 |

0 | sin(θ) | 0 | cos(θ) |

ZW rotation matrix

1 | 0 | 0 | 0 |

0 | 1 | 0 | 0 |

0 | 0 | cos(θ) | -sin(θ) |

0 | 0 | sin(θ) | cos(θ) |

It's not 13 Billion light years!

Effectively every popular science outlet, when describing cosmology, states or implies five things that make me crazy:

- The universe started from a single point.
- The universe is 13 billion years old.
- The light from distant objects has taken 13 billion years to get to us.
- Those objects are moving away from us at the speed of light.
- The universe is 13 billion light years in radius.

#### The answer

The answer involves a concept not discussed by science editors: The Hubble flow. The Hubble flow is essentially the constant creation of space between objects. It started with the Big Bang and it is why objects seem to be moving away from each other. So, how far are distant objects? Well distant objects appear as far away as the time it takes light to get to us from them. But if you ponder this definition closely, you'll see it is deeply counterintuitive. Normally when we see an object a hundred meters away, it is in fact 100m away. This is because we measure distance by means other than the time it took light to reach us. We use triangulation from the separation of our two eyes. Or we use parallax: distant mountains don't move very much when we walk from side to side but nearby trees do. An astronomical scales these techniques are unworkable so we might measure the distance to objects by timing how long it takes light to reach us. Of course if there were a way to slow down light on its way to us, the object would appear much further but we normally wouldn't imagine such a thing. But as it turns out there is such a thing and on the large-scale its ubiquitous. It's a consequence of the Hubble flow.#### Hubble flow

Let's supposed two objects are 1m apart shortly after the Big Bang. And that the Hubble flow is inserting huge amounts of space between them. A photon of light might leave one object heading towards the other. But as it travels additional space is added between the two objects. This delays its arrival by quite a bit. When the second object finally receives this photon it might think " this photon took quite a while to reach me therefore the other object is many meters away". If we scale this analogy up we can see that light from objects that were close when they emitted that light, may have taken billions of years to reach us and so those objects merely appear to have been billions of light years away when they emitted it. So the Hubble flow is different than normal movement in that important way: objects are not in fact moving apart, space is being created between them and therefore in front of any photons traveling from one to the other. The Hubble flow is different than normal movement in another important way: objects cannot move apart faster than the speed of light. But the Hubble flow could insert so much space per second that they seem to be moving apart at faster than the speed of light. This point leads to some other remarkable consequences for example there may be objects we will never see because more space is being inserted than a photon can ever traverse even at the speed of light. And there is yet another possibility: there may some objects from whom light has not reached us yet but will eventually. Also there's a possibility that the Hubble flow changes over time. That there was a time early in the universe when a great deal of space was inserted suddenly. This is called the inflationary theory but as it would make this already long winded explanation even longer, I'll leave its implications to the readers imagination.PID controller for a target orientation

Suppose we have an object we want to rotate to a goal direction. So we have our current direction and a goal direction and we want to apply torque to reach that goal.
Let us express or current and goal direction as points on a unit sphere.
Likewise angular acceleration is expressed as a 3D vector axis and an angle. The angle is the acceleration around that axis that the object is experiencing.
So those two points on the unit sphere (current direction and goal direction) lie on a geodesic on that unit sphere. That geodesic tell us the direction we wish to accelerate. And of course the axis of angular acceleration must be perpendicular to this.
So we can compute

Vector3 axis = curentDirection cossProduct goalDirection

Now we need the magnitude of that acceleration, ie the angle of the angular acceleration axis angle. This is obviously related to the angle between curentDirection

and goalDirection

.
The magnitude of the cross product relates to this angle. We know that the magnitude of the cross product equals the product of the two vector lengths times the sine of angle between them.
ie magnitude of the cross product of A and B = magnitude of A times magnitude of B times sin( angle between A and B ) or pseudo code:
magnitude( A cossProduct B ) = magnitude( A ) x magnitude( B ) x sin( angleBetween( A and B ) )

Since the current and goal direction lie on a unit sphere, they are normalized by definition. So the product of their magnitudes is one. This leaves magnitude of axis = 1 * 1 * sin( angle between current and goal) or pseudo code:
magnitude( A cossProduct B ) = 1 x 1 x sin( angleBetween( A and B ) )

Thus to get that angle we merely need the arc sin of the magnitude of axis.
ie angle = arcsin(magnitude of axis ) resulting in units of radians/sec^2
New we can find the needed change in angular velocity
Vector3 deltaAngleVel = axis.normalized * angle

To convert this delta velocity into a torque, we will need to take rotational inertia into account . If there is inertia ix

to rotation around the global x axis and iy

resistance around the global y axis and iz

resistance to the z axis then we can calculate the torque with
Vector3 torque = deltaAngleVel.x * ix, deltaAngleVel.y * iy, deltaAngleVel.z * iz

We can just add this torque to a non rotating rigid body to create a spin of angle per second along the geodesic connecting curentDirection to goalDirection.
And if the rigidbody is already spinning in a different direction, we merely need to subtract the existing spin from the desired spin. so
generalizedDeltaAngleVel = deltaAngleVel - alreadyExsitingSpin

Applying this straight forward subtraction at intervals will bring the rotation to a stop at exactly the desired direction at an optimal speed with no oscillation.
#### PID

The astute may have noticed that the torque above strongly resembles a PID spring force and infact we can treat it eactly like this. We easily add a spring constant by doing: generalizedDeltaAngleVel = springConstant * ( deltaAngleVel - alreadyExsitingSpin ) This will produce an oscillation between two directions which we can fix with a damping force. Normally a damping force be simply calculated as dampingForce = -dampingConstant*curentVelocity But because we are dealing with torque we need to integrate it into the desired velocity from which we will calculate the torque. This is very simple: generalizedDeltaAngleVel = (springConstant * ( deltaAngleVel - alreadyExsitingSpin ))+( -dampingConstant* alreadyExsitingSpin) And torque is calculated just as above Vector3 torque = generalizedDeltaAngleVel.x * ix, generalizedDeltaAngleVel.y * iy, generalizedDeltaAngleVel.z * iznon commutative dice

non commutative dice have an interesting property resembling rock paper scissors: die A may beat die B and B may beat C but remarkably C will beat A.
This example doesn't really spread the odds evenly between every dice pair, but it makes the explanation pretty easy to follow (and so requires less writing on my part).
The first die, die A, is the simplest, it has a three on every face. So obviously it will beat any role less than three and lose to any role greater than three.
The die B is pretty straightforward too; there are fours on four of the faces and ones on two of the faces. It's immediately obvious that when we roll die A and die B together, die B will win 2/3 of the time. B beats A.
The die C is similar. With four twos and two fives. So similarly when rolled with die A, die A will obviously always role three. Die C will loose to four out of six times when it roles twos . It will win 2 out of 6 times when it rolls fives. So the case where A vs C makes A wins 2/3 of the time. A beats C.
The final case is the most complicated where die B and die C are rolled together. 1/3 of the time, die C will roll a five. In this case dieC wins no mater what die B does. So we see die C wins at least 1/3 of the time. But what about the rest of the time?
Well, of the remain the remaing 2/3 of the time, B will then roll a four two thirds of the time and win. And one third of this time, die B will roll a one and loose.
So the remain the remaing 2/3 of the time, die C wins 1/3 * 2/3 of the time or 2/9 of the time and losses 2/3 *2/3 = 4/9 of the time. Adding this to its automaitic 1/3 of the time wins mentioned above, we have: C wins 1/3 + 2/ 9 = 5/9 of the time, more than half. So when die B is rolled with die C, C wins more than half the time. C beats B.
So A beats C, B beats A, C beats B: there is no best die.

The apparent location of distant objects moving near light speed.

Typically this question is reversed: given that I see an object at P that moves a constant velocity V, what is its current position A. And this is very easily answered: P = A + V*t which we solve for A: A = P - V*t where t is the time light took to get to us from P
But in a simulation we need the inverse of this question: given the current position A what position B does the object appear at for an observer.
Any obvious answer looks like it will require a combination of trigonometry and calculus to solve, both of which are software performance nightmares. But what if there were a solution that only required matrix transforms and arithmetic? Then not only would we have acceptable performance, we might even be able to offload the work to the gpu.

#### And there is such a solution

The first step is to simplify the observer's coordinates so we don't have to worry about them. So we have an observer, the object's current position, and its trajectory. Let's move all these such that the observer is at the origin. This is obviously simple: just subtract the observers position from everything. Now we don't have to worry about the observer anymore. The second step is to reduce the problem to just one coordinate instead of three. Let's rotate the current position and trajectory so they both lie on the Y coordinate. To do this we treat the velocity as a direction and create a quaternion that rotates from that direction to {0,1,0} and apply it to the objects position and velocity. We form a quaternion from these two directions d1 and d2 from the cross product and dot product as follows:
a = crossProdcut( d1, d2 )
q.x = a.x, q.y = a.y, q.z = a.z
q.w = 1 + dotProduct( normalized( d1 ), normalized( d2 ) )
q = normalized( q )

Now we use that quaternion q to rotate the velocity to yield rotated velocity V. Because of this rotation velocity V only varies along the y coordinate, the others are now 0. Let's call this v.
So now the rotated object moves only along the y axis. The actual position of the object at any arbitrary time t is D(t). The x coordinate is unvarying x0 and the z coordinate is unvarying z0 and the y coordinate is v*t so D(t)={ x0, v*t, z0 }

likewise B(t) = { x0, y(t), z0 }

#### What is y(t)

At time t, we see the object at B(t) The time it take light to reach the origin (the observer) from B(t) is t' This is related to the distance from B(t) to the origin sot' = magnitude( B(t) ) / c

The object has been moving for time t' since it was seen at B(t) so the object was at B(t) but t' ago so
B(t) = D(t-t')

Substituting from above we have
D(t-t') = D( t – magnitude( B(t) ) / c )

So B(t) = D( t – magnitude( B(t) ) / c )

Now we only need to solve for B(t)
And because we have rotated everything so that only y changes, we only need to solve for y(t) at time t.
y(t) = v( t - t' )
y(t) = v ( t - magnitude( B(t) ) / c )

recall that B(t) = { x0, y(t), z0 }

so magnitude( B(t) ) = sqrt( x0^2 + y(t)^2 + z0^2 )

substituting we get
y(t) = v * ( t - sqrt( x0^2 + y(t)^2 + z0^2 ) / c )

solving for y(t) we get
y(t) = 1/(1-v^2/c^2) (v*t - v/c * sqrt( v^2*t^2 + (x0^2 + z0^2) * (1-v^2/c^2))

recall that the current position of the object is D( t ) = { x0, v*t, z0}

and now we know the observed position B( t ) = { x0, y(t), z0 }

this then needs to be rotated by the inverse of the quaternion used above and we're done.
Proposal for a Universal time

Having recently been annoyed by daylight savings time, I have decided there should be a universal and practical way to measure time:
Take the average temperature of the background radiation and peg time to the way it (or the average temperature of the universe ) cools down.
This won't give us datetime numbers for the very early universe but time may be screwed up then anyway.

Dual boot with MBR intact

Linux Windows dual boot is normally done by installing Grub in place of the MBR then chainloading to windows. This is well documented on several boards.
I haven't seen the reverse described where the MBR is left intact and allowing a boot to GRUB on another part of the disk.
This is actually a pretty good idea because the windows configuration is left intact, you don’t need to change the windows boot sector, removing Linux is trivial, you can do radical actions such as fdisk /mbr without worry.
And actually it is necessary in several cases. For example, certain Dell computers rewrite the MBR seemingly randomly. (Dell uses a technology called DSR which does this)
The process is relatively straightforward but obviously messing with the MBR can render machine unbootable so don't try this unless you know what's going on
1. Install Linux from a CD
2. Install GRUB on the boot partition of Linux (eg./dev/sda5) instead of MBR
3 Boot to Linux using using a CD
4. Extracted the GRUB bootloader onto an external device by copying the first 512 bytes of the boot sector to a diskette.
eg sudo dd if=/dev/sda5 of=grub.bin bs=512 count=1
5. Place this copy on the windows drive at " c:\grub.bin"
6. Tell windows about it by editing c:/boot.ini and add the operating systems c:\grub.bin="Ubuntu 4.10 "

Sweat of the brow IP

In the myriad discussions of intellectual property law, there is a generally unnoticed point that I believe has profound importance: Patent grants generally reject a "sweat of the brow" doctrine. The idea is that something that results form a flash of insight should also be eligible for patent.
But this fundamentally misunderstands the nature of innovation in the present day. A new invention does not occur when someone has an idea in the shower. That is merely an idea and ideas are specifically not patentable. New non-obvious inventions are the result of intensive research and study. The purpose of patents today is to protect capital-intensive R&D and allow compensation for it. And so the justification for patents is to protect that investment by keeping companies from free-loading on another's expensive R&D
By rejecting the "sweat of the brow" doctrine, patent law has become divorced from its economic justification: to stop free-loading. And indeed, where no great effort was involved in the invention, it very strongly resembles simple rent-seeking.

Interstellar trade

If there ever is a greater interstellar human civilization there is at least one interesting problem:
The time dilation colonist experience will mean they will be permanently behind the Earth by a number of years equal to the distance they travel. And their knowledge will be double that obsolescence.
For example: they leave at the speed of light for a star 50 light years away. On Earth 50 years pass but no time passes for the colonist. So when they arrive they are technically where the earth was 50 years before. They set up shop and transmit their wonderful ideas and technology, which takes 50 more years to be received. Which means it's 100 years out of date.
Meanwhile the Earth beams its ideas to the colonists. Even on they day the land, the colonist can tune into their favorite sitcom without missing an episode (because of course they were traveling side by side with those transmissions and the subsequent fifty years were following behind in a lovely tail of vapid entertainment). So the Earth's information is still quite useful to them.
And this disparity in value can never be corrected. If the Earth moans about the uselessness of the colonist 100 year old ideas and sends suggestions, the colonist will have been listening to the Earth for 100 years anyway before they hear those suggestions and will already moved far beyond them. Yet still be 50 years behind.
So a trade in ideas also seems unworkable.

very old stuff

Several ancient toys using java 1.0 to demonstrate some home made quaternion and other 3d transform libraries can be found here