<jupyter_start><jupyter_text>Poles/Zeros/Gain from Look-Up Tables (Prototyping)KNK Feb 2021: Working off of BSM's original ipynb, factoring out numeric methods<jupyter_code>from mth5_test_data.util import MTH5_TEST_DATA_DIR
print(MTH5_TEST_DATA_DIR)
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd

from matplotlib.gridspec import GridSpec
from plot_helpers import plot_response
from scipy import signal

#from tf_helpers import Fit_ZerosPolesGain_toFrequencyResponse_LM
from tf_helpers import Fit_ZerosPolesGain_toFrequencyResponse_LLSQ
from mt_metadata.timeseries.filters.plotting_helpers import plot_response
#from plot_helpers import plot_response
from iris_mt_scratch.sandbox.plot_helpers import plot_complex_response
from vectfit_python import vectfit
from zen import zen_response

%matplotlib inline

DEG2RAD = np.pi/180<jupyter_output><empty_output><jupyter_text>Load the data from data repo<jupyter_code>calibration_data_dir = os.path.join(MTH5_TEST_DATA_DIR, 'calibration_files', 'LEMI_Calibration', 'Lemi-039_N131BxByBz')
calibration_file_path = os.path.join(calibration_data_dir, 'Lemi-039_N131.cal' )
df = pd.read_csv(calibration_file_path, skiprows=3, delimiter='\t')
column_labels = ['frequency', 'amplitude_x', 'phase_x', 'amplitude_y', 'phase_y', 'amplitude_z', 'phase_z']
columns_map = dict(zip(df.columns, column_labels))
df = df.rename(columns=columns_map)
print(df.columns)<jupyter_output>Index(['frequency', 'amplitude_x', 'phase_x', 'amplitude_y', 'phase_y',
       'amplitude_z', 'phase_z'],
      dtype='object')<jupyter_text>Preliminary Response Curve<jupyter_code>x_phase_radians = DEG2RAD*df['phase_x']
x_phasor = np.exp(1.J*x_phase_radians) #(np.cos( x_phase_radians)+ 1.j*np.sin(x_phase_radians))
x_complex_response = df['amplitude_x'] * x_phasor

plot_complex_response(df['frequency'], x_complex_response)
#plt.semilogx(df['frequency'], deg2rad*df['phase_x'])<jupyter_output>2021-03-08T18:16:30 [line 141] numexpr.utils._init_num_threads - INFO: NumExpr defaulting to 8 threads.<jupyter_text>Estimate the response amplitude in the pass band<jupyter_code>nominal_passband_response_amplitude = df['amplitude_x'].iloc[10].mean()
print('nominal_passband_response_amplitude', nominal_passband_response_amplitude)<jupyter_output>nominal_passband_response_amplitude 0.02847<jupyter_text>Express in units of deciBels and plot on semilog<jupyter_code>amplitude_in_db = 20*np.log10(df['amplitude_x']/nominal_passband_response_amplitude)
plt.semilogx(df['frequency'], amplitude_in_db)
plt.semilogx(df['frequency'], amplitude_in_db, 'r*')
plt.grid(which='both')
plt.title('Amplitude Response of Instrument in dB')<jupyter_output><empty_output><jupyter_text>Note that the decay is log-linear after about 5Hz... Extend the line a bit further in log space to get an estimate of dB/decade<jupyter_code>#The last 13 points are approximately a straight line on the log plot

frq_line = df['frequency'][-13:]
ampl_line = df['amplitude_x'][-13:]
db_line = amplitude_in_db[-13:]
log_f = np.array(np.log10(frq_line))
log_a = np.array(np.log10(ampl_line))
db = np.array(db_line)
#plt.loglog(frq_line, ampl_line)
plt.plot(log_f, db_line)
plt.xlabel('log$_{10}$ frequency')
plt.ylabel('dB')
#(log_a[-1]-log_a[0]) / (log_f[-1]-log_f[0])
db_per_decade = (db[-1]-db[0]) / (log_f[-1]-log_f[0])
print("db_per_decade", db_per_decade)#
#We need to plot this line in terms of dB, on a log plot where we clearly mark either decades or octaves
#This line may need to be extended forward in frequency, which we can do by extrapolating the line (lin logspace)<jupyter_output>db_per_decade -21.8456644503<jupyter_text>looks like approximately -20dB/decade, which theoretically corresponds to a single pole.  We obtain a decent fit with one pole, but a better fit with 3poles and one zero.If we want to optimize we need to add constraints however.  In particular, a regulariization that forces 1: complex poles to be conjugates, and the 'odd' poles to be pure real.  This means modifying the LLSQ process.  Here is the solutions from Ben's ipynb using Linear Least Squares<jupyter_code>n_zeros = 1; n_poles = 3;
#n_zeros = 0; n_poles = 1;

