This will finally resolve #147.

Merge and rewrite

• [x] LSIIR, invLSIIR and invLSIIR_unc based on Monte Carlo,
• [x] LSFIR, invLSFIR, invLSFIR_unc, invLSFIR_uncMC.
• [x] After reviewing all this, we should merge #251 into this one here and then merge all that into prepare_major_release_2.0.0

Checkout the docs as well to be sure they look nicely: https://pydynamic.readthedocs.io/en/fix-147_rewrite_inverse_methods_of_model_estimation

To avoid unambiguous method naming we will combine all methods from identification and devonvolution in one new module model_estimation including the fit_filter.py methods with same namings in both modules. That requires us to rename some methods out of devonvolution.fit_filter.py. Because of the following incompatibility with previous versions of PyDynamic we need to insert a deprecation warning into the next minor release and inform users about the upcoming change.

In cases of non-stationary uncertainty (king="diag"), the output of FIRuncFilter returns an uncertainty-result, that seems to be shifted.

Looking into the code, this is very likely caused by the equations to calculate UncCov: L.182 + L.186 + L.190 + L.192. Here, theta and Ulow are multiplied. However, for values in theta it holds theta[i+1] refers to an earlier point in time than theta[i] - while in Ulow it is reversed. Therefore, during the multiplication values are combined, that do not correspond in their time-order. To fix this, theta probably needs to be reversed, to achieve correct (time-ascending) order.

It seems like this misbehaviour is not covered by tests or examples so far. Therefore a suitable test/example should be introduced (Maybe against some MC-method). A comparison with the upcoming IIRuncFilter could also be an option.

Environment

• OS: Windows 10
• PyDynamic Version 1.6.0 (direct from master)

MWE showing comparison with MC-result

import matplotlib.pyplot as plt
import numpy as np

from PyDynamic.uncertainty.propagate_MonteCarlo import MC
from PyDynamic.uncertainty.propagate_filter import FIRuncFilter

n_signal = 100
n_filter = 50

