"""Functions for fitting head positions with (c)HPI coils.

``compute_head_pos`` can be used to:

1. Drop coils whose GOF are below ``gof_limit``. If fewer than 3 coils
   remain, abandon fitting for the chunk.
2. Fit dev_head_t quaternion (using ``_fit_chpi_quat_subset``),
   iteratively dropping coils (as long as 3 remain) to find the best GOF
   (using ``_fit_chpi_quat``).
3. If fewer than 3 coils meet the ``dist_limit`` criteria following
   projection of the fitted device coil locations into the head frame,
   abandon fitting for the chunk.

The function ``filter_chpi`` uses the same linear model to filter cHPI
and (optionally) line frequencies from the data.
"""

# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.

import copy
import itertools
from functools import partial

import numpy as np
from scipy.linalg import orth
from scipy.optimize import fmin_cobyla
from scipy.spatial.distance import cdist

from ._fiff.constants import FIFF
from ._fiff.meas_info import Info, _simplify_info
from ._fiff.pick import (
    _picks_to_idx,
    pick_channels,
    pick_channels_regexp,
    pick_info,
    pick_types,
)
from ._fiff.proj import Projection, setup_proj
from .channels.channels import _get_meg_system
from .cov import compute_whitener, make_ad_hoc_cov
from .dipole import _make_guesses
from .event import find_events
from .fixes import jit
from .forward import _concatenate_coils, _create_meg_coils, _magnetic_dipole_field_vec
from .io import BaseRaw
from .io.ctf.trans import _make_ctf_coord_trans_set
from .io.kit.constants import KIT
from .io.kit.kit import RawKIT as _RawKIT
from .preprocessing.maxwell import (
    _get_mf_picks_fix_mags,
    _prep_mf_coils,
    _regularize_out,
    _sss_basis,
)
from .transforms import (
    _angle_between_quats,
    _fit_matched_points,
    _quat_to_affine,
    als_ras_trans,
    apply_trans,
    invert_transform,
    quat_to_rot,
    rot_to_quat,
)
from .utils import (
    ProgressBar,
    _check_fname,
    _check_option,
    _on_missing,
    _pl,
    _validate_type,
    _verbose_safe_false,
    logger,
    use_log_level,
    verbose,
    warn,
)

# Eventually we should add:
#   hpicons
#   high-passing of data during fits
#   parsing cHPI coil information from acq pars, then to PSD if necessary


# ############################################################################
# Reading from text or FIF file


def read_head_pos(fname):
    """Read MaxFilter-formatted head position parameters.

    Parameters
    ----------
    fname : path-like
        The filename to read. This can be produced by e.g.,
        ``maxfilter -headpos <name>.pos``.

    Returns
    -------
    pos : array, shape (N, 10)
        The position and quaternion parameters from cHPI fitting.

    See Also
    --------
    write_head_pos
    head_pos_to_trans_rot_t

    Notes
    -----
    .. versionadded:: 0.12
    """
    _check_fname(fname, must_exist=True, overwrite="read")
    data = np.loadtxt(fname, skiprows=1)  # first line is header, skip it
    data.shape = (-1, 10)  # ensure it's the right size even if empty
    if np.isnan(data).any():  # make sure we didn't do something dumb
        raise RuntimeError(f"positions could not be read properly from {fname}")
    return data


def write_head_pos(fname, pos):
    """Write MaxFilter-formatted head position parameters.

    Parameters
    ----------
    fname : path-like
        The filename to write.
    pos : array, shape (N, 10)
        The position and quaternion parameters from cHPI fitting.

    See Also
    --------
    read_head_pos
    head_pos_to_trans_rot_t

    Notes
    -----
    .. versionadded:: 0.12
    """
    _check_fname(fname, overwrite=True)
    pos = np.array(pos, np.float64)
    if pos.ndim != 2 or pos.shape[1] != 10:
        raise ValueError("pos must be a 2D array of shape (N, 10)")
    with open(fname, "wb") as fid:
        fid.write(
            " Time       q1       q2       q3       q4       q5       "
            "q6       g-value  error    velocity\n".encode("ASCII")
        )
        for p in pos:
            fmts = ["% 9.3f"] + ["% 8.5f"] * 9
            fid.write(((" " + " ".join(fmts) + "\n") % tuple(p)).encode("ASCII"))


def head_pos_to_trans_rot_t(quats):
    """Convert Maxfilter-formatted head position quaternions.

    Parameters
    ----------
    quats : ndarray, shape (N, 10)
        MaxFilter-formatted position and quaternion parameters.

    Returns
    -------
    translation : ndarray, shape (N, 3)
        Translations at each time point.
    rotation : ndarray, shape (N, 3, 3)
        Rotations at each time point.
    t : ndarray, shape (N,)
        The time points.

    See Also
    --------
    read_head_pos
    write_head_pos
    """
    t = quats[..., 0].copy()
    rotation = quat_to_rot(quats[..., 1:4])
    translation = quats[..., 4:7].copy()
    return translation, rotation, t


@verbose
def extract_chpi_locs_ctf(raw, verbose=None):
    r"""Extract cHPI locations from CTF data.

    Parameters
    ----------
    raw : instance of Raw
        Raw data with CTF cHPI information.
    %(verbose)s

    Returns
    -------
    %(chpi_locs)s

    Notes
    -----
    CTF continuous head monitoring stores the x,y,z location (m) of each chpi
    coil as separate channels in the dataset:

    - ``HLC001[123]\\*`` - nasion
    - ``HLC002[123]\\*`` - lpa
    - ``HLC003[123]\\*`` - rpa

    This extracts these positions for use with
    :func:`~mne.chpi.compute_head_pos`.

    .. versionadded:: 0.20
    """
    # Pick channels corresponding to the cHPI positions
    hpi_picks = pick_channels_regexp(raw.info["ch_names"], "HLC00[123][123].*")

    # make sure we get 9 channels
    if len(hpi_picks) != 9:
        raise RuntimeError("Could not find all 9 cHPI channels")

    # get indices in alphabetical order
    sorted_picks = np.array(sorted(hpi_picks, key=lambda k: raw.info["ch_names"][k]))

    # make picks to match order of dig cardinial ident codes.
    # LPA (HPIC002[123]-*), NAS(HPIC001[123]-*), RPA(HPIC003[123]-*)
    hpi_picks = sorted_picks[[3, 4, 5, 0, 1, 2, 6, 7, 8]]
    del sorted_picks

    # process the entire run
    time_sl = slice(0, len(raw.times))
    chpi_data = raw[hpi_picks, time_sl][0]

    # transforms
    tmp_trans = _make_ctf_coord_trans_set(None, None)
    ctf_dev_dev_t = tmp_trans["t_ctf_dev_dev"]
    del tmp_trans

    # find indices where chpi locations change
    indices = [0]
    indices.extend(np.where(np.all(np.diff(chpi_data, axis=1), axis=0))[0] + 1)
    # data in channels are in ctf device coordinates (cm)
    rrs = chpi_data[:, indices].T.reshape(len(indices), 3, 3)  # m
    # map to mne device coords
    rrs = apply_trans(ctf_dev_dev_t, rrs)
    gofs = np.ones(rrs.shape[:2])  # not encoded, set all good
    moments = np.zeros(rrs.shape)  # not encoded, set all zero
    times = raw.times[indices] + raw._first_time
    return dict(rrs=rrs, gofs=gofs, times=times, moments=moments)


