1. You should probably go ahead and implement at least LU decomposition, QR factorization, and maybe some eigenvalue/eigenvector stuff if you want to claim "advanced".
2. You should use the recursive definition of the determinant in the >3x3 matrix case.
3. You should consider adding some vector methods of the form: ax+by (for scalar a & b, and vector x & y). This can be useful for doing linear interpolation and whatnot.
4. You should add cross products/outer products for vectors.
5. You should have a mechanism for generating a Gaussian or normal distribution, or at least supply some standard mappings from your (presumably uniform) random function to something more useful (Gaussian, Possoin, Weibull, whatever).
6. An FFT would be nice. :)
7. Quaternions would also be handy for graphics folks, or motion folks.
8. Complex number support would be cool.
9. Support for optimization, like integer programming, linear programming, or even least-squares or whatever would be nice.
~
Again, great work--I look forward to seeing how this project evolves!
This is a fantastic list. There are a few people who've expressed interest in working on eigenvalue work. I'll definitely be adding a couple of these in, in the next few days. I'll keep you posted.
Sparse and dense matrix math and eigensolvers would be really nice. Especially solvers that let you efficiently get the eigenvectors associated with the smallest eigenvalues as this comes up in graph and image partitioning.
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.
Agreed. I would have expected to see things like ICA, PCA, SVD, pseudorandom number generators for various discrete and continuous distributions, misc. multiplicative and arithmetic functions, misc. metrics (of the minkowski kind) and so on to label the library as advanced.
I've used http://www.numericjs.com/ to handle basic vector/matrix arithmetic and for numerically solving minimization problems to create constraint-based UIs.
Can you give an overview of how your library is different? Is it just a different set of math tools or is the architecture somehow different? Thanks!
Yeah, numeric.js is awesome and very fast compared to other libraries. Would be interesting to see a comparison in terms of performance between these two.
When you drag and drop the shapes, the drawing changes in non-trivial ways. For example, try dragging a shape deep down in the recursion hierarchy.
In a normal drag-and-drop application, it's fairly easy to compute how to, say, adjust the left/top of a div in response to the mouse movement events. However, with Recursive Drawing, I knew that there were certain properties I'd have to adjust during a drag operation, but because of all the nested transformations, it was tricky to figure out the math of how exactly to adjust them.
So instead, I set up a constraint problem. I know that the exact spot I've mousedown'ed on needs to stay under my mouse no matter what. That's the fundamental constraint of the drag-and-drop gesture. Then I set the properties that I knew I could change (e.g. positions or scaling/rotation of a specific shape) and set them as free variables in an energy minimization function. I solved the problem numerically using numeric.js's uncmin (unconstrained minimization, an algorithm originally written in FORTRAN, I believe).
I probably could have figured out how to do this without numerical methods, but this approach was a god-send when rapid prototyping the interactions.
Another good use of constraint systems for UIs is Ivan Sutherland's Sketchpad. His paper on it is long but really worth at least skimming through if you're interested in developing next-generation UIs!
http://www.cl.cam.ac.uk/techreports/UCAM-CL-TR-574.pdf
I also really like Rebecca Fiebrink's Wekinator, which is a framework for making musical instruments from arbitrary inputs. The approach is similar in that she uses numerical methods to solve for the constraints (the training data).
http://wekinator.cs.princeton.edu/
I think there's a lot of potential in constraint-based UIs and for me it's the most compelling reason to have good numerical libraries in javascript.
To hazard a guess, I would say it's setting up min/max/preferred width/height or margins on UI components and then dynamically computing the best way to lay them out at display time or when the window size changes. Java has a GridBagConstraints layout manager class that does something like this, and similar layout management shows up elsewhere. Apple's new AutoLayout is a collection of constraints too.
This is awesome! I started thinking about doing something similar a few months ago, but the project I needed it for kind of fell through and I moved on. Anyway, a few thoughts I had at the time:
-should support FFT/convolution
-could tie into http://sylvester.jcoglan.com/ for
vector/matrix stuff (excellent library!)
-I emailed jcoglan about the possibility of adding
complex number support to Sylvester, he encouraged me
to try but thought it might introduce too much function
call overhead (caching function results might help a
bit?)
-I didn't get to look at it too closely, but
https://github.com/jtobey/javascript-bignum looks
interesting, as does https://github.com/pr1001/MathPlus
-polynomial fit (least squares regression) and
interpolation stuff is in high demand
-it should be super easy to integrate output with d3.js
(like the equivalent of MATLAB plot function)
Honestly, I doubt I would have been able to pull it off, but I intend to contribute to this project as much as I can.
Great, it'd be sweet to have you involved! FFT, polynomial fit, and d3.js would be huge steps. For me, the most important thing is ensuring that all changes are application driven. Building a library that's designed to be used is top priority.
Great. I'm studying for my operations research exam, and I've wanted to "study" by implementing some of the algorithms we did by hand through code.
Is there interest in some more advanced matrix algebra functions - specifically Gaussian reduction, simplex, sensitivity analysis, and maybe some transport/network problems and integer programming?
My primary focus from here on out is statistics and matrix algebra focused on powering analytics. So yes, some of those functions would be of interest to me.
What's the use-case-scenario for this? Seems like if you had complex calculations you needed to perform, you would do it on a server side language like Python.
A couple of observations:
1. You should probably go ahead and implement at least LU decomposition, QR factorization, and maybe some eigenvalue/eigenvector stuff if you want to claim "advanced".
2. You should use the recursive definition of the determinant in the >3x3 matrix case.
3. You should consider adding some vector methods of the form: ax+by (for scalar a & b, and vector x & y). This can be useful for doing linear interpolation and whatnot.
4. You should add cross products/outer products for vectors.
5. You should have a mechanism for generating a Gaussian or normal distribution, or at least supply some standard mappings from your (presumably uniform) random function to something more useful (Gaussian, Possoin, Weibull, whatever).
6. An FFT would be nice. :)
7. Quaternions would also be handy for graphics folks, or motion folks.
8. Complex number support would be cool.
9. Support for optimization, like integer programming, linear programming, or even least-squares or whatever would be nice.
~
Again, great work--I look forward to seeing how this project evolves!