The Kitchin Research Group

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

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!

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

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

org-mode source

Org-mode version = 9.0.7

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:

# code

# below

# some code

# lots of code in large block

# Even more code

# The end of the long block

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

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

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

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

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

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

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

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

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

(setq org-babel-async-ipython t)

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

print(3 * a)

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

org-mode source

Org-mode version = 9.0.5

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.

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