@verbose
def extract_chpi_locs_kit(raw, stim_channel="MISC 064", *, verbose=None):
    """Extract cHPI locations from KIT data.

    Parameters
    ----------
    raw : instance of RawKIT
        Raw data with KIT cHPI information.
    stim_channel : str
        The stimulus channel that encodes HPI measurement intervals.
    %(verbose)s

    Returns
    -------
    %(chpi_locs)s

    Notes
    -----
    .. versionadded:: 0.23
    """
    _validate_type(raw, (_RawKIT,), "raw")
    stim_chs = [
        raw.info["ch_names"][pick]
        for pick in pick_types(raw.info, stim=True, misc=True, ref_meg=False)
    ]
    _validate_type(stim_channel, str, "stim_channel")
    _check_option("stim_channel", stim_channel, stim_chs)
    idx = raw.ch_names.index(stim_channel)
    safe_false = _verbose_safe_false()
    events_on = find_events(
        raw, stim_channel=raw.ch_names[idx], output="onset", verbose=safe_false
    )[:, 0]
    events_off = find_events(
        raw, stim_channel=raw.ch_names[idx], output="offset", verbose=safe_false
    )[:, 0]
    bad = False
    if len(events_on) == 0 or len(events_off) == 0:
        bad = True
    else:
        if events_on[-1] > events_off[-1]:
            events_on = events_on[:-1]
        if events_on.size != events_off.size or not (events_on < events_off).all():
            bad = True
    if bad:
        raise RuntimeError(
            f"Could not find appropriate cHPI intervals from {stim_channel}"
        )
    # use the midpoint for times
    times = (events_on + events_off) / (2 * raw.info["sfreq"])
    del events_on, events_off
    # XXX remove first two rows. It is unknown currently if there is a way to
    # determine from the con file the number of initial pulses that
    # indicate the start of reading. The number is shown by opening the con
    # file in MEG160, but I couldn't find the value in the .con file, so it
    # may just always be 2...
    times = times[2:]
    n_coils = 5  # KIT always has 5 (hard-coded in reader)
    header = raw._raw_extras[0]["dirs"][KIT.DIR_INDEX_CHPI_DATA]
    dtype = np.dtype([("good", "<u4"), ("data", "<f8", (4,))])
    assert dtype.itemsize == header["size"], (dtype.itemsize, header["size"])
    all_data = list()
    for fname in raw._filenames:
        with open(fname) as fid:
            fid.seek(header["offset"])
            all_data.append(
                np.fromfile(fid, dtype, count=header["count"]).reshape(-1, n_coils)
            )
    data = np.concatenate(all_data)
    extra = ""
    if len(times) < len(data):
        extra = f", truncating to {len(times)} based on events"
    logger.info(f"Found {len(data)} cHPI measurement{_pl(len(data))}{extra}")
    data = data[: len(times)]
    # good is not currently used, but keep this in case we want it later
    # good = data['good'] == 1
    data = data["data"]
    rrs, gofs = data[:, :, :3], data[:, :, 3]
    rrs = apply_trans(als_ras_trans, rrs)
    moments = np.zeros(rrs.shape)  # not encoded, set all zero
    return dict(rrs=rrs, gofs=gofs, times=times, moments=moments)


# ############################################################################
# Estimate positions from data


@verbose
def get_chpi_info(info, on_missing="raise", verbose=None):
    """Retrieve cHPI information from the data.

    Parameters
    ----------
    %(info_not_none)s
    %(on_missing_chpi)s
    %(verbose)s

    Returns
    -------
    hpi_freqs : array, shape (n_coils,)
        The frequency used for each individual cHPI coil.
    hpi_pick : int | None
        The index of the ``STIM`` channel containing information about when
        which cHPI coils were switched on.
    hpi_on : array, shape (n_coils,)
        The values coding for the "on" state of each individual cHPI coil.

    Notes
    -----
    .. versionadded:: 0.24
    """
    _validate_type(item=info, item_name="info", types=Info)
    _check_option(
        parameter="on_missing",
        value=on_missing,
        allowed_values=["ignore", "raise", "warn"],
    )

    if len(info["hpi_meas"]) == 0 or (
        "coil_freq" not in info["hpi_meas"][0]["hpi_coils"][0]
    ):
        _on_missing(
            on_missing,
            msg="No appropriate cHPI information found in "
            'info["hpi_meas"] and info["hpi_subsystem"]',
        )
        return np.empty(0), None, np.empty(0)

    hpi_coils = sorted(
        info["hpi_meas"][-1]["hpi_coils"], key=lambda x: x["number"]
    )  # ascending (info) order

    # get frequencies
    hpi_freqs = np.array([float(x["coil_freq"]) for x in hpi_coils])
    logger.info(
        f"Using {len(hpi_freqs)} HPI coils: {' '.join(str(int(s)) for s in hpi_freqs)} "
        "Hz"
    )

    # how cHPI active is indicated in the FIF file
    hpi_sub = info["hpi_subsystem"]
    hpi_pick = None  # there is no pick!
    if hpi_sub is not None:
        if "event_channel" in hpi_sub:
            hpi_pick = pick_channels(
                info["ch_names"], [hpi_sub["event_channel"]], ordered=False
            )
            hpi_pick = hpi_pick[0] if len(hpi_pick) > 0 else None
        # grab codes indicating a coil is active
        hpi_on = [coil["event_bits"][0] for coil in hpi_sub["hpi_coils"]]
        # not all HPI coils will actually be used
        hpi_on = np.array([hpi_on[hc["number"] - 1] for hc in hpi_coils])
        # mask for coils that may be active
        hpi_mask = np.array([event_bit != 0 for event_bit in hpi_on])
        hpi_on = hpi_on[hpi_mask]
        hpi_freqs = hpi_freqs[hpi_mask]
    else:
        hpi_on = np.zeros(len(hpi_freqs))

    return hpi_freqs, hpi_pick, hpi_on


