Add missing Matrix methods and documentation (#564)

* Add missing Matrix methods and documentation

* Rename adjoint to H, fix __repr__, pass to numpy

* Make H a property that calls getH

* Update docs

doc/docs/Python_User_Interface.md

@@ -577,6 +577,46 @@ geometry = [mp.Prism(vertices, height=1.5, center=mp.Vector3(), material=cSi)]
 
 ```
 
+### 3x3 Matrix
+
+**`Matrix`(c1 [`Vector3`], c2 [`Vector3`], c3 [`Vector3`])**
+The `Matrix` class represents a 3x3 matrix with c1, c2 and c3 as its columns.
+
+```
+m.transpose()
+m.getH() or m.H
+m.determinant()
+m.inverse()
+```
+
+Return the transpose, adjoint (conjugate transpose), determinant, or inverse of the given matrix.
+
+```
+m1 + m2
+m1 - m2
+m1 * m2
+```
+
+Return the sum, difference, or product of the given matrices.
+
+```
+v * m
+m * v
+```
+
+Returns the (3-vector) product of the matrix `m` by the vector `v`, with the vector multiplied on the left or the right respectively.
+
+```
+s * m
+m * s
+```
+
+Scales the matrix `m` by the number `s`.
+
+**`meep.get_rotation_matrix`(axis [`Vector3`], theta)**
+
+Like `Vector3.rotate`, except returns the (unitary) rotation matrix that performs the given rotation. i.e., `get_rotation_matrix(axis, theta) * v` produces the same result as `v.rotate(axis, theta)`.
+
 ### Symmetry
 
 This class is used for the `symmetries` input variable to specify symmetries which must preserve both the structure *and* the sources. Any number of symmetries can be exploited simultaneously but there is no point in specifying redundant symmetries: the computational cell can be reduced by at most a factor of 4 in 2d and 8 in 3d. See also [Exploiting Symmetry](Exploiting_Symmetry.md). This is the base class of the specific symmetries below, so normally you don't create it directly. However, it has two properties which are shared by all symmetries:

python/geom.py

@@ -396,8 +396,24 @@ def __mul__(self, m):
         else:
             raise TypeError("No operation known for 'Matrix * {}'".format(type(m)))
 
+    def __add__(self, m):
+        return Matrix(self.c1 + m.c1, self.c2 + m.c2, self.c3 + m.c3)
+
+    def __sub__(self, m):
+        return Matrix(self.c1 - m.c1, self.c2 - m.c2, self.c3 - m.c3)
+
     def __repr__(self):
-        return "<{}\n {}\n {}>".format(self.c1, self.c2, self.c3)
+        r0 = self.row(0)
+        r1 = self.row(1)
+        r2 = self.row(2)
+        return "<<{} {} {}>\n <{} {} {}>\n <{} {} {}>>".format(r0[0], r0[1], r0[2],
+                                                               r1[0], r1[1], r1[2],
+                                                               r2[0], r2[1], r2[2])
+
+    def __array__(self):
+        return np.array([self.row(0).__array__(),
+                         self.row(1).__array__(),
+                         self.row(2).__array__()])
 
     def row(self, i):
         return Vector3(self.c1[i], self.c2[i], self.c3[i])
@@ -434,9 +450,15 @@ def determinant(self):
         ])
         return sum1 - sum2
 
+    def conj(self):
+        return Matrix(self.c1.conj(), self.c2.conj(), self.c3.conj())
+
     def transpose(self):
         return Matrix(self.row(0), self.row(1), self.row(2))
 
+    def getH(self):
+        return self.transpose().conj()
+
     def inverse(self):
         v1x = self[1][1] * self[2][2] - self[1][2] * self[2][1]
         v1y = self[1][2] * self[2][0] - self[1][0] * self[2][2]
@@ -457,6 +479,8 @@ def inverse(self):
 
         return m.scale(1 / self.determinant())
 
+    H = property(getH, None)
+
 
 class Lattice(object):
 
@@ -707,3 +731,11 @@ def newton(x, a, b, dx):
 
     return newton(x_guess, pick_bound(lambda aa: aa < 0),
                   pick_bound(lambda aa: aa > 0), x_max - x_min)
+
+
+def get_rotation_matrix(axis, theta):
+    """ Returns the rotation matrix for rotating by theta around axis """
+
+    return Matrix(Vector3(x=1).rotate(axis, theta),
+                  Vector3(y=1).rotate(axis, theta),
+                  Vector3(z=1).rotate(axis, theta))

python/meep.i

@@ -1258,6 +1258,7 @@ kpoint_list get_eigenmode_coefficients_and_kpoints(meep::fields *f, meep::dft_fl
         cartesian_to_reciprocal,
         reciprocal_to_cartesian,
         find_root_deriv,
+        get_rotation_matrix,
     )
     from .simulation import (
         Absorber,

python/tests/geom.py

@@ -465,6 +465,19 @@ def test_mm_mult(self):
 
         self.matrix_eq(exp, res)
 
+    def test_add(self):
+        result = self.identity + self.identity
+        self.assertEqual(result.row(0), mp.Vector3(x=2))
+        self.assertEqual(result.row(1), mp.Vector3(y=2))
+        self.assertEqual(result.row(2), mp.Vector3(z=2))
+
+    def test_sub(self):
+        ones_matrix = mp.Matrix(ones(), ones(), ones())
+        result = ones_matrix - ones_matrix
+        self.assertEqual(result.row(0), zeros())
+        self.assertEqual(result.row(1), zeros())
+        self.assertEqual(result.row(2), zeros())
+
     def test_mv_mult(self):
         lattice = mp.Lattice(size=mp.Vector3(1, 7),
                              basis1=mp.Vector3(math.sqrt(3) / 2, 0.5),
@@ -519,6 +532,35 @@ def test_inverse(self):
 
         self.matrix_close(exp, res)
 
+    def test_get_rotation_matrix(self):
+        result = mp.get_rotation_matrix(ones(), 5)
+        self.assertTrue(result.c1.close(mp.Vector3(0.5224414569754843, -0.3148559165969717, 0.7924144596214877)))
+        self.assertTrue(result.c2.close(mp.Vector3(0.7924144596214877, 0.5224414569754843, -0.3148559165969717)))
+        self.assertTrue(result.c3.close(mp.Vector3(-0.3148559165969717, 0.7924144596214877, 0.5224414569754843)))
+
+    def test_conj(self):
+        m = mp.Matrix(mp.Vector3(x=1+1j), mp.Vector3(y=1+1j), mp.Vector3(z=1+1j))
+        result = m.conj()
+        self.assertEqual(result.c1, mp.Vector3(x=1-1j))
+        self.assertEqual(result.c2, mp.Vector3(y=1-1j))
+        self.assertEqual(result.c3, mp.Vector3(z=1-1j))
+
+    def test_adjoint(self):
+        m = mp.Matrix(mp.Vector3(1+1j), mp.Vector3(1+1j), mp.Vector3(1+1j))
+        getH_result = m.getH()
+        H_result = m.H
+        self.assertEqual(getH_result.c1, mp.Vector3(1-1j, 1-1j, 1-1j))
+        self.assertEqual(getH_result.c2, mp.Vector3())
+        self.assertEqual(getH_result.c3, mp.Vector3())
+        np.testing.assert_allclose(getH_result, H_result)
+
+    def test_to_numpy_array(self):
+        m = mp.Matrix(mp.Vector3(1+1j), mp.Vector3(1+1j), mp.Vector3(1+1j))
+        adjoint = m.H
+        m_arr = np.matrix(m)
+        np_adjoint = m_arr.H
+        np.testing.assert_allclose(adjoint, np_adjoint)
+
 
 if __name__ == '__main__':
     unittest.main()