from abc import ABC, abstractmethod
from typing import Dict, List

import numpy as np

from ase.cell import Cell
from ase.dft.kpoints import BandPath, parse_path_string, sc_special_points
from ase.utils import pbc2pbc

_degrees = np.pi / 180


class BravaisLattice(ABC):
    """Represent Bravais lattices and data related to the Brillouin zone.

    There are 14 3D Bravais classes: CUB, FCC, BCC, ..., and TRI, and
    five 2D classes.

    Each class stores basic static information:

    >>> from ase.lattice import FCC, MCL
    >>> FCC.name
    'FCC'
    >>> FCC.longname
    'face-centred cubic'
    >>> FCC.pearson_symbol
    'cF'
    >>> MCL.parameters
    ('a', 'b', 'c', 'alpha')

    Each class can be instantiated with the specific lattice parameters
    that apply to that lattice:

    >>> MCL(3, 4, 5, 80)
    MCL(a=3, b=4, c=5, alpha=80)

    """
    # These parameters can be set by the @bravais decorator for a subclass.
    # (We could also use metaclasses to do this, but that's more abstract)
    name = None  # e.g. 'CUB', 'BCT', 'ORCF', ...
    longname = None  # e.g. 'cubic', 'body-centred tetragonal', ...
    parameters = None  # e.g. ('a', 'c')
    variants = None  # e.g. {'BCT1': <variant object>,
    #                        'BCT2': <variant object>}
    ndim = None

    def __init__(self, **kwargs):
        p = {}
        eps = kwargs.pop('eps', 2e-4)
        for k, v in kwargs.items():
            p[k] = float(v)
        assert set(p) == set(self.parameters)
        self._parameters = p
        self._eps = eps

        if len(self.variants) == 1:
            # If there's only one it has the same name as the lattice type
            self._variant = self.variants[self.name]
        else:
            name = self._variant_name(**self._parameters)
            self._variant = self.variants[name]

    @property
    def variant(self) -> str:
        """Return name of lattice variant.

        >>> from ase.lattice import BCT

        >>> BCT(3, 5).variant
        'BCT2'
        """
        return self._variant.name

    def __getattr__(self, name: str):
        if name in self._parameters:
            return self._parameters[name]
        return self.__getattribute__(name)  # Raises error

    def vars(self) -> Dict[str, float]:
        """Get parameter names and values of this lattice as a dictionary."""
        return dict(self._parameters)

    def conventional(self) -> 'BravaisLattice':
        """Get the conventional cell corresponding to this lattice."""
        cls = bravais_lattices[self.conventional_cls]
        return cls(**self._parameters)

    def tocell(self) -> Cell:
        """Return this lattice as a :class:`~ase.cell.Cell` object."""
        cell = self._cell(**self._parameters)
        return Cell(cell)

    def cellpar(self) -> np.ndarray:
        """Get cell lengths and angles as array of length 6.

        See :func:`ase.geometry.Cell.cellpar`."""
        # (Just a brute-force implementation)
        cell = self.tocell()
        return cell.cellpar()

    @property
    def special_path(self) -> str:
        """Get default special k-point path for this lattice as a string.

        >>> BCT(3, 5).special_path
        'GXYSGZS1NPY1Z,XP'
        """
        return self._variant.special_path

    @property
    def special_point_names(self) -> List[str]:
        """Return all special point names as a list of strings.

        >>> from ase.lattice import BCT

        >>> BCT(3, 5).special_point_names
        ['G', 'N', 'P', 'S', 'S1', 'X', 'Y', 'Y1', 'Z']
        """
        labels = parse_path_string(self._variant.special_point_names)
        assert len(labels) == 1  # list of lists
        return labels[0]

    def get_special_points_array(self) -> np.ndarray:
        """Return all special points for this lattice as an array.

        Ordering is consistent with special_point_names."""
        if self._variant.special_points is not None:
            # Fixed dictionary of special points
            d = self.get_special_points()
            labels = self.special_point_names
            assert len(d) == len(labels)
            points = np.empty((len(d), 3))
            for i, label in enumerate(labels):
                points[i] = d[label]
            return points

        # Special points depend on lattice parameters:
        points = self._special_points(variant=self._variant,
                                      **self._parameters)
        assert len(points) == len(self.special_point_names)
        return np.array(points)

    def get_special_points(self) -> Dict[str, np.ndarray]:
        """Return a dictionary of named special k-points for this lattice."""
        if self._variant.special_points is not None:
            return self._variant.special_points

        labels = self.special_point_names
        points = self.get_special_points_array()

        return dict(zip(labels, points))

    def plot_bz(self, path=None, special_points=None, **plotkwargs):
        """Plot the reciprocal cell and default bandpath."""
        # Create a generic bandpath (no interpolated kpoints):
        bandpath = self.bandpath(path=path, special_points=special_points,
                                 npoints=0)
        return bandpath.plot(dimension=self.ndim, **plotkwargs)

    def bandpath(self, path=None, npoints=None, special_points=None,
                 density=None) -> BandPath:
        """Return a :class:`~ase.dft.kpoints.BandPath` for this lattice.

        See :meth:`ase.cell.Cell.bandpath` for description of parameters.

        >>> from ase.lattice import BCT

        >>> BCT(3, 5).bandpath()
        BandPath(path='GXYSGZS1NPY1Z,XP', cell=[3x3], \
special_points={GNPSS1XYY1Z}, kpts=[51x3])

        .. note:: This produces the standard band path following AFlow
           conventions.  If your cell does not follow this convention,
           you will need :meth:`ase.cell.Cell.bandpath` instead or
           the kpoints may not correspond to your particular cell.

        """
        if special_points is None:
            special_points = self.get_special_points()

        if path is None:
            path = self._variant.special_path
        elif not isinstance(path, str):
            from ase.dft.kpoints import resolve_custom_points
            path, special_points = resolve_custom_points(path,
                                                         special_points,
                                                         self._eps)

        cell = self.tocell()
        bandpath = BandPath(cell=cell, path=path,
                            special_points=special_points)
        return bandpath.interpolate(npoints=npoints, density=density)

    @abstractmethod
    def _cell(self, **kwargs):
        """Return a Cell object from this Bravais lattice.

        Arguments are the dictionary of Bravais parameters."""

    def _special_points(self, **kwargs):
        """Return the special point coordinates as an npoints x 3 sequence.

        Subclasses typically return a nested list.

        Ordering must be same as kpoint labels.

        Arguments are the dictionary of Bravais parameters and the variant."""
        raise NotImplementedError

    def _variant_name(self, **kwargs):
        """Return the name (e.g. 'ORCF3') of variant.

        Arguments will be the dictionary of Bravais parameters."""
        raise NotImplementedError

    def __format__(self, spec):
        tokens = []
        if not spec:
            spec = '.6g'
        template = f'{{}}={{:{spec}}}'

        for name in self.parameters:
            value = self._parameters[name]
            tokens.append(template.format(name, value))
        return '{}({})'.format(self.name, ', '.join(tokens))

    def __str__(self) -> str:
        return self.__format__('')

    def __repr__(self) -> str:
        return self.__format__('.20g')

    def description(self) -> str:
        """Return complete description of lattice and Brillouin zone."""
        points = self.get_special_points()
        labels = self.special_point_names

        coordstring = '\n'.join(['    {:2s} {:7.4f} {:7.4f} {:7.4f}'
                                 .format(label, *points[label])
                                 for label in labels])

        string = """\
{repr}
  {variant}
  Special point coordinates:
{special_points}
""".format(repr=str(self),
           variant=self._variant,
           special_points=coordstring)
        return string

    @classmethod
    def type_description(cls):
        """Return complete description of this Bravais lattice type."""
        desc = """\
Lattice name: {name}
  Long name: {longname}
  Parameters: {parameters}
""".format(**vars(cls))

        chunks = [desc]
        for name in cls.variant_names:
            var = cls.variants[name]
            txt = str(var)
            lines = ['  ' + L for L in txt.splitlines()]
            lines.append('')
            chunks.extend(lines)

        return '\n'.join(chunks)


