def test_lif(Simulator):
    """Test that the dynamic model approximately matches the rates"""
    rng = np.random.RandomState(85243)

    dt = 0.001
    n = 5000
    x = 0.5
    max_rates = rng.uniform(low=10, high=200, size=n)
    intercepts = rng.uniform(low=-1, high=1, size=n)

    m = nengo.Network()
    with m:
        ins = nengo.Node(x)
        ens = nengo.Ensemble(
            nengo.LIF(n), 1, max_rates=max_rates, intercepts=intercepts)
        nengo.Connection(ins, ens.neurons, transform=np.ones((n, 1)))
        spike_probe = nengo.Probe(ens.neurons, "output")

    sim = Simulator(m, dt=dt)

    t_final = 1.0
    sim.run(t_final)
    spikes = sim.data[spike_probe].sum(0)

    math_rates = ens.neurons.rates(
        x, *ens.neurons.gain_bias(max_rates, intercepts))
    sim_rates = spikes / t_final
    logger.debug("ME = %f", (sim_rates - math_rates).mean())
    logger.debug("RMSE = %f",
                 rms(sim_rates - math_rates) / (rms(math_rates) + 1e-20))
    assert np.sum(math_rates > 0) > 0.5 * n, (
        "At least 50% of neurons must fire")
    assert np.allclose(sim_rates, math_rates, atol=1, rtol=0.02)

def test_decoder_solver(Solver, plt, rng):
    if isinstance(Solver, tuple):
        Solver, args, kwargs = Solver
    else:
        args, kwargs = (), {}

    dims = 1
    n_neurons = 100
    n_points = 500

    rates = get_rate_function(n_neurons, dims, rng=rng)
    E = get_encoders(n_neurons, dims, rng=rng)

    train = get_eval_points(n_points, dims, rng=rng)
    Atrain = rates(np.dot(train, E))

    D, _ = Solver(*args, **kwargs)(Atrain, train, rng=rng)

    test = get_eval_points(n_points, dims, rng=rng, sort=True)
    Atest = rates(np.dot(test, E))
    est = np.dot(Atest, D)
    rel_rmse = rms(est - test) / rms(test)

    plt.plot(test, np.zeros_like(test), 'k--')
    plt.plot(test, test - est)
    plt.title("relative RMSE: %0.2e" % rel_rmse)

    atol = 3.5e-2 if Solver is LstsqNoise else 1.5e-2
    assert np.allclose(test, est, atol=atol, rtol=1e-3)
    assert rel_rmse < 0.02

def test_decoder_solver(Solver):
    rng = np.random.RandomState(39408)

    dims = 1
    n_neurons = 100
    n_points = 500

    rates = get_rate_function(n_neurons, dims, rng=rng)
    E = get_encoders(n_neurons, dims, rng=rng)

    train = get_eval_points(n_points, dims, rng=rng)
    Atrain = rates(np.dot(train, E))

    D, _ = Solver()(Atrain, train, rng=rng)

    test = get_eval_points(n_points, dims, rng=rng, sort=True)
    Atest = rates(np.dot(test, E))
    est = np.dot(Atest, D)
    rel_rmse = rms(est - test) / rms(test)

    with Plotter(nengo.Simulator) as plt:
        plt.plot(test, np.zeros_like(test), 'k--')
        plt.plot(test, test - est)
        plt.title("relative RMSE: %0.2e" % rel_rmse)
        plt.savefig('test_decoders.test_decoder_solver.%s.pdf'
                    % Solver.__name__)
        plt.close()

    assert np.allclose(test, est, atol=3e-2, rtol=1e-3)
    assert rel_rmse < 0.02

