(LearnerND) add advanced usage example

[basnijholt]

(original issue on GitLab)

opened by Bas Nijholt (@basnijholt) at 2018-10-20T12:51:05.772Z

We should add the following "real world usage" code to the tutorial as an "Advanced example" once we merged gitlab:!127 and gitlab:!124.

This as a downloadable file kwant_functions.py

from functools import lru_cache
import numpy as np
import scipy.linalg
import scipy.spatial
import kwant


@lru_cache()
def create_syst(unit_cell):
    lat = kwant.lattice.Polyatomic(unit_cell, [(0, 0, 0)])
    syst = kwant.Builder(kwant.TranslationalSymmetry(*lat.prim_vecs))
    syst[lat.shape(lambda _: True, (0, 0, 0))] = 6
    syst[lat.neighbors()] = -1
    return kwant.wraparound.wraparound(syst).finalized()


def get_brillouin_zone(unit_cell):
    syst = create_syst(unit_cell)
    A = get_A(syst)
    neighbours = kwant.linalg.lll.voronoi(A)
    lattice_points = np.concatenate(([[0, 0, 0]], neighbours))
    lattice_points = 2 * np.pi * (lattice_points @ A.T)
    vor = scipy.spatial.Voronoi(lattice_points)
    brillouin_zone = vor.vertices[vor.regions[vor.point_region[0]]]
    return scipy.spatial.ConvexHull(brillouin_zone)


def momentum_to_lattice(k, syst):
    A = get_A(syst)
    k, residuals = scipy.linalg.lstsq(A, k)[:2]
    if np.any(abs(residuals) > 1e-7):
        raise RuntimeError("Requested momentum doesn't correspond"
                           " to any lattice momentum.")
    return k


def get_A(syst):
    B = np.asarray(syst._wrapped_symmetry.periods).T
    return np.linalg.pinv(B).T


def energies(k, unit_cell):
    syst = create_syst(unit_cell)
    k_x, k_y, k_z = momentum_to_lattice(k, syst)
    params = {'k_x': k_x, 'k_y': k_y, 'k_z': k_z}
    H = syst.hamiltonian_submatrix(params=params)
    energies = np.linalg.eigvalsh(H)
    return min(energies)

This in the tutorial:

from functools import partial

from ipywidgets import interact_manual
import numpy as np

import adaptive
from kwant_functions import get_brillouin_zone, energies

adaptive.notebook_extension()

# Define the lattice vectors of some common unit cells
lattices = dict(
    hexegonal=(
        (0, 1, 0),
        (np.cos(-np.pi / 6), np.sin(-np.pi / 6), 0),
        (0, 0, 1)
    ),
    simple_cubic=(
        (1, 0, 0),
        (0, 1, 0),
        (0, 0, 1)
    ),
    fcc=(
        (0, .5, .5),
        (.5, .5, 0),
        (.5, 0, .5)
    ),
    bcc=(
        (-.5, .5, .5),
        (.5, -.5, .5),
        (.5, .5, -.5)
    )
)

learners = []
fnames = []
for name, unit_cell in lattices.items():
    hull = get_brillouin_zone(unit_cell)
    learner = adaptive.LearnerND(
        partial(energies, unit_cell=unit_cell), hull, loss_per_simplex=None
    )
    fnames.append(name)
    learners.append(learner)

learner = adaptive.BalancingLearner(learners, strategy="npoints")
adaptive.runner.simple(learner, goal=lambda l: l.learners[0].npoints > 1000)

# XXX: maybe this could even be a `BalancingLearner` method.
def select(name, learners, fnames):
    return learners[fnames.index(name)]


def iso(unit_cell, level=8.5):
    l = select(unit_cell, learners, fnames)
    return l.plot_isosurface(level=level)


def plot_tri(unit_cell):
    l = select(unit_cell, learners, fnames)
    return l.plot_3D()

[basnijholt]

originally posted by Jorn Hoofwijk (@Jorn) at 2018-10-25T09:30:09.244Z on GitLab

# XXX: maybe this could even be a `BalancingLearner` method.

Maybe the api could be something like

bl = BalancingLearner(learners)
learner_a = bl.get("abc")

[basnijholt]

To focus on a isoline use:

def isoline_loss_function(y_iso, sigma, priority=1):
    from adaptive.learner.learnerND import default_loss

    def gaussian(x, mu, sigma):
        return np.exp(-(x - mu) ** 2 / sigma ** 2 / 2)

    def loss(simplex, values, value_scale):
        distance = np.mean([abs(y_iso * value_scale - y) for y in values])
        return priority * gaussian(distance, 0, sigma) + default_loss(simplex, values, value_scale)

    return loss


level = 0.1
loss = isoline_loss_function(level, sigma=0.4, priority=0.5)