A partial symbolic numeric solver in emacs-lisp

| categories: emacs, emacs-lisp, 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

Collecting entries from files in a directory

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

I am running a class where students will be generating files that contain their answers. I want to quickly get a list of counts of all the answers. For example, six students will create a file called animal.dat in a directory called example/<studentid>, and that file will contain their favorite animal. I made some example files to test this idea out. Here are the contents.

cat example/*/animal.dat
dog
cat
dog
bird
dog
bird

You can see there are three dogs, two birds and a cat. I want code to do this counting, because in my real application there will be 58 of these files, and lots of times I need to aggregate them. Let us start with a simple example that counts the elements in a list.

(let ((animals '(dog cat dog bird dog bird))
      (counts '())
      place)
  (dolist (animal animals)
    (setq place (assoc animal counts))
    (message "place = %s" place)
    (if place
        (setf (cdr place) (+ 1 (cdr place)))
      (setq counts (cons `(,animal . 1) counts))))
  counts)

((bird . 2) (cat . 1) (dog . 3))

Let us turn that into a function.

(defun counts (list)
  (let ((counts '())
        place)
    (dolist (el list)   
      (setq place (assoc el  counts))
    (message "place = %s" place)
    (if place
        (setf (cdr place) (+ 1 (cdr place)))
      (setq counts (cons `(,el . 1) counts))))
    counts))

(counts '(dog cat dog bird dog bird))

((bird . 2) (cat . 1) (dog . 3))

Nice. Now we need a simple way to get that list. We need to glob the files to find them, then open them and read the value. Here is a way to get the files.

(f-entries "example"
           (lambda (f)
             (string= (file-name-nondirectory f) "animal.dat"))
           t)
/Users/jkitchin/blogofile-jkitchin.github.com/blog/collect-entries/example/s1/animal.dat /Users/jkitchin/blogofile-jkitchin.github.com/blog/collect-entries/example/s2/animal.dat /Users/jkitchin/blogofile-jkitchin.github.com/blog/collect-entries/example/s3/animal.dat /Users/jkitchin/blogofile-jkitchin.github.com/blog/collect-entries/example/s4/animal.dat /Users/jkitchin/blogofile-jkitchin.github.com/blog/collect-entries/example/s5/animal.dat /Users/jkitchin/blogofile-jkitchin.github.com/blog/collect-entries/example/s6/animal.dat

Now we just need to run a mapcar over this list of files.

(mapcar
 (lambda (f)
   (with-temp-buffer
     (insert-file-contents f)
     (s-trim (buffer-string))))
 (f-entries "example"
            (lambda (f)
              (string= (file-name-nondirectory f) "animal.dat"))
            t))
dog cat dog bird dog bird

Finally, putting this together, we have some code that maps over the files, and counts the entries.

(counts
 (mapcar
  (lambda (f)
    (with-temp-buffer
      (insert-file-contents f)
      (s-trim (buffer-string))))
  (f-entries "example"
             (lambda (f)
               (string= (file-name-nondirectory f) "animal.dat"))
             t)))

((bird . 2) (cat . 1) (dog . 3))

This will be helpful in dealing with 58 entries during my class!

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

org-mode source

Org-mode version = 8.2.7c

Read and Post Comments

Autogenerating functions in emacs-lisp

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

I have a need to generate a lot of similar functions, and I do not want to cut and paste the code. I want to generate the functions with code. This seems to be what macros are for in emacs lisp.

As a prototype example, we will make functions that raise a number to a power. We want functions like power-3 and power-4 that raise numbers to the third and fourth powers. We will define functions like this for the numbers 0-9.

Here we define the macro. i do not want to get into the nitty gritty details of macro definitions here.

(defmacro make-power-n (n)
 `(defun ,(intern (format "power-%s" n)) (arg) (expt arg ,n)))

(make-power-n 4)

(power-4 4)
256

Now we use the macro and mapcar on it onto a list of numbers. We have to eval the macro in the mapcar lambda function.

(defmacro make-power-n (n)
 `(defun ,(intern (format "power-%s" n)) (arg) (expt arg ,n)))

(mapcar (lambda (x) (eval `(make-power-n ,x))) '(0 1 2 3 4 5 6 7 8 9))
 
;; example of a few functions
(list (power-0 3) (power-1 3) (power-2 3))
1 3 9

It works! We created 10 functions in a little bit of code.

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

org-mode source

Org-mode version = 8.2.6

Read and Post Comments

Language specific default headers for code blocks in org-mode

| categories: org-mode, emacs-lisp | tags: | View Comments

I use code blocks in org-mode a lot. I usually code in Python, and in Python I usually write code that prints output which I want to see. So I almost always want the code blocks to return the output, and not the value of the last function. I have set my default header args like this:

org-babel-default-header-args
(:exports . both) (:results . replace output) (:session . none) (:cache . no) (:noweb . no) (:hlines . no) (:tangle . no)

However, I would prefer that emacs-lisp blocks use value for the results. I know I can get that by putting :results value in the code block header, but that annoys me. I learned today from http://orgmode.org/worg/org-contrib/babel/header-args.html that you can make language specific default headers!

This code in my init file sets emacs-lisp specific default headers:

(setq org-babel-default-header-args:emacs-lisp 
      (cons '(:results . "value")
            (assq-delete-all :results org-babel-default-header-args)))

That way I do not have type :results value at the top of every elisp block. Of course, if I want the output I could specify :results output in the block.

org-babel-default-header-args:emacs-lisp
(:results . value) (:exports . both) (:session . none) (:cache . no) (:noweb . no) (:hlines . no) (:tangle . no)

Problem solved!

On a related note, I find I write so many blocks of python and elisp I added these templates:

;; add <p for python expansion
(add-to-list 'org-structure-template-alist
             '("p" "#+BEGIN_SRC python\n?\n#+END_SRC" "<src lang=\"python\">\n?\n</src>"))

;; add <el for emacs-lisp expansion
(add-to-list 'org-structure-template-alist
             '("el" "#+BEGIN_SRC emacs-lisp\n?\n#+END_SRC" "<src lang=\"emacs-lisp\">\n?\n</src>"))

I probably could have also coded the :results header into those too. They add a tidbit of convenience so I do not have to type python or emacs-lisp after expanding a source block with <s.

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

org-mode source

Org-mode version = 8.2.5g

Read and Post Comments

Another alternative to string templates

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

In the last post I explored a way to expand a string template that was more readable than the usual format. Today I look at another approach where I use sexp expansions to accomplish the same thing. The idea is to embed lisp expressions and replace them by what they evaluate to.

In emacs-lisp, if we have a command in a string, we can "read" it, and then eval it.

Here we get the user-full-name:

(eval (read "user-full-name"))
John Kitchin

We can use this on variables too.

(setq some-variable "test")
(eval (read "some-variable"))
test

So, if we use a syntax to identify what to replace, we can substitute in the values. Let us try %() as the syntax.

(defun expand-template (s)
  "expand a template containing %() with the eval of its contents"
  (replace-regexp-in-string "%(\\([^)]+\\))"
                            (lambda (arg)
                              (format "%s" (eval (read (substring arg 2 -1))))) s))


(let ((key "kitchin-2014")
      (author "Kitchin, J. R.")
      (journal "HACS")
      (year "2014")
      (volume "1")
      (pages "1--10")
      (doi "10.1.1.109/hacs.1.10")
      (url "http://hacs.org/10.1.1.109/hacs.1.10")
      (pdf-dir "/home/jkitchin/pdfs")
      (template "
 :PROPERTIES:
  :Custom_ID: %(key)
  :AUTHOR: %(author
  :JOURNAL: %(journal)
  :YEAR: %(year)
  :VOLUME: %(volume)
  :PAGES: %(pages)
  :DOI: %(doi)
  :URL: %(url)
 :END:
[[cite:%(key)]] [[file:%(pdf-dir)/%(key).pdf][pdf]]\n\n"))

(expand-template template))
 :PROPERTIES:
  :Custom_ID: kitchin-2014
  :AUTHOR: Kitchin, J. R.
  :YEAR: 2014
  :VOLUME: 1
  :PAGES: 1--10
  :DOI: 10.1.1.109/hacs.1.10
  :URL: http://hacs.org/10.1.1.109/hacs.1.10
 :END:
[[cite:kitchin-2014]] [[file:/home/jkitchin/pdfs/kitchin-2014.pdf][pdf]]

That is pretty nice. I like it better than the plist expansion I used before. Presumably these variables would already be defined somewhere in your code.

I thought of trying this on a more complex expansion, and discovered a weakness in the regexp that finds the expansion values. It turns out to be simpler to use %{} as the delimiter than %(), because you may want nested parentheses. The regexp above does not correctly match sets of parentheses.

(defun expand-template (s)
  "expand a template containing %{} with the eval of its contents"
  (replace-regexp-in-string "%{\\([^}]+\\)}"
                            (lambda (arg)
                              (let ((sexp (substring arg 2 -1)))
                                (format "%s" (eval (read sexp))))) s))

(expand-template "2 * 2 = %{(* 2 2)}")
2 * 2 = 4

I am not sure this is a desirable way to make a template, with multiline code to be expanded, but at least this works!

(defun expand-template (s)
  "expand a template containing %{} with the eval of its contents"
  (replace-regexp-in-string "%{\\([^}]+\\)}"
                            (lambda (arg)
                              (let ((sexp (substring arg 2 -1)))
                                (format "%s" (eval (read sexp))))) s))

(expand-template "The result is %{(progn
  (if (> 4 3)
      'true
    'false))}")
The result is true

The regexp used in the expansion is not very robust. In particular if there is a } in the code, it will probably fail because the regexp does not match closing } correctly. Fixing that is beyond me right now!

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

org-mode source

Org-mode version = 8.2.5g

Read and Post Comments

Next Page ยป