def test_fastlif(plt):
    """Test that the dynamic model approximately matches the rates."""
    # Based nengo.tests.test_neurons.test_lif
    rng = np.random.RandomState(10)

    dt = 1e-3
    n = 5000
    x = 0.5
    encoders = np.ones((n, 1))
    max_rates = rng.uniform(low=10, high=200, size=n)
    intercepts = rng.uniform(low=-1, high=1, size=n)

    m = nengo.Network()
    with m:
        ins = nengo.Node(x)
        ens = nengo.Ensemble(n, dimensions=1,
                             neuron_type=FastLIF(),
                             encoders=encoders,
                             max_rates=max_rates,
                             intercepts=intercepts)
        nengo.Connection(
            ins, ens.neurons, transform=np.ones((n, 1)), synapse=None)
        spike_probe = nengo.Probe(ens.neurons)
        voltage_probe = nengo.Probe(ens.neurons, 'voltage')
        ref_probe = nengo.Probe(ens.neurons, 'refractory_time')

    t_final = 1.0
    with nengo.Simulator(m, dt=dt) as sim:
        sim.run(t_final)

    i = 3
    plt.subplot(311)
    plt.plot(sim.trange(), sim.data[spike_probe][:, :i])
    plt.subplot(312)
    plt.plot(sim.trange(), sim.data[voltage_probe][:, :i])
    plt.subplot(313)
    plt.plot(sim.trange(), sim.data[ref_probe][:, :i])
    plt.ylim([-dt, ens.neuron_type.tau_ref + dt])

    # check rates against analytic rates
    math_rates = ens.neuron_type.rates(
        x, *ens.neuron_type.gain_bias(max_rates, intercepts))
    spikes = sim.data[spike_probe]
    sim_rates = (spikes > 0).sum(0) / t_final
    print("ME = %f" % (sim_rates - math_rates).mean())
    print("RMSE = %f" % (
        rms(sim_rates - math_rates) / (rms(math_rates) + 1e-20)))
    assert np.sum(math_rates > 0) > 0.5 * n, (
        "At least 50% of neurons must fire")
    assert np.allclose(sim_rates, math_rates, atol=1, rtol=0.02)

    # if voltage and ref time are non-constant, the probe is doing something
    assert np.abs(np.diff(sim.data[voltage_probe])).sum() > 1
    assert np.abs(np.diff(sim.data[ref_probe])).sum() > 1

    # compute spike counts after each timestep
    actual_counts = (spikes > 0).cumsum(axis=0)
    expected_counts = np.outer(sim.trange(), math_rates)
    assert (abs(actual_counts - expected_counts) < 1).all()

def test_subsolvers(Solver, seed, rng, tol=1e-2):
    get_rng = lambda: np.random.RandomState(seed)

    A, b = get_system(500, 100, 5, rng=rng)
    x0, _ = Solver(solver=cholesky)(A, b, rng=get_rng())

    subsolvers = [conjgrad, block_conjgrad]
    for subsolver in subsolvers:
        x, info = Solver(solver=subsolver, tol=tol)(A, b, rng=get_rng())
        rel_rmse = rms(x - x0) / rms(x0)
        assert rel_rmse < 4 * tol

def _test_rates(Simulator, rates, name=None):
    if name is None:
        name = rates.__name__

    n = 100
    max_rates = 50 * np.ones(n)
    # max_rates = 200 * np.ones(n)
    intercepts = np.linspace(-0.99, 0.99, n)
    encoders = np.ones((n, 1))

    model = nengo.Network()
    with model:
        model.config[nengo.Ensemble].max_rates = max_rates
        model.config[nengo.Ensemble].intercepts = intercepts
        model.config[nengo.Ensemble].encoders = encoders
        u = nengo.Node(output=whitenoise(1, 5, seed=8393))
        a = nengo.Ensemble(n, 1, neuron_type=nengo.LIFRate())
        b = nengo.Ensemble(n, 1, neuron_type=nengo.LIF())
        nengo.Connection(u, a, synapse=0)
        nengo.Connection(u, b, synapse=0)
        up = nengo.Probe(u)
        ap = nengo.Probe(a.neurons)
        bp = nengo.Probe(b.neurons)

    dt = 1e-3
    sim = Simulator(model, dt=dt)
    sim.run(2.)

    t = sim.trange()
    x = sim.data[up]
    a_rates = sim.data[ap] / dt
    spikes = sim.data[bp]
    b_rates = rates(t, spikes)

    with Plotter(Simulator) as plt:
        ax = plt.subplot(411)
        plt.plot(t, x)
        ax = plt.subplot(412)
        implot(plt, t, intercepts, a_rates.T, ax=ax)
        ax.set_ylabel('intercept')
        ax = plt.subplot(413)
        implot(plt, t, intercepts, b_rates.T, ax=ax)
        ax.set_ylabel('intercept')
        ax = plt.subplot(414)
        implot(plt, t, intercepts, (b_rates - a_rates).T, ax=ax)
        ax.set_xlabel('time [s]')
        ax.set_ylabel('intercept')
        plt.savefig('utils.test_neurons.test_rates.%s.pdf' % name)
        plt.close()

    tmask = (t > 0.1) & (t < 1.9)
    relative_rmse = rms(b_rates[tmask] - a_rates[tmask]) / rms(a_rates[tmask])
    return relative_rmse