class Variant:
    variant_desc = """\
Variant name: {name}
  Special point names: {special_point_names}
  Default path: {special_path}
"""

    def __init__(self, name, special_point_names, special_path,
                 special_points=None):
        self.name = name
        self.special_point_names = special_point_names
        self.special_path = special_path
        if special_points is not None:
            special_points = dict(special_points)
            for key, value in special_points.items():
                special_points[key] = np.array(value)
        self.special_points = special_points
        # XXX Should make special_points available as a single array as well
        # (easier to transform that way)

    def __str__(self) -> str:
        return self.variant_desc.format(**vars(self))


bravais_names = []
bravais_lattices = {}
bravais_classes = {}


def bravaisclass(longname, crystal_family, lattice_system, pearson_symbol,
                 parameters, variants, ndim=3):
    """Decorator for Bravais lattice classes.

    This sets a number of class variables and processes the information
    about different variants into a list of Variant objects."""

    def decorate(cls):
        btype = cls.__name__
        cls.name = btype
        cls.longname = longname
        cls.crystal_family = crystal_family
        cls.lattice_system = lattice_system
        cls.pearson_symbol = pearson_symbol
        cls.parameters = tuple(parameters)
        cls.variant_names = []
        cls.variants = {}
        cls.ndim = ndim

        for [name, special_point_names, special_path,
             special_points] in variants:
            cls.variant_names.append(name)
            cls.variants[name] = Variant(name, special_point_names,
                                         special_path, special_points)

        # Register in global list and dictionary
        bravais_names.append(btype)
        bravais_lattices[btype] = cls
        bravais_classes[pearson_symbol] = cls
        return cls

    return decorate


# Common mappings of primitive to conventional cells:
_identity = np.identity(3, int)
_fcc_map = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
_bcc_map = np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]])


class UnconventionalLattice(ValueError):
    pass


class Cubic(BravaisLattice):
    """Abstract class for cubic lattices."""
    conventional_cls = 'CUB'

    def __init__(self, a):
        super().__init__(a=a)


@bravaisclass('primitive cubic', 'cubic', 'cubic', 'cP', 'a',
              [['CUB', 'GXRM', 'GXMGRX,MR', sc_special_points['cubic']]])
class CUB(Cubic):
    conventional_cellmap = _identity

    def _cell(self, a):
        return a * np.eye(3)


@bravaisclass('face-centred cubic', 'cubic', 'cubic', 'cF', 'a',
              [['FCC', 'GKLUWX', 'GXWKGLUWLK,UX', sc_special_points['fcc']]])
class FCC(Cubic):
    conventional_cellmap = _bcc_map

    def _cell(self, a):
        return 0.5 * np.array([[0., a, a], [a, 0, a], [a, a, 0]])


@bravaisclass('body-centred cubic', 'cubic', 'cubic', 'cI', 'a',
              [['BCC', 'GHPN', 'GHNGPH,PN', sc_special_points['bcc']]])
class BCC(Cubic):
    conventional_cellmap = _fcc_map

    def _cell(self, a):
        return 0.5 * np.array([[-a, a, a], [a, -a, a], [a, a, -a]])


@bravaisclass('primitive tetragonal', 'tetragonal', 'tetragonal', 'tP', 'ac',
              [['TET', 'GAMRXZ', 'GXMGZRAZ,XR,MA',
                sc_special_points['tetragonal']]])
class TET(BravaisLattice):
    conventional_cls = 'TET'
    conventional_cellmap = _identity

    def __init__(self, a, c):
        super().__init__(a=a, c=c)

    def _cell(self, a, c):
        return np.diag(np.array([a, a, c]))


# XXX in BCT2 we use S for Sigma.
# Also in other places I think
@bravaisclass('body-centred tetragonal', 'tetragonal', 'tetragonal', 'tI',
              'ac',
              [['BCT1', 'GMNPXZZ1', 'GXMGZPNZ1M,XP', None],
               ['BCT2', 'GNPSS1XYY1Z', 'GXYSGZS1NPY1Z,XP', None]])
class BCT(BravaisLattice):
    conventional_cls = 'TET'
    conventional_cellmap = _fcc_map

    def __init__(self, a, c):
        super().__init__(a=a, c=c)

    def _cell(self, a, c):
        return 0.5 * np.array([[-a, a, c], [a, -a, c], [a, a, -c]])

    def _variant_name(self, a, c):
        return 'BCT1' if c < a else 'BCT2'

    def _special_points(self, a, c, variant):
        a2 = a * a
        c2 = c * c

        assert variant.name in self.variants

        if variant.name == 'BCT1':
            eta = .25 * (1 + c2 / a2)
            points = [[0, 0, 0],
                      [-.5, .5, .5],
                      [0., .5, 0.],
                      [.25, .25, .25],
                      [0., 0., .5],
                      [eta, eta, -eta],
                      [-eta, 1 - eta, eta]]
        else:
            eta = .25 * (1 + a2 / c2)  # Not same eta as BCT1!
            zeta = 0.5 * a2 / c2
            points = [[0., .0, 0.],
                      [0., .5, 0.],
                      [.25, .25, .25],
                      [-eta, eta, eta],
                      [eta, 1 - eta, -eta],
                      [0., 0., .5],
                      [-zeta, zeta, .5],
                      [.5, .5, -zeta],
                      [.5, .5, -.5]]
        return points