@verbose
def _get_hpi_initial_fit(info, adjust=False, verbose=None):
    """Get HPI fit locations from raw."""
    if info["hpi_results"] is None or len(info["hpi_results"]) == 0:
        raise RuntimeError("no initial cHPI head localization performed")

    hpi_result = info["hpi_results"][-1]
    hpi_dig = sorted(
        [d for d in info["dig"] if d["kind"] == FIFF.FIFFV_POINT_HPI],
        key=lambda x: x["ident"],
    )  # ascending (dig) order
    if len(hpi_dig) == 0:  # CTF data, probably
        msg = "HPIFIT: No HPI dig points, using hpifit result"
        hpi_dig = sorted(hpi_result["dig_points"], key=lambda x: x["ident"])
        if all(
            d["coord_frame"] in (FIFF.FIFFV_COORD_DEVICE, FIFF.FIFFV_COORD_UNKNOWN)
            for d in hpi_dig
        ):
            # Do not modify in place!
            hpi_dig = copy.deepcopy(hpi_dig)
            msg += " transformed to head coords"
            for dig in hpi_dig:
                dig.update(
                    r=apply_trans(info["dev_head_t"], dig["r"]),
                    coord_frame=FIFF.FIFFV_COORD_HEAD,
                )
        logger.debug(msg)

    # zero-based indexing, dig->info
    # CTF does not populate some entries so we use .get here
    pos_order = hpi_result.get("order", np.arange(1, len(hpi_dig) + 1)) - 1
    used = hpi_result.get("used", np.arange(len(hpi_dig)))
    dist_limit = hpi_result.get("dist_limit", 0.005)
    good_limit = hpi_result.get("good_limit", 0.98)
    goodness = hpi_result.get("goodness", np.ones(len(hpi_dig)))

    # this shouldn't happen, eventually we could add the transforms
    # necessary to put it in head coords
    if not all(d["coord_frame"] == FIFF.FIFFV_COORD_HEAD for d in hpi_dig):
        raise RuntimeError("cHPI coordinate frame incorrect")
    # Give the user some info
    logger.info(
        f"HPIFIT: {len(pos_order)} coils digitized in order "
        f"{' '.join(str(o + 1) for o in pos_order)}"
    )
    logger.debug(
        f"HPIFIT: {len(used)} coils accepted: {' '.join(str(h) for h in used)}"
    )
    hpi_rrs = np.array([d["r"] for d in hpi_dig])[pos_order]
    assert len(hpi_rrs) >= 3

    # Fitting errors
    hpi_rrs_fit = sorted(
        [d for d in info["hpi_results"][-1]["dig_points"]], key=lambda x: x["ident"]
    )
    hpi_rrs_fit = np.array([d["r"] for d in hpi_rrs_fit])
    # hpi_result['dig_points'] are in FIFFV_COORD_UNKNOWN coords, but this
    # is probably a misnomer because it should be FIFFV_COORD_DEVICE for this
    # to work
    assert hpi_result["coord_trans"]["to"] == FIFF.FIFFV_COORD_HEAD
    hpi_rrs_fit = apply_trans(hpi_result["coord_trans"]["trans"], hpi_rrs_fit)
    if "moments" in hpi_result:
        logger.debug(f"Hpi coil moments {hpi_result['moments'].shape[::-1]}:")
        for moment in hpi_result["moments"]:
            logger.debug(f"{moment[0]:g} {moment[1]:g} {moment[2]:g}")
    errors = np.linalg.norm(hpi_rrs - hpi_rrs_fit, axis=1)
    logger.debug(f"HPIFIT errors:  {', '.join(f'{1000 * e:0.1f}' for e in errors)} mm.")
    if errors.sum() < len(errors) * dist_limit:
        logger.info("HPI consistency of isotrak and hpifit is OK.")
    elif not adjust and (len(used) == len(hpi_dig)):
        warn("HPI consistency of isotrak and hpifit is poor.")
    else:
        # adjust HPI coil locations using the hpifit transformation
        for hi, (err, r_fit) in enumerate(zip(errors, hpi_rrs_fit)):
            # transform to head frame
            d = 1000 * err
            if not adjust:
                if err >= dist_limit:
                    warn(
                        f"Discrepancy of HPI coil {hi + 1} isotrak and hpifit is "
                        f"{d:.1f} mm!"
                    )
            elif hi + 1 not in used:
                if goodness[hi] >= good_limit:
                    logger.info(
                        f"Note: HPI coil {hi + 1} isotrak is adjusted by {d:.1f} mm!"
                    )
                    hpi_rrs[hi] = r_fit
                else:
                    warn(
                        f"Discrepancy of HPI coil {hi + 1} isotrak and hpifit of "
                        f"{d:.1f} mm was not adjusted!"
                    )
    logger.debug(
        f"HP fitting limits: err = {1000 * dist_limit:.1f} mm, gval = {good_limit:.3f}."
    )

    return hpi_rrs.astype(float)


def _magnetic_dipole_objective(
    x, B, B2, coils, whitener, too_close, return_moment=False
):
    """Project data onto right eigenvectors of whitened forward."""
    fwd = _magnetic_dipole_field_vec(x[np.newaxis], coils, too_close)
    out, u, s, one = _magnetic_dipole_delta(fwd, whitener, B, B2)
    if return_moment:
        one /= s
        Q = np.dot(one, u.T)
        out = (out, Q)
    return out


@jit()
def _magnetic_dipole_delta(fwd, whitener, B, B2):
    # Here we use .T to get whitener to Fortran order, which speeds things up
    fwd = np.dot(fwd, whitener.T)
    u, s, v = np.linalg.svd(fwd, full_matrices=False)
    one = np.dot(v, B)
    Bm2 = np.dot(one, one)
    return B2 - Bm2, u, s, one


def _magnetic_dipole_delta_multi(whitened_fwd_svd, B, B2):
    # Here we use .T to get whitener to Fortran order, which speeds things up
    one = np.matmul(whitened_fwd_svd, B)
    Bm2 = np.sum(one * one, axis=1)
    return B2 - Bm2


def _fit_magnetic_dipole(B_orig, x0, too_close, whitener, coils, guesses):
    """Fit a single bit of data (x0 = pos)."""
    B = np.dot(whitener, B_orig)
    B2 = np.dot(B, B)
    objective = partial(
        _magnetic_dipole_objective,
        B=B,
        B2=B2,
        coils=coils,
        whitener=whitener,
        too_close=too_close,
    )
    if guesses is not None:
        res0 = objective(x0)
        res = _magnetic_dipole_delta_multi(guesses["whitened_fwd_svd"], B, B2)
        assert res.shape == (guesses["rr"].shape[0],)
        idx = np.argmin(res)
        if res[idx] < res0:
            x0 = guesses["rr"][idx]
    x = fmin_cobyla(objective, x0, (), rhobeg=1e-3, rhoend=1e-5, disp=False)
    gof, moment = objective(x, return_moment=True)
    gof = 1.0 - gof / B2
    return x, gof, moment


@jit()
def _chpi_objective(x, coil_dev_rrs, coil_head_rrs):
    """Compute objective function."""
    d = np.dot(coil_dev_rrs, quat_to_rot(x[:3]).T)
    d += x[3:]
    d -= coil_head_rrs
    d *= d
    return d.sum()


def _fit_chpi_quat(coil_dev_rrs, coil_head_rrs):
    """Fit rotation and translation (quaternion) parameters for cHPI coils."""
    denom = np.linalg.norm(coil_head_rrs - np.mean(coil_head_rrs, axis=0))
    denom *= denom
    # We could try to solve it the analytic way:
    # XXX someday we could choose to weight these points by their goodness
    # of fit somehow.
    quat = _fit_matched_points(coil_dev_rrs, coil_head_rrs)[0]
    gof = 1.0 - _chpi_objective(quat, coil_dev_rrs, coil_head_rrs) / denom
    return quat, gof


def _fit_coil_order_dev_head_trans(dev_pnts, head_pnts, bias=True):
    """Compute Device to Head transform allowing for permutiatons of points."""
    id_quat = np.zeros(6)
    best_order = None
    best_g = -999
    best_quat = id_quat
    for this_order in itertools.permutations(np.arange(len(head_pnts))):
        head_pnts_tmp = head_pnts[np.array(this_order)]
        this_quat, g = _fit_chpi_quat(dev_pnts, head_pnts_tmp)
        assert np.linalg.det(quat_to_rot(this_quat[:3])) > 0.9999
        if bias:
            # For symmetrical arrangements, flips can produce roughly
            # equivalent g values. To avoid this, heavily penalize
            # large rotations.
            rotation = _angle_between_quats(this_quat[:3], np.zeros(3))
            check_g = g * max(1.0 - rotation / np.pi, 0) ** 0.25
        else:
            check_g = g
        if check_g > best_g:
            out_g = g
            best_g = check_g
            best_order = np.array(this_order)
            best_quat = this_quat

    # Convert Quaterion to transform
    dev_head_t = _quat_to_affine(best_quat)
    return dev_head_t, best_order, out_g


