Adding numerical methods to emacs with dynamic modules

| categories: emacs | tags: | View Comments

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.

1 Using the GSL to calculate a Bessel function value

A canonical example of using the GSL is given at https://www.gnu.org/software/gsl/manual/html_node/An-Example-Program.html. Here it is for reference. It just calculates a value for a Bessel function. We save this program in a file called example.c.

#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;
}

We have to compile and run this program. Here are the commands to do that.

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

That is a lot of code and steps to get one number.

What I would like to do instead is this:

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

So, enter the dynamic module!

2 A GSL dynamic module for a Bessel function

To create the dynamic module we need a small c file that wraps the GSL function and adds it to the Emacs environment. Here is the smallest example I could come up with for this one function. Basically, we create a function that takes the emacs environment and arguments, extract what we want from them, and use that to calculate what we need and return it to the environment. Then we register that function and define what the module provides in emacs_module_init.

#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;
}

Now we compile it into a shared library.

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

That creates our shared library in gsl-sf-bessel.so.

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

That is the same answer as we got before. Here is the documentation we defined. It could use some improvement, e.g. to note that the argument has to be a float, and that only one argument is allowed. I am not sure why the signature doesn't show a single argument.

(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)

I am not a very skilled C-programmer yet, so I don't know how hard it would be to make this function accept integers as well, or to vectorize it so you could have an arbitrary number of args to it and return a list.

3 Summary

