Equation of a plane through three points

| categories: python | tags:

We are given three points, and we seek the equation of the plane that goes through them. The method is straight forward. A plane is defined by the equation:

\(a x + b y + c z = d\)

and we just need the coefficients. The \(a, b, c\) coefficients are obtained from a vector normal to the plane, and \(d\) is calculated separately. We get the normal vector from the cross-product of two vectors connecting the points, and we get \(d\) from the dot product of the normal vector with any one of the point position vectors.

Finally, given the equation, we want to generate a mesh that samples the plane, and plot the mesh and original points to verify the plane goes through the points. Here is the implementation.

import numpy as np

p1 = np.array([1, 2, 3])
p2 = np.array([4, 6, 9])
p3 = np.array([12, 11, 9])

# These two vectors are in the plane
v1 = p3 - p1
v2 = p2 - p1

# the cross product is a vector normal to the plane
cp = np.cross(v1, v2)
a, b, c = cp

# This evaluates a * x3 + b * y3 + c * z3 which equals d
d = np.dot(cp, p3)

print('The equation is {0}x + {1}y + {2}z = {3}'.format(a, b, c, d))

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

x = np.linspace(-2, 14, 5)
y = np.linspace(-2, 14, 5)
X, Y = np.meshgrid(x, y)

Z = (d - a * X - b * Y) / c

# plot the mesh. Each array is 2D, so we flatten them to 1D arrays
ax.plot(X.flatten(),
        Y.flatten(),
        Z.flatten(), 'bo ')

# plot the original points. We use zip to get 1D lists of x, y and z
# coordinates.
ax.plot(*zip(p1, p2, p3), color='r', linestyle=' ', marker='o')

# adjust the view so we can see the point/plane alignment
ax.view_init(0, 22)
plt.tight_layout()
plt.savefig('images/plane.png')
plt.show()
The equation is 30x + -48y + 17z = -15

It looks like the blue points form a plane that contains the red points.

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

org-mode source

Org-mode version = 8.2.10

Discuss on Twitter

Building atomic clusters in ase

| categories: ase, python | tags:

I was perusing the ase codebase, and came across the cluster module . This does not seem to be documented in the main docs, so here are some examples of using it. The module provides some functions to make atomic clusters for simulations.

Below I show some typical usage. First, we look at an icosahedron with three shells.

from ase.cluster.icosahedron import Icosahedron
from ase.io import write

atoms = Icosahedron('Au', noshells=3)
print atoms

write('images/Au-icosahedron-3.png', atoms)
Atoms(symbols='Au55', positions=..., tags=..., cell=[9.816495585723144, 9.816495585723144, 9.816495585723144], pbc=[False, False, False])

Even with only three shells, there are already 55 atoms in this cluster!

How about a decahedron? There are more parameters to set here. I am not sure what the depth of the Marks re-entrance is.

from ase.cluster.decahedron import Decahedron
from ase.io import write

atoms = Decahedron('Pt',
                   p=5,  # natoms on 100 face normal to 5-fold axis
                   q=2,  # natoms 0n 100 parallel to 5-fold axis
                   r=0)  # depth of the Marks re-entrance?

print('#atoms = {}'.format(len(atoms)))

write('images/decahedron.png', atoms)
#atoms = 156

You can see the 5-fold symmetry here. We can make octahedra too. Here, the length is the number of atoms on the edge.

from ase.cluster.octahedron import Octahedron
from ase.io import write

atoms = Octahedron('Cu', length=5)
print atoms
write('images/octahedron.png', atoms)
Cluster(symbols='Cu85', positions=..., cell=[14.44, 14.44, 14.44], pbc=[False, False, False])

Finally, we can make particles based on a Wulff construction! We provide a list of surfaces, and their surface energies, with an approximate size we want, the structure to make the particle in, and what to do if there is not an exact number of atoms matching our size. We choose to round below here, so that the particle is not bigger than our size. In this example I totally made up the surface energies, with (100) as the lowest, so the particle comes out looking like a cube.

from ase.cluster.wulff import wulff_construction
from ase.io import write

atoms = wulff_construction('Pd',
                           surfaces=[(1, 0, 0),
                                     (1, 1, 1),
                                     (1, 1, 0)],
                           energies=[0.1, 0.5, 0.15],
                           size=100,
                           structure='fcc',
                           rounding='below')