def check_orc(a, b, c):
    if not a < b < c:
        raise UnconventionalLattice('Expected a < b < c, got {}, {}, {}'
                                    .format(a, b, c))


class Orthorhombic(BravaisLattice):
    """Abstract class for orthorhombic types."""

    def __init__(self, a, b, c):
        check_orc(a, b, c)
        super().__init__(a=a, b=b, c=c)


@bravaisclass('primitive orthorhombic', 'orthorhombic', 'orthorhombic', 'oP',
              'abc',
              [['ORC', 'GRSTUXYZ', 'GXSYGZURTZ,YT,UX,SR',
                sc_special_points['orthorhombic']]])
class ORC(Orthorhombic):
    conventional_cls = 'ORC'
    conventional_cellmap = _identity

    def _cell(self, a, b, c):
        return np.diag([a, b, c]).astype(float)


@bravaisclass('face-centred orthorhombic', 'orthorhombic', 'orthorhombic',
              'oF', 'abc',
              [['ORCF1', 'GAA1LTXX1YZ', 'GYTZGXA1Y,TX1,XAZ,LG', None],
               ['ORCF2', 'GCC1DD1LHH1XYZ', 'GYCDXGZD1HC,C1Z,XH1,HY,LG', None],
               ['ORCF3', 'GAA1LTXX1YZ', 'GYTZGXA1Y,XAZ,LG', None]])
class ORCF(Orthorhombic):
    conventional_cls = 'ORC'
    conventional_cellmap = _bcc_map

    def _cell(self, a, b, c):
        return 0.5 * np.array([[0, b, c], [a, 0, c], [a, b, 0]])

    def _special_points(self, a, b, c, variant):
        a2 = a * a
        b2 = b * b
        c2 = c * c
        xminus = 0.25 * (1 + a2 / b2 - a2 / c2)
        xplus = 0.25 * (1 + a2 / b2 + a2 / c2)

        if variant.name == 'ORCF1' or variant.name == 'ORCF3':
            zeta = xminus
            eta = xplus

            points = [[0, 0, 0],
                      [.5, .5 + zeta, zeta],
                      [.5, .5 - zeta, 1 - zeta],
                      [.5, .5, .5],
                      [1., .5, .5],
                      [0., eta, eta],
                      [1., 1 - eta, 1 - eta],
                      [.5, 0, .5],
                      [.5, .5, 0]]
        else:
            assert variant.name == 'ORCF2'
            phi = 0.25 * (1 + c2 / b2 - c2 / a2)
            delta = 0.25 * (1 + b2 / a2 - b2 / c2)
            eta = xminus

            points = [[0, 0, 0],
                      [.5, .5 - eta, 1 - eta],
                      [.5, .5 + eta, eta],
                      [.5 - delta, .5, 1 - delta],
                      [.5 + delta, .5, delta],
                      [.5, .5, .5],
                      [1 - phi, .5 - phi, .5],
                      [phi, .5 + phi, .5],
                      [0., .5, .5],
                      [.5, 0., .5],
                      [.5, .5, 0.]]

        return points

    def _variant_name(self, a, b, c):
        diff = 1.0 / (a * a) - 1.0 / (b * b) - 1.0 / (c * c)
        if abs(diff) < self._eps:
            return 'ORCF3'
        return 'ORCF1' if diff > 0 else 'ORCF2'


@bravaisclass('body-centred orthorhombic', 'orthorhombic', 'orthorhombic',
              'oI', 'abc',
              [['ORCI', 'GLL1L2RSTWXX1YY1Z', 'GXLTWRX1ZGYSW,L1Y,Y1Z', None]])
class ORCI(Orthorhombic):
    conventional_cls = 'ORC'
    conventional_cellmap = _fcc_map

    def _cell(self, a, b, c):
        return 0.5 * np.array([[-a, b, c], [a, -b, c], [a, b, -c]])

    def _special_points(self, a, b, c, variant):
        a2 = a**2
        b2 = b**2
        c2 = c**2

        zeta = .25 * (1 + a2 / c2)
        eta = .25 * (1 + b2 / c2)
        delta = .25 * (b2 - a2) / c2
        mu = .25 * (a2 + b2) / c2

        points = [[0., 0., 0.],
                  [-mu, mu, .5 - delta],
                  [mu, -mu, .5 + delta],
                  [.5 - delta, .5 + delta, -mu],
                  [0, .5, 0],
                  [.5, 0, 0],
                  [0., 0., .5],
                  [.25, .25, .25],
                  [-zeta, zeta, zeta],
                  [zeta, 1 - zeta, -zeta],
                  [eta, -eta, eta],
                  [1 - eta, eta, -eta],
                  [.5, .5, -.5]]
        return points


@bravaisclass('base-centred orthorhombic', 'orthorhombic', 'orthorhombic',
              'oC', 'abc',
              [['ORCC', 'GAA1RSTXX1YZ', 'GXSRAZGYX1A1TY,ZT', None]])
class ORCC(BravaisLattice):
    conventional_cls = 'ORC'
    conventional_cellmap = np.array([[1, 1, 0], [-1, 1, 0], [0, 0, 1]])

    def __init__(self, a, b, c):
        # ORCC is the only ORCx lattice with a<b and not a<b<c
        if a >= b:
            raise UnconventionalLattice(f'Expected a < b, got a={a}, b={b}')
        super().__init__(a=a, b=b, c=c)

    def _cell(self, a, b, c):
        return np.array([[0.5 * a, -0.5 * b, 0], [0.5 * a, 0.5 * b, 0],
                         [0, 0, c]])

    def _special_points(self, a, b, c, variant):
        zeta = .25 * (1 + a * a / (b * b))
        points = [[0, 0, 0],
                  [zeta, zeta, .5],
                  [-zeta, 1 - zeta, .5],
                  [0, .5, .5],
                  [0, .5, 0],
                  [-.5, .5, .5],
                  [zeta, zeta, 0],
                  [-zeta, 1 - zeta, 0],
                  [-.5, .5, 0],
                  [0, 0, .5]]
        return points


@bravaisclass('primitive hexagonal', 'hexagonal', 'hexagonal', 'hP',
              'ac',
              [['HEX', 'GMKALH', 'GMKGALHA,LM,KH',
                sc_special_points['hexagonal']]])
class HEX(BravaisLattice):
    conventional_cls = 'HEX'
    conventional_cellmap = _identity

    def __init__(self, a, c):
        super().__init__(a=a, c=c)

    def _cell(self, a, c):
        x = 0.5 * np.sqrt(3)
        return np.array([[0.5 * a, -x * a, 0], [0.5 * a, x * a, 0],
                         [0., 0., c]])