Dynamic modules look promising to extend Emacs with. This example is about the simplest function from the GSL there is. There are many more (https://www.gnu.org/software/gsl/doc/html/index.html) functions that do linear algebra on arrays, integration or optimization of functions, interpolation of data, etc. I don't have a sense yet of how easy it will be to integrate these into a module.

It looks like you are not limited to writing these in C. There is an example of a plugin written in Rust here, and a framework to write them in Go. Maybe any language that can make a shared library with the required plugin_is_GPL_compatible symbol and emacs_module_init function would work. Those examples do not look significantly easier to write than the C versions though since I am not that fluent in those languages either.

There are some challenges to figure out in developing and using dynamic modules. Here are a few:

  1. The documentation on what is possible is not that great yet, so there is a lot of exploring to do. There are a fair number of examples out there though to learn from (https://github.com/emacs-pe/emacs-modules). The official example shows a lot of the functionality.
  2. I guess it will be tricky to distribute these. I don't know how easy it would be to build all the libraries for each platform for distribution on MELPA for example. I don't think there is a standard way to incorporate a compile step in elisp package installation. Also, you need an Emacs version of at least 25 with the dynamic module feature compiled in. It is not yet a default enabled option. The required emacs-module.h should be gotten from the emacs build, so people with binaries might not be able to build it anyway.
  3. Users will need the libraries the dynamic module uses. In this example, they will need libgsl.
  4. Once you require the module, it does not seem possible to modify it, rebuild it, and reload it. It appears you have to close Emacs and reload it. That is tedious.

It would be nice to have a more generic foreign function interface that would allow you to develop more on the elisp side. One effort in that direction is https://github.com/tromey/emacs-ffi. It looks like it might be a lot simpler to use than creating a dynamic module. Once it is installed, it looks like you can write elisp code to wrap the library functions. I will write about this on another day.

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

org-mode source

Org-mode version = 9.0.7

Read and Post Comments

Adding keymaps to src blocks via org-font-lock-hook

| categories: emacs, orgmode | tags: | View Comments

Table of Contents

I had an idea to use custom keymaps in src-blocks. For example, you could then use lispy directly in your org-files without entering org-special-edit, or the elpy key-bindings in python blocks. There are other solutions I have seen, e.g. polymode, that claim to do this. You might guess that if they worked, I would not be writing this! There was some nice discussion about this idea on the org-mode mailing list, and Nicolas Goaziou pointed out this might be accomplished with the org-font-lock-hook.

You can check out the video here:

It was relatively easy to figure out how to do this. Keymaps can be added to regions during font-lock, so I just had to hook into the org-mode font lock system with a function to find the src blocks and add the keymap as a text-property. That took three steps:

  1. Define the keymaps to use. I use an a-list of (language . map) for this.
  2. Define the font-lock function. This will add the keymap properties to src-blocks.
  3. Define a minor mode to toggle this feature on and off.

Here is the definition of the keymaps. Generally I just copy the mode-map I want and then add some things to them. For example sometimes it is still a good idea to jump into the org-special-edit mode. For example, if you try to use a command in a Python block to send the buffer to the repl while in org-mode you are sure to get an error! You might also want to add the C-c C-e export command if you use that a lot. An alternative approach, of course, is to copy the org-map and add additional bindings to it. The choice is up to you.

(require 'lispy)
(require 'elpy)

(setq scimax-src-block-keymaps
      `(("ipython" . ,(let ((map (make-composed-keymap
                                  `(,elpy-mode-map ,python-mode-map ,pyvenv-mode-map)
                                  org-mode-map)))
                        ;; In org-mode I define RET so we f
                        (define-key map (kbd "<return>") 'newline)
                        (define-key map (kbd "C-c C-c") 'org-ctrl-c-ctrl-c)
                        map))
        ("python" . ,(let ((map (make-composed-keymap
                                 `(,elpy-mode-map ,python-mode-map ,pyvenv-mode-map)
                                 org-mode-map)))
                       ;; In org-mode I define RET so we f
                       (define-key map (kbd "<return>") 'newline)
                       (define-key map (kbd "C-c C-c") 'org-ctrl-c-ctrl-c)
                       map))
        ("emacs-lisp" . ,(let ((map (make-composed-keymap `(,lispy-mode-map
                                                            ,emacs-lisp-mode-map
                                                            ,outline-minor-mode-map)
                                                          org-mode-map)))
                           (define-key map (kbd "C-c C-c") 'org-ctrl-c-ctrl-c)
                           map))))

Next we define the function that will apply the keymap to each src block. The keymaps are only applied when they are defined in the variable above. This function is derived from org-fontify-meta-lines-and-blocks-1.

(defun scimax-add-keymap-to-src-blocks (limit)
  "Add keymaps to src-blocks defined in `scimax-src-block-keymaps'."
  (let ((case-fold-search t)
        lang)
    (while (re-search-forward org-babel-src-block-regexp limit t)
      (let ((lang (match-string 2))
            (beg (match-beginning 0))
            (end (match-end 0)))
        (if (assoc (org-no-properties lang) scimax-src-block-keymaps)
            (progn
              (add-text-properties
               beg end `(local-map ,(cdr (assoc
                                          (org-no-properties lang)
                                          scimax-src-block-keymaps))))
              (add-text-properties
               beg end `(cursor-sensor-functions
                         ((lambda (win prev-pos sym)
                            ;; This simulates a mouse click and makes a menu change
                            (org-mouse-down-mouse nil)))))))))))

Here we create an advice to trick any functions that need to know the major mode. We only apply the spoof if we are in org-mode and in a src block though. Otherwise we call the original function. So far lispy–eval is the only function I have needed it for. This might be a general strategy though to do other things like narrow to the src-block, or even go into special edit mode temporarily if there are commands that require it.

(defun scimax-spoof-mode (orig-func &rest args)
  "Advice function to spoof commands in org-mode src blocks.
It is for commands that depend on the major mode. One example is
`lispy--eval'."
  (if (org-in-src-block-p)
      (let ((major-mode (intern (format "%s-mode" (first (org-babel-get-src-block-info))))))
        (apply orig-func args))
    (apply orig-func args)))

We define a minor mode so we can toggle this on and off. Here we add the function to the org-font-lock-hook and advise the lispy–eval function. I had to add the font-lock-function to the end of the org-font-lock hook for some reason, and also add local-map as an extra-managed property so it would be removed when we toggle it off.

(define-minor-mode scimax-src-keymap-mode
  "Minor mode to add mode keymaps to src-blocks."
  :init-value nil
  (if scimax-src-keymap-mode
      (progn
        (add-hook 'org-font-lock-hook #'scimax-add-keymap-to-src-blocks t)
        (add-to-list 'font-lock-extra-managed-props 'local-map)
        (add-to-list 'font-lock-extra-managed-props 'cursor-sensor-functions)
        (advice-add 'lispy--eval :around 'scimax-spoof-mode)
        (cursor-sensor-mode +1))
    (remove-hook 'org-font-lock-hook #'scimax-add-keymap-to-src-blocks)
    (advice-remove 'lispy--eval 'scimax-spoof-mode)
    (cursor-sensor-mode -1))
  (font-lock-fontify-buffer))

(add-hook 'org-mode-hook (lambda ()
                           (scimax-src-keymap-mode +1)))

That is it! I am pretty sure this is a good idea. It helps a lot when you are writing a lot of short code blocks and near equal amounts of text (like in this blog post). It also helps write the code since many things like indentation, parentheses, etc. are automatically handled. That is what I used to go into special-edit mode all the time for!

I have not used this long enough to know if it causes any other surprises. If you try it and find any, leave a comment!

1 Update

It turns out you can have the best of all the worlds by combining keymaps. The make-composed-keymap creates a new keymap that combines a keymaps and falls through to a parent keymap. So here we use that to combine several keymaps, falling through to org-mode. The only subtlety I have come across is that I remapped <return> in orgmode to scimax/org-return, and not all modes define it, so I redefine it in some places to just be newline. Also to keep C-c C-c for executing the block, I add that back too.

I use a few maps here, and some of them seem to just add menus that are only active when your cursor is in the block. Pretty handy!

(setq scimax-src-block-keymaps
      `(("ipython" . ,(let ((map (make-composed-keymap
                                  `(,elpy-mode-map ,python-mode-map ,pyvenv-mode-map)
                                  org-mode-map)))
                        ;; In org-mode I define RET so we f
                        (define-key map (kbd "<return>") 'newline)
                        (define-key map (kbd "C-c C-c") 'org-ctrl-c-ctrl-c)
                        map))
        ("python" . ,(let ((map (make-composed-keymap
                                 `(,elpy-mode-map ,python-mode-map ,pyvenv-mode-map)
                                 org-mode-map)))
                       ;; In org-mode I define RET so we f
                       (define-key map (kbd "<return>") 'newline)
                       (define-key map (kbd "C-c C-c") 'org-ctrl-c-ctrl-c)
                       map))
        ("emacs-lisp" . ,(let ((map (make-composed-keymap `(,lispy-mode-map
                                                            ,emacs-lisp-mode-map
                                                            ,outline-minor-mode-map)
                                                          org-mode-map)))
                           (define-key map (kbd "C-c C-c") 'org-ctrl-c-ctrl-c)
                           map))))

2 Update #2

The previous version had some issues where it would only add a keymap to the first block. The code in this post now addresses that and uses cursor-sensor-functions to make sure we change key map on entering and leaving blocks. That might mean you need an emacs of at least version 25 to use this. I guess it will work with an earlier version, but the cursor-sensor-functions might get ignored. You might have to comment out the cursor-sensor-mode line

Thanks to those brave people alpha-testing this and helping refine the idea!

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

org-mode source

Org-mode version = 9.0.7

Read and Post Comments

Org-mode and ipython enhancements in scimax

| categories: emacs, orgmode, ipython | tags: | View Comments

We have made some improvements to using Ipython in org-mode in the past including:

  1. Inline figures
  2. Export to Jupyter notebooks

Today I will talk about a few new features and improvements I have introduced to scimax for using org-mode and Ipython together.

The video for this post might be more obvious than the post:

1 Some convenience functions

There are a few nice shortcuts in the Jupyter notebook. Now we have some convenient commands in scimax to mimic those. My favorites are adding cells above or below the current cell. You can insert a new src block above the current one with (M-x org-babel-insert-block). You can use a prefix arg to insert it below the current block.

# code
# below
# some code

I am particularly fond of splitting a large block into two smaller blocks. Use (M-x org-babel-split-src-block) to do that and leave the point in the upper block. Use a prefix arg to leave the point in the lower block.

# lots of code in large block
# Even more code
# The end of the long block

You can execute all the blocks up to the current point with (M-x org-babel-execute-to-point).

2 ob-ipython-inspect works

In the original ob-ipython I found that ob-ipython-inspect did not work unless you were in special edit mode. That is too inconvenient. I modified a few functions to work directly from the org-buffer. I bind this to M-. in org-mode.

%matplotlib inline
import numpy as np

import matplotlib.pyplot as plt

# Compute areas and colors
N = 150
r = 2 * np.random.rand(N)
theta = 2 * np.pi * np.random.rand(N)
area = 200 * r**2
colors = theta

ax = plt.subplot(111, projection='polar')
c = ax.scatter(theta, r, c=colors, s=area, cmap='hsv', alpha=0.75)

<matplotlib.figure.Figure at 0x114ded710>

3 Getting selective output from Ipython

Out of the box Ipython returns a lot of results. This block, for example returns a plain text, image and latex result as output.

from sympy import *
# commenting out init_printing() results in no output
init_printing()

var('x y')
x**2 + y

2 x + y

We can select which one we want with a new header argument :ob-ipython-results. For this block you can give it the value of text/plain, text/latex or image/png.

var('x y')
x**2 + y

2 x + y

Or to get the image:

var('x y')
x**2 + y

This shows up with pandas too. This block creates a table of data and then shows the first 5 rows. Ipython returns both plain text and html here.

import pandas as pd
import numpy as np
import datetime as dt

def makeSim(nHosps, nPatients):
    df = pd.DataFrame()
    df['patientid'] = range(nPatients)
    df['hospid'] = np.random.randint(0, nHosps, nPatients)
    df['sex'] = np.random.randint(0, 2, nPatients)
    df['age'] = np.random.normal(65,18, nPatients)
    df['race'] = np.random.randint(0, 4, nPatients)
    df['cptCode'] = np.random.randint(1, 100, nPatients)
    df['rdm30d'] = np.random.uniform(0, 1, nPatients) < 0.1
    df['mort30d'] = np.random.uniform(0, 1, nPatients) < 0.2
    df['los'] = np.random.normal(8, 2, nPatients)
    return df

discharges = makeSim(50, 10000)
discharges.head()

patientid hospid sex age race cptCode rdm30d mort30d los 0 0 10 1 64.311947 0 8 False False 8.036793 1 1 6 0 82.951484 1 73 True False 7.996024 2 2 27 1 53.064501 3 95 False False 9.015144 3 3 37 0 64.799128 0 93 False False 10.099032 4 4 46 0 99.111394 2 25 False False 11.711427

patientid hospid sex age race cptCode rdm30d mort30d los
0 0 10 1 64.311947 0 8 False False 8.036793
1 1 6 0 82.951484 1 73 True False 7.996024
2 2 27 1 53.064501 3 95 False False 9.015144
3 3 37 0 64.799128 0 93 False False 10.099032
4 4 46 0 99.111394 2 25 False False 11.711427

We can use the header to select only the plain text output!

import pandas as pd
import numpy as np
import datetime as dt

def makeSim(nHosps, nPatients):
    df = pd.DataFrame()
    df['patientid'] = range(nPatients)
    df['hospid'] = np.random.randint(0, nHosps, nPatients)
    df['sex'] = np.random.randint(0, 2, nPatients)
    df['age'] = np.random.normal(65,18, nPatients)
    df['race'] = np.random.randint(0, 4, nPatients)
    df['cptCode'] = np.random.randint(1, 100, nPatients)
    df['rdm30d'] = np.random.uniform(0, 1, nPatients) < 0.1
    df['mort30d'] = np.random.uniform(0, 1, nPatients) < 0.2
    df['los'] = np.random.normal(8, 2, nPatients)
    return df

discharges = makeSim(50, 10000)
discharges.head()

patientid hospid sex age race cptCode rdm30d mort30d los 0 0 21 0 73.633836 1 38 False False 7.144019 1 1 16 1 67.518804 3 23 False False 3.340534 2 2 15 0 44.139033 0 8 False False 9.258706 3 3 29 1 45.510276 2 5 False False 10.590245 4 4 7 0 52.974924 2 4 False True 5.811064

4 Where was that error?

A somewhat annoying feature of running cells in org-mode is when there is an exception there has not been a good way to jump to the line that caused the error to edit it. The lines in the src block are not numbered, so in a large block it can be tedious to find the line. In scimax, when you get an exception it will number the lines in the src block, and when you press q in the exception traceback buffer it will jump to the line in the block where the error occurred.

print(1)
#raise Exception('Here')
print(2)

1 2

If you don't like the numbers add this to your init file:

(setq ob-ipython-number-on-exception nil)

5 Asynchronous Ipython

I have made a few improvements to the asynchronous workflow in Ipython. We now have a calculation queue, so you can use C-c C-c to execute several blocks in a row, and they will run asynchronously in the order you ran them. While they are running you can continue using Emacs, e.g. writing that paper, reading email, checking RSS feeds, tetris, … This also lets you run all the blocks up to the current point (M-x org-babel-execute-ipython-buffer-to-point-async) or the whole buffer (of Ipython) blocks asynchronously (M-x org-babel-execute-ipython-buffer-async).

To turn this on by default put this in your init file:

(setq org-babel-async-ipython t)

This requires all src blocks to have a name, and running the block will give it a name if you have not named the block. By default we use human-readable names. While the block is running, there will be a link indicating it is running. You can click on the link to cancel it. Running subsequent blocks will queue them to be run when the first block is done.

Here is an example:

import time
time.sleep(5)
a = 5
print('done')
print(3 * a)

15

Occasionally you will run into an issue. You can clear the queue with org-babel-async-ipython-clear-queue.

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

org-mode source

Org-mode version = 9.0.5

Read and Post Comments

A partial symbolic numeric solver in emacs-lisp

| categories: emacs-lisp, emacs, math | tags: | View Comments

I have been exploring ways to use emacs-lisp to express scientific ideas. In this post, we explore a partial symbolic numeric solver in Emacs-lisp. This involves some syntactic developments to more clearly identify something we want to solve for and to then generate the code required to solve it.

In section The Newton solver you can find a simple implementation of a Newton solver in emacs-lisp. This function allows you to numerically solve equations that can be written in the form \(f(x) = 0\) for \(x\) given an initial guess. You write a function for \(f(x)\) and pass the function to the solver. This is a standard approach used in Python with fsolve, for example. Here is an example of solving a trivial problem: \(x - 4 = 0\) just to check that it works. We use a lambda function for \(f(x) = x - 4 = 0\). The answer is \(x=4\).

(newton-f (lambda (x) (- x 4)) 2)

That syntax is not too bad, but we have the whole lambda expression in there, and some repetition of what we want to solve for as an argument and in the function. It would be interesting if we could just have an expression that gets solved, e.g. (newton-f (- x? 4) 2) where x? indicates the thing to solve for.

We can do that! We can take an expression, flatten it and find the variable names that end with ?. We should check that there is only one, but for now we don't. Here is an example that does that. I use a nested expression here just to illustrate that the code finds the x? variable correctly.

(require 'dash)

(let ((body '((* (- x? 4) 1))))
  (loop for item in (-flatten body)
        if (and (symbolp item) (s-ends-with? "?" (symbol-name item)))
        collect item))

So, given an expression we can identify the unknown that should be the argument to a lambda function. So, we create a macro that takes that expression and constructs a function to solve it, then calls newton-f on it. The macro is syntactically useful here because we do not have to quote the expression. Here is that macro.

(defmacro solve (expression guess)
  `(newton-f
    (lambda ,(loop for item in (-flatten expression)
                   if (and (symbolp item) (s-ends-with? "?" (symbol-name item)))
                   collect item)
      ,expression)
    ,guess))

I call this a partial symbolic solver because we do some introspection symbolically to identify what to solve for, and then construct the code required to solve it. Here is that trivial example (x? - 4 = 0). It just shows we can have some nesting and it still works. I am not so thrilled with the initial guess, but this is an iterative solver, so you either need an initial guess, or a solution range.

(solve (* (- x? 4) 1) 3)

Here is what that expands into:

(macroexpand '(solve (* (- x? 4) 1) 3))

It expands into what we would have written in the first place. The benefit to us is less typing, and a simpler syntax. Both of those reduce the opportunity to make errors!

A more realistic problem might be: Reactant A flows into a continuously stirred tank reactor at a rate of \(F_{A0} = 1\) mol/min with a volumetric flow of \(v_0 = 1\) L/min.. The reactor achieves 50% conversion (\(X\)) of A to products. The reaction rate law is known to be \(-r_A = k C_A\) with \(k = 0.1\) 1/min. Estimate the volume of the reactor. If you have taken my class in reaction engineering, you know the following facts:

  • The exit molar flow is defined by \(F_A = F_{A0} (1 - X)\)
  • The exit concentration is \(C_A = F_A / v_0\)
  • The mole balance is defined by \(0 = F_{A0} - F_A + r_A V\)

That is all we need; we can solve for \(V\) from the last equation. This is simple enough you might do the algebra to get: \(V = \frac{F_{A0} - F_A}{-r_A}\) which can be simply evaluated. We use our solver here and compare it to the evaluation.

Here is the solver:

(let* ((Fa0 1)
       (X 0.5)
       (Fa (* Fa0 (- 1 X)))
       (k 0.1)
       (v0 1)
       (Ca (/ Fa v0))
       (r (* k Ca))
       (ra (* r -1)))
  (solve (+ Fa0 (* Fa -1) (* ra V?)) 2))

It is pretty hard to imagine doing something like this in Python! It would probably involve parsing a string.

Here is the evaluation from our algebra:

(let* ((Fa0 1)
       (X 0.5)
       (Fa (* Fa0 (- 1 X)))
       (k 0.1)
       (v0 1)
       (Ca (/ Fa v0))
       (r (* k Ca))
       (ra (* r -1)))
  (/ (- Fa0 Fa) (* -1 ra)))

Within the tolerance specified in newton-f, these are the same.

This is just the tip of the iceberg though. You may have noticed that none of the variables in the let* had any descriptions. Sure, you could put some comments after them, but those are not really part of the code.

Also, we had to define the variables in advance of the expression. That is a limitation of how computers work, they cannot evaluate undefined variables. It constrains how we can express the idea. What if we could instead specify the equation first, then the data? That way we are clear what we are trying to do at a higher level, and fill in the details later. Suppose we wanted a syntax like the block below instead. Here we emphasize the equation we are solving first, and then define the variables and quantities used in the equation, and finally the guess that we use to find the solution.

(solve1
 (eqn (+ Fa0 (* -1 Fa) (* ra V?)))
 (data ((k 0.1 "rate constant 1/min")
        (Ca0 1.0 "feed concentration")
        (v0 1 "volumetric flow L/min")
        (Fa0 (* v0 Ca0) "Inlet molar flow")
        (X 0.5 "Desired conversion")
        (Fa (* Fa0 (- 1 X)) "Exit molar flow")
        (Ca (/ Fa v0) "exit concentration")
        (ra (* -1 k Ca) "rate in the reactor")))
 (guess 8))

That is achievable with the solve1 macro below! It too has some limitations, mostly the order of the data block still has to be correct, e.g. you cannot use a variable before it is defined. It would take some serious macro-fu to sort these so that everything is defined in the right order! Still, it allows you to express an executable idea in the order we defined. The strings in this syntax for documenting the variables are ignored, but they could be used in the macro to print useful information or something else you could imagine. You could also make them mandatory to encourage documentation.

(defmacro solve1 (&rest body)
  (let ((expression (second (assoc 'eqn body)))
        (data (loop for d in (second (assoc 'data body))
                    collect (list (first d) (second d))))
        (guess (second (assoc 'guess body))))
    `(let* ,data
       (newton-f
        (lambda ,(loop for item in (-flatten expression)
                       if (and (symbolp item) (s-ends-with? "?" (symbol-name item)))
                       collect item)
          ,expression)
        ,guess))))

To summarize, lisp macros allow us to rewrite the syntax of code before it is evaluated. This gives us the opportunity to inspect it, and generate new code, e.g. functions with arguments based on the contents of expressions, to save us typing. It also allows us to define ideas in a near arbitrary order that make sense to us, and then rearrange them so they make sense to the computer. See, for example, this post for an example of changing how functions are defined.

This seems to be heading in the domain specific language direction. I think it would be very helpful in engineering problem solving to build up tools like this. They could start out simple for new students to use. They never need to see the macro parts of this, just to learn how to use them for problem solving. These beginner tools would be limited in what they could do to minimize how much lisp is required to be learned so students can focus on the problem solving. Eventually they might outgrow them, and in the process transition to having the full lisp language at their disposal for problem solving.

1 The Newton solver in emacs-lisp

This is an emacs-lisp implementation of Newton's method. It is a simple implementation for a single variable. The tolerance and step-size are hard-coded for this post because we focus on the partial symbolic solver, not the best solver methods.

;; See https://en.wikipedia.org/wiki/Newton%27s_method for the method

(defun newton-f (func x0)
  "Solve the equation FUNC(x)=0 using Newton's method.
X0 is an initial guess."
  (let* ((tolerance 1e-6)
         (x x0)
         (dx 1e-6)
         fx fpx)
    (while (> (abs (funcall func x)) tolerance)
      (setq fx (funcall func x)
            fpx (/ (- (funcall func (+ x dx)) (funcall func (- x dx))) (* 2 dx))
            x (- x (/ fx fpx))))
    x))

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

org-mode source

Org-mode version = 9.0.5

Read and Post Comments

A new and improved Emacs gnuplot DSL

| categories: emacs, plotting, lisp | tags: | View Comments

A significant limitation of the previous DSL I wrote is that all the plotting commands have to go in one macro. It would be nice to accumulate them in different forms, and when you want to run them all. A classic way to do that in Emacs lisp is to make a global variable, e.g. *gnuplot-cmds* and append commands to it. Then when you want to, run the commands.

A more modern approach is to use a closure to encapsulate the commands. Here is a "let over lambda" that defines a few functions that encapsulate an enclosed variable gnuplot-commands. We define one function to add commands to the list of commands, one to clear the commands, one to generate the gnuplot script as a string, and one to run the program. The enclosed variable gnuplot-commands is basically only accessible by these functions. It is encapsulated, similar to if we defined a class in Python then made an instance of it with an attribute that was accessible only be instance methods. On one hand, this "protects" the variable, and keeps it out of the global namespace. On the other hand, we lose the documentation that would have come with a defvar, and we have to define a function to access the contents of that variable.

(let ((gnuplot-commands '("set terminal qt")))

  (defun gnuplot-add-cmd (s)
    "Append the command S to gnuplot-cmds."
    (setq gnuplot-commands (append gnuplot-commands (list s))))

  (defun gnuplot-clear ()
    (setq gnuplot-commands '("set terminal qt")))

  (defun gnuplot-script ()
    (s-join "\n" gnuplot-commands)))
gnuplot-script

To run the commands, we define this function. It does not need to be in the closure because it only accesses the commands through functions we defined in the closure.

(defun gnuplot-show ()
    (let* ((temporary-file-directory ".")
           (cmdfile (make-temp-file "gnuplot-cmds-" nil ".gpl"))
           (shellcmd (format "gnuplot --persist -c \"%s\"" cmdfile))
           (cmds (gnuplot-script)))
      (with-temp-file cmdfile
        (insert cmds))
      (shell-command shellcmd)
      (delete-file cmdfile)
      cmds))
gnuplot-show

Last time I noted I had a new idea for the DSL syntax that would give us more flexibility to inject variables and code into the DSL. The idea is to use keywords, symbols that start with :, to indicate they should be replaced by the value of the non-keyword symbol in the environment, and for any form that starts with : to evaluate that form. So, (: - 5 4) would get replaced by 1. Here is the new macro for that.

(defmacro kargs (&rest args)
  "Convert symbols to strings, quote strings, and (expr) to what they evaluate to."
  `(s-join " " (list ,@(cl-mapcan
                        (lambda (s)
                          (list
                           (cond
                            ((keywordp s)
                             (format "%s"
                                     (symbol-value (intern (substring (symbol-name s) 1)))))
                            ((symbolp s)
                             (symbol-name s))
                            ((stringp s)
                             (format "\"%s\"" s))
                            ((and (listp s) (eq : (car s)))
                             `(with-output-to-string
                                (princ ,(cdr s))))
                            (t
                             (format "%s" s)))))
                        args))))
kargs

Now, our gnuplot macro is simpler, since all it does is add commands to the list. If the form is a string, we add it as is, if the form starts with (: stuff) we evaluate the cdr of the form, and otherwise, we pass the form contents to the kargs macro for processing.

(defmacro gnuplot (&rest forms)
  `(loop for s in (list ,@(mapcar (lambda (x)
                                    (cond
                                     ((stringp x)
                                      x)
                                     ((and (listp x) (eq : (car x)))
                                      `,(cdr x))
                                     (t
                                      `(kargs ,@x))))
                                  forms))
         do (gnuplot-add-cmd s)))
gnuplot

What did that gain us? First, we can break up a script so we can talk about it, maybe do some calculations, etc… Let's look at the example at http://gnuplot.sourceforge.net/demo/linkedaxes.html.

We can start with the basic settings.

(gnuplot-clear)

(gnuplot
 (set terminal png)
 (set output "linkedaxes.png")
 (set encoding utf8)
 (set key outside Left)
 (set bmargin 5)
 (set tmargin 6)
 (set style data lines)
 (set tics in)
 (set ticslevel 0.5)
 (set xlabel  "X-ray energy in eV")

 (set format y  \'%5.1fe\')
 (set title " Anomalous scattering factors ")
 (set xrange  [9000:14400])
 (set offset 0\,0\,1.0\,0)
 (set xtics nomirror)
 (set link x via 12398./x inverse 12398./x)

 (set x2label  "X-ray wavelength in Å")
 (set x2tics 0.1  format "%.1f Å" nomirror))

We need to download some data files. We can do that, and add another line to the gnuplot script. The escaping on the quotes and commas is especially tedious in this one ;) but, we don't need those pesky line-continuations here.

(shell-command "wget http://www.bmsc.washington.edu/scatter/data/Br.dat")
(shell-command "wget http://www.bmsc.washington.edu/scatter/data/Ta.dat")


(gnuplot
 (plot "Br.dat" volatile using 1:3 title \'Br f\"\'  lt 1 lw 3\, \'\' volatile using 1:2 title "Br f'"  lt 1 lw 1\,
       "Ta.dat" volatile using 1:3 title \'Ta f\"\' lt 2 lw 3\, \'\' volatile using 1:2 title \"Ta f\'\"  lt 2 lw 1))

(gnuplot-script)
set terminal qt
set terminal png
set output "linkedaxes.png"
set encoding utf8
set key outside Left
set bmargin 5
set tmargin 6
set style data lines
set tics in
set ticslevel 0.5
set xlabel "X-ray energy in eV"
set format y '%5.1fe'
set title " Anomalous scattering factors "
set xrange [9000:14400]
set offset 0,0,1.0,0
set xtics nomirror
set link x via 12398./x inverse 12398./x
set x2label "X-ray wavelength in Å"
set x2tics 0.1 format "%.1f Å" nomirror
plot "Br.dat" volatile using 1:3 title 'Br f"' lt 1 lw 3, '' volatile using 1:2 title "Br f'" lt 1 lw 1, "Ta.dat" volatile using 1:3 title 'Ta f"' lt 2 lw 3, '' volatile using 1:2 title "Ta f'" lt 2 lw 1

Finally, we can set the output to png, and run our program.

(gnuplot-show)

Looks good.

What about the fancy keyword formatting? Here is an example of that in action. :term gets replaced by the term variable, :png gets replaced by the filename, and :x is replaced by 4.

(gnuplot-clear)
(let ((x 4)
      (term "png")
      (png "\"polar.png\""))
  (gnuplot
   (set terminal :term)
   (set output :png)
   (set polar)
   (set dummy t)
   (plot sin\( :x *t\) \,cos\( :x *t\))
   (set offset 0\,0\,0\,0)))

(gnuplot-show)
set terminal qt
set terminal png
set output "polar.png"
set polar
set dummy t
plot sin( 4 *t) ,cos( 4 *t)
set offset 0,0,0,0

There are a few nuances I didn't expect. First, you have to escape the parentheses in this case because otherwise it looks like a form that will be ignored. Second, you have to quote the string to get quotes into the gnuplot script. Third, there has to be a space before and after the keywords for emacs to parse it correctly and do the substitution.

Let's look at one last example that uses the (: form). We reproduce a figure from http://gnuplot.sourceforge.net/demo/transparent_solids.html here.

(gnuplot-clear)
(gnuplot
 (set terminal pngcairo  background "#ffffff" enhanced font "arial,9" fontscale 1.0 size 512\, 384 )
 (set output "transparent-solids.png")
 ;; construct the title
 (set title (: format "\"%s\"" (concat "Interlocking Tori - PM3D surface" "with depth sorting and transparency")))

 ;; use lisp code to create a gnuplot command
 (: concat "unset" " " "border")

 (unset key)
 (set object 1 rect from screen 0\, 0\, 0 to screen 1\, 1\, 0 behind)
 (set object 1 rect fc  rgb \"gray\"  fillstyle solid 1.0  border -1)
 (set view 64\, 345\, 1.24375\, 0.995902)
 (set isosamples 50\, 20)
 (unset xtics)
 (unset ytics)
 (unset ztics)
 (set dummy u\,v)
 (set parametric)
 (set urange [ -pi : pi ])
 (set vrange [ -pi : pi ])

 (set style fill  transparent solid 0.30 border)
 (set pm3d depthorder)
 (set palette rgbformulae 8\, 9\, 7)
 (set pm3d interpolate 1\,1 flush begin noftriangles border lt black linewidth 0.500 dashtype solid corners2color mean)
 (set colorbox vertical origin screen 0.9\, 0.2\, 0 size screen 0.05\, 0.6\, 0 front  noinvert bdefault)

 (splot (: concat "cos(u)+.5*cos(u)*cos(v),sin(u)+.5*sin(u)*cos(v),.5*sin(v) with pm3d,"
           "1+cos(u)+.5*cos(u)*cos(v),.5*sin(v),sin(u)+.5*sin(u)*cos(v) with pm3d")))
(gnuplot-show)
set terminal qt
set terminal pngcairo background "#ffffff" enhanced font "arial,9" fontscale 1.0 size 512, 384
set output "transparent-solids.png"
set title "Interlocking Tori - PM3D surfacewith depth sorting and transparency"
unset border
unset key
set object 1 rect from screen 0, 0, 0 to screen 1, 1, 0 behind
set object 1 rect fc rgb "gray" fillstyle solid 1.0 border -1
set view 64, 345, 1.24375, 0.995902
set isosamples 50, 20
unset xtics
unset ytics
unset ztics
set dummy u,v
set parametric
set urange [-pi : pi]
set vrange [-pi : pi]
set style fill transparent solid 0.3 border
set pm3d depthorder
set palette rgbformulae 8, 9, 7
set pm3d interpolate 1,1 flush begin noftriangles border lt black linewidth 0.5 dashtype solid corners2color mean
set colorbox vertical origin screen 0.9, 0.2, 0 size screen 0.05, 0.6, 0 front noinvert bdefault
splot cos(u)+.5*cos(u)*cos(v),sin(u)+.5*sin(u)*cos(v),.5*sin(v) with pm3d,1+cos(u)+.5*cos(u)*cos(v),.5*sin(v),sin(u)+.5*sin(u)*cos(v) with pm3d

Overall this seems like an improvement to the DSL. I didn't invent the idea of reusing keywords this way out of the blue. In On Lisp, Paul graham uses "special" variable names in Chapter 18, where he shows how to use gensyms for special purposes, and also variables with special names like ?x. Even Emacs is using a variation of this idea. Check out this new let-alist macro:

(let-alist '((x . 5))
  (+ 1 .x))
6

There is a special variable inside the body that is a dot-name. The macro expands to provide a value for that symbol. I wonder if I should have tried to use an approach like this instead. Maybe another day. After I read and study the four defuns and single defmacro that make this possible!

You can see here what happens:

(macroexpand '(let-alist '((x . 5))
  (+ 1 .x)))
(let
    ((alist
      '((x . 5))))
  (let ((\.x (cdr (assq 'x alist))))
    (+ 1 \.x)))

The macro builds up an internal alist for the dot-names.

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

org-mode source

Org-mode version = 9.1.6

Read and Post Comments

« Previous Page -- Next Page »