print atoms
write('images/wulff.png', atoms)
Cluster(symbols='Pd63', positions=..., cell=[7.779999999999999, 7.779999999999999, 7.779999999999999], pbc=[False, False, False])

This is one handy module, if you need to make clusters for some kind of simulation!

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

org-mode source

Org-mode version = 8.2.10

Discuss on Twitter

Capturing stderr from Python in org-mode - take 2

| categories: python, orgmode, emacs | tags:

In a previous post I wrote about a sandbox module to help capture stderr in Python code blocks in org-mode. That module worked, but ran as a script.

stderr is not captured in the output of a code block in org-mode. For example:

import sys
print >>sys.stdout, "message on stdout"
print >>sys.stderr, "testing stderr"
message on stdout

The messages to stderr just disappears. Not good for code like this:

from scipy.integrate import odeint

def ode(y, x):
    return -k * x

xspan = [0, 1]
y0 = 1

sol = odeint(ode, y0, xspan)
print sol
[[ 1.]
 [ 1.]]

There is an error in that code, k is not defined. If you run that as a script, you get this output:

>>> Traceback (most recent call last):
  File "/var/folders/5q/lllv2yf95hg_n6h6kjttbmdw0000gn/T//python-69413hLF.py", line 4, in ode
    return -k * x
NameError: global name 'k' is not defined
Traceback (most recent call last):
  File "/var/folders/5q/lllv2yf95hg_n6h6kjttbmdw0000gn/T//python-69413hLF.py", line 4, in ode
    return -k * x
NameError: global name 'k' is not defined
Traceback (most recent call last):
  File "/var/folders/5q/lllv2yf95hg_n6h6kjttbmdw0000gn/T//python-69413hLF.py", line 4, in ode
    return -k * x
NameError: global name 'k' is not defined
Traceback (most recent call last):
  File "/var/folders/5q/lllv2yf95hg_n6h6kjttbmdw0000gn/T//python-69413hLF.py", line 4, in ode
    return -k * x
NameError: global name 'k' is not defined

But, that is evidently going to stderr, and not getting captured in org-mode. Boo. A silent error that returns a value! This behavior of odeint may be fixed in scipy 0.15, but it is a general deficiency of org-mode babel code blocks. So, today I am looking back into a way to fix it. We try something as mundane as just redefining stderr in Python at runtime.

import sys
sys.stderr = sys.stdout

print >>sys.stdout, "message on stdout"
print >>sys.stderr, "testing stderr"
message on stdout
testing stderr

That works fine. Let us test it with the other block.

import sys
sys.stderr = sys.stdout

from scipy.integrate import odeint

def ode(y, x):
    return -k * x

xspan = [0, 1]
y0 = 1

sol = odeint(ode, y0, xspan)
print sol
Traceback (most recent call last):
  File "<stdin>", line 6, in ode
NameError: global name 'k' is not defined
Traceback (most recent call last):
  File "<stdin>", line 6, in ode
NameError: global name 'k' is not defined
Traceback (most recent call last):
  File "<stdin>", line 6, in ode
NameError: global name 'k' is not defined
Traceback (most recent call last):
  File "<stdin>", line 6, in ode
NameError: global name 'k' is not defined
[[ 1.]
 [ 1.]]

Sweet, we get the errors. We still get the returned value, but it is immediately obvious something is wrong. I have wrapped that little tidbit into a Python module in pycse.orgmode , which you can import to get the same effect.

import pycse.orgmode

from scipy.integrate import odeint

def ode(y, x):
    return -k * x

xspan = [0, 1]
y0 = 1

sol = odeint(ode, y0, xspan)
print sol
Traceback (most recent call last):
  File "<stdin>", line 2, in ode
NameError: global name 'k' is not defined
Traceback (most recent call last):
  File "<stdin>", line 2, in ode
NameError: global name 'k' is not defined
Traceback (most recent call last):
  File "<stdin>", line 2, in ode
NameError: global name 'k' is not defined
Traceback (most recent call last):
  File "<stdin>", line 2, in ode
NameError: global name 'k' is not defined