@verbose
def _setup_hpi_amplitude_fitting(
    info, t_window, remove_aliased=False, ext_order=1, allow_empty=False, verbose=None
):
    """Generate HPI structure for HPI localization."""
    # grab basic info.
    on_missing = "raise" if not allow_empty else "ignore"
    hpi_freqs, hpi_pick, hpi_ons = get_chpi_info(info, on_missing=on_missing)

    # check for maxwell filtering
    for ent in info["proc_history"]:
        for key in ("sss_info", "max_st"):
            if len(ent["max_info"]["sss_info"]) > 0:
                warn(
                    "Fitting cHPI amplitudes after Maxwell filtering may not work, "
                    "consider fitting on the original data."
                )
                break

    _validate_type(t_window, (str, "numeric"), "t_window")
    if info["line_freq"] is not None:
        line_freqs = np.arange(
            info["line_freq"], info["sfreq"] / 3.0, info["line_freq"]
        )
    else:
        line_freqs = np.zeros([0])
    lfs = " ".join(f"{lf}" for lf in line_freqs)
    logger.info(f"Line interference frequencies: {lfs} Hz")
    # worry about resampled/filtered data.
    # What to do e.g. if Raw has been resampled and some of our
    # HPI freqs would now be aliased
    highest = info.get("lowpass")
    highest = info["sfreq"] / 2.0 if highest is None else highest
    keepers = hpi_freqs <= highest
    if remove_aliased:
        hpi_freqs = hpi_freqs[keepers]
        hpi_ons = hpi_ons[keepers]
    elif not keepers.all():
        raise RuntimeError(
            f"Found HPI frequencies {hpi_freqs[~keepers].tolist()} above the lowpass ("
            f"or Nyquist) frequency {highest:0.1f}"
        )
    # calculate optimal window length.
    if isinstance(t_window, str):
        _check_option("t_window", t_window, ("auto",), extra="if a string")
        if len(hpi_freqs):
            all_freqs = np.concatenate((hpi_freqs, line_freqs))
            delta_freqs = np.diff(np.unique(all_freqs))
            t_window = max(5.0 / all_freqs.min(), 1.0 / delta_freqs.min())
        else:
            t_window = 0.2
    t_window = float(t_window)
    if t_window <= 0:
        raise ValueError(f"t_window ({t_window}) must be > 0")
    logger.info(f"Using time window: {1000 * t_window:0.1f} ms")
    window_nsamp = np.rint(t_window * info["sfreq"]).astype(int)
    model = _setup_hpi_glm(hpi_freqs, line_freqs, info["sfreq"], window_nsamp)
    inv_model = np.linalg.pinv(model)
    inv_model_reord = _reorder_inv_model(inv_model, len(hpi_freqs))
    proj, proj_op, meg_picks = _setup_ext_proj(info, ext_order)
    # include mag and grad picks separately, for SNR computations
    mag_subpicks = _picks_to_idx(info, "mag", allow_empty=True)
    mag_subpicks = np.searchsorted(meg_picks, mag_subpicks)
    grad_subpicks = _picks_to_idx(info, "grad", allow_empty=True)
    grad_subpicks = np.searchsorted(meg_picks, grad_subpicks)
    # Set up magnetic dipole fits
    hpi = dict(
        meg_picks=meg_picks,
        mag_subpicks=mag_subpicks,
        grad_subpicks=grad_subpicks,
        hpi_pick=hpi_pick,
        model=model,
        inv_model=inv_model,
        t_window=t_window,
        inv_model_reord=inv_model_reord,
        on=hpi_ons,
        n_window=window_nsamp,
        proj=proj,
        proj_op=proj_op,
        freqs=hpi_freqs,
        line_freqs=line_freqs,
    )
    return hpi


def _setup_hpi_glm(hpi_freqs, line_freqs, sfreq, window_nsamp):
    """Initialize a general linear model for HPI amplitude estimation."""
    slope = np.linspace(-0.5, 0.5, window_nsamp)[:, np.newaxis]
    radians_per_sec = 2 * np.pi * np.arange(window_nsamp, dtype=float) / sfreq
    f_t = hpi_freqs[np.newaxis, :] * radians_per_sec[:, np.newaxis]
    l_t = line_freqs[np.newaxis, :] * radians_per_sec[:, np.newaxis]
    model = [
        np.sin(f_t),
        np.cos(f_t),  # hpi freqs
        np.sin(l_t),
        np.cos(l_t),  # line freqs
        slope,
        np.ones_like(slope),
    ]  # drift, DC
    return np.hstack(model)


@jit()
def _reorder_inv_model(inv_model, n_freqs):
    # Reorder for faster computation
    idx = np.arange(2 * n_freqs).reshape(2, n_freqs).T.ravel()
    return inv_model[idx]


def _setup_ext_proj(info, ext_order):
    meg_picks = pick_types(info, meg=True, eeg=False, exclude="bads")
    info = pick_info(_simplify_info(info), meg_picks)  # makes a copy
    _, _, _, _, mag_or_fine = _get_mf_picks_fix_mags(
        info, int_order=0, ext_order=ext_order, ignore_ref=True, verbose="error"
    )
    mf_coils = _prep_mf_coils(info, verbose="error")
    ext = _sss_basis(
        dict(origin=(0.0, 0.0, 0.0), int_order=0, ext_order=ext_order), mf_coils
    ).T
    out_removes = _regularize_out(0, 1, mag_or_fine, [])
    ext = ext[~np.isin(np.arange(len(ext)), out_removes)]
    ext = orth(ext.T).T
    assert ext.shape[1] == len(meg_picks)
    proj = Projection(
        kind=FIFF.FIFFV_PROJ_ITEM_HOMOG_FIELD,
        desc="SSS",
        active=False,
        data=dict(
            data=ext, ncol=info["nchan"], col_names=info["ch_names"], nrow=len(ext)
        ),
    )
    with info._unlock():
        info["projs"] = [proj]
    proj_op, _ = setup_proj(
        info, add_eeg_ref=False, activate=False, verbose=_verbose_safe_false()
    )
    assert proj_op.shape == (len(meg_picks),) * 2
    return proj, proj_op, meg_picks


def _time_prefix(fit_time):
    """Format log messages."""
    return (f"    t={fit_time:0.3f}:").ljust(17)


def _fit_chpi_amplitudes(raw, time_sl, hpi, snr=False):
    """Fit amplitudes for each channel from each of the N cHPI sinusoids.

    Returns
    -------
    sin_fit : ndarray, shape (n_freqs, n_channels)
        The sin amplitudes matching each cHPI frequency.
        Will be all nan if this time window should be skipped.
    snr : ndarray, shape (n_freqs, 2)
        Estimated SNR for this window, separately for mag and grad channels.
    """
    # No need to detrend the data because our model has a DC term
    with use_log_level(False):
        # loads good channels
        this_data = raw[hpi["meg_picks"], time_sl][0]

    # which HPI coils to use
    if hpi["hpi_pick"] is not None:
        with use_log_level(False):
            # loads hpi_stim channel
            chpi_data = raw[hpi["hpi_pick"], time_sl][0]

        ons = (np.round(chpi_data).astype(np.int64) & hpi["on"][:, np.newaxis]).astype(
            bool
        )
        n_on = ons.all(axis=-1).sum(axis=0)
        if not (n_on >= 3).all():
            return None
    if snr:
        return _fast_fit_snr(
            this_data,
            len(hpi["freqs"]),
            hpi["model"],
            hpi["inv_model"],
            hpi["mag_subpicks"],
            hpi["grad_subpicks"],
        )
    return _fast_fit(
        this_data,
        hpi["proj_op"],
        len(hpi["freqs"]),
        hpi["model"],
        hpi["inv_model_reord"],
    )


@jit()
def _fast_fit(this_data, proj, n_freqs, model, inv_model_reord):
    # first or last window
    if this_data.shape[1] != model.shape[0]:
        model = model[: this_data.shape[1]]
        inv_model_reord = _reorder_inv_model(np.linalg.pinv(model), n_freqs)
    proj_data = proj @ this_data
    X = inv_model_reord @ proj_data.T

    sin_fit = np.zeros((n_freqs, X.shape[1]))
    for fi in range(n_freqs):
        # use SVD across all sensors to estimate the sinusoid phase
        u, s, vt = np.linalg.svd(X[2 * fi : 2 * fi + 2], full_matrices=False)
        # the first component holds the predominant phase direction
        # (so ignore the second, effectively doing s[1] = 0):
        sin_fit[fi] = vt[0] * s[0]
    return sin_fit