@bravaisclass('primitive rhombohedral', 'hexagonal', 'rhombohedral', 'hR',
              ('a', 'alpha'),
              [['RHL1', 'GBB1FLL1PP1P2QXZ', 'GLB1,BZGX,QFP1Z,LP', None],
               ['RHL2', 'GFLPP1QQ1Z', 'GPZQGFP1Q1LZ', None]])
class RHL(BravaisLattice):
    conventional_cls = 'RHL'
    conventional_cellmap = _identity

    def __init__(self, a, alpha):
        if alpha >= 120:
            raise UnconventionalLattice('Need alpha < 120 degrees, got {}'
                                        .format(alpha))
        super().__init__(a=a, alpha=alpha)

    def _cell(self, a, alpha):
        alpha *= np.pi / 180
        acosa = a * np.cos(alpha)
        acosa2 = a * np.cos(0.5 * alpha)
        asina2 = a * np.sin(0.5 * alpha)
        acosfrac = acosa / acosa2
        xx = (1 - acosfrac**2)
        assert xx > 0.0
        return np.array([[acosa2, -asina2, 0], [acosa2, asina2, 0],
                         [a * acosfrac, 0, a * xx**0.5]])

    def _variant_name(self, a, alpha):
        return 'RHL1' if alpha < 90 else 'RHL2'

    def _special_points(self, a, alpha, variant):
        if variant.name == 'RHL1':
            cosa = np.cos(alpha * _degrees)
            eta = (1 + 4 * cosa) / (2 + 4 * cosa)
            nu = .75 - 0.5 * eta
            points = [[0, 0, 0],
                      [eta, .5, 1 - eta],
                      [.5, 1 - eta, eta - 1],
                      [.5, .5, 0],
                      [.5, 0, 0],
                      [0, 0, -.5],
                      [eta, nu, nu],
                      [1 - nu, 1 - nu, 1 - eta],
                      [nu, nu, eta - 1],
                      [1 - nu, nu, 0],
                      [nu, 0, -nu],
                      [.5, .5, .5]]
        else:
            eta = 1 / (2 * np.tan(alpha * _degrees / 2)**2)
            nu = .75 - 0.5 * eta
            points = [[0, 0, 0],
                      [.5, -.5, 0],
                      [.5, 0, 0],
                      [1 - nu, -nu, 1 - nu],
                      [nu, nu - 1, nu - 1],
                      [eta, eta, eta],
                      [1 - eta, -eta, -eta],
                      [.5, -.5, .5]]
        return points


def check_mcl(a, b, c, alpha):
    if b > c or alpha >= 90:
        raise UnconventionalLattice('Expected b <= c, alpha < 90; '
                                    'got a={}, b={}, c={}, alpha={}'
                                    .format(a, b, c, alpha))


@bravaisclass('primitive monoclinic', 'monoclinic', 'monoclinic', 'mP',
              ('a', 'b', 'c', 'alpha'),
              [['MCL', 'GACDD1EHH1H2MM1M2XYY1Z', 'GYHCEM1AXH1,MDZ,YD', None]])
class MCL(BravaisLattice):
    conventional_cls = 'MCL'
    conventional_cellmap = _identity

    def __init__(self, a, b, c, alpha):
        check_mcl(a, b, c, alpha)
        super().__init__(a=a, b=b, c=c, alpha=alpha)

    def _cell(self, a, b, c, alpha):
        alpha *= _degrees
        return np.array([[a, 0, 0], [0, b, 0],
                         [0, c * np.cos(alpha), c * np.sin(alpha)]])

    def _special_points(self, a, b, c, alpha, variant):
        cosa = np.cos(alpha * _degrees)
        eta = (1 - b * cosa / c) / (2 * np.sin(alpha * _degrees)**2)
        nu = .5 - eta * c * cosa / b

        points = [[0, 0, 0],
                  [.5, .5, 0],
                  [0, .5, .5],
                  [.5, 0, .5],
                  [.5, 0, -.5],
                  [.5, .5, .5],
                  [0, eta, 1 - nu],
                  [0, 1 - eta, nu],
                  [0, eta, -nu],
                  [.5, eta, 1 - nu],
                  [.5, 1 - eta, nu],
                  [.5, eta, -nu],
                  [0, .5, 0],
                  [0, 0, .5],
                  [0, 0, -.5],
                  [.5, 0, 0]]
        return points

    def _variant_name(self, a, b, c, alpha):
        check_mcl(a, b, c, alpha)
        return 'MCL'


@bravaisclass('base-centred monoclinic', 'monoclinic', 'monoclinic', 'mC',
              ('a', 'b', 'c', 'alpha'),
              [['MCLC1', 'GNN1FF1F2F3II1LMXX1X2YY1Z',
                'GYFLI,I1ZF1,YX1,XGN,MG', None],
               ['MCLC2', 'GNN1FF1F2F3II1LMXX1X2YY1Z',
                'GYFLI,I1ZF1,NGM', None],
               ['MCLC3', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',
                'GYFHZIF1,H1Y1XGN,MG', None],
               ['MCLC4', 'GFF1F2HH1H2IMNN1XYY1Y2Y3Z',
                'GYFHZI,H1Y1XGN,MG', None],
               ['MCLC5', 'GFF1F2HH1H2II1LMNN1XYY1Y2Y3Z',
                'GYFLI,I1ZHF1,H1Y1XGN,MG', None]])