print("PARAMETERIZATION: {} zeros, {} poles".format(n_zeros, n_poles))
w_obs = 2*np.pi*df['frequency']
resp_obs =  x_complex_response#df['amplitude_x']
LEMI_LM_ZPK = Fit_ZerosPolesGain_toFrequencyResponse_LLSQ(w_obs, 
                                                        resp_obs, 
                                                        n_zeros, n_poles, useSKiter=False, 
                                                         regularize=False)
print(LEMI_LM_ZPK)
plot_response(resp_obs=resp_obs, w_obs=w_obs, x_units='frequency', 
              zpk_pred=LEMI_LM_ZPK, w_values=w_obs)<jupyter_output>PARAMETERIZATION: 1 zeros, 3 poles
ZerosPolesGainContinuous(
array([-32.6930273+2.42939461j]),
array([-1032.30977561-43.56871821j,   -19.41778864-14.74673119j,
         -18.02197649+16.03155314j]),
(-526.39801470067096-20.959507362663881j),
dt: None
)
[ 179.79540659  179.75234024  179.70208872  179.64345864  179.57507195
  179.49527966  179.40217388  179.29356278  179.16684673  179.01898005
  178.8464565   178.64517016  178.410237    178.13603411  177.81581382
  177.44202253  177.0053953   176.49513373  175.89830217  175.19965376
  174.38086096  173.41946116  172.28806507  170.95272126  169.37066644
  167.48772299  165.23414341  162.52155145  159.23920498  155.25505894
  150.43229023  144.6666272   137.9610296   130.52651898  122.81176697
  115.41563443  108.85673225  103.40970581   99.0809631    95.70961046
   93.07178329   90.95164762   89.16554533   87.57102514   86.05923484
   84.54805184   82.97347641   81.28337588   79.43222352   77.3808073 ]<jupyter_text>We are getting complex valued zpk elements  Lets try a purely numerical approach using scipy optimize and some regularization that penalizes deviations from the fundamental theorem of algebra solutions If there is an odd number of poles/zeros there must be at least one purely real rootAll other roots must come in complex-conjugate pairsThis means that the sum of the complex components must be zero for both the poles, the zeros and their sum.Similarly the so-called gain factor (k) should be purely real valued.  Set up the  problem using the classic d=G(m) framework    -We are given an empirical resposne curve, which is a complex-valued function of frequency.  This corresponds to our data **d**    - The model **m** is composed of three elements: first there is a tuple of complex-numbers, one per zero, zecond is a tuple of complex-numbers, one per pole, and third is a scalar multiplier.    - G is the function linking **d** and **m**    - Ben cast this as a linear problem but in this approach I am going to just search the model space using a direct evaluation of the zpk formulaTo define the model vector we need to first define the solution _order_.  The order is defined by two integers  (n_zeros, n_poles).  Generally n_poles >= n_zeros.The entire system is parameterized by the Laplace variable _s_  = *j* $\omega$Ben parameterized the system like this:$ \displaystyle h(s) = \frac{a_{m}s^m + a_{m-1}s^{m-1} + ... + a_{1}s^{1} + a_{0}}{b_{n}s^{n} + b_{n-1}s^{n-1} + ... + b_{1}s^{1} + b_{0}} $We are going to use$ \displaystyle d (s) = k  \frac{(s-z_1)(s-z_2)...(s-z_m) }{(s-p_1)(s-p_2)...(s-p_n)} $<jupyter_code>def add_regularization_term(zeros, poles, k, misfit, aa=0.1, verbose=False):
    if verbose:
        print('misfit in', misfit)
    #punish unbalanced imaginary components
    misfit +=  aa * np.abs(np.mean(np.sum( np.imag(zeros) )) )
    misfit += np.abs(np.mean(np.sum(np.imag(poles))))
    misfit +=  aa * np.abs(np.imag(k))
    #punish poles on rhs of plane
    for pole in poles:
        if np.real(pole)>0:
            misfit += aa * np.abs(np.real(pole) )
    if verbose:
        print('misfit out', misfit)
    return misfit
X_COMPLEX_RESPONSE = x_complex_response
f = df['frequency'].to_numpy()
W = 2*np.pi*f
S = 1.j*W
N_ZEROS = 1
N_POLES = 3
z0 = np.zeros(N_ZEROS, dtype=np.complex128)
p0 = np.zeros(N_POLES, dtype=np.complex128)
k0 = 1.0
MODEL_VECTOR_0 = np.hstack((z0, p0, k0))
print('MODEL_VECTOR_0',MODEL_VECTOR_0)
#Now we need a function that takes the MODEL_VECTOR as input and creates a lambda function for h(s)