def test_subsolvers(solver, tol=1e-2):
    rng = np.random.RandomState(89)
    get_rng = lambda: np.random.RandomState(87)

    A, b = get_system(500, 100, 5, rng=rng)
    x0, _ = solver(A, b, rng=get_rng(), solver=_cholesky)

    subsolvers = [_conjgrad, _block_conjgrad]
    for subsolver in subsolvers:
        x, info = solver(A, b, rng=get_rng(), solver=subsolver, tol=tol)
        rel_rmse = rms(x - x0) / rms(x0)
        assert rel_rmse < 3 * tol

def test_lif(Simulator):
    """Test that the dynamic model approximately matches the rates"""
    rng = np.random.RandomState(85243)

    dt = 0.001
    n = 5000
    x = 0.5
    encoders = np.ones((n, 1))
    max_rates = rng.uniform(low=10, high=200, size=n)
    intercepts = rng.uniform(low=-1, high=1, size=n)

    m = nengo.Network()
    with m:
        ins = nengo.Node(x)
        ens = nengo.Ensemble(
            n, dimensions=1, neuron_type=nengo.LIF(),
            encoders=encoders, max_rates=max_rates, intercepts=intercepts)
        nengo.Connection(ins, ens.neurons, transform=np.ones((n, 1)))
        spike_probe = nengo.Probe(ens.neurons)
        voltage_probe = nengo.Probe(ens.neurons, 'voltage')
        ref_probe = nengo.Probe(ens.neurons, 'refractory_time')

    sim = Simulator(m, dt=dt)

    t_final = 1.0
    sim.run(t_final)

    with Plotter(Simulator) as plt:
        i = 3
        plt.subplot(311)
        plt.plot(sim.trange(), sim.data[spike_probe][:, :i])
        plt.subplot(312)
        plt.plot(sim.trange(), sim.data[voltage_probe][:, :i])
        plt.subplot(313)
        plt.plot(sim.trange(), sim.data[ref_probe][:, :i])
        plt.ylim([-dt, ens.neuron_type.tau_ref + dt])
        plt.savefig('test_neurons.test_lif.pdf')
        plt.close()

    # check rates against analytic rates
    math_rates = ens.neuron_type.rates(
        x, *ens.neuron_type.gain_bias(max_rates, intercepts))
    sim_rates = sim.data[spike_probe].sum(0) / t_final
    logger.debug("ME = %f", (sim_rates - math_rates).mean())
    logger.debug("RMSE = %f",
                 rms(sim_rates - math_rates) / (rms(math_rates) + 1e-20))
    assert np.sum(math_rates > 0) > 0.5 * n, (
        "At least 50% of neurons must fire")
    assert np.allclose(sim_rates, math_rates, atol=1, rtol=0.02)

    # if voltage and ref time are non-constant, the probe is doing something
    assert np.abs(np.diff(sim.data[voltage_probe])).sum() > 1
    assert np.abs(np.diff(sim.data[ref_probe])).sum() > 1

