The Kitchin Research Group

I do a lot of scientific programming, and it is one of the reasons I have been learning to extend Emacs with dynamic modules. They have allowed me to add physical constants, numerical integration, root finding and linear algebra from established c-libraries to Emacs. Today I am taking a break from that and finally getting to another one of the reasons I started playing around with dynamic modules: zeromq. Zeromq is a messaging library that Jupyter uses to communicate with kernels. I thought we might get a smoother integration with Emacs and Jupyter if we could use zeromq directly to communicate between org-mode and the kernel. Currently we have to run a web server that does the communication for us via http requests. We won't solve the Jupyter problem today, but we will look at communicating with a Zeromq server from Emacs.

This might have lots of other useful applications. Suppose Emacs could communicate directly with other zeromq servers to retrieve data from, perhaps scientific data. It might even be possible for Emacs to run its own zeromq server, and other instances of Emacs could then communicate with it. Collaborative editing anyone?

Here we just implement the "Hello world" client example in the zeromq guide. The code for the server, a c-client, the mod-zmq library, a makefile and tests can be found at https://github.com/jkitchin/emacs-modules/tree/master/zeromq. All the server does is receive a string, and then send a response (in this case just the string "World") back to the client.

To run this, make sure to run the hwserver executable in a terminal. I wrapped the zeromq commands required to implement the client into a dynamic module. Since this example focuses on strings, the module returns strings to Emacs. I am not sure if that is always the right thing to do, as zeromq more generically uses bytes, but I will have to wait until I know more about zeromq to know if this is an issue.

This dynamic module uses a new feature that none of the previous posts used, and that is the user_ptr. These allow you to essentially return a reference pointer back to emacs that you can pass back to another function. That way they stay alive between function calls. For example, here we have to create a context and socket and pass these items to functions like zmq_send and zmq_recv.

The directory this library is in is not on my path, so we load it like this:

(add-to-list 'load-path (expand-file-name "."))
(require 'mod-zmq)

Here are the functions and their signatures that have been implemented so far. I only implemented the ones I needed for the client. The signatures may change in the future; this is just a proof of concept for now for the purpose of building the client.

(apropos-command "zmq*" t)
(with-current-buffer "*Apropos*" (buffer-string))

Type RET on a type label to view its full documentation.

zmq-close
  Function: (zmq-close SOCKET)
zmq-connect
  Function: (zmq-connect SOCKET ENDPOINT)
zmq-ctx-destroy
  Function: (zmq-ctx-destroy CONTEXT)
zmq-ctx-new
  Function: (zmq-ctx-new)
zmq-recv
  Function: (zmq-recv SOCKET LEN FLAGS)
zmq-send
  Function: (zmq-send SOCKET BUF FLAGS)
zmq-socket
  Function: (zmq-socket CONTEXT TYPE)

You can see the c code for the client here: hwclient.c . Here is a simple elisp version of the hwclient that basically does the same thing! The main difference is I added a while loop around the zmq-recv because sometimes it returns -1 and no result. So, here we loop until the return result is not -1. That seems to do the right thing.

(let* ((context (zmq-ctx-new))
       (socket (zmq-socket context ZMQ-REQ))
       (recv-ret -1)
       (result))

  (zmq-connect socket "tcp://localhost:5555")
  (zmq-send socket "Hello" 0)

  (while (= recv-ret -1)
    (setq result (zmq-recv socket 10 0)
          recv-ret (second result)))

  (print result)

  (zmq-close socket)
  (zmq-ctx-destroy context))

("World" 5)

Basically this creates the context, then the socket, and connects to it on port 5555 of the localhost where the server is running. Then we send the string "Hello". The server returns the string "World" and tells us it sent 5 bytes. Then we close the socket and destroy the context. There is a lot of code in the module to make this happen. A lot of it is converting args in emacs functions to things we can use in c, running a few lines of zmq commands, and then code to convert those results back to emacs values. Finally, there is code to register each function and define docstrings for them. I am not totally convinced this is the best way to do this, but it does work! An alternative might be emacs-ffi, which might enable most of this to be developed in just elisp.

Copyright (C) 2017 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 9.0.7

The goal of this post is to be able to solve equations like this one:

\[\left(\begin{array}{cccc} 0.18& 0.60& 0.57& 0.96 \\ 0.41& 0.24& 0.99& 0.58 \\ 0.14& 0.30& 0.97& 0.66 \\ 0.51& 0.13& 0.19& 0.85 \end{array} \right ) \left ( \begin{array}{c} x_0 \\ x_1 \\ x_2 \\ x_3 \end{array} \right ) = \left ( \begin{array}{c} 1.0 \\ 2.0 \\ 3.0 \\ 4.0 \end{array} \right ) \]

The answer is given as

\[x = \left ( \begin{array}{c} -4.05205 \\ -12.6056 \\ 1.66091 \\ 8.69377 \end{array} \right ) \]

The syntax we want to use is shown below, and we want it to return a vector containing the solution:

(let ((A [[0.18 0.60 0.57 0.96]
          [0.41 0.24 0.99 0.58]
          [0.14 0.30 0.97 0.66]
          [0.51 0.13 0.19 0.85]])
      (b [1.0 2.0 3.0 4.0]))
  (gsl-linalg-LU-solve A b))

Rather than put all the code in here like I have for the past several posts, I started a git repo at https://github.com/jkitchin/emacs-modules that contains this code.

The module for this post can be found here: https://github.com/jkitchin/emacs-modules/blob/master/gsl-linalg.c. There are a few notable features in it. First, I started writing/collecting some helper functions to make these modules simpler to write. For example, look how nice this looks to declare the functions and provide the feature.

int emacs_module_init(struct emacs_runtime *ert)
{
  emacs_env *env = ert->get_environment(ert);
  
  DEFUN("gsl-linalg-LU-solve", Fgsl_linalg_LU_solve, 2, 2,
        "(gsl-linalg-LU-solve A b).\n" \
        "Solve A x = b for x.\n" \
        "Returns a vector containing the solution x.",
        NULL);
  provide(env, "gsl-linalg");
  
  return 0;
}

The DEFUN and provide function are defined in https://github.com/jkitchin/emacs-modules/blob/master/emacs-module-helpers.c.

Within the module itself, we have to loop over the inputs to create the arrays that GSL wants to solve this problem. Second, after the solution is obtained, we have to build up a vector to return. The solution is in a gsl_vector, and we need to create an array of emacs_value elements containing each element of the gsl_vector as a float, and then create a vector to return to emacs. I use vectors here because it was easy to get the size of the b vector, which is also related to the size of the A matrix.

The repo has a Makefile in it, so we can build this module with:

make gsl-linalg.so

Once it is compiled, we load it like this. In this post we are in the emacs-modules directory where the gsl-linalg.so library is, and it is not on my load-path, so I add it here.

(add-to-list 'load-path (expand-file-name "."))
(require 'gsl-linalg)

gsl-linalg

Here is one function in the module:

(describe-function 'gsl-linalg-LU-solve)

gsl-linalg-LU-solve is a Lisp function.

(gsl-linalg-LU-solve &rest ARGS)

For more information check the manuals.

(gsl-linalg-LU-solve A b).
Solve A x = b for x.
Returns a vector containing the solution x.

Now, we can solve linear equations like this:

(gsl-linalg-LU-solve
 [[0.18 0.60 0.57 0.96]
  [0.41 0.24 0.99 0.58]
  [0.14 0.30 0.97 0.66]
  [0.51 0.13 0.19 0.85]]
 [1.0 2.0 3.0 4.0])

[-4.052050229573973 -12.605611395906903 1.6609116267088417 8.693766928795227]

We have a limited ability to confirm this answer. I have written a function that uses blas for multiplication of 2d vectors. You can see from this:

(gsl-blas-dgemm [[0.18 0.60 0.57 0.96]
                 [0.41 0.24 0.99 0.58]
                 [0.14 0.30 0.97 0.66]
                 [0.51 0.13 0.19 0.85]]
                [[-4.052050229573973]
                 [-12.605611395906903]
                 [1.6609116267088417]
                 [8.693766928795227]])

[[1.0] [1.9999999999999991] [2.9999999999999996] [4.0]]

That within float that indeed \(A x = b\).

The main limitation of this module at the moment is that you have to use vectors; you cannot put in a list of numbers. It is possible to make it take lists and vectors, but for now I am leaving it at vectors. Also, it only produces solutions of float numbers (not integers).

The module does not handle 1d vectors well,, e.g. in gsl-linalg-LU-solve example, the right hand side is implied to be a column vector, and we don't have the array broadcasting features of Python yet. Those are doable things for some future day perhaps. For now I am happy to have figured out how to handle arrays!

Copyright (C) 2017 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 9.0.7

The GNU Scientific Library defines a lot of physical constants. Since we are exploring how to make Emacs a more scientific environment to work in, it would be nice to import these constants to elisp. We do that through a dynamic module. This turned out to be tricky. I thought we could just use a funcall to defconst or defvar, but these are special forms and you cannot funcall them. @polizta on Stackoverflow pointed me to the path that led to success: You make a list like '(defconst sym val doc) and then eval it. That can be funcall'd, and it works nicely in the module below. It is a growing theme that it takes much hacking around to figure out how to do things like this.

The only other notable feature of this module is that I created a defconst function to make adding multiple constants less verbose. Here I only add two constants. There are about 408 constants defined in gsl_const_*.h, so brevity might be a good idea! Here is the module.

#include <gsl/gsl_const_mksa.h>
#include <string.h>
#include "emacs-module.h"

int plugin_is_GPL_compatible;

// I assume here that all values will be double. I can't think of any that would
// be ints
static void defconst (emacs_env *env, const char *name, double value, const char *doc)
{
  // These are functions we will call
  emacs_value eval = env->intern(env, "eval");  
  emacs_value list = env->intern(env, "list");

  // These will make up the list we will eventally eval
  emacs_value fdefconst = env->intern(env, "defconst");
  emacs_value sym = env->intern(env, name);
  emacs_value val = env->make_float(env, value);
  emacs_value sdoc = env->make_string(env, doc, strlen(doc));

  // make a list of (defconst sym val doc)
  emacs_value largs[] = {fdefconst, sym, val, sdoc};
  emacs_value qlist = env->funcall(env, list, 4, &largs);   

  // now eval the list of symbols
  emacs_value args[] = { qlist };  
  env->funcall(env, eval, 1, &args);
}

int emacs_module_init(struct emacs_runtime *ert)
{
  emacs_env *env = ert->get_environment(ert);

  defconst(env, "GSL-CONST-MKSA-SPEED-OF-LIGHT",
           GSL_CONST_MKSA_SPEED_OF_LIGHT,
           "Speed of light in vacuum (m/s).");
  
  defconst(env, "GSL-CONST-MKSA-PLANCKS-CONSTANT-H",
           GSL_CONST_MKSA_PLANCKS_CONSTANT_H,
           "Plank's constant, h");

  // This is what allows the shared library to provide a feature 
  emacs_value provide = env->intern(env, "provide");
  emacs_value provide_args[] = { env->intern(env, "gsl-constants") };
  env->funcall(env, provide, 1, provide_args);
  
  return 0;
}

Regular gcc will work to compile this module.

rm -f gsl-constants.so gsl-constants.o
gcc -Wall -I/usr/local/include -fPIC -c gsl-constants.c
gcc -shared -L/usr/local/include -lgsl -o gsl-constants.so gsl-constants.o

Here is in action.

(add-to-list 'load-path "/Users/jkitchin/vc/blogofile-jkitchin.github.com/_blog/dynamic-module/")
(require 'gsl-constants)
GSL-CONST-MKSA-SPEED-OF-LIGHT

299792458.0

We can see there is a docstring on that constant:

(describe-variable 'GSL-CONST-MKSA-SPEED-OF-LIGHT)

GSL-CONST-MKSA-SPEED-OF-LIGHT's value is 299792458.0


  This variable can be risky when used as a file-local variable.

Documentation:
Speed of light in vacuum (m/s).

For more information check the manuals.

It is worth thinking about what we accomplished here. The value of each constant in GSL is stored in a header file. The units are stored in a comment next to the value, and the documentation is in an html page somewhere. It is not easy to introspect that! Getting it all into an Emacs variable makes that more introspectable, and findable. That means while typing elisp code you will get completion on these variables. Check this out:

(apropos-variable "GSL-*")
(with-current-buffer "*Apropos*" (buffer-string))

Type RET on a type label to view its full documentation.

GSL-CONST-MKSA-PLANCKS-CONSTANT-H
  Variable: Plank's constant, h
GSL-CONST-MKSA-SPEED-OF-LIGHT
  Variable: Speed of light in vacuum (m/s).

It seems like it might be possible to partially automate creation of this module by parsing the gsl_const*.h files. There is no automating adding the doc strings though, I am pretty sure that will have to be done by hand ;(

Copyright (C) 2017 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 9.0.7

Here we work out how to run this program: https://www.gnu.org/software/gsl/doc/html/integration.html#adaptive-integration-example in a dynamic module in emacs. The goal is to be able to evaluate \(\int_0^1 x^{-1/2} \log(x) dx\). According to the example page the answer is -4. We will define an integration function that takes at least a function and integration bounds as arguments, and several optional arguments to specify tolerances and limits. In other words we want to evaluate integrals of the form:

\(\int_a^b f(x; params) dx\)

I want that to happen in an elisp function with a signature like:

(gsl-integration-qags (lambda (x params) body) a b &optional params epsabs epsrel limit)

And that function will return a list containing (result error-estimate). Here is the C-code that makes this happen. It is more complex that the last example, and only compiles with gcc that allows nested functions. I don't know how to write this without that feature. This is more complex also because you have to create a workspace to do the integration inside the function that does the integration. The C-module also has extra code in it to allow for optional arguments.

#include <gsl/gsl_integration.h>
#include "emacs-module.h"

int plugin_is_GPL_compatible;

static emacs_value F_gsl_integrate (emacs_env *env, ptrdiff_t nargs, emacs_value args[], void *data)
{
  // nested function - only supported as an extension in gcc
  double f (double x, void *params) 
  {
    emacs_value fn = args[0];  // function we will integrate
    emacs_value x2[] = { env->make_float(env, x), params };
    emacs_value y = env->funcall(env, fn, 2, &x2);   
    
    return env->extract_float (env, y);
  }

  double a = env->extract_float (env, args[1]);
  double b = env->extract_float (env, args[2]);

  // default values for optional arguments
  double epsabs = 0.0;
  double epsrel = 1e-7;
  size_t limit = 1000;
  double result, error; 

  // Here is how I handle the optional arguments
  // (gsl-integrate func a b params epsabs epsrel limit)
  gsl_function F;
  F.function = &f;
  if (nargs >= 4) {F.params = args[3];}
  if (nargs >= 5 && env->is_not_nil(env, args[4])) {epsabs = env->extract_float(env, args[4]);}
  if (nargs >= 6 && env->is_not_nil(env, args[5])) {epsrel = env->extract_float(env, args[5]);}
  if (nargs >= 7 && env->is_not_nil(env, args[6])) {limit = env->extract_integer(env, args[6]);}

  gsl_integration_workspace * w = gsl_integration_workspace_alloc (limit);

  gsl_integration_qags (&F, // gsl_function pointer
                        a, // lower integration bound
                        b, // upper integration bound
                        epsabs, // absolute error tolerance
                        epsrel, // relative error tolerance
                        limit, // max number of subintervals for integration
                        w, // the workspace
                        // pointers to put results and error in
                        &result, &error);

  gsl_integration_workspace_free (w);
    
  // make a list of (result error) to return
  emacs_value Qlist = env->intern(env, "list");
  emacs_value Qresult = env->make_float (env, result);
  emacs_value Qerror = env->make_float (env, error);
  emacs_value list_args[] = { Qresult, Qerror };
  return env->funcall(env, Qlist, 2, list_args);
}

int emacs_module_init(struct emacs_runtime *ert)
{
  emacs_env *env = ert->get_environment(ert);
  
  // Here we create the function.
  emacs_value fset = env->intern(env, "fset");
  emacs_value args[2];
  args[0] = env->intern(env, "gsl-integration-qags"); // symbol to create for function
  // The function we set that symbol to.
  args[1] = env->make_function(env,
                               3, // min nargs
                               7, // max nargs
                               F_gsl_integrate,
                               "(gsl-integration-qags F A B &optional PARAMS EPSABS EPSREL LIMIT)\n" \
                               "Integrate F(x; params) from A to B.\n" \
                               "F is a function of a single variable and parameters.\n" \
                               "A is the lower bound of integration\n"  \
                               "B is the upper bound of integration.\n" \
                               "Optional parameters:\n"\
                               "PARAMS is a list of params to pass to F.\n" \
                               "EPSABS is a float (default 0.0) and is the absolute error tolerance.\n" \
                               "EPSREL is a float (default 1e-7) and is the relative error tolerance.\n" \
                               "LIMIT is the maximum number of subintervals for the integration (default 1000).\n" \
                               "Returns (list result error-estimate).\n" \
                               "See https://www.gnu.org/software/gsl/manual/html_node/QAGS-adaptive-integration-with-singularities.html.",
                               0);
  // This is basically (fset 'gsl-integration-qags (lambda func))
  env->funcall(env, fset, 2, args);
  
  // This is what allows the shared library to provide a feature 
  emacs_value provide = env->intern(env, "provide");
  emacs_value provide_args[] = { env->intern(env, "gsl-integration") };
  env->funcall(env, provide, 1, provide_args);
  
  return 0;
}

Building this was moderately tricky. It appears the first gcc on my path uses clang which does not support nested functions in C. I don't know enough C to figure out how to do this without a nested function though, since the function has to be defined at run-time based on the emacs env and args. gcc does support inline functions, so the code below uses a gcc that does compile it.

rm -f gsl-integration.so gsl-integration.o
/usr/local/Cellar/gcc/6.1.0/bin/gcc-6 -Wall -I/usr/local/include -fPIC -c gsl-integration.c
/usr/local/Cellar/gcc/6.1.0/bin/gcc-6  -shared -L/usr/local/include -lgsl -o gsl-integration.so gsl-integration.o

Now we add this directory to our path since it is not on it and require our new module.

(add-to-list 'load-path "/Users/jkitchin/vc/blogofile-jkitchin.github.com/_blog/dynamic-module/")
(require 'gsl-integration)

gsl-integration

Let us see our new function in action. We evaluate \(\int_0^1 x^{-1/2} \log(x) dx\). According to the example page the answer is -4. Here is an example where we ignore the parameters. You have to be careful; Emacs sometimes segfaults and crashes if you use an integer or float argument when it expects the other type.

(gsl-integration-qags (lambda (x params) (/ (log x) (sqrt x))) 0.0 1.0)

Here are some optional arguments.

(gsl-integration-qags (lambda (x params) (/ (log x) (sqrt x))) 0.0 1.0 nil nil 0.01)

Nice, with a larger epsrel argument we get a larger error. Note the arguments are positional, so we have to include them all just to set the epsrel argument. How about an easier example with parameters that we actually use. Here we integrate a constant, and set the value of the constant from the params arg. The integral should be the area of a rectangle of length 1 and width of the param we use.

(list
 (gsl-integration-qags (lambda (x params) (first params)) 0.0 1.0 '(1.0))
 (gsl-integration-qags (lambda (x params) (first params)) 0.0 1.0 '(0.5)))

Wow! It actually works!!! That was harder won success than usual for me. I am claiming victory for now and leaving the following notes to future me:

1. It would be nice to have optional keyword arguments. This would take some handling of the arguments beyond what I know how to do for now, unless it is possible to pull in something like plist-get the way we pull in fset, provide and list in this example.
2. Error checking on types would be helpful. It is not good for Emacs to crash because 0 is not 0.0!
3. In numpy there is often a feature to get full_output. Here, the workspace created in the function has more information available in a struct that might be helpful to have access to at times. It seems like it might be possible to get that here too.

Copyright (C) 2017 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 9.0.7

There is a relatively new feature in Emacs 25 that allows you to extend Emacs using compiled libraries (http://diobla.info/blog-archive/modules-tut.html). This could be very helpful in a few ways:

1. To add functionality that exists in other libraries, e.g.
   1. libyaml
   2. libmemcached
   3. Embedding Ruby in Emacs
2. Interface Emacs with hardware, e.g. a joystick, or ejecting a CD.
3. To speed up slow elisp functions
   1. A c implementation of a fibonacci function is 150 times faster than an elisp version here.
   2. This json parser is up to 4 times faster than the json.el library for some operations.

I am interested in this in particular to bring numerical methods into Emacs. It is fair to ask why. Even I think the numpy/scipy/matplotlib Python stack is currently unparalleled in functionality for scientific programming. But I like writing elisp code so much more! So, we will take a look today at a simple example of integrating a function from the GNU Scientific Library into Emacs.

#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>

int
main (void)
{
  double x = 5.0;
  double y = gsl_sf_bessel_J0 (x);
  printf ("J0(%g) = %.18e\n", x, y);
  return 0;
}

gcc -Wall -I/usr/local/include/ -c example.c
gcc -L/usr/local/include -lgsl example.o 
./a.out

(require 'gsl-sf-bessel)
(gsl-sf-bessel-J0 5)

#include <assert.h>
#include <gsl/gsl_sf_bessel.h>
#include "emacs-module.h"

int plugin_is_GPL_compatible;

static emacs_value
F_gsl_sf_bessel_J0 (emacs_env *env, ptrdiff_t nargs, emacs_value args[], void *data)
{
  assert (nargs == 1);
  double x = env->extract_float (env, args[0]);
  return env->make_float (env, gsl_sf_bessel_J0 (x));
}

int
emacs_module_init(struct emacs_runtime *ert)
{
        emacs_env *env = ert->get_environment(ert);

        emacs_value gsl_sf_bessel_J0_fn = env->make_function(env, 1, 1, F_gsl_sf_bessel_J0, "Regular cylindrical Bessel function of zeroth order, J_0(x)", NULL);

        emacs_value Qfset = env->intern(env, "fset");
        emacs_value Q_gsl_sf_bessel_J0 = env->intern(env, "gsl-sf-bessel-J0");
        emacs_value fset_args[] = { Q_gsl_sf_bessel_J0, gsl_sf_bessel_J0_fn };
        env->funcall(env, Qfset, 2, fset_args);

        emacs_value Qprovide = env->intern(env, "provide");
        emacs_value Q_gsl_sf_bessel = env->intern(env, "gsl-sf-bessel");
        emacs_value provide_args[] = { Q_gsl_sf_bessel };
        env->funcall(env, Qprovide, 1, provide_args);

        return 0;
}

gcc -Wall -I/usr/local/include -fPIC -c gsl-sf-bessel.c
gcc -shared -L/usr/local/include -lgsl -o gsl-sf-bessel.so gsl-sf-bessel.o

ls *.so

(add-to-list 'load-path "/Users/jkitchin/vc/blogofile-jkitchin.github.com/_blog/dynamic-module/")
(require 'gsl-sf-bessel)
(gsl-sf-bessel-J0 5.0)

-0.17759677131433826

(describe-function 'gsl-sf-bessel-J0)

gsl-sf-bessel-J0 is a Lisp function.

(gsl-sf-bessel-J0 &rest ARGS)

For more information check the manuals.

Regular cylindrical Bessel function of zeroth order, J_0(x)

Copyright (C) 2017 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 9.0.7