[[ 1.]
 [ 1.]]

Finally, you can avoid the import by setting your org-babel Python command like this:

(setq org-babel-python-command "python -i -c \"import pycse.orgmode\"")
python -i -c "import pycse.orgmode"

Now, we run our faulty block again:

from scipy.integrate import odeint

def ode(y, x):
    return -k * x

xspan = [0, 1]
y0 = 1

sol = odeint(ode, y0, xspan)
print sol
Traceback (most recent call last):
  File "<stdin>", line 2, in ode
NameError: global name 'k' is not defined
Traceback (most recent call last):
  File "<stdin>", line 2, in ode
NameError: global name 'k' is not defined
Traceback (most recent call last):
  File "<stdin>", line 2, in ode
NameError: global name 'k' is not defined
Traceback (most recent call last):
  File "<stdin>", line 2, in ode
NameError: global name 'k' is not defined

[[ 1.]
 [ 1.]]

Excellent. The stderr is captured.

And we get basically the same output as before for regular code blocks. There is an extra line before and after the output for some reason. I can live with that!

print 6 + 7
13

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

org-mode source

Org-mode version = 8.2.10

Discuss on Twitter

A new mode for Python documentation

| categories: python, emacs | tags:

[2014-12-22 Mon] update: See this in action at http://www.youtube.com/watch?v=G_r7wTcVK54 , and see the latest source at https://github.com/jkitchin/jmax/blob/master/pydoc.el .

The emacs-lisp documentation in Emacs is inspiring. It is interlinked, you can click on links to open source files, other commands, etc… Python documentation is not that nice. It should be.

I wrote a little pydoc function:

(defun pydoc (name)
  "Display pydoc information for NAME in a buffer named *pydoc*."
  (interactive "sName of function or module: ")
  (switch-to-buffer-other-window "*pydoc*")
  (erase-buffer)
  (insert (shell-command-to-string (format "python -m pydoc %s" name)))
  (goto-char (point-min)))

which at least accesses python documentation in emacs. It looks like this:

But, this lacks functionality. I want there to be useful links in this, so I can click on the filename to open the source, or click on the packages to get their documentation. Below, we walk through a few functions that will operate on the buffer and put text properties on different pieces.

First, let us make the source file clickable so it opens the source.