class MCLC(BravaisLattice):
    conventional_cls = 'MCL'
    conventional_cellmap = np.array([[1, -1, 0], [1, 1, 0], [0, 0, 1]])

    def __init__(self, a, b, c, alpha):
        check_mcl(a, b, c, alpha)
        super().__init__(a=a, b=b, c=c, alpha=alpha)

    def _cell(self, a, b, c, alpha):
        alpha *= np.pi / 180
        return np.array([[0.5 * a, 0.5 * b, 0], [-0.5 * a, 0.5 * b, 0],
                         [0, c * np.cos(alpha), c * np.sin(alpha)]])

    def _variant_name(self, a, b, c, alpha):
        # from ase.geometry.cell import mclc
        # okay, this is a bit hacky

        # We need the same parameters here as when determining the points.
        # Right now we just repeat the code:
        check_mcl(a, b, c, alpha)

        a2 = a * a
        b2 = b * b
        cosa = np.cos(alpha * _degrees)
        sina = np.sin(alpha * _degrees)
        sina2 = sina**2

        cell = self.tocell()
        lengths_angles = Cell(cell.reciprocal()).cellpar()

        kgamma = lengths_angles[-1]

        eps = self._eps
        # We should not compare angles in degrees versus lengths with
        # the same precision.
        if abs(kgamma - 90) < eps:
            variant = 2
        elif kgamma > 90:
            variant = 1
        elif kgamma < 90:
            num = b * cosa / c + b2 * sina2 / a2
            if abs(num - 1) < eps:
                variant = 4
            elif num < 1:
                variant = 3
            else:
                variant = 5
        variant = 'MCLC' + str(variant)
        return variant

    def _special_points(self, a, b, c, alpha, variant):
        variant = int(variant.name[-1])

        a2 = a * a
        b2 = b * b
        # c2 = c * c
        cosa = np.cos(alpha * _degrees)
        sina = np.sin(alpha * _degrees)
        sina2 = sina**2

        if variant == 1 or variant == 2:
            zeta = (2 - b * cosa / c) / (4 * sina2)
            eta = 0.5 + 2 * zeta * c * cosa / b
            psi = .75 - a2 / (4 * b2 * sina * sina)
            phi = psi + (.75 - psi) * b * cosa / c

            points = [[0, 0, 0],
                      [.5, 0, 0],
                      [0, -.5, 0],
                      [1 - zeta, 1 - zeta, 1 - eta],
                      [zeta, zeta, eta],
                      [-zeta, -zeta, 1 - eta],
                      [1 - zeta, -zeta, 1 - eta],
                      [phi, 1 - phi, .5],
                      [1 - phi, phi - 1, .5],
                      [.5, .5, .5],
                      [.5, 0, .5],
                      [1 - psi, psi - 1, 0],
                      [psi, 1 - psi, 0],
                      [psi - 1, -psi, 0],
                      [.5, .5, 0],
                      [-.5, -.5, 0],
                      [0, 0, .5]]
        elif variant == 3 or variant == 4:
            mu = .25 * (1 + b2 / a2)
            delta = b * c * cosa / (2 * a2)
            zeta = mu - 0.25 + (1 - b * cosa / c) / (4 * sina2)
            eta = 0.5 + 2 * zeta * c * cosa / b
            phi = 1 + zeta - 2 * mu
            psi = eta - 2 * delta

            points = [[0, 0, 0],
                      [1 - phi, 1 - phi, 1 - psi],
                      [phi, phi - 1, psi],
                      [1 - phi, -phi, 1 - psi],
                      [zeta, zeta, eta],
                      [1 - zeta, -zeta, 1 - eta],
                      [-zeta, -zeta, 1 - eta],
                      [.5, -.5, .5],
                      [.5, 0, .5],
                      [.5, 0, 0],
                      [0, -.5, 0],
                      [.5, -.5, 0],
                      [mu, mu, delta],
                      [1 - mu, -mu, -delta],
                      [-mu, -mu, -delta],
                      [mu, mu - 1, delta],
                      [0, 0, .5]]
        elif variant == 5:
            zeta = .25 * (b2 / a2 + (1 - b * cosa / c) / sina2)
            eta = 0.5 + 2 * zeta * c * cosa / b
            mu = .5 * eta + b2 / (4 * a2) - b * c * cosa / (2 * a2)
            nu = 2 * mu - zeta
            omega = (4 * nu - 1 - b2 * sina2 / a2) * c / (2 * b * cosa)
            delta = zeta * c * cosa / b + omega / 2 - .25
            rho = 1 - zeta * a2 / b2

            points = [[0, 0, 0],
                      [nu, nu, omega],
                      [1 - nu, 1 - nu, 1 - omega],
                      [nu, nu - 1, omega],
                      [zeta, zeta, eta],
                      [1 - zeta, -zeta, 1 - eta],
                      [-zeta, -zeta, 1 - eta],
                      [rho, 1 - rho, .5],
                      [1 - rho, rho - 1, .5],
                      [.5, .5, .5],
                      [.5, 0, .5],
                      [.5, 0, 0],
                      [0, -.5, 0],
                      [.5, -.5, 0],
                      [mu, mu, delta],
                      [1 - mu, -mu, -delta],
                      [-mu, -mu, -delta],
                      [mu, mu - 1, delta],
                      [0, 0, .5]]

        return points


tri_angles_explanation = """\
Angles kalpha, kbeta and kgamma of TRI lattice must be 1) all greater \
than 90 degrees with kgamma being the smallest, or 2) all smaller than \
90 with kgamma being the largest, or 3) kgamma=90 being the \
smallest of the three, or 4) kgamma=90 being the largest of the three.  \
Angles of reciprocal lattice are kalpha={}, kbeta={}, kgamma={}.  \
If you don't care, please use Cell.fromcellpar() instead."""

# XXX labels, paths, are all the same.


@bravaisclass('primitive triclinic', 'triclinic', 'triclinic', 'aP',
              ('a', 'b', 'c', 'alpha', 'beta', 'gamma'),
              [['TRI1a', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None],
               ['TRI2a', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None],
               ['TRI1b', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None],
               ['TRI2b', 'GLMNRXYZ', 'XGY,LGZ,NGM,RG', None]])
class TRI(BravaisLattice):
    conventional_cls = 'TRI'
    conventional_cellmap = _identity

    def __init__(self, a, b, c, alpha, beta, gamma):
        super().__init__(a=a, b=b, c=c, alpha=alpha, beta=beta,
                         gamma=gamma)

    def _cell(self, a, b, c, alpha, beta, gamma):
        alpha, beta, gamma = np.array([alpha, beta, gamma])
        singamma = np.sin(gamma * _degrees)
        cosgamma = np.cos(gamma * _degrees)
        cosbeta = np.cos(beta * _degrees)
        cosalpha = np.cos(alpha * _degrees)
        a3x = c * cosbeta
        a3y = c / singamma * (cosalpha - cosbeta * cosgamma)
        a3z = c / singamma * np.sqrt(singamma**2 - cosalpha**2 - cosbeta**2
                                     + 2 * cosalpha * cosbeta * cosgamma)
        return np.array([[a, 0, 0], [b * cosgamma, b * singamma, 0],
                         [a3x, a3y, a3z]])

    def _variant_name(self, a, b, c, alpha, beta, gamma):
        cell = Cell.new([a, b, c, alpha, beta, gamma])
        icellpar = Cell(cell.reciprocal()).cellpar()
        kangles = kalpha, kbeta, kgamma = icellpar[3:]

        def raise_unconventional():
            raise UnconventionalLattice(tri_angles_explanation
                                        .format(*kangles))

        eps = self._eps
        if abs(kgamma - 90) < eps:
            if kalpha > 90 and kbeta > 90:
                var = '2a'
            elif kalpha < 90 and kbeta < 90:
                var = '2b'
            else:
                # Is this possible?  Maybe due to epsilon
                raise_unconventional()
        elif all(kangles > 90):
            if kgamma > min(kangles):
                raise_unconventional()
            var = '1a'
        elif all(kangles < 90):  # and kgamma > max(kalpha, kbeta):
            if kgamma < max(kangles):
                raise_unconventional()
            var = '1b'
        else:
            raise_unconventional()

        return 'TRI' + var

    def _special_points(self, a, b, c, alpha, beta, gamma, variant):
        # (None of the points actually depend on any parameters)
        # (We should store the points openly on the variant objects)
        if variant.name == 'TRI1a' or variant.name == 'TRI2a':
            points = [[0., 0., 0.],
                      [.5, .5, 0],
                      [0, .5, .5],
                      [.5, 0, .5],
                      [.5, .5, .5],
                      [.5, 0, 0],
                      [0, .5, 0],
                      [0, 0, .5]]
        else:
            points = [[0, 0, 0],
                      [.5, -.5, 0],
                      [0, 0, .5],
                      [-.5, -.5, .5],
                      [0, -.5, .5],
                      [0, -0.5, 0],
                      [.5, 0, 0],
                      [-.5, 0, .5]]

        return points


