This is pretty cool. I've been using CoffeeScript to post-process MD simulation results because it's such an easy language to get things done in. Tools like this will make it even easier.
Although, could someone clarify for me the floating point issues mentioned by Steve? I wasn't aware there were issues...
It is the same as trying to store 1/3 = 0.333... as a finite decimal. Basically, a floating point number only has a finite space to put the component after the decimal point, so it gets truncated and made inaccurate, just like truncating 1/3 to 0.3333 will mean that calculations become wrong: 3*0.3333 = 0.9999 ≠ 1.
Ah, ok. Thank you. This isn't a problem for me then. It appears to be just typical floating-point math. I was thinking perhaps JavaScript did something else weird that I didn't know about.
Out of interest, what sort of MD/post-processing do you need to do? I've been in a similar position recently, but using CoffeeScript wouldn't have occurred to me.
By the by, for typical MD processing, as long as you're using sensible units you shouldn't need to worry about floating point precision. Just don't use SI units, because then you might end up with some very very small numbers and that's when it'll start to bite you!
For instance, I've been simulating some graphitic crystallites and part of the analysis includes analyzing the motion of the individual crystallite planes.
This is how easy that is:
(Note: it appears the horizontal scroll bar for this code is hidden by default for OS X browsers.)
eigenvalues = require 'eigenvalues.js'
regressCoordinatesToNormal = (coords) ->
min = (array) ->
array.reduce ((x, y, i) -> if x.val < y then x else {index: i, val: y}), {index: null, val: Infinity}
centerOfMass = (coords) ->
coords.reduce ((x, y) -> (x[i] + (y[i] / coords.length) for v, i in y)), [0, 0, 0]
normalUnitVector = (coords) ->
c = centerOfMass coords
matrix = for row in [0..2]
for col in [0..2]
((coord[row] - c[row]) * (coord[col] - c[col]) for coord in coords).reduce (x, y) -> x + y
results = eigenvalues.calculateEigenvalues matrix
results.eigenvectors[min(results.real).index]
normalUnitVector coords
This code just takes in a list of coordinates and performs a linear plane regression on them to return the normal vector of the plane. Then if I want the vertical distance between two planes of atoms (these are circular planes that kind of wobble around), then it's just:
normals = (regressCoordinatesToNormal coords for coords in planes)
planeToPlaneVector = distanceVector planeCenterOfMass[0], planeCenterOfMass[1]
projectionDistance = dotProduct planeToPlaneVector, normals[0]
I left out some details, but you can get the gist of it.
I can't really comment on that actually. For some reason, of all the programming languages I know, I never bothered to learn Python. If you need to make sure your results are 100% accurate then I would recommend a library that has been designed and tested by others, so NumPy would probably be better than creating your own routines. I do really like the clean syntax and support for functional(-ish) programming that CoffeeScript provides though.
Glad you like it! I'd love to know how it goes if you make use of it.
In your browser console, run:
0.3 + 0.544
You'll get back
0.8440000000000001
In practice (especially if in the browser), the difference is either irrelevant (in CSS 1px == 1.1px) or not significant. But in accumulation, it could be problematic.
That's why you don't compare floating point numbers with the = operator.
You use something like:
fabs( a-b ) < epsilon
where epsilon is usually defined by the language as a very small quantity (example values in C are 1E-8 for a float, 1E-15 for double), or you can just use a hard coded value appropriate to the values you are comparing.
That "where epsilon is usually defined by the language as a very small quantity" is BAD advice. double.epsilon and float.epsilon (std::numeric_limits::epsilon in C++) are "machine epsilons": numbers equal to the difference between 1 and the next representable value. In other words: 1 and 1+epsilon can be represented exactly, but there are no representable numbers between the two. The distance between representable numbers never goes down when you move away from zero. That means, that, for x,y >= 1 there is no difference between
x == y
and
|x - y| < epsilon
As I said, things get worse if your comparison is between larger numbers. For example, the smallest double larger than 1024 is 1024+1024 epsilon (IIRC; I am too lazy to double-check that now)
The epsilon used in numerical algorithms is an entirely different beast than the machine epsilon.
That epsilon you should pick as follows: make a numerical analysis of your problem, choose a good algorithm, and determine the desired accuracy of your answer. From those, derive a (relative, absolute, whatever) error you can live with.
Alternatively, pick a reasonable value from thin air and hope for the best/test your code to get confidence that it will return good values (do not do this when programming flight control software, pacemakers, etc). Oftentimes, it is not really hard to produce a reasonable value. For example, an iterative procedure that produces pixel coordinate likely can stop once the absolute error is less than .01 pixel, and possibly a lot earlier.
Thanks for the clarification, I did say they were examples (they were the first results I found in a quick google) and that you should use epsilon values appropriate to the values you are comparing.
Although, could someone clarify for me the floating point issues mentioned by Steve? I wasn't aware there were issues...