@jit()
def _fast_fit_snr(this_data, n_freqs, model, inv_model, mag_picks, grad_picks):
    # first or last window
    if this_data.shape[1] != model.shape[0]:
        model = model[: this_data.shape[1]]
        inv_model = np.linalg.pinv(model)
    coefs = np.ascontiguousarray(inv_model) @ np.ascontiguousarray(this_data.T)
    # average sin & cos terms (special property of sinusoids: power=A²/2)
    hpi_power = (coefs[:n_freqs] ** 2 + coefs[n_freqs : (2 * n_freqs)] ** 2) / 2
    resid = this_data - np.ascontiguousarray((model @ coefs).T)
    # can't use np.var(..., axis=1) with Numba, so do it manually:
    resid_mean = np.atleast_2d(resid.sum(axis=1) / resid.shape[1]).T
    squared_devs = np.abs(resid - resid_mean) ** 2
    resid_var = squared_devs.sum(axis=1) / squared_devs.shape[1]
    # output array will be (n_freqs, 3 * n_ch_types). The 3 columns for each
    # channel type are the SNR, the mean cHPI power and the residual variance
    # (which gets tiled to shape (n_freqs,) because it's a scalar).
    snrs = np.empty((n_freqs, 0))
    # average power & compute residual variance separately for each ch type
    for _picks in (mag_picks, grad_picks):
        if len(_picks):
            avg_power = hpi_power[:, _picks].sum(axis=1) / len(_picks)
            avg_resid = resid_var[_picks].mean() * np.ones(n_freqs)
            snr = 10 * np.log10(avg_power / avg_resid)
            snrs = np.hstack((snrs, np.stack((snr, avg_power, avg_resid), 1)))
    return snrs


def _check_chpi_param(chpi_, name):
    if name == "chpi_locs":
        want_ndims = dict(times=1, rrs=3, moments=3, gofs=2)
        extra_keys = list()
    else:
        assert name == "chpi_amplitudes"
        want_ndims = dict(times=1, slopes=3)
        extra_keys = ["proj"]

    _validate_type(chpi_, dict, name)
    want_keys = list(want_ndims.keys()) + extra_keys
    if set(want_keys).symmetric_difference(chpi_):
        raise ValueError(
            f"{name} must be a dict with entries {want_keys}, got "
            f"{sorted(chpi_.keys())}"
        )
    n_times = None
    for key, want_ndim in want_ndims.items():
        key_str = f"{name}[{key}]"
        val = chpi_[key]
        _validate_type(val, np.ndarray, key_str)
        shape = val.shape
        if val.ndim != want_ndim:
            raise ValueError(f"{key_str} must have ndim={want_ndim}, got {val.ndim}")
        if n_times is None and key != "proj":
            n_times = shape[0]
        if n_times != shape[0] and key != "proj":
            raise ValueError(
                f"{name} have inconsistent number of time points in {want_keys}"
            )
    if name == "chpi_locs":
        n_coils = chpi_["rrs"].shape[1]
        for key in ("gofs", "moments"):
            val = chpi_[key]
            if val.shape[1] != n_coils:
                raise ValueError(
                    f'chpi_locs["rrs"] had values for {n_coils} coils but '
                    f'chpi_locs["{key}"] had values for {val.shape[1]} coils'
                )
        for key in ("rrs", "moments"):
            val = chpi_[key]
            if val.shape[2] != 3:
                raise ValueError(
                    f'chpi_locs["{key}"].shape[2] must be 3, got shape {shape}'
                )
    else:
        assert name == "chpi_amplitudes"
        slopes, proj = chpi_["slopes"], chpi_["proj"]
        _validate_type(proj, Projection, 'chpi_amplitudes["proj"]')
        n_ch = len(proj["data"]["col_names"])
        if slopes.shape[0] != n_times or slopes.shape[2] != n_ch:
            raise ValueError(
                f"slopes must have shape[0]=={n_times} and shape[2]=={n_ch}, got shape "
                f"{slopes.shape}"
            )


@verbose
def compute_head_pos(
    info, chpi_locs, dist_limit=0.005, gof_limit=0.98, adjust_dig=False, verbose=None
):
    """Compute time-varying head positions.

    Parameters
    ----------
    %(info_not_none)s
    %(chpi_locs)s
        Typically obtained by :func:`~mne.chpi.compute_chpi_locs` or
        :func:`~mne.chpi.extract_chpi_locs_ctf`.
    dist_limit : float
        Minimum distance (m) to accept for coil position fitting.
    gof_limit : float
        Minimum goodness of fit to accept for each coil.
    %(adjust_dig_chpi)s
    %(verbose)s

    Returns
    -------
    quats : ndarray, shape (n_pos, 10)
        The ``[t, q1, q2, q3, x, y, z, gof, err, v]`` for each fit.

    See Also
    --------
    compute_chpi_locs
    extract_chpi_locs_ctf
    read_head_pos
    write_head_pos

    Notes
    -----
    .. versionadded:: 0.20
    """
    _check_chpi_param(chpi_locs, "chpi_locs")
    _validate_type(info, Info, "info")
    hpi_dig_head_rrs = _get_hpi_initial_fit(info, adjust=adjust_dig, verbose="error")
    n_coils = len(hpi_dig_head_rrs)
    coil_dev_rrs = apply_trans(invert_transform(info["dev_head_t"]), hpi_dig_head_rrs)
    dev_head_t = info["dev_head_t"]["trans"]
    pos_0 = dev_head_t[:3, 3]
    last = dict(
        quat_fit_time=-0.1,
        coil_dev_rrs=coil_dev_rrs,
        quat=np.concatenate([rot_to_quat(dev_head_t[:3, :3]), dev_head_t[:3, 3]]),
    )
    del coil_dev_rrs
    quats = []
    for fit_time, this_coil_dev_rrs, g_coils in zip(
        *(chpi_locs[key] for key in ("times", "rrs", "gofs"))
    ):
        use_idx = np.where(g_coils >= gof_limit)[0]

        #
        # 1. Check number of good ones
        #
        if len(use_idx) < 3:
            gofs = ", ".join(f"{g:0.2f}" for g in g_coils)
            warn(
                f"{_time_prefix(fit_time)}{len(use_idx)}/{n_coils} "
                "good HPI fits, cannot determine the transformation "
                f"({gofs} GOF)!"
            )
            continue

        #
        # 2. Fit the head translation and rotation params (minimize error
        #    between coil positions and the head coil digitization
        #    positions) iteratively using different sets of coils.
        #
        this_quat, g, use_idx = _fit_chpi_quat_subset(
            this_coil_dev_rrs, hpi_dig_head_rrs, use_idx
        )

        #
        # 3. Stop if < 3 good
        #

        # Convert quaterion to transform
        this_dev_head_t = _quat_to_affine(this_quat)
        est_coil_head_rrs = apply_trans(this_dev_head_t, this_coil_dev_rrs)
        errs = np.linalg.norm(hpi_dig_head_rrs - est_coil_head_rrs, axis=1)
        n_good = ((g_coils >= gof_limit) & (errs < dist_limit)).sum()
        if n_good < 3:
            warn_str = ", ".join(
                f"{1000 * e:0.1f}::{g:0.2f}" for e, g in zip(errs, g_coils)
            )
            warn(
                f"{_time_prefix(fit_time)}{n_good}/{n_coils} good HPI fits, cannot "
                f"determine the transformation ({warn_str} mm/GOF)!"
            )
            continue

        # velocities, in device coords, of HPI coils
        dt = fit_time - last["quat_fit_time"]
        vs = tuple(
            1000.0
            * np.linalg.norm(last["coil_dev_rrs"] - this_coil_dev_rrs, axis=1)
            / dt
        )
        logger.info(
            _time_prefix(fit_time)
            + (
                "%s/%s good HPI fits, movements [mm/s] = "
                + " / ".join(["% 8.1f"] * n_coils)
            )
            % ((n_good, n_coils) + vs)
        )

        # Log results
        # MaxFilter averages over a 200 ms window for display, but we don't
        for ii in range(n_coils):
            if ii in use_idx:
                start, end = " ", "/"
            else:
                start, end = "(", ")"
            log_str = (
                "    "
                + start
                + "{0:6.1f} {1:6.1f} {2:6.1f} / "
                + "{3:6.1f} {4:6.1f} {5:6.1f} / "
                + "g = {6:0.3f} err = {7:4.1f} "
                + end
            )
            vals = np.concatenate(
                (
                    1000 * hpi_dig_head_rrs[ii],
                    1000 * est_coil_head_rrs[ii],
                    [g_coils[ii], 1000 * errs[ii]],
                )
            )
            if len(use_idx) >= 3:
                if ii <= 2:
                    log_str += "{8:6.3f} {9:6.3f} {10:6.3f}"
                    vals = np.concatenate((vals, this_dev_head_t[ii, :3]))
                elif ii == 3:
                    log_str += "{8:6.1f} {9:6.1f} {10:6.1f}"
                    vals = np.concatenate((vals, this_dev_head_t[:3, 3] * 1000.0))
            logger.debug(log_str.format(*vals))

        # resulting errors in head coil positions
        d = np.linalg.norm(last["quat"][3:] - this_quat[3:])  # m
        r = _angle_between_quats(last["quat"][:3], this_quat[:3]) / dt
        v = d / dt  # m/s
        d = 100 * np.linalg.norm(this_quat[3:] - pos_0)  # dis from 1st
        logger.debug(
            f"    #t = {fit_time:0.3f}, #e = {100 * errs.mean():0.2f} cm, #g = {g:0.3f}"
            f", #v = {100 * v:0.2f} cm/s, #r = {r:0.2f} rad/s, #d = {d:0.2f} cm"
        )
        q_rep = " ".join(f"{qq:8.5f}" for qq in this_quat)
        logger.debug(f"    #t = {fit_time:0.3f}, #q = {q_rep}")

        quats.append(
            np.concatenate(([fit_time], this_quat, [g], [errs[use_idx].mean()], [v]))
        )
        last["quat_fit_time"] = fit_time
        last["quat"] = this_quat
        last["coil_dev_rrs"] = this_coil_dev_rrs
    quats = np.array(quats, np.float64)
    quats = np.zeros((0, 10)) if quats.size == 0 else quats
    return quats