def get_subset_points(names, points):
    newpoints = {name: points[name] for name in names}
    return newpoints


@bravaisclass('primitive oblique', 'monoclinic', None, 'mp',
              ('a', 'b', 'alpha'), [['OBL', 'GYHCH1X', 'GYHCH1XG', None]],
              ndim=2)
class OBL(BravaisLattice):
    def __init__(self, a, b, alpha, **kwargs):
        check_rect(a, b)
        if alpha >= 90:
            raise UnconventionalLattice(
                f'Expected alpha < 90, got alpha={alpha}')
        super().__init__(a=a, b=b, alpha=alpha, **kwargs)

    def _cell(self, a, b, alpha):
        cosa = np.cos(alpha * _degrees)
        sina = np.sin(alpha * _degrees)

        return np.array([[a, 0, 0],
                         [b * cosa, b * sina, 0],
                         [0., 0., 0.]])

    def _special_points(self, a, b, alpha, variant):
        cosa = np.cos(alpha * _degrees)
        eta = (1 - a * cosa / b) / (2 * np.sin(alpha * _degrees)**2)
        nu = .5 - eta * b * cosa / a

        points = [[0, 0, 0],
                  [0, 0.5, 0],
                  [eta, 1 - nu, 0],
                  [.5, .5, 0],
                  [1 - eta, nu, 0],
                  [.5, 0, 0]]

        return points


@bravaisclass('primitive hexagonal', 'hexagonal', None, 'hp', 'a',
              [['HEX2D', 'GMK', 'GMKG',
                get_subset_points('GMK',
                                  sc_special_points['hexagonal'])]],
              ndim=2)
class HEX2D(BravaisLattice):
    def __init__(self, a, **kwargs):
        super().__init__(a=a, **kwargs)

    def _cell(self, a):
        x = 0.5 * np.sqrt(3)
        return np.array([[a, 0, 0],
                         [-0.5 * a, x * a, 0],
                         [0., 0., 0.]])


def check_rect(a, b):
    if a >= b:
        raise UnconventionalLattice(f'Expected a < b, got a={a}, b={b}')


@bravaisclass('primitive rectangular', 'orthorhombic', None, 'op', 'ab',
              [['RECT', 'GXSY', 'GXSYGS',
                get_subset_points('GXSY',
                                  sc_special_points['orthorhombic'])]],
              ndim=2)
class RECT(BravaisLattice):
    def __init__(self, a, b, **kwargs):
        check_rect(a, b)
        super().__init__(a=a, b=b, **kwargs)

    def _cell(self, a, b):
        return np.array([[a, 0, 0],
                         [0, b, 0],
                         [0, 0, 0.]])


@bravaisclass('centred rectangular', 'orthorhombic', None, 'oc',
              ('a', 'alpha'), [['CRECT', 'GXA1Y', 'GXA1YG', None]], ndim=2)
class CRECT(BravaisLattice):
    def __init__(self, a, alpha, **kwargs):
        # It would probably be better to define the CRECT cell
        # by (a, b) rather than (a, alpha).  Then we can require a < b
        # like in ordinary RECT.
        #
        # In 3D, all lattices in the same family generally take
        # identical parameters.
        if alpha >= 90:
            raise UnconventionalLattice(
                f'Expected alpha < 90.  Got alpha={alpha}')
        super().__init__(a=a, alpha=alpha, **kwargs)

    def _cell(self, a, alpha):
        x = np.cos(alpha * _degrees)
        y = np.sin(alpha * _degrees)
        return np.array([[a, 0, 0],
                         [a * x, a * y, 0],
                         [0, 0, 0.]])

    def _special_points(self, a, alpha, variant):
        sina2 = np.sin(alpha / 2 * _degrees)**2
        sina = np.sin(alpha * _degrees)**2
        eta = sina2 / sina
        cosa = np.cos(alpha * _degrees)
        xi = eta * cosa

        points = [[0, 0, 0],
                  [eta, - eta, 0],
                  [0.5 + xi, 0.5 - xi, 0],
                  [0.5, 0.5, 0]]

        return points


@bravaisclass('primitive square', 'tetragonal', None, 'tp', ('a',),
              [['SQR', 'GMX', 'MGXM',
                get_subset_points('GMX', sc_special_points['tetragonal'])]],
              ndim=2)
class SQR(BravaisLattice):
    def __init__(self, a, **kwargs):
        super().__init__(a=a, **kwargs)

    def _cell(self, a):
        return np.array([[a, 0, 0],
                         [0, a, 0],
                         [0, 0, 0.]])


@bravaisclass('primitive line', 'line', None, '?', ('a',),
              [['LINE', 'GX', 'GX', {'G': [0, 0, 0], 'X': [0.5, 0, 0]}]],
              ndim=1)
class LINE(BravaisLattice):
    def __init__(self, a, **kwargs):
        super().__init__(a=a, **kwargs)

    def _cell(self, a):
        return np.array([[a, 0.0, 0.0],
                         [0.0, 0.0, 0.0],
                         [0.0, 0.0, 0.0]])