def objective_function(model_vector, verbose=False):
    #unpack the model vector
    zeros = model_vector[:N_ZEROS]
    poles = model_vector[N_ZEROS:N_ZEROS+N_POLES]
    k = model_vector[-1]
    if verbose:
        print('zeros', zeros)
        print('poles', poles)
        print('k',k)
    num, den = signal.zpk2tf(zeros, poles, k)
    w, h = signal.freqs_zpk(zeros, poles, k, worN=W)
    residual = X_COMPLEX_RESPONSE - h
    misfit = np.sqrt(np.mean(np.abs(residual**2)))
    misfit = add_regularization_term(zeros, poles, k, misfit, aa=0.1)
    print('misfit',misfit)
    return misfit

objective_function(MODEL_VECTOR_0, verbose=True)
from scipy.optimize import minimize
method = 'nelder-mead'
#method = 'Powell'
#method = 'CG'
#method = 'Newton-CG'
res = minimize(objective_function, MODEL_VECTOR_0, method=method, options={'xatol': 1e-10, 'disp': True, })
res.x<jupyter_output><empty_output><jupyter_text>Looks like we are not getting back complex values here... thich could<jupyter_code>zeros = res.x[:N_ZEROS]
poles = res.x[N_ZEROS:N_ZEROS+N_POLES]
k = res.x[-1]
def make_zpk_from_model_vector(model_vector):
    zeros_real = model_vector[:2*N_ZEROS]
    poles_real = model_vector[2*N_ZEROS:2*N_ZEROS+2*N_POLES]
    k = model_vector[-1]
    zeros = np.zeros(N_ZEROS, dtype=np.complex128)
    for i_z in range(N_ZEROS):
        zeros[i_z] = np.complex(zeros_real[0+2*i_z], zeros_real[1+2*i_z])
    poles = np.zeros(N_POLES, dtype=np.complex128)
    for i_p in range(N_POLES):
        poles[i_p] = np.complex(poles_real[0+2*i_p], poles_real[1+2*i_p])
    return zeros, poles, k
        
X_COMPLEX_RESPONSE = x_complex_response
f = df['frequency'].to_numpy()
W = 2*np.pi*f
S = 1.j*W
N_ZEROS = 1
N_POLES = 3
z0 = np.zeros(2*N_ZEROS)#, dtype=np.complex128)
p0 = np.zeros(2*N_POLES)#, dtype=np.complex128)
k0 = 1.0
MODEL_VECTOR_0 = np.hstack((z0, p0, k0))
print(MODEL_VECTOR_0)
print(MODEL_VECTOR_0.dtype)
#Now we need a function that takes the MODEL_VECTOR as input and creates a lambda function for h(s)

def objective_function(model_vector):
    zeros, poles, k = make_zpk_from_model_vector(model_vector)
#     zeros_real = model_vector[:2*N_ZEROS]
#     poles_real = model_vector[2*N_ZEROS:2*N_ZEROS+2*N_POLES]
#     #print('poles_real', poles_real)
#     k = model_vector[-1]
#     #print('ok')
#     zeros = np.zeros(N_ZEROS, dtype=np.complex128)
#     for i_z in range(N_ZEROS):
#         zeros[i_z] = np.complex(zeros_real[0+2*i_z], zeros_real[1+2*i_z])
#     poles = np.zeros(N_POLES, dtype=np.complex128)
#     for i_p in range(N_POLES):
#     #    print(i_p)
#         poles[i_p] = np.complex(poles_real[0+2*i_p], poles_real[1+2*i_p])
        
        
#     #print('zeros',zeros)
#     #print(poles)
#     #print(k)
    num, den = signal.zpk2tf(zeros, poles, k)
    w, h = signal.freqs_zpk(zeros, poles, k, worN=W)
    #print(h)
    #print(X_COMPLEX_RESPONSE)
    residual = X_COMPLEX_RESPONSE - h
    #print(residual)
    misfit = np.sqrt(np.mean(np.abs(residual**2)))
    misfit = add_regularization_term(zeros, poles, k, misfit, aa=100.0)
    
    print(misfit)
    return misfit