def test_alif(Simulator):
    """Test ALIF and ALIFRate by comparing them to each other"""

    n = 100
    max_rates = 50 * np.ones(n)
    intercepts = np.linspace(-0.99, 0.99, n)
    encoders = np.ones((n, 1))
    nparams = dict(tau_n=1, inc_n=10e-3)
    eparams = dict(n_neurons=n, max_rates=max_rates,
                   intercepts=intercepts, encoders=encoders)

    model = nengo.Network()
    with model:
        u = nengo.Node(output=0.5)
        a = nengo.Ensemble(neuron_type=nengo.AdaptiveLIFRate(**nparams),
                           dimensions=1,
                           **eparams)
        b = nengo.Ensemble(neuron_type=nengo.AdaptiveLIF(**nparams),
                           dimensions=1,
                           **eparams)
        nengo.Connection(u, a, synapse=0)
        nengo.Connection(u, b, synapse=0)
        ap = nengo.Probe(a, "spikes", synapse=0)
        bp = nengo.Probe(b, "spikes", synapse=0)

    dt = 1e-3
    sim = Simulator(model, dt=dt)
    sim.run(2.)

    t = sim.trange()
    a_rates = sim.data[ap] / dt
    spikes = sim.data[bp]
    b_rates = rates_kernel(t, spikes)

    tmask = (t > 0.1) & (t < 1.7)
    rel_rmse = rms(b_rates[tmask] - a_rates[tmask]) / rms(a_rates[tmask])

    with Plotter(Simulator) as plt:
        ax = plt.subplot(311)
        implot(plt, t, intercepts[::-1], a_rates.T, ax=ax)
        ax.set_ylabel('input')
        ax = plt.subplot(312)
        implot(plt, t, intercepts[::-1], b_rates.T, ax=ax)
        ax.set_ylabel('input')
        ax = plt.subplot(313)
        implot(plt, t, intercepts[::-1], (b_rates - a_rates)[tmask].T, ax=ax)
        ax.set_xlabel('time [s]')
        ax.set_ylabel('input')
        plt.savefig('test_neurons.test_alif.pdf')
        plt.close()

    assert rel_rmse < 0.07

def test_alif(Simulator, plt):
    """Test ALIF and ALIFRate by comparing them to each other"""

    n = 100
    max_rates = 50 * np.ones(n)
    intercepts = np.linspace(-0.99, 0.99, n)
    encoders = np.ones((n, 1))
    nparams = dict(tau_n=1, inc_n=10e-3)
    eparams = dict(n_neurons=n, max_rates=max_rates,
                   intercepts=intercepts, encoders=encoders)

    model = nengo.Network()
    with model:
        u = nengo.Node(output=0.5)
        a = nengo.Ensemble(neuron_type=AdaptiveLIFRate(**nparams),
                           dimensions=1,
                           **eparams)
        b = nengo.Ensemble(neuron_type=AdaptiveLIF(**nparams),
                           dimensions=1,
                           **eparams)
        nengo.Connection(u, a, synapse=0)
        nengo.Connection(u, b, synapse=0)
        ap = nengo.Probe(a.neurons)
        bp = nengo.Probe(b.neurons)

    with Simulator(model) as sim:
        sim.run(2.)

    t = sim.trange()
    a_rates = sim.data[ap]
    spikes = sim.data[bp]
    b_rates = nengo.Lowpass(0.04).filtfilt(spikes)

    tmask = (t > 0.1) & (t < 1.7)
    rel_rmse = rms(b_rates[tmask] - a_rates[tmask]) / rms(a_rates[tmask])

    ax = plt.subplot(311)
    implot(plt, t, intercepts[::-1], a_rates.T, ax=ax)
    ax.set_ylabel('input')
    ax = plt.subplot(312)
    implot(plt, t, intercepts[::-1], b_rates.T, ax=ax)
    ax.set_ylabel('input')
    ax = plt.subplot(313)
    implot(plt, t, intercepts[::-1], (b_rates - a_rates)[tmask].T, ax=ax)
    ax.set_xlabel('time [s]')
    ax.set_ylabel('input')

    assert rel_rmse < 0.07

