All the Maths

[mikeal]

From a few conversations with @mikolalysenko I learned that languages like R and libraries like scipy which are known for amazing math are just binding to the same old Fortran libraries (and in a few cases C). That said, I don't think a binding layer will work for JavaScript, especially if we want to use them in the browser, but I do think we could use emscripten to transpile them in to JS, hopefully asm.js.

I got as far as writing most of a script to pull down calgo, which is laid out incredibly horribly, extract them all and generate package.json files which we could eventually checkin to git.

My attempts at getting emscripten to work ended at getting dragonegg to find the right llvm :(

Anyway, my instinct is that the best strategy is to:

• Pull down all these ancient libraries in to a repo and check them in to git (they aren't in any source control as is).
• Write scripts which can automate creating a package.json and individual npm package for each algorithm.
• Write benchmarks to get the VM engineers competing for performance, hopefully ending in superior performance to the original Fortran libraries.

Thoughts?

I know that @groundwater is interested in this as I imagine @substack and @maxogden are as well. Additionally, I saw @rwaldron throw out some links on twitter which I think were related to crazy math in JS.

[gaborcsardi]

@derekm

R does generate SSE/SIMD instructions depending on the implementation of specific libraries or routines. But, primarily, R achieves parallel vectorization by means of CPU multiprocessing, and does so inherently any time R's native array processing syntax is used.

I am not sure what exactly you mean here, but to me all of R's array processing seems very sequential. See e.g. double binary operators here: https://github.com/wch/r-source/blob/170a2a4a030e84a2437aee1a8d15c310a3b22527/src/main/arithmetic.c#L904 and the macros here: https://github.com/wch/r-source/blob/fbf5cdf29d923395b537a9893f46af1aa75e38f3/src/include/R_ext/Itermacros.h#L32

Am I missing something? How is this inherently parallel?

[tkelman]

@eush77

Operator overloading is rarely useful for regular use cases, and bringing more complexity to a single language because of some edge cases doesn't make any sense.

Overloading is very useful for numerical work. It's common to want to define a custom type that behaves like a number syntactically but with special behavior - arbitrary-precision arithmetic, polynomials, dual numbers for automatic differentiation, probability distributions, placeholder variables for modeling optimization problems, etc. But it's not very pleasant to do when overloading is tacked-on to a language that is otherwise not designed for it. Using JavaScript as a compilation target for a better-suited language seems like a better approach when it comes to numerical computing in the browser. It would likely save some work to pick an LLVM-based language to start from - Julia, Rust, or maybe even Swift would be reasonable choices. I've heard GHCJS works very well which would be another option, somewhat lacking in the array libraries or ease-of-use categories though.

@derekm

Python also has a threading model that nicely supports non-I/O bound numerical calculations.

Huh? Python's threading model for numerical calculations is "do threading (and all of the actual work) in C."

@gaborcsardi

all of R's array processing seems very sequential

Any time R calls out to BLAS, you get multithreading (and SIMD) by using OpenBLAS or MKL. You can run into issues with excessive fork-join, but I understand packages like plyr are designed to decompose problems at a higher level to mitigate this.

[gaborcsardi]

@tkelman:

Any time R calls out to BLAS, you get multithreading (and SIMD) by using OpenBLAS or MKL.

Right. But R does not use BLAS for array operations. AFAIK, these are the only operations it uses BLAS/LAPACK for: https://github.com/wch/r-source/blob/d75f39d532819ccc8251f93b8ab10d5b83aac89a/src/main/names.c#L926 (the do_lapack stuff). This is SVD, a linear solver, plus some decompositions, and nothing else.

but I understand packages like plyr are designed to decompose problems at a higher level to mitigate this.

Yes, some packages work around sequential R, but you need to do extra work to use them. It's not that you just load the package and then array operations are parallel.

IMO there is nothing inherently parallel in R or its packages. (One could actually override the default array operations to make them SIMD parallel, but I don't know of any package that does this. It is too much work, I guess.)

[tkelman]

But R does not use BLAS for array operations. AFAIK, these are the only operations it uses BLAS/LAPACK for

Ah, good to know, you know much more about the internals of R than I do. I'm a little surprised to not see gemv or gemm (matrix-vector or matrix-matrix multiplication, for those not familiar with the Fortran acronyms http://www.netlib.org/blas/blasqr.pdf) there.

One could actually override the default array operations to make them SIMD parallel, but I don't know of any package that does this.

Some packages use Rcpp to write extensions in C++ where you can do whatever you want in your own code, right? But maybe it's considered bad form to have your R package substantially change default behaviors.

[gaborcsardi]

Ah, good to know, you know much more about the internals of R than I do. I'm a little surprised to not see gemv or gemm (matrix-vector or matrix-matrix multiplication,

You are actually right, the above ones are only the ones called directly from R code. I missed a couple more, they are here:

~/works/r-source/src/main (trunk)$ grep F77_CALL *
array.c:        F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
array.c:    F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
array.c:    F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);
array.c:    F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
array.c:    F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
array.c:    F77_CALL(dsyrk)(uplo, trans, &nr, &nc, &one, x, &nr, &zero, z, &nr);
array.c:    F77_CALL(dgemm)(transa, transb, &nrx, &nry, &ncx, &one,
array.c:    F77_CALL(zgemm)(transa, transb, &nrx, &nry, &ncx, &one,
array.c:    F77_CALL(dtrsm)("L", upper ? "U" : "L", trans ? "T" : "N", "N",
main.c:    F77_CALL(intpr)("dummy", &i, &i, &i);

So there is indeed gemm as least (if your matrices do not contain NA). Btw. as most of BLAS/LAPACK is included in R, you can call BLAS/LAPACK from packages as well, and some packages do that, e.g. there are calls from the base stats package, but this again just a handful of functions.

So the situation is slightly better than what I said above, but not much, and array operations are certainly sequential.

Some packages use Rcpp to write extensions in C++ where you can do whatever you want in your own code, right?

Right. With some limits. (Btw. you don't need Rcpp for C/C++ code in R package, base R has an API for it. Rcpp just makes it easier.)

[derekm]

Okay... I only recently began prototyping stuff in R, it is hard not to notice the differences between vector operations and for loops, either in the literature or in your own code.

The literature is a little misleading claiming things in R having vectorizations or being vector operations or being a vectorized array processing language. When they say vectorization they must mean in syntax only and not in operation or processing, their claims having nothing to do with vector processors.

Performance differences must be solely interpreted R on for-loop syntax vs. compiled C loops on so-called "vector" syntax.

On Thursday, December 11, 2014, Gabor Csardi notifications@github.com wrote:

Ah, good to know, you know much more about the internals of R than I do. I'm a little surprised to not see gemv or gemm (matrix-vector or matrix-matrix multiplication,

You are actually right, the above ones are only the ones called directly from R code. I missed a couple more, they are here:

~/works/r-source/src/main (trunk)$ grep F77_CALL * array.c: F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one, array.c: F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one, array.c: F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc); array.c: F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one, array.c: F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one, array.c: F77_CALL(dsyrk)(uplo, trans, &nr, &nc, &one, x, &nr, &zero, z, &nr); array.c: F77_CALL(dgemm)(transa, transb, &nrx, &nry, &ncx, &one, array.c: F77_CALL(zgemm)(transa, transb, &nrx, &nry, &ncx, &one, array.c: F77_CALL(dtrsm)("L", upper ? "U" : "L", trans ? "T" : "N", "N", main.c: F77_CALL(intpr)("dummy", &i, &i, &i);

So there is indeed gemm as least (if your matrices do not contain NA). Btw. as most of BLAS/LAPACK is included in R, you can call BLAS/LAPACK from packages as well, and some packages do that, e.g. there are calls from the base stats package, but this again just a handful of functions.

So the situation is slightly better than what I said above, but not much, and array operations are certainly sequential.

Some packages use Rcpp to write extensions in C++ where you can do whatever you want in your own code, right?

Right. With some limits. (Btw. you don't need Rcpp for C/C++ code in R package, base R has an API for it. Rcpp just makes it easier.)

— Reply to this email directly or view it on GitHub https://github.com/node-forward/discussions/issues/1#issuecomment-66725503 .

[gaborcsardi]

Performance differences must be solely interpreted R on for-loop syntax vs. compiled C loops on so-called "vector" syntax.

I believe so. Plus there is the potential, that they will be parallel some time in the future.

[Planeshifter]

Great to see people working to make JavaScript better for numerical computing!

I do not see any compelling argument why R, a language which I am quite familiar with as an aspiring statistician, would be better suited for numerical work than JavaScript. However, the biggest disadvantage of JS at the moment is the lack of appropriate, well maintained libraries.

For example, it would be great to have a well-implemented linear-algebra package. I do not know of any such package currently. There are a couple node modules binding to C or C++ libraries (such as the mentioned eigenjs or armadillo, which I have been working on), but they are all unfinished. And the pure JS implementations often lack a lot of functionality, e.g. matrix decompositions.

Concerning operator-overloading, one could use sweet.js macros to implement it, at the cost of having to compile one's code.

Besides for node.js, it would be nice to have access to such libraries in the brower, too. For example, I recently wanted to create a d3 visualization plotting some simulation of logistic regression results. The ability to easily carry out the model fitting inside the browser would have been of great help in this case.

So porting some of the ancient math libraries to JS via emscripten sounds like a good idea to me. So far, I played around a bit with it and tried to port parts of the GNU Scientific Library (GSL) to JS via emscripten. While the process is trivial for functions accepting or returning standard types such as floats or doubles, the documentation is lacking in how to properly convert functions accepting or returning custom C structs (the Embind and WebIDL facilities are only tailored at C++ code).

For example, consider the following artificial example:

#include <math.h>

struct Vector3d{
  double x;
  double y;
  double z;
};

double length(struct Vector3d vec){
  double vec_sum = vec.x*vec.x + vec.y*vec.y + vec.z*vec.z;
  return sqrt(vec_sum);
}

struct Vector3d createVec(double x, double y, double z){
  struct Vector3d ret;
  ret.x = x;
  ret.y = y;
  ret.z = z;
  return ret;
}

I am very new to emscripten, so I have no idea how I could port such code as to make it callable from JS, in the sense that I want to end up with a constructor function than can create 3d vectors and a length function to calculate their length.

Has anyone figured out how to do this?

[smikes]

From the very small amount of time I've spent playing around with emscripten, it seems that it will compile a complete program + runtime environment into a javascript "program" that can be run in a browser. It does not have a facility for separately compiling one C/C++ function idiomatically into javascript.

One of my goals with femscripten (see above) is to make it possible (in the special case of reasonably pure functions) to produce javascript functions from C/C++/FORTRAN source, without having to link in the entire runtime+memory model that emscripten uses in order to emulate the complete program environment.

[Planeshifter]

The nice thing is that emscripten does provide facilities to expose functions directly to JS. I created a git repo with the gsl code and my shell script which exports a large number of functions for use in JavaScript. Here is the link: https://github.com/Planeshifter/gsl-js It is also immediately possible to require the created *.js file in node, as in var gsl = require("./output.js"). This will give you access to all the exported gsl functions listed in https://github.com/Planeshifter/gsl-js/blob/master/exported_functions

I have not investigated the performance hit incurred because of the overhead of having the emscripten runtime in the back, but the claim on their website is that it should be reasonably fast.

[smikes]

Oh, awesome. I hadn't found that yet. I will have to try that with the FORTRAN frontend, it may save me a lot of work. Thanks for pointing it out.

[Planeshifter]

Quick update: I asked on the emscripten repo about how to deal with structs, and meanwhile figured out how to export the C example I posted above via Embind: https://github.com/kripken/emscripten/issues/3083. So I do not think that there are any fundamental barriers preventing us from exposing entire C libraries to JS via emscripten.

[Planeshifter]

I played around a bit and ended up transpiling all of the complex number functionality of the GSL library into a node package. The GitHub repo with some documentation is here: https://github.com/Planeshifter/gsl-complex-js. Installation via npm: npm install gsl-complex. All functions are there, and using emscripten to export them ended up being straightforward (see the bottom of the complex.cpp file for the Embind bindings I had to create). I only had to create getter and setter functions for the real and imaginary part of a complex number in the original C code.