arvoelke / nengolib Goto Github PK

We usually want to connect these things together under the assumption that our input is not filtered, and we must do the output filtering ourselves.

class MultiProcess(Process):

    def __init__(self, *processes):
        self._processes = processes
        super(MultiProcess, self).__init__()

    def make_step(self, *args, **kwargs):
        steps = [p.make_step(*args, **kwargs) for p in self._processes]
        def step(t):
            return np.sum([s(t) for s in steps], axis=0)
        return step

Too lazy to put as a PR atm. This is a reminder to self.

• eig
• svd

Note: matrix_rank is not in scipy.

Adding the synapse (1. - ~z) / dt to a connection is an elegant / cheap way to implement differentiation.

Reference licenses from used libraries:

• nengo
• numpy
• scipy (especially for discrete2cont)
• sobol_seq.py

Probably also:

• ipython
• sphinx
• pytest
• seaborn

Due to how the solver is set up, building the model multiple times (via nengo.Simulator, without redefining the network) will raise an error.

• discrete2cont (2c7b6ec)
• tfestimate (scipy/scipy#5737)
• minreal (#19)
• balreal
  • continuous (#19)
  • discrete
• modred
  • continuous (#19)
  • discrete
• PEM (ARX):
  • tfest
  • ssest

Also include example of how to scale the radii based on the H2-norms for arbitrarily filtered white-noise input. Related to #6.

• Documentation
  • #1 Generate Documentation
  • #3 Pretty README
  • Example/demonstration notebooks for all key features
    • #6 Discretized dynamics
    • Network using number-theoretic method
    • Heterogeneous synapses
    • LinearFilter efficient simulation
    • LinearSystem manipulation
  • Docstrings should reference related notebooks and vice-versa
• #4 License
• Key Functionality
  • #8 Matlab routines for model reduction and system identification (Partial)
  • #9 Leech lattice
  • #13 Manipulating LTI systems
  • #22 LinearNetwork class
• Testing
  • 100% Test Coverage
  • #46 Target release of Nengo 2.1.0
  • #17 Python 3 support

Seems to be some numerical issues with PadeDelay of higher order. Requires investigation. Using HankelNorm as default normalizer for now.

Pure delay, and a box filter.

This may be useful for analysis. A work-in-progress:

import numpy as np
from scipy.misc import factorial, pade

# Taylor series
c = 1. / factorial(np.arange(10))

# Array representing the Toeplitz matrix
a = c[-2::-1]
assert len(a) % 2 == 1  # odd

# Form the Toeplitz matrix explicitly
N = len(a) // 2 + 1  # even
T = np.empty((N, N))
for i in range(N):
    T[i, :] = a[N-i-1:2*N-i-1]

# Solution vector
y = -c[N:]
assert len(y) == N

# Levinson algorithm using the underlying a vector
f = np.zeros(N)
b = np.zeros(N)

f[0] = b[0] = 1. / a[N-1]
x = y[0] * b

for n in range(1, N):
    ef = a[N-n-1:N-1].dot(f[:n])
    eb = a[N:N+n].dot(b[:n])
    ex = a[N-n-1:N-1].dot(x[:n])

    fn = np.append(f[:n], [0])
    bn = np.append([0], b[:n])
    
    f[:n+1] = 1. / (1 - eb*ef) * fn - ef / (1 - eb*ef) * bn
    b[:n+1] = -eb / (1 - eb*ef) * fn + 1. / (1 - eb*ef) * bn
    x[:n+1] += (y[n] - ex) * b[:n+1]

assert np.allclose(np.dot(T, x), y)

# Relationship to Pade approximants
assert np.allclose(pade(c, N)[1], np.append(x[::-1], [1]))

References:

• A new algorithm for computing Pade approximants (2011)
• https://en.wikipedia.org/wiki/Levinson_recursion

Also note that computing d = D^{m+1}[q(x)*f(x) - p(x)] eliminates p and reveals the higher-order error O(x^n) in terms of q and f:

from scipy.special import binom

m = 2
n = 3

# Exponential taylor series
a = 1. / factorial(np.arange(m + n + 1))

assert m + n + 1 == len(a)
p, q = pade(a, n)
assert len(p) == m
assert len(q) == n

a = np.append(a, [0]*(m+1))  # padding for convenience
def x_pow(n): return np.poly1d([1., 0]) ** n

d = np.poly1d([0.])
for l in range(m+n+1):
    for i in range(l, l+m+2):
        d += a[i] * binom(m+1, i-l) * factorial(i) / factorial(l) * np.polymul(
            x_pow(l), np.polyder(q, m+1-(i-l)))
print d

and the derivatives can be grouped together into a single linear differential operator for each index l.

Derived using:

• https://en.wikipedia.org/wiki/General_Leibniz_rule
• https://www-sop.inria.fr/apics/anap03/PadeTalk.pdf (generalizing the first part of page 14)

By adding an analog attribute that cont2discrete and discrete2cont will flip. This should enable some better error checking in the few places that assume a synapse is in one form or another (apply_filter, ss2sim, NengoLinearFilterMixin, balreal). Then make a discrete q operator. Also replace all instances of scipy.signal.cont2discrete with nengolib.signal.cont2discrete.

Currently some methods that operate in the state-space return 4-tuples. Previously this seemed like the clear thing to do, but now it feels inconsistent and slightly difficult to remember since the line is sometimes blurry.

• L1-Norm (#43)
• H2-Norm (#103)
• Balanced Realization (#85)
• Identity realization (#103)
• Controllable
• Observable
• Mechanism for pretending the input is something + arbitrary filter
• Document in LinearNetwork notebook

Would be good to expose this to the user in case they have derivatives they can supply, e.g.,

from numpy.linalg import matrix_power

def deriv_mapping(sys, H):  # both analog or both discretized
    gain = np.sum(H.num)
    c = H.den / gain
    k = len(H)
    A, B, C, D = sys.ss
    powA = [matrix_power(A, i) for i in range(k + 1)]
    AH = np.sum([c[i] * powA[i] for i in range(k + 1)], axis=0)
    BH = [np.dot(np.sum([c[i] * powA[i-j-1] for i in range(j+1, k+1)], axis=0), B)
          for j in range(k)]
    return AH, BH

Now that the release has been made: https://github.com/nengo/nengo/releases

• Create a new ensemble array node and pass the solver through there (until nengo/nengo#1040).
• See if there's a way to bypass problem from nengo/nengo#1044 in Nengo 2.1.0.
• Update reservoir notebook to use better regularization for the structured example.

Show connection between delay <-> deconvolution <-> acausal filtering. Be explicit about sources of error with respect to input's Laplace transform.

Not sure if I agree with the default way that nengo's filt does things. Messes up the impulse response and such. Related to #70.

Add a network to subclass LinearNetwork with a delay system and reasonable defaults. Provide an easy mechanism for decoding functions of the input history via,

def delay_readout(q, thetap, theta):
    # Normalized C matrix
    c = np.zeros(q)
    for i in range(q):
        j = np.arange(i+1, dtype=np.float64)
        c[i] += 1 / binom(q, i) * np.sum(
            binom(q, j) * binom(2*q - 1 - j, i - j) *
            (-thetap/theta)**(i - j))
    return c

Depends on #42, since changes to the state-space need to be reflected in the above C matrix (via the corresponding similarity transform).

Also necessary for the representation to be distributed (i.e., not an ensemble array) in order to get the most general class of nonlinear functions.

To do principle 3 in the continuous case.

Use pytest. Look into difficulty of integrating with Travis-CI and/or Codecov.

Lazily-compute and store them once in the nengo cache directory. Then lazy-load when needed.

We shouldn't mutate the user's parameters. This leads to difficult bugs when the radii is iteratively updated according to some scheme.

i.e., F(H(s)^{-1}) as a LinearSystem

• Support non-linear functions of the dynamical state.
• Ensure (with example/test) that MIMO state-space works as expected.

And the same for tau close to 0.

See gain-control branch.

In the time-domain: (tau1 * nengolib.Lowpass(tau1) - tau2 * nengolib.Lowpass(tau2)) / (tau1 - tau2)
In the frequency-domain: 1 / (tau1 * s + 1) / (tau2 * s + 1)

Include derivation that these two definitions are equivalent (with the latter being well-defined when tau1 == tau2).

Note: unit tests have already been made.

When all of the synapses are single order, the block diagonal A matrix is actually most diagonal, and so we can do an elementwise multiply rather than a matrix multiply.

ValidationError: LinearSystem.default_dt: Unconfigurable parameters have no defaults. Please ensure the value of the parameter is set before trying to access it.

This is because I intentionally don't invoke the superclass of LinearSystem, which is where default_dt gets set (even if it never actually gets used by the underlying filter).

Add badges for build/coverage, as well as some nice images.

Rather than using decompose_states and then simulating it multiple times, which is inefficient.

This is needed by the state normalizers for instance.

[Deleted]

Create a class that accepts a system in any number of forms:

• state-space
• transfer function
• zero-pole
• nengo.synapses.LinearFilter

Then functions such as impulse, is_exp_stable, scale_state, and state_norm must accept these objects. The objects may also support multiplication, operation, and comparisons.

Use Sphinx to model after scipy.

Sometimes ss2sim (via LinearSystem) will destabilize sys when discretizing it with the given dt. This occurs for larger order delays that are not properly balanced / normalized. A workaround is to supply dt=None. It would be better to gracefully detect and handle this problem automatically, since it's not always apparent that something is going to go wrong.

The normalizer is the only reason this isn't done atm. But still don't have a use case for this anyways. Might be good for when the recurrent function is discrete sample data, but that will be more important for nonlinear systems. Putting this here as a reminder to self.