def celldiff(cell1, cell2):
    """Return a unitless measure of the difference between two cells."""
    cell1 = Cell.ascell(cell1).complete()
    cell2 = Cell.ascell(cell2).complete()
    v1v2 = cell1.volume * cell2.volume
    if v1v2 < 1e-10:
        # (Proposed cell may be linearly dependent)
        return np.inf

    scale = v1v2**(-1. / 3.)  # --> 1/Ang^2
    x1 = cell1 @ cell1.T
    x2 = cell2 @ cell2.T
    dev = scale * np.abs(x2 - x1).max()
    return dev


def get_lattice_from_canonical_cell(cell, eps=2e-4):
    """Return a Bravais lattice representing the given cell.

    This works only for cells that are derived from the standard form
    (as generated by lat.tocell()) or rotations thereof.

    If the given cell does not resemble the known form of a Bravais
    lattice, raise RuntimeError."""
    return LatticeChecker(cell, eps).match()


def identify_lattice(cell, eps=2e-4, *, pbc=True):
    """Find Bravais lattice representing this cell.

    Returns Bravais lattice object representing the cell along with
    an operation that, applied to the cell, yields the same lengths
    and angles as the Bravais lattice object."""
    from ase.geometry.bravais_type_engine import niggli_op_table

    pbc = cell.any(1) & pbc2pbc(pbc)
    npbc = sum(pbc)

    cell = cell.uncomplete(pbc)
    rcell, reduction_op = cell.niggli_reduce(eps=eps)

    # We tabulate the cell's Niggli-mapped versions so we don't need to
    # redo any work when the same Niggli-operation appears multiple times
    # in the table:
    memory = {}

    # We loop through the most symmetric kinds (CUB etc.) and return
    # the first one we find:
    for latname in LatticeChecker.check_orders[npbc]:
        # There may be multiple Niggli operations that produce valid
        # lattices, at least for MCL.  In that case we will pick the
        # one whose angle is closest to 90, but it means we cannot
        # just return the first one we find so we must remember then:
        matching_lattices = []

        for op_key in niggli_op_table[latname]:
            checker_and_op = memory.get(op_key)
            if checker_and_op is None:
                normalization_op = np.array(op_key).reshape(3, 3)
                candidate = Cell(np.linalg.inv(normalization_op.T) @ rcell)
                checker = LatticeChecker(candidate, eps=eps)
                memory[op_key] = (checker, normalization_op)
            else:
                checker, normalization_op = checker_and_op

            lat = checker.query(latname)
            if lat is not None:
                op = normalization_op @ np.linalg.inv(reduction_op)
                matching_lattices.append((lat, op))

        if not matching_lattices:
            continue  # Move to next Bravais lattice

        lat, op = pick_best_lattice(matching_lattices)

        if npbc == 2 and op[2, 2] < 0:
            op = flip_2d_handedness(op)

        return lat, op

    raise RuntimeError('Failed to recognize lattice')


def flip_2d_handedness(op):
    # The 3x3 operation may flip the z axis, but then the x/y
    # components are necessarily also left-handed which
    # means a defacto left-handed 2D bandpath.
    #
    # We repair this by applying an operation that unflips the
    # z axis and interchanges x/y:
    repair_op = np.array([[0, 1, 0], [1, 0, 0], [0, 0, -1]])
    return repair_op @ op


def pick_best_lattice(matching_lattices):
    """Return (lat, op) with lowest orthogonality defect."""
    best = None
    best_defect = np.inf
    for lat, op in matching_lattices:
        cell = lat.tocell().complete()
        orthogonality_defect = np.prod(cell.lengths()) / cell.volume
        if orthogonality_defect < best_defect:
            best = lat, op
            best_defect = orthogonality_defect
    return best