def _fit_chpi_quat_subset(coil_dev_rrs, coil_head_rrs, use_idx):
    quat, g = _fit_chpi_quat(coil_dev_rrs[use_idx], coil_head_rrs[use_idx])
    out_idx = use_idx.copy()
    if len(use_idx) > 3:  # try dropping one (recursively)
        for di in range(len(use_idx)):
            this_use_idx = list(use_idx[:di]) + list(use_idx[di + 1 :])
            this_quat, this_g, this_use_idx = _fit_chpi_quat_subset(
                coil_dev_rrs, coil_head_rrs, this_use_idx
            )
            if this_g > g:
                quat, g, out_idx = this_quat, this_g, this_use_idx
    return quat, g, np.array(out_idx, int)


@verbose
def compute_chpi_snr(
    raw, t_step_min=0.01, t_window="auto", ext_order=1, tmin=0, tmax=None, verbose=None
):
    """Compute time-varying estimates of cHPI SNR.

    Parameters
    ----------
    raw : instance of Raw
        Raw data with cHPI information.
    t_step_min : float
        Minimum time step to use.
    %(t_window_chpi_t)s
    %(ext_order_chpi)s
    %(tmin_raw)s
    %(tmax_raw)s
    %(verbose)s

    Returns
    -------
    chpi_snrs : dict
        The time-varying cHPI SNR estimates, with entries "times", "freqs",
        "snr_mag", "power_mag", and "resid_mag" (and/or "snr_grad",
        "power_grad", and "resid_grad", depending on which channel types are
        present in ``raw``).

    See Also
    --------
    mne.chpi.compute_chpi_locs, mne.chpi.compute_chpi_amplitudes

    Notes
    -----
    .. versionadded:: 0.24
    """
    return _compute_chpi_amp_or_snr(
        raw, t_step_min, t_window, ext_order, tmin, tmax, verbose, snr=True
    )


@verbose
def compute_chpi_amplitudes(
    raw, t_step_min=0.01, t_window="auto", ext_order=1, tmin=0, tmax=None, verbose=None
):
    """Compute time-varying cHPI amplitudes.

    Parameters
    ----------
    raw : instance of Raw
        Raw data with cHPI information.
    t_step_min : float
        Minimum time step to use.
    %(t_window_chpi_t)s
    %(ext_order_chpi)s
    %(tmin_raw)s
    %(tmax_raw)s
    %(verbose)s

    Returns
    -------
    %(chpi_amplitudes)s

    See Also
    --------
    mne.chpi.compute_chpi_locs, mne.chpi.compute_chpi_snr

    Notes
    -----
    This function will:

    1. Get HPI frequencies,  HPI status channel, HPI status bits,
       and digitization order using ``_setup_hpi_amplitude_fitting``.
    2. Window data using ``t_window`` (half before and half after ``t``) and
       ``t_step_min``.
    3. Use a linear model (DC + linear slope + sin + cos terms) to fit
       sinusoidal amplitudes to MEG channels.
       It uses SVD to determine the phase/amplitude of the sinusoids.

    In "auto" mode, ``t_window`` will be set to the longer of:

    1. Five cycles of the lowest HPI or line frequency.
          Ensures that the frequency estimate is stable.
    2. The reciprocal of the smallest difference between HPI and line freqs.
          Ensures that neighboring frequencies can be disambiguated.

    The output is meant to be used with :func:`~mne.chpi.compute_chpi_locs`.

    .. versionadded:: 0.20
    """
    return _compute_chpi_amp_or_snr(
        raw, t_step_min, t_window, ext_order, tmin, tmax, verbose
    )