def randomized_svd(A, Y, sigma, rng=np.random,
                   n_components=60, n_oversamples=10, **kwargs):
    """Solve the least-squares system using a randomized (partial) SVD.

    Parameters
    ----------
    n_components : int (default is 50)
        The number of SVD components to compute. A small survey of activity
        matrices suggests that the first 50 components capture almost all
        the variance.
    n_oversamples: int (default is 10)
        The number of additional samples on the range of A.
    n_iter : int (default is 0)
        The number of power iterations to perform (can help with noisy data).

    See also
    --------
    ``sklearn.utils.extmath.randomized_svd`` for details about the parameters.
    """
    from sklearn.utils.extmath import randomized_svd as sklearn_randomized_svd

    Y, m, n, _, matrix_in = _format_system(A, Y)
    if min(m, n) <= n_components + n_oversamples:
        # more efficient to do a full SVD
        return svd(A, Y, sigma, rng=rng)

    U, s, V = sklearn_randomized_svd(
        A, n_components, random_state=rng, **kwargs)
    si = s / (s**2 + m * sigma**2)
    X = np.dot(V.T, si[:, None] * np.dot(U.T, Y))
    info = {'rmses': npext.rms(Y - np.dot(A, X), axis=0)}
    return X if matrix_in else X.flatten(), info

def test_regularization(Simulator, nl_nodirect, plt):

    # TODO: multiple trials per parameter set, with different seeds

    Solvers = [LstsqL2, LstsqL2nz]
    neurons = np.array([10, 20, 50, 100])
    regs = np.linspace(0.01, 0.3, 16)
    filters = np.linspace(0, 0.03, 11)

    buf = 0.2  # buffer for initial transients
    dt = 1e-3
    tfinal = 3 + buf

    def input_function(t):
        return np.interp(t, [1, 3], [-1, 1], left=-1, right=1)

    model = nengo.Network('test_regularization')
    with model:
        model.config[nengo.Ensemble].neuron_type = nl_nodirect()
        u = nengo.Node(output=input_function)
        up = nengo.Probe(u)

        probes = np.zeros(
            (len(Solvers), len(neurons), len(regs), len(filters)),
            dtype='object')

        for j, n_neurons in enumerate(neurons):
            a = nengo.Ensemble(n_neurons, dimensions=1)
            nengo.Connection(u, a)

            for i, Solver in enumerate(Solvers):
                for k, reg in enumerate(regs):
                    for l, synapse in enumerate(filters):
                        probes[i, j, k, l] = nengo.Probe(
                            a, solver=Solver(reg=reg), synapse=synapse)

    sim = Simulator(model, dt=dt)
    sim.run(tfinal)
    t = sim.trange()

    ref = sim.data[up]
    rmse_buf = lambda a, b: rms(a[t > buf] - b[t > buf])
    rmses = np.zeros(probes.shape)
    for i, probe in enumerate(probes.flat):
        rmses.flat[i] = rmse_buf(sim.data[probe], ref)
    rmses = rmses - rmses[:, :, [0], :]

    plt.figure(figsize=(8, 12))
    X, Y = np.meshgrid(filters, regs)

    for i, Solver in enumerate(Solvers):
        for j, n_neurons in enumerate(neurons):
            plt.subplot(len(neurons), len(Solvers), len(Solvers)*j + i + 1)
            Z = rmses[i, j, :, :]
            plt.contourf(X, Y, Z, levels=np.linspace(Z.min(), Z.max(), 21))
            plt.xlabel('filter')
            plt.ylabel('reg')
            plt.title("%s (N=%d)" % (Solver.__name__, n_neurons))

    plt.tight_layout()