x = np.append(np.arange(n_signal//2), np.zeros(n_signal//2))
ux = 0.1 * np.abs(x)

theta = np.array([0.5,0.5] + [0] * (n_filter - 2))
utheta = np.zeros((n_filter, n_filter))


# run filter
y, uy = FIRuncFilter(
    x, ux, theta, utheta, kind="diag"
)  # apply uncertain FIR filter (GUM formula)
yMC, uyMC = MC(
    x, ux, theta, [1.0], utheta, runs=1000
)  # apply uncertain FIR filter (Monte Carlo)


# compare FIR and MC results
fig, ax = plt.subplots(nrows=2)

ax[0].plot(x, label="x")
ax[0].plot(theta, label="theta")
ax[0].plot(y, label="y")
ax[0].plot(yMC, label="y (MC)")
ax[0].legend()

ax[1].plot(ux, label="ux")
ax[1].plot(np.diag(utheta), label="utheta")
ax[1].plot(uy, label="uy")
ax[1].plot(np.sqrt(np.diag(uyMC)), label="uy (MC)")
ax[1].legend()

plt.show()

Output of MWE Current situation:

Expected result:

Note The effect becomes especially obvious for long filter lengths and if the signal gets appended by many zeros.

Requires detailled review of implemented new methods

• check normalization of created ideal autocorrelation
• check normalization of generated colored noise from white noise
• check execution of propagate-filter (compare sigma as float or np.ndarray)
• review tests

This pull request addresses issue #175 in some way.

By introducing a new function _fir_filter we can propagate full covariance information into the output of an FIR-filter. Based on this function a wrapper FIRuncFilter_2 is introduced, mimicing the behaviour of the existing FIRuncFilter. (And thus preparing a potential replacement lateron.)

Benefits of the new function(s):

• input can now have full covariance information as well
• return full covariance of output
• reduced complexity of the code
• no more python-for loops (all is done using 2D-convolution)
• reduce computations if Utheta == None or Ux is None (this was already done for the FIRuncFilter, but in a much more complex way)
• control over initial conditions of the FIR filter (at least limited)

A visual comparison with Monte Carlo covariance result show good agreement.

TODO: A quick runtime comparison of both methods should be done to evaluate potential performance gains/losses.

I just worked with the sine and multi_sine methods and the freq-argument gets used in different ways, than I would expect it:

In misc.testsignals.sine line 185 should be replaced by x = amp * np.sin(2 * np.pi * freq * time) (all multiplication). Alternativly, the freq parameter could be renamed to period.

In misc.testsignals.multi_sine line 213 should be replaced by x += amp * np.sin(2 * np.pi * freq * time) (scale sine by 2pi). Alternativly, the freq parameter could be renamed to angular_frequency.

To have a matching signature of both sine and multi_sine, I suggest to keep using the freq-argument, but adjust their use in the code as proposed.

Some references for the suggested changes can be found here: https://de.wikipedia.org/wiki/Sinuston https://en.wikipedia.org/wiki/Angular_frequency

Some digital sensors generate their own sample clock. This clock fluctuates considerably which leads to not equidistant sample distances. FFT, DFT and wavelet transformations need equidistant samples, so it would be nice if PyDynamic could implement interpolation methods with error poropagation.

Beste regards Benedikt

Apparently one of the tolerances in our test suite for checking ARMA() is too strict. See some of the most recent failed instances:

• https://app.circleci.com/pipelines/github/PTB-M4D/PyDynamic/2711/workflows/666c9563-ff15-4b60-892a-9d4380271275/jobs/14973/tests#failed-test-0
• https://app.circleci.com/pipelines/github/PTB-M4D/PyDynamic/2711/workflows/666c9563-ff15-4b60-892a-9d4380271275/jobs/14975/tests#failed-test-0
• https://app.circleci.com/pipelines/github/PTB-M4D/PyDynamic/2709/workflows/c4768521-6f49-4295-addc-1ec2fd4162b1/jobs/14961/tests#failed-test-0
• https://app.circleci.com/pipelines/github/PTB-M4D/PyDynamic/2705/workflows/9c5f9554-7840-449b-8ab1-85d7f2570dde/jobs/14938/tests#failed-test-0

It looks like we should increase to a value of atol = 0.22. The line above that assertion fails less often, but the following failed instance suggests to raise the tolerance for that as well:

• https://app.circleci.com/pipelines/github/PTB-M4D/PyDynamic/2704/workflows/164c507d-354e-441d-a5ea-d940d5e1d907/jobs/14929/tests#failed-test-0

Seems as if rtol = 0.22 might be a good choice as well.

For convenience, this implements the use of 1D-arrays for U1 and U2 in PyDynamic.uncertainty.convolve_unc to specify the standard uncertainties instead of full covariance matrices.

Specification of uncertainties using full covariance matrices remains possible and the output is still returned as full covariance matrix. Hence, no breaking change is introduced.

This PR makes required changes on:

• relevant test strategy
• docstring
• actual code change

I just updated to version v1.10.0. When I tried to check and verify, if I have new version (I also checked with same command that previous installed version was v1.9.0) I got an error message:

import PyDynamic PyDynamic.__version__

AttributeError Traceback (most recent call last) in 1 import PyDynamic ----> 2 PyDynamic.__version__

AttributeError: module 'PyDynamic' has no attribute '__version__'

Other attributes like "__doc__" or "__package__" are available.

Here we collect all preparations of the next major release with all those BREAKING CHANGES. After merging this we should have:

• [x] resolved #166
• [x] resolved #163
• [x] resolved #154
• [x] resolved #147 (which is handled by #189)
• [x] resolved #135 (which was handled by #188)
• [x] resolved #134 (which was handled by #180)
• [x] resolved #133 (which was handled by #158)
• [x] resolved #109 (which was handled by #143)
• [x] resolved #86 (which was handled by #142)
• [x] resolved #66 (which was handled by #187)
• [x] resolved #49
• [x] resolved #31 (which was handled by #181).

Describe the bug UMC_generic has an error in the calculation of the covariance matrix after the first block. The covariance therefore has a major mismatch to a basic MC implementation if only one (or a few) blocks are evaluated. If the number of blocks is big, this error wears off due to averaging and inherent random fluctuations.

To Reproduce

from PyDynamic.uncertainty.propagate_MonteCarlo import UMC_generic
import numpy as np
import matplotlib.pyplot as plt

# init
x = np.linspace(0,5,num=11)
ux = np.full_like(x, 0.5)
n_runs = 2000

draw_samples = lambda size: np.random.normal(x, ux, size=(size,x.size))
evaluate = lambda X: X + 1

# PyDynamic
y, Uy, _, output_shape = UMC_generic(draw_samples, evaluate=evaluate, runs=n_runs, blocksize=2000, n_cpu=1) 

# naiv Monte Carlo
results = []
for sample in draw_samples(n_runs):
    result = evaluate(sample)
    results.append(result)
y_mc = np.mean(results, axis=0)
Uy_mc = np.cov(results, rowvar=False)

# comparison
fig, ax = plt.subplots(1, 3)
ax[0].plot(y - y_mc)
ax[1].plot(np.diag(Uy) - np.diag(Uy_mc))
ax[2].imshow(Uy - Uy_mc)
plt.show()

# comparison
print("median y - y_mc   (ideal: 0.0): ", np.median(y - y_mc))
print("median Uy - Uy_mc (ideal: 0.0): ", np.median(Uy - Uy_mc))
print("median diag(Uy) / np.diag(Uy_mc) (ideal: 1.0): ", np.median(np.diag(Uy) / np.diag(Uy_mc)))

Exemplary Output:

median y - y_mc   (ideal: 0.0):  -0.014300805445340181
median Uy - Uy_mc (ideal: 0.0):  0.8981734069108699
median diag(Uy) / np.diag(Uy_mc) (ideal: 1.0):  1987.8304933280958

As can be seen, if only one block (of 2000 samples) is executed, the elements of the main diagonal are roughly factor 2000 apart from the basic MC implementation.

Potential Solution Here, the code should normalize by the number of samples in the current block.

Environment (please complete the following information):

• OS: Windows 10
• PyDynamic Version 2.3.0

At the moment we provide no shapes in the docstrings of this module. This is inconsistent with how we provide documentation for modules like interpolate.

Is your feature request related to a problem? Please describe. At the moment the signal class if not provided by the user, computes the sampling interval length as kind of a weighted sum of the different occuring interval legths in the time vector. The stored interval length is taken as the arithmetic mean of all unique interval lengths in the time vector. The sampling frequency if not provided then is taken as the reciprocal of that. This can result in kind of inconsistent data in an instances attributes, e.g. in case there is one single interval being much larger than all other possibly very short intervals. If you then only look at the frequency and the values, you might draw wrong conclusions. We do not see clearly though, what are the use cases of the two attributes.

Describe the solution you'd like We wish for a profound reasoning behind the way of computing the interval length and sampling frequency, which is then well documented in the corresponding part of the docs.

Describe alternatives you've considered One way would be to find the most reasonable mean (unweighted arithmetic/geometric, harmonic, median) of interval length and compute that for the provided time vector. One could as well compute the frequency and interval length only when needed (when is that?) and maybe even interpolate the values at the resulting time instances for these cases?!

Additional context We are grateful for any input on that from any expert user just as a comment in this issue to start with.

Meanwhile we have a Jupyter Notebook for the FIR example and an rst file containing some more formulas. We should combine all info in the notebook and delete the other file.

For instance in the module model_estimation.fit_filter we are converting back and forth those two formats and could probably improve performance by sticking to one form. This requires thorough checks of the formulas for the arithmetic operation though. Other modules are affected as well probably, e.g. the GUM2DFT related stuff.