#h = lambda:s
objective_function(MODEL_VECTOR_0)
from scipy.optimize import minimize
method = 'nelder-mead'
#method = 'Powell'
#method = 'CG'
#method = 'Newton-CG'
res = minimize(objective_function, MODEL_VECTOR_0, method=method, options={'xatol': 1e-10, 'disp': True, })
res.x
zeros, poles, k = make_zpk_from_model_vector(res.x)
print('zeros', zeros)
print('poles', poles)
print('k', k)
plot_response(resp_obs=resp_obs, w_obs=w_obs, x_units='frequency', 
              zpk_pred=LEMI_LM_ZPK, w_values=w_obs)
NIMS_LM_ZPK.poles =np.array([-1032.9,   -18.71-15.37j,
         -18.71+15.37j])
NIMS_LM_ZPK.zeros[0] = -32.78<jupyter_output><empty_output><jupyter_text>Non-Linear Least Squares FittingUsing nonlinear least squares with Levenberg-Marquardt step length damping to fit poles and zeros representation to look-up table of complex system response. Nominally, the target function is of the form             $ \displaystyle zpk(s) = k \frac{(s-z_{1})(s-z_{2})(s-z_{3})...(s-z_{m})}{(s-p_{1})(s-p_{2})(s-p_{3})...(s-p_{n})}$where s is complex frequency ($s = i\omega$, where $\omega$ is angular frequency), $k$ is the overall system gain, $z$ represents the $m$ system zeros, and $p$ represents the $n$ system poles.   HOWEVER. We're actually going to fit the look-up table with the transfer function form of the system response function:        $ \displaystyle h(s) = \frac{a_{m}s^m + a_{m-1}s^{m-1} + ... + a_{1}s^{1} + a_{0}}{b_{n}s^{n} + b_{n-1}s^{n-1} + ... + b_{1}s^{1} + b_{0}} $where $ a_{x} $ and $ b_{x} $ are the (real) coefficients on the numerator and denominator polynomials, $m$ is the number of system zeros (also the degree of the numerator polynomial), $n$ is the number of system poles (also the degree of the denominator polynomial), and $s$ is complex frequency ($s = i\omega$, where $\omega$ in angular frequency).This form has several benefits. First of all, because the coefficients are real, the roots of these polynomials (and therefore the poles/zeros) will occur as complex conjugate pairs, which is what we expect/require for a real linear system operating on a real time series (I think...). Having real coefficients in this polynomial form also makes the math a bit easier.Then, to find the poles/zeros, we just have to factor the numerator and denominator of this rational function. Operationally, we'll just use Python tools to do so.For the nonlinear least squares fitting, the Jacobian is fairly simple to calculate:$ \displaystyle \frac{\partial h(s)}{\partial a_{m}} = \frac{s^{m}}{b_{n}s^{n} + b_{n-1}s^{n-1} + ... + b_{1}s^{1} + b_{0}} = \frac{s^{m}}{d(s)} $(Other partial derivatives with respect to the numerator terms are similar.)$ \displaystyle \frac{\partial h(s)}{\partial b_{n}} = -\frac{a_{m}s^m + a_{m-1}s^{m-1} + ... + a_{1}s^{1} + a_{0}}{(b_{n}s^{n} + b_{n-1}s^{n-1} + ... + b_{1}s^{1} + b_{0})^{2}} s^{n} = -\frac{n(s)}{d(s)^{2}}s^{n} $(Other partial derivatives with respect to the denominator terms are similar.)Here, $n(s)$ and $d(s)$ are the numerator and denominator polynomials, respectively.In the matrices here, all the numerator terms will be listed first, then all the denominator terms. **Nonlinear LM fit to NIMS B response**<jupyter_code>print("PARAMETERIZATION: 0 zeros, 3 poles")
test_period = np.logspace(-1, 5, num=100)
test_w = 2. * np.pi / test_period
NIMS_w, NIMS_Bresp = signal.freqresp(NIMS_Magnetic_3PoleLowpass, 
                                     w=test_w)   # w is in rad/sec
NIMS_LM_ZPK = Fit_ZerosPolesGain_toFrequencyResponse_LM(NIMS_w, 
                                                        NIMS_Bresp, 
                                                        0, 3)
plot_response(zpk_obs=NIMS_Magnetic_3PoleLowpass, zpk_pred=NIMS_LM_ZPK, 
              w_values=test_w)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12203791469194313, 'fraction_numerical': 0.061940494997367035, 'mean_word_length': 4.334620786516854, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '407615', 'license': 'MIT', 'n_tokens_mistral': 6269, 'n_tokens_neox': 5622, 'n_words': 1393, 'path': 'iris_mt_scratch/sandbox/poles_and_zeros/lemi_fitting_backup_20210308.ipynb', 'repo_name': 'simpeg-research/iris-mt-scratch', 'size': '15192'}