deprecate

README.md

@@ -1,5 +1,11 @@
 # fdfd_tools
 
+** DEPRECATED **
+
+The functionality in this module is now provided by [meanas](https://mpxd.net/code/jan/meanas).
+
+-----------------------
+
 **fdfd_tools** is a python package containing utilities for
 creating and analyzing 2D and 3D finite-difference frequency-domain (FDFD)
 electromagnetic simulations.

examples/bloch.py

@@ -30,11 +30,11 @@ g2.shifts = numpy.zeros((6,3))
 g2.grids = [numpy.zeros(g.shape) for _ in range(6)]
 
 epsilon = [g.grids[0],] * 3
-reciprocal_lattice = numpy.diag(1000/numpy.array([x_period, y_period, z_period])) #cols are vectors
+reciprocal_lattice = numpy.diag(1e6/numpy.array([x_period, y_period, z_period])) #cols are vectors
 
 #print('Finding k at 1550nm')
-#k, f = bloch.find_k(frequency=1000/1550,
-#                    tolerance=(1000 * (1/1550 - 1/1551)),
+#k, f = bloch.find_k(frequency=1/1550,
+#                    tolerance=(1/1550 - 1/1551),
 #                    direction=[1, 0, 0],
 #                    G_matrix=reciprocal_lattice,
 #                    epsilon=epsilon,
@@ -47,10 +47,10 @@ for k0x in [.25]:
     k0 = numpy.array([k0x, 0, 0])
 
     kmag = norm(reciprocal_lattice @ k0)
-    tolerance = (1000/1550) * 1e-4/1.5  # df = f * dn_eff / n
+    tolerance = (1e6/1550) * 1e-4/1.5  # df = f * dn_eff / n
     logger.info('tolerance {}'.format(tolerance))
 
-    n, v = bloch.eigsolve(4, k0, G_matrix=reciprocal_lattice, epsilon=epsilon, tolerance=tolerance**2)
+    n, v = bloch.eigsolve(4, k0, G_matrix=reciprocal_lattice, epsilon=epsilon, tolerance=tolerance)
     v2e = bloch.hmn_2_exyz(k0, G_matrix=reciprocal_lattice, epsilon=epsilon)
     v2h = bloch.hmn_2_hxyz(k0, G_matrix=reciprocal_lattice, epsilon=epsilon)
     ki = bloch.generate_kmn(k0, reciprocal_lattice, g.shape)

fdfd_tools/bloch.py

@@ -73,44 +73,21 @@ This module contains functions for generating and solving the
 
 '''
 
-from typing import Tuple, Callable
+from typing import List, Tuple, Callable, Dict
 import logging
 import numpy
-from numpy import pi, real, trace
-from numpy.fft import fftfreq
+from numpy.fft import fftn, ifftn, fftfreq
 import scipy
 import scipy.optimize
 from scipy.linalg import norm
 import scipy.sparse.linalg as spalg
 
+from .eigensolvers import rayleigh_quotient_iteration
 from . import field_t
 
 logger = logging.getLogger(__name__)
 
 
-try:
-    import pyfftw.interfaces.numpy_fft
-    import pyfftw.interfaces
-    import multiprocessing
-
-    pyfftw.interfaces.cache.enable()
-    pyfftw.interfaces.cache.set_keepalive_time(3600)
-    fftw_args = {
-        'threads': multiprocessing.cpu_count(),
-        'overwrite_input': True,
-        'planner_effort': 'FFTW_EXHAUSTIVE',
-        }
-
-    def fftn(*args, **kwargs):
-        return pyfftw.interfaces.numpy_fft.fftn(*args, **kwargs, **fftw_args)
-
-    def ifftn(*args, **kwargs):
-        return pyfftw.interfaces.numpy_fft.ifftn(*args, **kwargs, **fftw_args)
-
-except ImportError:
-    from numpy.fft import fftn, ifftn
-
-
 def generate_kmn(k0: numpy.ndarray,
                  G_matrix: numpy.ndarray,
                  shape: numpy.ndarray
@@ -278,7 +255,7 @@ def hmn_2_hxyz(k0: numpy.ndarray,
     :return: Function for converting h_mn into H_xyz
     """
     shape = epsilon[0].shape + (1,)
-    _k_mag, m, n = generate_kmn(k0, G_matrix, shape)
+    k_mag, m, n = generate_kmn(k0, G_matrix, shape)
 
     def operator(h: numpy.ndarray):
         hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
@@ -352,14 +329,147 @@ def inverse_maxwell_operator_approx(k0: numpy.ndarray,
         d_xyz = fftn(ifftn(e_xyz, axes=range(3)) * epsilon, axes=range(3))
 
         # cross product and transform into mn basis   crossinv_t2c
-        h_m = numpy.sum(d_xyz * n, axis=3)[:, :, :, None] / +k_mag
-        h_n = numpy.sum(d_xyz * m, axis=3)[:, :, :, None] / -k_mag
+        h_m = numpy.sum(e_xyz * n, axis=3)[:, :, :, None] / +k_mag
+        h_n = numpy.sum(e_xyz * m, axis=3)[:, :, :, None] / -k_mag
 
         return numpy.hstack((h_m.ravel(), h_n.ravel()))
 
     return operator
 
 
+def eigsolve(num_modes: int,
+             k0: numpy.ndarray,
+             G_matrix: numpy.ndarray,
+             epsilon: field_t,
+             mu: field_t = None,
+             tolerance = 1e-8,
+             ) -> Tuple[numpy.ndarray, numpy.ndarray]:
+    """
+    Find the first (lowest-frequency) num_modes eigenmodes with Bloch wavevector
+     k0 of the specified structure.
+
+    :param k0: Bloch wavevector, [k0x, k0y, k0z].
+    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
+    :param epsilon: Dielectric constant distribution for the simulation.
+        All fields are sampled at cell centers (i.e., NOT Yee-gridded)
+    :param mu: Magnetic permability distribution for the simulation.
+        Default None (1 everywhere).
+    :return: (eigenvalues, eigenvectors) where eigenvalues[i] corresponds to the
+        vector eigenvectors[i, :]
+    """
+    h_size = 2 * epsilon[0].size
+
+    kmag = norm(G_matrix @ k0)
+
+    '''
+    Generate the operators
+    '''
+    mop = maxwell_operator(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu)
+    imop = inverse_maxwell_operator_approx(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu)
+
+    scipy_op = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=mop)
+    scipy_iop = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=imop)
+
+    y_shape = (h_size, num_modes)
+
+    def rayleigh_quotient(Z: numpy.ndarray, approx_grad: bool = True):
+        """
+        Absolute value of the block Rayleigh quotient, and the associated gradient.
+
+        See Johnson and Joannopoulos, Opt. Expr. 8, 3 (2001) for details (full
+         citation in module docstring).
+
+        ===
+
+        Notes on my understanding of the procedure:
+
+        Minimize f(Y) = |trace((Y.H @ A @ Y)|, making use of Y = Z @ inv(Z.H @ Z)^(1/2)
+         (a polar orthogonalization of Y). This gives f(Z) = |trace(Z.H @ A @ Z @ U)|,
+         where U = inv(Z.H @ Z). We minimize the absolute value to find the eigenvalues
+         with smallest magnitude.
+
+        The gradient is P @ (A @ Z @ U), where P = (1 - Z @ U @ Z.H) is a projection
+         onto the space orthonormal to Z. If approx_grad is True, the approximate
+         inverse of the maxwell operator is used to precondition the gradient.
+        """
+        z = Z.view(dtype=complex).reshape(y_shape)
+        U = numpy.linalg.inv(z.conj().T @ z)
+        zU = z @ U
+        AzU = scipy_op @ zU
+        zTAzU = z.conj().T @ AzU
+        f = numpy.real(numpy.trace(zTAzU))
+        if approx_grad:
+            df_dy = scipy_iop @ (AzU - zU @ zTAzU)
+        else:
+            df_dy = (AzU - zU @ zTAzU)
+        
+        df_dy_flat = df_dy.view(dtype=float).ravel()
+        return numpy.abs(f), numpy.sign(f) * df_dy_flat
+
+    '''
+    Use the conjugate gradient method and the approximate gradient calculation to
+     quickly find approximate eigenvectors.
+    '''
+    result = scipy.optimize.minimize(rayleigh_quotient,
+                                     numpy.random.rand(*y_shape, 2),
+                                     jac=True,
+                                     method='L-BFGS-B',
+                                     tol=1e-20,
+                                     options={'maxiter': 2000, 'gtol':0, 'ftol':1e-20 , 'disp':True})#, 'maxls':80, 'm':30})
+
+
+    result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, True),
+                                     result.x,
+                                     jac=True,
+                                     method='L-BFGS-B',
+                                     tol=1e-20,
+                                     options={'maxiter': 2000, 'gtol':0, 'disp':True})
+
+    result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False),
+                                     result.x,
+                                     jac=True,
+                                     method='L-BFGS-B',
+                                     tol=1e-20,
+                                     options={'maxiter': 2000, 'gtol':0, 'disp':True})
+
+    for i in range(20):
+        result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False),
+                                     result.x,
+                                     jac=True,
+                                     method='L-BFGS-B',
+                                     tol=1e-20,
+                                     options={'maxiter': 70, 'gtol':0, 'disp':True})
+        if result.nit == 0:
+            # We took 0 steps, so re-running won't help
+            break
+
+
+    z = result.x.view(dtype=complex).reshape(y_shape)
+
+    '''
+    Recover eigenvectors from Z
+    '''
+    U = numpy.linalg.inv(z.conj().T @ z)
+    y = z @ scipy.linalg.sqrtm(U)
+    w = y.conj().T @ (scipy_op @ y)
+
+    eigvals, w_eigvecs = numpy.linalg.eig(w)
+    eigvecs = y @ w_eigvecs
+
+    for i in range(len(eigvals)):
+        v = eigvecs[:, i]
+        n = eigvals[i]
+        v /= norm(v)
+        eigness = norm(scipy_op @ v - (v.conj() @ (scipy_op @ v)) * v )
+        f = numpy.sqrt(-numpy.real(n))
+        df = numpy.sqrt(-numpy.real(n + eigness))
+        neff_err = kmag * (1/df - 1/f)
+        logger.info('eigness {}: {}\n neff_err: {}'.format(i, eigness, neff_err))
+
+    order = numpy.argsort(numpy.abs(eigvals))
+    return eigvals[order], eigvecs.T[order]
+
+
 def find_k(frequency: float,
            tolerance: float,
            direction: numpy.ndarray,
@@ -389,7 +499,7 @@ def find_k(frequency: float,
 
     def get_f(k0_mag: float, band: int = 0):
         k0 = direction * k0_mag
-        n, _v = eigsolve(band + 1, k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu)
+        n, _v = eigsolve(band + 1, k0, G_matrix=G_matrix, epsilon=epsilon)
         f = numpy.sqrt(numpy.abs(numpy.real(n[band])))
         return f
 
@@ -401,244 +511,3 @@ def find_k(frequency: float,
     return res.x * direction, res.fun + frequency
 
 
-def eigsolve(num_modes: int,
-             k0: numpy.ndarray,
-             G_matrix: numpy.ndarray,
-             epsilon: field_t,
-             mu: field_t = None,
-             tolerance: float = 1e-20,
-             max_iters: int = 10000,
-             reset_iters: int = 100,
-             ) -> Tuple[numpy.ndarray, numpy.ndarray]:
-    """
-    Find the first (lowest-frequency) num_modes eigenmodes with Bloch wavevector
-     k0 of the specified structure.
-
-    :param k0: Bloch wavevector, [k0x, k0y, k0z].
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        All fields are sampled at cell centers (i.e., NOT Yee-gridded)
-    :param mu: Magnetic permability distribution for the simulation.
-        Default None (1 everywhere).
-    :param tolerance: Solver stops when fractional change in the objective
-        trace(Z.H @ A @ Z @ inv(Z Z.H)) is smaller than the tolerance
-    :return: (eigenvalues, eigenvectors) where eigenvalues[i] corresponds to the
-        vector eigenvectors[i, :]
-    """
-    h_size = 2 * epsilon[0].size
-
-    kmag = norm(G_matrix @ k0)
-
-    '''
-    Generate the operators
-    '''
-    mop = maxwell_operator(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu)
-    imop = inverse_maxwell_operator_approx(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu)
-
-    scipy_op = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=mop)
-    scipy_iop = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=imop)
-
-    y_shape = (h_size, num_modes)
-
-    prev_E = 0
-    d_scale = 1
-    prev_traceGtKG = 0
-    #prev_theta = 0.5
-    D = numpy.zeros(shape=y_shape, dtype=complex)
-
-    y0 = None
-    if y0 is None:
-        Z = numpy.random.rand(*y_shape) + 1j * numpy.random.rand(*y_shape)
-    else:
-        Z = y0
-
-    while True:
-        Z *= num_modes / norm(Z)
-        ZtZ = Z.conj().T @ Z
-        try:
-            U = numpy.linalg.inv(ZtZ)
-        except numpy.linalg.LinAlgError:
-            Z = numpy.random.rand(*y_shape) + 1j * numpy.random.rand(*y_shape)
-            continue
-
-        trace_U = real(trace(U))
-        if trace_U > 1e8 * num_modes:
-            Z = Z @ scipy.linalg.sqrtm(U).conj().T
-            prev_traceGtKG = 0
-            continue
-        break
-
-    for i in range(max_iters):
-        ZtZ = Z.conj().T @ Z
-        U = numpy.linalg.inv(ZtZ)
-        AZ = scipy_op @ Z
-        AZU = AZ @ U
-        ZtAZU = Z.conj().T @ AZU
-        E_signed = real(trace(ZtAZU))
-        sgn = numpy.sign(E_signed)
-        E = numpy.abs(E_signed)
-        G = (AZU - Z @ U @ ZtAZU) * sgn
-
-        if i > 0 and abs(E - prev_E) < tolerance * 0.5 * (E + prev_E + 1e-7):
-            logging.info('Optimization succeded: {} - 5e-8 < {} * {} / 2'.format(abs(E - prev_E), tolerance, E + prev_E))
-            break
-
-        KG = scipy_iop @ G
-        traceGtKG = _rtrace_AtB(G, KG)
-
-        if prev_traceGtKG == 0 or i % reset_iters == 0:
-            logger.info('CG reset')
-            gamma = 0
-        else:
-            gamma = traceGtKG / prev_traceGtKG
-
-        D = gamma / d_scale * D + KG
-        d_scale = num_modes / norm(D)
-        D *= d_scale
-
-        ZtAZ = Z.conj().T @ AZ
-
-        AD = scipy_op @ D
-        DtD = D.conj().T @ D
-        DtAD = D.conj().T @ AD
-
-        symZtD = _symmetrize(Z.conj().T @ D)
-        symZtAD = _symmetrize(Z.conj().T @ AD)
-
-        Qi_memo = [None, None]
-        def Qi_func(theta):
-            nonlocal Qi_memo
-            if Qi_memo[0] == theta:
-                return Qi_memo[1]
-
-            c = numpy.cos(theta)
-            s = numpy.sin(theta)
-            Q = c*c * ZtZ + s*s * DtD + 2*s*c * symZtD
-            try:
-                Qi = numpy.linalg.inv(Q)
-            except numpy.linalg.LinAlgError:
-                logger.info('taylor Qi')
-                # if c or s small, taylor expand
-                if c < 1e-4 * s and c != 0:
-                    DtDi = numpy.linalg.inv(DtD)
-                    Qi = DtDi / (s*s) - 2*c/(s*s*s) * (DtDi @ (DtDi @ symZtD).conj().T)
-                elif s < 1e-4 * c and s != 0:
-                    ZtZi = numpy.linalg.inv(ZtZ)
-                    Qi = ZtZi / (c*c) - 2*s/(c*c*c) * (ZtZi @ (ZtZi @ symZtD).conj().T)
-                else:
-                    raise Exception('Inexplicable singularity in trace_func')
-            Qi_memo[0] = theta
-            Qi_memo[1] = Qi
-            return Qi
-
-        def trace_func(theta):
-            c = numpy.cos(theta)
-            s = numpy.sin(theta)
-            Qi = Qi_func(theta)
-            R = c*c * ZtAZ + s*s * DtAD + 2*s*c * symZtAD
-            trace = _rtrace_AtB(R, Qi)
-            return numpy.abs(trace)
-
-        '''
-        def trace_deriv(theta):
-            Qi = Qi_func(theta)
-            c2 = numpy.cos(2 * theta)
-            s2 = numpy.sin(2 * theta)
-            F = -0.5*s2 *  (ZtAZ - DtAD) + c2 * symZtAD
-            trace_deriv = _rtrace_AtB(Qi, F)
-
-            G = Qi @ F.conj().T @ Qi.conj().T
-            H = -0.5*s2 * (ZtZ - DtD) + c2 * symZtD
-            trace_deriv -= _rtrace_AtB(G, H)
-
-            trace_deriv *= 2
-            return trace_deriv * sgn
-
-        U_sZtD = U @ symZtD
-
-        dE = 2.0 * (_rtrace_AtB(U, symZtAD) -
-                    _rtrace_AtB(ZtAZU, U_sZtD))
-
-        d2E = 2 * (_rtrace_AtB(U, DtAD) -
-                   _rtrace_AtB(ZtAZU, U @ (DtD - 4 * symZtD @ U_sZtD)) -
-               4 * _rtrace_AtB(U, symZtAD @ U_sZtD))
-
-        # Newton-Raphson to find a root of the first derivative:
-        theta = -dE/d2E
-
-        if d2E < 0 or abs(theta) >= pi:
-            theta = -abs(prev_theta) * numpy.sign(dE)
-
-        # theta, new_E, new_dE = linmin(theta, E, dE, 0.1, min(tolerance, 1e-6), 1e-14, 0, -numpy.sign(dE) * K_PI, trace_func)
-        theta, n, _, new_E, _, _new_dE = scipy.optimize.line_search(trace_func, trace_deriv, xk=theta, pk=numpy.ones((1,1)), gfk=dE, old_fval=E, c1=min(tolerance, 1e-6), c2=0.1, amax=pi)
-        '''
-        result = scipy.optimize.minimize_scalar(trace_func, bounds=(0, pi), tol=tolerance)
-        new_E = result.fun
-        theta = result.x
-
-        improvement = numpy.abs(E - new_E) * 2 / numpy.abs(E + new_E)
-        logger.info('linmin improvement {}'.format(improvement))
-        Z *= numpy.cos(theta)
-        Z += D * numpy.sin(theta)
-
-        prev_traceGtKG = traceGtKG
-        #prev_theta = theta
-        prev_E = E
-
-    '''
-    Recover eigenvectors from Z
-    '''
-    U = numpy.linalg.inv(Z.conj().T @ Z)
-    Y = Z @ scipy.linalg.sqrtm(U)
-    W = Y.conj().T @ (scipy_op @ Y)
-
-    eigvals, W_eigvecs = numpy.linalg.eig(W)
-    eigvecs = Y @ W_eigvecs
-
-    for i in range(len(eigvals)):
-        v = eigvecs[:, i]
-        n = eigvals[i]
-        v /= norm(v)
-        eigness = norm(scipy_op @ v - (v.conj() @ (scipy_op @ v)) * v )
-        f = numpy.sqrt(-numpy.real(n))
-        df = numpy.sqrt(-numpy.real(n + eigness))
-        neff_err = kmag * (1/df - 1/f)
-        logger.info('eigness {}: {}\n neff_err: {}'.format(i, eigness, neff_err))
-
-    order = numpy.argsort(numpy.abs(eigvals))
-    return eigvals[order], eigvecs.T[order]
-
-'''
-def linmin(x_guess, f0, df0, x_max, f_tol=0.1, df_tol=min(tolerance, 1e-6), x_tol=1e-14, x_min=0, linmin_func):
-    if df0 > 0:
-        x0, f0, df0 = linmin(-x_guess, f0, -df0, -x_max, f_tol, df_tol, x_tol, -x_min, lambda q, dq: -linmin_func(q, dq))
-        return -x0, f0, -df0
-    elif df0 == 0:
-        return 0, f0, df0
-    else:
-        x = x_guess
-        fx = f0
-        dfx = df0
-
-        isave = numpy.zeros((2,), numpy.intc)
-        dsave = numpy.zeros((13,), float)
-
-        x, fx, dfx, task = minpack2.dsrch(x, fx, dfx, f_tol, df_tol, x_tol, task,
-                                          x_min, x_max, isave, dsave)
-        for i in range(int(1e6)):
-            if task != 'F':
-                logging.info('search converged in {} iterations'.format(i))
-                break
-            fx = f(x, dfx)
-            x, fx, dfx, task = minpack2.dsrch(x, fx, dfx, f_tol, df_tol, x_tol, task,
-                                              x_min, x_max, isave, dsave)
-
-        return x, fx, dfx
-'''
-
-def _rtrace_AtB(A, B):
-    return real(numpy.sum(A.conj() * B))
-
-def _symmetrize(A):
-    return (A + A.conj().T) * 0.5
-