def _compute_chpi_amp_or_snr(
    raw,
    t_step_min=0.01,
    t_window="auto",
    ext_order=1,
    tmin=0,
    tmax=None,
    verbose=None,
    snr=False,
):
    """Compute cHPI amplitude or SNR.

    See compute_chpi_amplitudes for parameter descriptions. One additional
    boolean parameter ``snr`` signals whether to return SNR instead of
    amplitude.
    """
    hpi = _setup_hpi_amplitude_fitting(raw.info, t_window, ext_order=ext_order)
    tmin, tmax = raw._tmin_tmax_to_start_stop(tmin, tmax)
    tmin = tmin / raw.info["sfreq"]
    tmax = tmax / raw.info["sfreq"]
    need_win = hpi["t_window"] / 2.0
    fit_idxs = raw.time_as_index(
        np.arange(tmin + need_win, tmax, t_step_min), use_rounding=True
    )
    logger.info(
        f"Fitting {len(hpi['freqs'])} HPI coil locations at up to "
        f"{len(fit_idxs)} time points ({tmax - tmin:.1f} s duration)"
    )
    del tmin, tmax
    sin_fits = dict()
    sin_fits["proj"] = hpi["proj"]
    sin_fits["times"] = (
        np.round(fit_idxs + raw.first_samp - hpi["n_window"] / 2.0) / raw.info["sfreq"]
    )
    n_times = len(sin_fits["times"])
    n_freqs = len(hpi["freqs"])
    n_chans = len(sin_fits["proj"]["data"]["col_names"])
    if snr:
        del sin_fits["proj"]
        sin_fits["freqs"] = hpi["freqs"]
        ch_types = raw.get_channel_types()
        grad_offset = 3 if "mag" in ch_types else 0
        for ch_type in ("mag", "grad"):
            if ch_type in ch_types:
                for key in ("snr", "power", "resid"):
                    cols = 1 if key == "resid" else n_freqs
                    sin_fits[f"{ch_type}_{key}"] = np.empty((n_times, cols))
    else:
        sin_fits["slopes"] = np.empty((n_times, n_freqs, n_chans))
    message = f"cHPI {'SNRs' if snr else 'amplitudes'}"
    for mi, midpt in enumerate(ProgressBar(fit_idxs, mesg=message)):
        #
        # 0. determine samples to fit.
        #
        time_sl = midpt - hpi["n_window"] // 2
        time_sl = slice(max(time_sl, 0), min(time_sl + hpi["n_window"], len(raw.times)))

        #
        # 1. Fit amplitudes for each channel from each of the N sinusoids
        #
        amps_or_snrs = _fit_chpi_amplitudes(raw, time_sl, hpi, snr)
        if snr:
            if amps_or_snrs is None:
                amps_or_snrs = np.full((n_freqs, grad_offset + 3), np.nan)
            # unpack the SNR estimates. mag & grad are returned in one array
            # (because of Numba) so take care with which column is which.
            # note that mean residual is a scalar (same for all HPI freqs) but
            # is returned as a (tiled) vector (again, because Numba) so that's
            # why below we take amps_or_snrs[0, 2] instead of [:, 2]
            ch_types = raw.get_channel_types()
            if "mag" in ch_types:
                sin_fits["mag_snr"][mi] = amps_or_snrs[:, 0]  # SNR
                sin_fits["mag_power"][mi] = amps_or_snrs[:, 1]  # mean power
                sin_fits["mag_resid"][mi] = amps_or_snrs[0, 2]  # mean resid
            if "grad" in ch_types:
                sin_fits["grad_snr"][mi] = amps_or_snrs[:, grad_offset]
                sin_fits["grad_power"][mi] = amps_or_snrs[:, grad_offset + 1]
                sin_fits["grad_resid"][mi] = amps_or_snrs[0, grad_offset + 2]
        else:
            sin_fits["slopes"][mi] = amps_or_snrs
    return sin_fits


@verbose
def compute_chpi_locs(
    info,
    chpi_amplitudes,
    t_step_max=1.0,
    too_close="raise",
    adjust_dig=False,
    verbose=None,
):
    """Compute locations of each cHPI coils over time.

    Parameters
    ----------
    %(info_not_none)s
    %(chpi_amplitudes)s
        Typically obtained by :func:`mne.chpi.compute_chpi_amplitudes`.
    t_step_max : float
        Maximum time step to use.
    too_close : str
        How to handle HPI positions too close to the sensors,
        can be ``'raise'`` (default), ``'warning'``, or ``'info'``.
    %(adjust_dig_chpi)s
    %(verbose)s

    Returns
    -------
    %(chpi_locs)s

    See Also
    --------
    compute_chpi_amplitudes
    compute_head_pos
    read_head_pos
    write_head_pos
    extract_chpi_locs_ctf

    Notes
    -----
    This function is designed to take the output of
    :func:`mne.chpi.compute_chpi_amplitudes` and:

    1. Get HPI coil locations (as digitized in ``info['dig']``) in head coords.
    2. If the amplitudes are 98%% correlated with last position
       (and Δt < t_step_max), skip fitting.
    3. Fit magnetic dipoles using the amplitudes for each coil frequency.

    The number of fitted points ``n_pos`` will depend on the velocity of head
    movements as well as ``t_step_max`` (and ``t_step_min`` from
    :func:`mne.chpi.compute_chpi_amplitudes`).

    .. versionadded:: 0.20
    """
    # Set up magnetic dipole fits
    _check_option("too_close", too_close, ["raise", "warning", "info"])
    _check_chpi_param(chpi_amplitudes, "chpi_amplitudes")
    _validate_type(info, Info, "info")
    sin_fits = chpi_amplitudes  # use the old name below
    del chpi_amplitudes
    proj = sin_fits["proj"]
    meg_picks = pick_channels(info["ch_names"], proj["data"]["col_names"], ordered=True)
    info = pick_info(info, meg_picks)  # makes a copy
    with info._unlock():
        info["projs"] = [proj]
    del meg_picks, proj
    meg_coils = _concatenate_coils(_create_meg_coils(info["chs"], "accurate"))

    # Set up external model for interference suppression
    safe_false = _verbose_safe_false()
    cov = make_ad_hoc_cov(info, verbose=safe_false)
    whitener, _ = compute_whitener(cov, info, verbose=safe_false)

    # Make some location guesses (1 cm grid)
    R = np.linalg.norm(meg_coils[0], axis=1).min()
    guesses = _make_guesses(
        dict(R=R, r0=np.zeros(3)), 0.01, 0.0, 0.005, verbose=safe_false
    )[0]["rr"]
    logger.info(
        f"Computing {len(guesses)} HPI location guesses "
        f"(1 cm grid in a {R * 100:.1f} cm sphere)"
    )
    fwd = _magnetic_dipole_field_vec(guesses, meg_coils, too_close)
    fwd = np.dot(fwd, whitener.T)
    fwd.shape = (guesses.shape[0], 3, -1)
    fwd = np.linalg.svd(fwd, full_matrices=False)[2]
    guesses = dict(rr=guesses, whitened_fwd_svd=fwd)
    del fwd, R

    iter_ = list(zip(sin_fits["times"], sin_fits["slopes"]))
    chpi_locs = dict(times=[], rrs=[], gofs=[], moments=[])
    # setup last iteration structure
    hpi_dig_dev_rrs = apply_trans(
        invert_transform(info["dev_head_t"])["trans"],
        _get_hpi_initial_fit(info, adjust=adjust_dig),
    )
    last = dict(
        sin_fit=None,
        coil_fit_time=sin_fits["times"][0] - 1,
        coil_dev_rrs=hpi_dig_dev_rrs,
    )
    n_hpi = len(hpi_dig_dev_rrs)
    del hpi_dig_dev_rrs
    for fit_time, sin_fit in ProgressBar(iter_, mesg="cHPI locations "):
        # skip this window if bad
        if not np.isfinite(sin_fit).all():
            continue

        # check if data has sufficiently changed
        if last["sin_fit"] is not None:  # first iteration
            corrs = np.array(
                [np.corrcoef(s, lst)[0, 1] for s, lst in zip(sin_fit, last["sin_fit"])]
            )
            corrs *= corrs
            # check to see if we need to continue
            if (
                fit_time - last["coil_fit_time"] <= t_step_max - 1e-7
                and (corrs > 0.98).sum() >= 3
            ):
                # don't need to refit data
                continue

        # update 'last' sin_fit *before* inplace sign mult
        last["sin_fit"] = sin_fit.copy()

        #
        # 2. Fit magnetic dipole for each coil to obtain coil positions
        #    in device coordinates
        #
        coil_fits = [
            _fit_magnetic_dipole(f, x0, too_close, whitener, meg_coils, guesses)
            for f, x0 in zip(sin_fit, last["coil_dev_rrs"])
        ]
        rrs, gofs, moments = zip(*coil_fits)
        chpi_locs["times"].append(fit_time)
        chpi_locs["rrs"].append(rrs)
        chpi_locs["gofs"].append(gofs)
        chpi_locs["moments"].append(moments)
        last["coil_fit_time"] = fit_time
        last["coil_dev_rrs"] = rrs
    n_times = len(chpi_locs["times"])
    shapes = dict(
        times=(n_times,),
        rrs=(n_times, n_hpi, 3),
        gofs=(n_times, n_hpi),
        moments=(n_times, n_hpi, 3),
    )
    for key, val in chpi_locs.items():
        chpi_locs[key] = np.array(val, float).reshape(shapes[key])
    return chpi_locs