def lstsq_drop(A, Y, rng, E=None, noise_amp=0.1, drop=0.25, solver=lstsq_L2nz):
    """Find sparser decoders/weights by dropping small values.

    This solver first solves for coefficients (decoders/weights) with
    L2 regularization, drops those nearest to zero, and retrains remaining.
    """
    Y, m, n, d, matrix_in = _format_system(A, Y)

    # solve for coefficients using standard solver
    X, info0 = solver(A, Y, rng=rng, noise_amp=noise_amp)
    X = np.dot(X, E) if E is not None else X

    # drop weights close to zero, based on `drop` ratio
    Xabs = np.sort(np.abs(X.flat))
    threshold = Xabs[int(np.round(drop * Xabs.size))]
    X[np.abs(X) < threshold] = 0

    # retrain nonzero weights
    Y = np.dot(Y, E) if E is not None else Y
    for i in range(X.shape[1]):
        nonzero = X[:, i] != 0
        if nonzero.sum() > 0:
            X[nonzero, i], info1 = solver(
                A[:, nonzero], Y[:, i], rng=rng, noise_amp=0.1 * noise_amp)

    info = {'rmses': npext.rms(Y - np.dot(A, X), axis=0),
            'info0': info0, 'info1': info1}
    return X if matrix_in else X.flatten(), info

 def __call__(self, A, Y, sigma, rng=None):
     Y, m, _, _, matrix_in = format_system(A, Y)
     U, s, V = np.linalg.svd(A, full_matrices=0)
     si = s / (s**2 + m * sigma**2)
     X = np.dot(V.T, si[:, None] * np.dot(U.T, Y))
     info = {'rmses': npext.rms(Y - np.dot(A, X), axis=0)}
     return X if matrix_in else X.flatten(), info

    def make_step(self, size_in, size_out, dt, rng):
        assert size_in[0] == 0
        assert size_out[0] == 1

        rate = 1. / dt

        orig_rate, orig = readwav(self.path)
        new_size = int(orig.size * (rate / orig_rate))
        wave = resample(orig, new_size)
        wave -= wave.mean()

        # Normalize wave to desired rms
        wave_rms = npext.rms(wave)
        wave *= (self.rms / wave_rms)

        if self.at_end == 'loop':

            def step_wavfileloop(t):
                idx = int(t * rate) % wave.size
                return wave[idx]
            return step_wavfileloop

        elif self.at_end == 'stop':

            def step_wavfilestop(t):
                idx = int(t * rate)
                if idx > wave.size:
                    return 0.
                else:
                    return wave[idx]
            return step_wavfilestop

 def __call__(self, A, Y, rng=None, E=None):
     Y = self.mul_encoders(Y, E)
     X, residuals2, rank, s = np.linalg.lstsq(A, Y, rcond=self.rcond)
     return X, {'rmses': npext.rms(Y - np.dot(A, X), axis=0),
                'residuals': np.sqrt(residuals2),
                'rank': rank,
                'singular_values': s}

def conjgrad_scipy(A, Y, sigma, tol=1e-4):
    """Solve the least-squares system using Scipy's conjugate gradient."""
    import scipy.sparse.linalg
    Y, m, n, d, matrix_in = _format_system(A, Y)

    damp = m * sigma**2
    calcAA = lambda x: np.dot(A.T, np.dot(A, x)) + damp * x
    G = scipy.sparse.linalg.LinearOperator(
        (n, n), matvec=calcAA, matmat=calcAA, dtype=A.dtype)
    B = np.dot(A.T, Y)

    X = np.zeros((n, d), dtype=B.dtype)
    infos = np.zeros(d, dtype='int')
    itns = np.zeros(d, dtype='int')
    for i in range(d):
        def callback(x):
            itns[i] += 1  # use the callback to count the number of iterations

        X[:, i], infos[i] = scipy.sparse.linalg.cg(
            G, B[:, i], tol=tol, callback=callback)

    info = {'rmses': npext.rms(Y - np.dot(A, X), axis=0),
            'iterations': itns,
            'info': infos}
    return X if matrix_in else X.flatten(), info