(defun pydoc-make-file-link ()
  "Find FILE in a pydoc buffer and make it a clickable link"
  (goto-char (point-min))
  (when (re-search-forward "^FILE
    \\(.*\\)$" nil t)

    (let ((map (make-sparse-keymap))
          (start (match-beginning 1))
          (end (match-end 1))
          (source-file (match-string 1)))
      
      ;; set file to be clickable to open the source
      (define-key map [mouse-1]
        `(lambda ()
          (interactive)
          (find-file ,source-file)))

      (set-text-properties
       start end
       `(local-map, map
                   font-lock-face (:foreground "blue"  :underline t)
                   mouse-face highlight
                   help-echo "mouse-1: click to open")))))
pydoc-make-file-link

Next, sometimes there are URLs in the python documentation. These should all open up in a browser when you click on them. Here we propertize anything we recognize as a URL to make it open when clicked on.

(defun pydoc-make-url-links ()
  (goto-char (point-min))
  (while (re-search-forward "\\(http\\(s\\)?://.*$\\)" nil t)
    (let ((map (make-sparse-keymap))
          (start (match-beginning 1))
          (end (match-end 1)))
        
      (define-key map [mouse-1]
        `(lambda ()
          (interactive)
          (browse-url ,(buffer-substring start end))))
        
      (set-text-properties
       start end
       `(local-map ,map
                   font-lock-face (:foreground "blue"  :underline t)
                   mouse-face highlight
                   help-echo (format "mouse-1: click to open"))))))

When we get documentation for a package, we should make each entry of the package clickable, so we can get to the documentation for that package easily. We store the name of the current package so we can construct the path to the subpackage.

(defun pydoc-get-name ()
  "get NAME and store locally"
  (make-variable-buffer-local 'pydoc-name)
  (goto-char (point-min))
  (when (re-search-forward "^NAME
\\s-*\\([^-][a-zA-Z]*\\)" nil t)
    (setq pydoc-name (match-string 1))))


(defun pydoc-make-package-links ()
  "make links in PACKAGE CONTENTS"
  (goto-char (point-min))
  (when (re-search-forward "^PACKAGE CONTENTS" nil t)
    (forward-line)

    (while (string-match
            "^    \\([a-zA-Z0-9_]*\\)[ ]?\\((package)\\)?"
            (buffer-substring
             (line-beginning-position)
             (line-end-position)))
                
      (let ((map (make-sparse-keymap))
            (start (match-beginning 1))
            (end (match-end 1))
            (package (concat
                      pydoc-name "."
                      (match-string 1
                                    (buffer-substring
                                     (line-beginning-position)
                                     (line-end-position))))))
        
        (define-key map [mouse-1]
          `(lambda ()
            (interactive)
            (pydoc ,package)))
          
        (set-text-properties
         (+ (line-beginning-position) start)
         (+ (line-beginning-position) end)
         `(local-map, map
                      font-lock-face (:foreground "blue"  :underline t)
                      mouse-face highlight
                      help-echo (format "mouse-1: click to open %s" ,package))))
      (forward-line))))

Next, we put some eye candy on function names and arguments. This won't do anything functionally, but it breaks up the monotony of all black text.

(defun pydoc-colorize-functions ()
  "Change color of function names and args."
  (goto-char (point-min))
  (when (re-search-forward "^Functions" nil t)  
    (while (re-search-forward "    \\([a-zA-z0-9-]+\\)(\\([^)]*\\))" nil t)
      (set-text-properties
       (match-beginning 1)
       (match-end 1)
       '(font-lock-face (:foreground "brown")))

      (set-text-properties
       (match-beginning 2)
       (match-end 2)
       '(font-lock-face (:foreground "red"))))))

I have gotten used to the [back] link in emacs-lisp documentation, so we try to emulate it here.

(defun pydoc-insert-back-link ()
  "Insert link to previous buffer"
  (goto-char (point-max)) 
  (insert "
[back]")
  (let ((map (make-sparse-keymap)))
    
    ;; set file to be clickable to open the source
    (define-key map [mouse-1]
      (lambda ()
        (interactive)
        (pydoc *pydoc-last*)))

      (set-text-properties
       (line-beginning-position)
       (line-end-position)
       `(local-map, map
                    font-lock-face (:foreground "blue"  :underline t)
                    mouse-face highlight
                    help-echo "mouse-1: click to return"))))
pydoc-insert-back-link

Ok, finally we remake the pydoc function.

(defvar *pydoc-current* nil
 "Stores current pydoc command")

(defvar *pydoc-last* nil
 "Stores the last pydoc command")

(defun pydoc (name)
  "Display pydoc information for NAME in a buffer named *pydoc*."
  (interactive "sName of function or module: ")

  (switch-to-buffer-other-window "*pydoc*")
  (setq buffer-read-only nil)
  (erase-buffer)
  (insert (shell-command-to-string (format "python -m pydoc %s" name)))
  (goto-char (point-min))

  ;; save 
  (when *pydoc-current*
      (setq *pydoc-last* *pydoc-current*))
  (setq *pydoc-current* name)


  (save-excursion
    (pydoc-get-name)
    (pydoc-make-url-links)
    (pydoc-make-file-link)
    (pydoc-make-package-links)
    (pydoc-colorize-functions)
    (pydoc-insert-back-link))

  ;; make read-only and press q to quit
  (setq buffer-read-only t)
  (use-local-map (copy-keymap org-mode-map))
  (local-set-key "q" #'(lambda () (interactive) (kill-buffer)))

  (font-lock-mode))
pydoc

Now, we get a much more functional pydoc:

Figure 2: Annotated screenshot

and with the colorized function names:

Admittedly, there seems to be a lot of boilerplate code for propertizing the strings, but it doesn't seem too bad. I will probably use this documentation tool this spring, so maybe I will think of new functionality to add to pydoc. Any ideas?

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

org-mode source

Org-mode version = 8.2.10

Discuss on Twitter

Using Pymacs to integrate Python into Emacs

| categories: emacs, python | tags:

Pymacs is a project that aims to integrate Python into Emacs, and vice versa. In this post, I am going to examine the Python into Emacs integration. I cloned the git repository, ran make install, and setup my init.el file like this, as suggested in the manual.

(add-to-list 'load-path (expand-file-name "Pymacs" starter-kit-dir))
(require 'pymacs)
(autoload 'pymacs-apply "pymacs")
(autoload 'pymacs-call "pymacs")
(autoload 'pymacs-eval "pymacs" nil t)
(autoload 'pymacs-exec "pymacs" nil t)
(autoload 'pymacs-load "pymacs" nil t)
(autoload 'pymacs-autoload "pymacs")

Pymacs provides some mapping of Python modules to emacs-lisp functions. You load modules in emacs-lisp, and then a dash-mangled version of the Python functions are available, in emacs lisp. Here is an example. We will load numpy, and find the maximum element of an array. For comparison, here is the Python script.

import numpy as np
print np.max(np.array([[1, 1], [2, 4]]))
4

Now, the corresponding emacs version using Pymacs.

(pymacs-load "numpy" "np-")
(np-max (np-array '((1 1) (2 4))))
4

Neat! The dot notation is basically replaced with dash notation, and we use a lisp list as the argument instead of an array. Otherwise, this looks almost identical. Now, let us consider something more complicated, and get the determinant of the array. We add a PREFIX to the load statement for numpy.linalg similar to what we would do in Python:

import numpy as np
import numpy.linalg as la
print la.det(np.array([[1, 1], [2, 4]]))
2.0

And in emacs-lisp:

(pymacs-load "numpy" "np-")
(pymacs-load "numpy.linalg" "la-")
(la-det (np-array '((1 1) (2 4))))
2.0

We can call functions from matplotlib to make a figure. For example:

(pymacs-load "matplotlib.pyplot" "plt-")
(let* ((x  '(1 2 3 4))
       (y  (mapcar (lambda (z) (* z z)) x)))
  (plt-plot x y)
  (plt-xlabel "x values")
  (plt-ylabel "x$^2$")
  (plt-savefig "plt-pymacs.png"))

This was a little subtle. It was necessary to save the lists as variables, and use the variables in the plot command.

I am not sure what this offers over just having a Python block present in org-mode though. Maybe it is more useful in emacs-lisp libraries where you want to bring in some numerical analysis. Or if you have some custom library of Python you would like to use in elisp. Here is a highly contrived example. Suppose we have a Python module with this special function that converts an argument to "J":

def special_func(x):
    return "J"

In Python, we might use it like this:

import my_python as mp
print [mp.special_func(x) for x in [1, 2, 3]]
['J', 'J', 'J']

We can import the module, and use the function in emacs-lisp too. The underscore in the function name is turned into a dash, which is a little confusing, but it works otherwise.

(pymacs-load "my_python" "mp-")
(mapcar 'mp-special-func '(1 2 3))
J J J

It does not seem possible to do everything though. For example, It is not clear how to pass functions through either side. For example, this does not work for fsolve, although it seems like it should.

(pymacs-load "scipy.optimize" "so-")

(defun objective (x)
  (- x 5))

(so-fsolve 'objective 3)

I get an error like this:

Pymacs loading scipy.optimize...done
pymacs-report-error: Python: Emacs: "(wrong-type-argument number-or-marker-p (pymacs-python . 47))"

The Python equivalent is here:

from scipy.optimize import fsolve
def objective(x):
    return x - 5

print fsolve(objective, 3)
[ 5.]

There is an open question on StackOverflow here on this issue. Overall, I find the project very interesting. It would be awesome if you could extend emacs more easily in other languages, especially scripting languages such as Python that have numerical and plotting capabilities. Right now, this is possible in limited ways. For example, Xah Lee describes an approach where an arbitrary script can take data on stdin, process it, and output the results to stdout. Emacs can capture this and use it to modify the buffer. This uses the shell-command features in Emacs. These scripts could be written in Python, Perl, Ruby, etc… This seems like a simpler and more flexible approach, except that it requires creating the shell commands and putting them on the executable path (as opposed to having Python modules on a PYTHONPATH). These lack the deep integration of documentation you get with emacs-lisp and Python functions.

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

org-mode source

Org-mode version = 8.2.7c

Discuss on Twitter
« Previous Page -- Next Page »