def _chpi_locs_to_times_dig(chpi_locs):
    """Reformat chpi_locs as list of dig (dict)."""
    dig = list()
    for rrs, gofs in zip(*(chpi_locs[key] for key in ("rrs", "gofs"))):
        dig.append(
            [
                {
                    "r": rr,
                    "ident": idx,
                    "gof": gof,
                    "kind": FIFF.FIFFV_POINT_HPI,
                    "coord_frame": FIFF.FIFFV_COORD_DEVICE,
                }
                for idx, (rr, gof) in enumerate(zip(rrs, gofs), 1)
            ]
        )
    return chpi_locs["times"], dig


@verbose
def filter_chpi(
    raw,
    include_line=True,
    t_step=0.01,
    t_window="auto",
    ext_order=1,
    allow_line_only=False,
    verbose=None,
):
    """Remove cHPI and line noise from data.

    .. note:: This function will only work properly if cHPI was on
              during the recording.

    Parameters
    ----------
    raw : instance of Raw
        Raw data with cHPI information. Must be preloaded. Operates in-place.
    include_line : bool
        If True, also filter line noise.
    t_step : float
        Time step to use for estimation, default is 0.01 (10 ms).
    %(t_window_chpi_t)s
    %(ext_order_chpi)s
    allow_line_only : bool
        If True, allow filtering line noise only. The default is False,
        which only allows the function to run when cHPI information is present.

        .. versionadded:: 0.20
    %(verbose)s

    Returns
    -------
    raw : instance of Raw
        The raw data.

    Notes
    -----
    cHPI signals are in general not stationary, because head movements act
    like amplitude modulators on cHPI signals. Thus it is recommended to
    use this procedure, which uses an iterative fitting method, to
    remove cHPI signals, as opposed to notch filtering.

    .. versionadded:: 0.12
    """
    _validate_type(raw, BaseRaw, "raw")
    if not raw.preload:
        raise RuntimeError("raw data must be preloaded")
    t_step = float(t_step)
    if t_step <= 0:
        raise ValueError(f"t_step ({t_step}) must be > 0")
    n_step = int(np.ceil(t_step * raw.info["sfreq"]))
    if include_line and raw.info["line_freq"] is None:
        raise RuntimeError(
            'include_line=True but raw.info["line_freq"] is '
            "None, consider setting it to the line frequency"
        )
    hpi = _setup_hpi_amplitude_fitting(
        raw.info,
        t_window,
        remove_aliased=True,
        ext_order=ext_order,
        allow_empty=allow_line_only,
        verbose=_verbose_safe_false(),
    )

    fit_idxs = np.arange(0, len(raw.times) + hpi["n_window"] // 2, n_step)
    n_freqs = len(hpi["freqs"])
    n_remove = 2 * n_freqs
    meg_picks = pick_types(raw.info, meg=True, exclude=())  # filter all chs
    n_times = len(raw.times)

    msg = f"Removing {n_freqs} cHPI"
    if include_line:
        n_remove += 2 * len(hpi["line_freqs"])
        msg += f" and {len(hpi['line_freqs'])} line harmonic"
    msg += f" frequencies from {len(meg_picks)} MEG channels"

    recon = np.dot(hpi["model"][:, :n_remove], hpi["inv_model"][:n_remove]).T
    logger.info(msg)
    chunks = list()  # the chunks to subtract
    last_endpt = 0
    pb = ProgressBar(fit_idxs, mesg="Filtering")
    for ii, midpt in enumerate(pb):
        left_edge = midpt - hpi["n_window"] // 2
        time_sl = slice(
            max(left_edge, 0), min(left_edge + hpi["n_window"], len(raw.times))
        )
        this_len = time_sl.stop - time_sl.start
        if this_len == hpi["n_window"]:
            this_recon = recon
        else:  # first or last window
            model = hpi["model"][:this_len]
            inv_model = np.linalg.pinv(model)
            this_recon = np.dot(model[:, :n_remove], inv_model[:n_remove]).T
        this_data = raw._data[meg_picks, time_sl]
        subt_pt = min(midpt + n_step, n_times)
        if last_endpt != subt_pt:
            fit_left_edge = left_edge - time_sl.start + hpi["n_window"] // 2
            fit_sl = slice(fit_left_edge, fit_left_edge + (subt_pt - last_endpt))
            chunks.append((subt_pt, np.dot(this_data, this_recon[:, fit_sl])))
        last_endpt = subt_pt

        # Consume (trailing) chunks that are now safe to remove because
        # our windows will no longer touch them
        if ii < len(fit_idxs) - 1:
            next_left_edge = fit_idxs[ii + 1] - hpi["n_window"] // 2
        else:
            next_left_edge = np.inf
        while len(chunks) > 0 and chunks[0][0] <= next_left_edge:
            right_edge, chunk = chunks.pop(0)
            raw._data[meg_picks, right_edge - chunk.shape[1] : right_edge] -= chunk
    return raw


def _compute_good_distances(hpi_coil_dists, new_pos, dist_limit=0.005):
    """Compute good coils based on distances."""
    these_dists = cdist(new_pos, new_pos)
    these_dists = np.abs(hpi_coil_dists - these_dists)
    # there is probably a better algorithm for finding the bad ones...
    good = False
    use_mask = np.ones(len(hpi_coil_dists), bool)
    while not good:
        d = these_dists[use_mask][:, use_mask]
        d_bad = d > dist_limit
        good = not d_bad.any()
        if not good:
            if use_mask.sum() == 2:
                use_mask[:] = False
                break  # failure
            # exclude next worst point
            badness = (d * d_bad).sum(axis=0)
            exclude_coils = np.where(use_mask)[0][np.argmax(badness)]
            use_mask[exclude_coils] = False
    return use_mask, these_dists


@verbose
def get_active_chpi(raw, *, on_missing="raise", verbose=None):
    """Determine how many HPI coils were active for a time point.

    Parameters
    ----------
    raw : instance of Raw
        Raw data with cHPI information.
    %(on_missing_chpi)s
    %(verbose)s

    Returns
    -------
    n_active : array, shape (n_times)
        The number of active cHPIs for every timepoint in raw.

    Notes
    -----
    .. versionadded:: 1.2
    """
    # get meg system
    system, _ = _get_meg_system(raw.info)

    # check whether we have a neuromag system
    if system not in ["122m", "306m"]:
        raise NotImplementedError(
            "Identifying active HPI channels is not implemented for other systems than "
            "neuromag."
        )
    # extract hpi info
    chpi_info = get_chpi_info(raw.info, on_missing=on_missing)
    if (len(chpi_info[2]) == 0) or (chpi_info[1] is None):
        return np.zeros_like(raw.times)

    # extract hpi time series and infer which one was on
    chpi_ts = raw[chpi_info[1]][0].astype(int)
    chpi_active = (chpi_ts & chpi_info[2][:, np.newaxis]).astype(bool)
    return chpi_active.sum(axis=0)