def lstsq_L1(A, Y, rng=np.random, E=None, l1=1e-4, l2=1e-6):
    """Least-squares with L1 and L2 regularization (elastic net).

    This method is well suited for creating sparse decoders or weight matrices.
    """
    import sklearn.linear_model

    # TODO: play around with these regularization constants (I just guessed).
    #   Do we need to scale regularization by number of neurons, to get same
    #   level of sparsity? esp. with weights? Currently, setting l1=1e-3 works
    #   well with weights when connecting 1D populations with 100 neurons each.
    a = l1 * A.max()      # L1 regularization
    b = l2 * A.max()**2   # L2 regularization
    alpha = a + b
    l1_ratio = a / (a + b)

    # --- solve least-squares A * X = Y
    if E is not None:
        Y = np.dot(Y, E)

    model = sklearn.linear_model.ElasticNet(
        alpha=alpha, l1_ratio=l1_ratio, fit_intercept=False, max_iter=1000)
    model.fit(A, Y)
    X = model.coef_.T
    X.shape = (A.shape[1], Y.shape[1]) if Y.ndim > 1 else (A.shape[1],)
    infos = {'rmses': npext.rms(Y - np.dot(A, X), axis=0)}
    return X, infos

def cholesky(A, y, sigma, transpose=None):
    """Solve the least-squares system using the Cholesky decomposition."""
    m, n = A.shape
    if transpose is None:
        # transpose if matrix is fat, but not if we have sigmas for each neuron
        transpose = m < n and sigma.size == 1

    if transpose:
        # substitution: x = A'*xbar, G*xbar = b where G = A*A' + lambda*I
        G = np.dot(A, A.T)
        b = y
    else:
        # multiplication by A': G*x = A'*b where G = A'*A + lambda*I
        G = np.dot(A.T, A)
        b = np.dot(A.T, y)

    # add L2 regularization term 'lambda' = m * sigma**2
    np.fill_diagonal(G, G.diagonal() + m * sigma**2)

    try:
        import scipy.linalg
        factor = scipy.linalg.cho_factor(G, overwrite_a=True)
        x = scipy.linalg.cho_solve(factor, b)
    except ImportError:
        L = np.linalg.cholesky(G)
        L = np.linalg.inv(L.T)
        x = np.dot(L, np.dot(L.T, b))

    x = np.dot(A.T, x) if transpose else x
    info = {'rmses': npext.rms(y - np.dot(A, x), axis=0)}
    return x, info

def test_eval_point_decoding(Simulator, nl_nodirect, plt, seed):
    with nengo.Network(seed=seed) as model:
        model.config[nengo.Ensemble].neuron_type = nl_nodirect()
        a = nengo.Ensemble(200, 2)
        b = nengo.Ensemble(100, 1)
        c = nengo.Connection(a, b, function=lambda x: x[0] * x[1])

    sim = Simulator(model)
    eval_points, targets, decoded = eval_point_decoding(c, sim)

    def contour(xy, z):
        xi = np.linspace(-1, 1, 101)
        yi = np.linspace(-1, 1, 101)
        zi = griddata(xy[:, 0], xy[:, 1], z.ravel(), xi, yi, interp="linear")
        plt.contourf(xi, yi, zi, cmap=plt.cm.seismic)
        plt.colorbar()

    plt.figure(figsize=(15, 5))
    plt.subplot(131)
    contour(eval_points, targets)
    plt.title("Target (desired decoding)")
    plt.subplot(132)
    plt.title("Actual decoding")
    contour(eval_points, decoded)
    plt.subplot(133)
    plt.title("Difference between actual and desired")
    contour(eval_points, decoded - targets)

    # Generous error check, just to make sure it's in the right ballpark.
    # Also make sure error is above zero, i.e. y != z
    error = rms(decoded - targets, axis=1).mean()
    assert error < 0.1 and error > 1e-8