class LatticeChecker:
    # The check order is slightly different than elsewhere listed order
    # as we need to check HEX/RHL before the ORCx family.
    check_orders = {
        1: ['LINE'],
        2: ['SQR', 'RECT', 'HEX2D', 'CRECT', 'OBL'],
        3: ['CUB', 'FCC', 'BCC', 'TET', 'BCT', 'HEX', 'RHL',
            'ORC', 'ORCF', 'ORCI', 'ORCC', 'MCL', 'MCLC', 'TRI']}

    def __init__(self, cell, eps=2e-4):
        """Generate Bravais lattices that look (or not) like the given cell.

        The cell must be reduced to canonical form, i.e., it must
        be possible to produce a cell with the same lengths and angles
        by directly through one of the Bravais lattice classes.

        Generally for internal use (this module).

        For each of the 14 Bravais lattices, this object can produce
        a lattice object which represents the same cell, or None if
        the tolerance eps is not met."""
        self.cell = cell
        self.eps = eps

        self.cellpar = cell.cellpar()
        self.lengths = self.A, self.B, self.C = self.cellpar[:3]
        self.angles = self.cellpar[3:]

        # Use a 'neutral' length for checking cubic lattices
        self.A0 = self.lengths.mean()

        # Vector of the diagonal and then off-diagonal dot products:
        #   [a1 · a1, a2 · a2, a3 · a3, a2 · a3, a3 · a1, a1 · a2]
        self.prods = (cell @ cell.T).flat[[0, 4, 8, 5, 2, 1]]

    def _check(self, latcls, *args):
        if any(arg <= 0 for arg in args):
            return None
        try:
            lat = latcls(*args)
        except UnconventionalLattice:
            return None

        newcell = lat.tocell()
        err = celldiff(self.cell, newcell)
        if err < self.eps:
            return lat

    def match(self):
        """Match cell against all lattices, returning most symmetric match.

        Returns the lattice object.  Raises RuntimeError on failure."""
        for name in self.check_orders[self.cell.rank]:
            lat = self.query(name)
            if lat:
                return lat
        raise RuntimeError('Could not find lattice type for cell '
                           'with lengths and angles {}'
                           .format(self.cell.cellpar().tolist()))

    def query(self, latname):
        """Match cell against named Bravais lattice.

        Return lattice object on success, None on failure."""
        meth = getattr(self, latname)
        lat = meth()
        return lat

    def LINE(self):
        return self._check(LINE, self.lengths[0])

    def SQR(self):
        return self._check(SQR, self.lengths[0])

    def RECT(self):
        return self._check(RECT, *self.lengths[:2])

    def CRECT(self):
        return self._check(CRECT, self.lengths[0], self.angles[2])

    def HEX2D(self):
        return self._check(HEX2D, self.lengths[0])

    def OBL(self):
        return self._check(OBL, *self.lengths[:2], self.angles[2])

    def CUB(self):
        # These methods (CUB, FCC, ...) all return a lattice object if
        # it matches, else None.
        return self._check(CUB, self.A0)

    def FCC(self):
        return self._check(FCC, np.sqrt(2) * self.A0)

    def BCC(self):
        return self._check(BCC, 2.0 * self.A0 / np.sqrt(3))

    def TET(self):
        return self._check(TET, self.A, self.C)

    def _bct_orci_lengths(self):
        # Coordinate-system independent relation for BCT and ORCI
        # standard cells:
        #   a1 · a1 + a2 · a3 == a² / 2
        #   a2 · a2 + a3 · a1 == a² / 2 (BCT)
        #                     == b² / 2 (ORCI)
        #   a3 · a3 + a1 · a2 == c² / 2
        # We use these to get a, b, and c in those cases.
        prods = self.prods
        lengthsqr = 2.0 * (prods[:3] + prods[3:])
        if any(lengthsqr < 0):
            return None
        return np.sqrt(lengthsqr)

    def BCT(self):
        lengths = self._bct_orci_lengths()
        if lengths is None:
            return None
        return self._check(BCT, lengths[0], lengths[2])

    def HEX(self):
        return self._check(HEX, self.A, self.C)

    def RHL(self):
        return self._check(RHL, self.A, self.angles[0])

    def ORC(self):
        return self._check(ORC, *self.lengths)

    def ORCF(self):
        # ORCF standard cell:
        #   a2 · a3 = a²/4
        #   a3 · a1 = b²/4
        #   a1 · a2 = c²/4
        prods = self.prods
        if all(prods[3:] > 0):
            orcf_abc = 2 * np.sqrt(prods[3:])
            return self._check(ORCF, *orcf_abc)

    def ORCI(self):
        lengths = self._bct_orci_lengths()
        if lengths is None:
            return None
        return self._check(ORCI, *lengths)

    def _orcc_ab(self):
        # ORCC: a1 · a1 + a2 · a3 = a²/2
        #       a2 · a2 - a2 · a3 = b²/2
        prods = self.prods
        orcc_sqr_ab = np.empty(2)
        orcc_sqr_ab[0] = 2.0 * (prods[0] + prods[5])
        orcc_sqr_ab[1] = 2.0 * (prods[1] - prods[5])
        if all(orcc_sqr_ab > 0):
            return np.sqrt(orcc_sqr_ab)

    def ORCC(self):
        orcc_lengths_ab = self._orcc_ab()
        if orcc_lengths_ab is None:
            return None
        return self._check(ORCC, *orcc_lengths_ab, self.C)

    def MCL(self):
        return self._check(MCL, *self.lengths, self.angles[0])

    def MCLC(self):
        # MCLC is similar to ORCC:
        orcc_ab = self._orcc_ab()
        if orcc_ab is None:
            return None

        prods = self.prods
        C = self.C
        mclc_a, mclc_b = orcc_ab[::-1]  # a, b reversed wrt. ORCC
        mclc_cosa = 2.0 * prods[3] / (mclc_b * C)
        if -1 < mclc_cosa < 1:
            mclc_alpha = np.arccos(mclc_cosa) * 180 / np.pi
            if mclc_b > C:
                # XXX Temporary fix for certain otherwise
                # unrecognizable lattices.
                #
                # This error could happen if the input lattice maps to
                # something just outside the domain of conventional
                # lattices (less than the tolerance).  Our solution is to
                # propose a nearby conventional lattice instead, which
                # will then be accepted if it's close enough.
                mclc_b = 0.5 * (mclc_b + C)
                C = mclc_b
            return self._check(MCLC, mclc_a, mclc_b, C, mclc_alpha)

    def TRI(self):
        return self._check(TRI, *self.cellpar)


def all_variants():
    """For testing and examples; yield all variants of all lattices."""
    a, b, c = 3., 4., 5.
    alpha = 55.0
    yield CUB(a)
    yield FCC(a)
    yield BCC(a)
    yield TET(a, c)
    bct1 = BCT(2 * a, c)
    bct2 = BCT(a, c)
    assert bct1.variant == 'BCT1'
    assert bct2.variant == 'BCT2'

    yield bct1
    yield bct2

    yield ORC(a, b, c)

    a0 = np.sqrt(1.0 / (1 / b**2 + 1 / c**2))
    orcf1 = ORCF(0.5 * a0, b, c)
    orcf2 = ORCF(1.2 * a0, b, c)
    orcf3 = ORCF(a0, b, c)
    assert orcf1.variant == 'ORCF1'
    assert orcf2.variant == 'ORCF2'
    assert orcf3.variant == 'ORCF3'
    yield orcf1
    yield orcf2
    yield orcf3

    yield ORCI(a, b, c)
    yield ORCC(a, b, c)

    yield HEX(a, c)

    rhl1 = RHL(a, alpha=55.0)
    assert rhl1.variant == 'RHL1'
    yield rhl1

    rhl2 = RHL(a, alpha=105.0)
    assert rhl2.variant == 'RHL2'
    yield rhl2

    # With these lengths, alpha < 65 (or so) would result in a lattice that
    # could also be represented with alpha > 65, which is more conventional.
    yield MCL(a, b, c, alpha=70.0)

    mclc1 = MCLC(a, b, c, 80)
    assert mclc1.variant == 'MCLC1'
    yield mclc1
    # mclc2 has same special points as mclc1

    mclc3 = MCLC(1.8 * a, b, c * 2, 80)
    assert mclc3.variant == 'MCLC3'
    yield mclc3
    # mclc4 has same special points as mclc3

    # XXX We should add MCLC2 and MCLC4 as well.

    mclc5 = MCLC(b, b, 1.1 * b, 70)
    assert mclc5.variant == 'MCLC5'
    yield mclc5

    def get_tri(kcellpar):
        # We build the TRI lattices from cellpars of reciprocal cell
        icell = Cell.fromcellpar(kcellpar)
        cellpar = Cell(4 * icell.reciprocal()).cellpar()
        return TRI(*cellpar)

    tri1a = get_tri([1., 1.2, 1.4, 120., 110., 100.])
    assert tri1a.variant == 'TRI1a'
    yield tri1a

    tri1b = get_tri([1., 1.2, 1.4, 50., 60., 70.])
    assert tri1b.variant == 'TRI1b'
    yield tri1b

    tri2a = get_tri([1., 1.2, 1.4, 120., 110., 90.])
    assert tri2a.variant == 'TRI2a'
    yield tri2a
    tri2b = get_tri([1., 1.2, 1.4, 50., 60., 90.])
    assert tri2b.variant == 'TRI2b'
    yield tri2b

    # Choose an OBL lattice that round-trip-converts to itself.
    # The default a/b/alpha parameters result in another representation
    # of the same lattice.
    yield OBL(a=3.0, b=3.35, alpha=77.85)
    yield RECT(a, b)
    yield CRECT(a, alpha=alpha)
    yield HEX2D(a)
    yield SQR(a)
    yield LINE(a)