<commit_before>"""
=========================================
VectorFitting (:mod:`skrf.vectorFitting`)
=========================================

.. autoclass:: VectorFitting
    :members:
    :member-order: groupwise
    :special-members: __init__

"""

import numpy as np
import os
from functools import wraps
try:
    from . import plotting    # will perform the correct setup for matplotlib before it is called below
    import matplotlib.pyplot as mplt
    from matplotlib.ticker import EngFormatter
except ImportError:
    mplt = None

import logging


def check_plotting(func):
    """
    This decorator checks if matplotlib.pyplot is available under the name mplt.
    If not, raise an RuntimeError.

    Raises
    ------
    RuntimeError
        When trying to run the decorated function without matplotlib
    """

    @wraps(func)
    def wrapper(*args, **kwargs):
        if mplt is None:
            raise RuntimeError('Plotting is not available')
        func(*args, **kwargs)

    return wrapper


class VectorFitting:
    """
    =========================================================
    VectorFitting (:class:`skrf.vectorFitting.VectorFitting`)
    =========================================================

    This class provides a Python implementation of the Vector Fitting algorithm and various functions for the fit
    analysis and export of SPICE equivalent circuits.

    Notes
    -----
    The fitting code is based on the original algorithm [1]_ and on two improvements for relaxed pole relocation [2]_
    and efficient (fast) solving [3]_. See also the Vector Fitting website [4]_ for further information and download of
    the papers listed below. A Matlab implementation is also available there for reference.

    References
    ----------
    .. [1] <NAME>, <NAME>, "Rational Approximation of Frequency Domain Responses by Vector Fitting", IEEE
        Transactions on Power Delivery, vol. 14, no. 3, pp. 1052-1061, July 1999, DOI: https://doi.org/10.1109/61.772353

    .. [2] <NAME>, "Improving the Pole Relocating Properties of Vector Fitting", IEEE Transactions on Power
        Delivery, vol. 21, no. 3, pp. 1587-1592, July 2006, DOI: https://doi.org/10.1109/TPWRD.2005.860281

    .. [3] <NAME>, <NAME>, <NAME>, <NAME>, "Marcomodeling of Multiport Systems Using a Fast
        Implementation of the Vector Fitting Method", IEEE Microwave and Wireless Components Letters, vol. 18, no. 6,
        pp. 383-385, June 2008, DOI: https://doi.org/10.1109/LMWC.2008.922585

    .. [4] Vector Fitting website: https://www.sintef.no/projectweb/vectorfitting/
    """

    def __init__(self, network):
        """
        Creates a VectorFitting instance based on a supplied :class:`skrf.network.Network` containing the frequency
        responses of the N-port.

        Parameters
        ----------
        network : :class:`skrf.network.Network`
            Network instance of the N-port holding the S-matrix to be fitted.
        """

        self.network = network
        self.initial_poles = None

        self.poles = None
        """ Instance variable holding the list of fitted poles. Will be initialized by :func:`vector_fit`. """

        self.zeros = None
        """ Instance variable holding the list of fitted zeros. Will be initialized by :func:`vector_fit`. """

        self.proportional_coeff = None
        """ Instance variable holding the list of fitted proportional coefficients. Will be initialized by 
        :func:`vector_fit`. """

        self.constant_coeff = None
        """ Instance variable holding the list of fitted constants. Will be initialized by :func:`vector_fit`. """

        self.max_iterations = 100
        """ Instance variable specifying the maximum number of iterations for the fitting process. To be changed by the
         user before calling :func:`vector_fit`. """

        self.max_tol = 1e-6
        """ Instance variable specifying the convergence criterion in terms of relative tolerance. To be changed by the
         user before calling :func:`vector_fit`. """

        self.d_res_history = []
        self.delta_max_history = []
        self.history_max_sigma = []

    def vector_fit(self, n_poles_real=2, n_poles_cmplx=2, init_pole_spacing='lin', parameter_type='S',
                   fit_constant=True, fit_proportional=False):
        """
        Main work routine performing the vector fit. The results will be stored in the class variables
        :attr:`poles`, :attr:`zeros`, :attr:`proportional_coeff` and :attr:`constant_coeff`.

        Parameters
        ----------
        n_poles_real : int, optional
            Number of initial real poles. See notes.

        n_poles_cmplx : int, optional
            Number of initial complex conjugate poles. See notes.

        init_pole_spacing : str, optional
            Type of initial pole spacing across the frequency interval of the S-matrix. Either linear (lin) or
            logarithmic (log).

        parameter_type : str, optional
            Representation type of the frequency responses to be fitted. Either *scattering* (:attr:`s` or :attr:`S`),
            *impedance* (:attr:`z` or :attr:`Z`) or *admittance* (:attr:`y` or :attr:`Y`). As scikit-rf can currently
            only read S parameters from a Touchstone file, the fit should also be performed on the original S
            parameters. Otherwise, scikit-rf will convert the responses from S to Z or Y, which might work for the fit
            but can cause other issues.

        fit_constant : bool, optional
            Include a constant term **d** in the fit.

        fit_proportional : bool, optional
            Include a proportional term **e** in the fit.

        Returns
        -------
        None

        Notes
        -----
        The required number of real or complex conjugate starting poles depends on the behaviour of the frequency
        responses. To fit a smooth response such as a low-pass characteristic, 1-3 real poles and no complex conjugate
        poles is usually sufficient. If resonances or other types of peaks are present in some or all of the responses,
        a similar number of complex conjugate poles is required. Be careful not to use too many poles, as excessive
        poles will not only increase the computation workload during the fitting and the subsequent use of the model,
        but they can also introduce unwanted resonances at frequencies well outside the fit interval.
        """

        # create initial poles and space them across the frequencies in the provided Touchstone file
        # use normalized frequencies during the iterations (seems to be more stable during least-squares fit)
        norm = np.average(self.network.f)
        freqs_norm = self.network.f / norm

        fmin = np.amin(freqs_norm)
        fmax = np.amax(freqs_norm)
        weight_regular = 1.0
        if init_pole_spacing == 'log':
            pole_freqs_real = np.geomspace(fmin, fmax, n_poles_real)
            pole_freqs_cmplx = np.geomspace(fmin, fmax, n_poles_cmplx)
        elif init_pole_spacing == 'lin':
            pole_freqs_real = np.linspace(fmin, fmax, n_poles_real)
            pole_freqs_cmplx = np.linspace(fmin, fmax, n_poles_cmplx)
        else:
            logging.warning('Invalid choice of initial pole spacing; proceeding with linear spacing')
            pole_freqs_real = np.linspace(fmin, fmax, n_poles_real)
            pole_freqs_cmplx = np.linspace(fmin, fmax, n_poles_cmplx)
        poles = []

        # add real poles
        for i, f in enumerate(pole_freqs_real):
            omega = 2 * np.pi * f
            poles.append((-1 / 100 + 0j) * omega)
        # add complex-conjugate poles (store only positive imaginary parts)
        for i, f in enumerate(pole_freqs_cmplx):
            omega = 2 * np.pi * f
            poles.append((-1 / 100 + 1j) * omega)
        poles = np.array(poles)

        # save initial poles (un-normalize first)
        self.initial_poles = poles * norm
        max_singular = 1

        logging.info('### Starting pole relocation process.\n')

        # stack frequency responses as a single vector
        # stacking order:
        # s11, s12, s13, ..., s21, s22, s23, ...
        freq_responses = []
        for i in range(self.network.nports):
            for j in range(self.network.nports):
                if parameter_type.lower() == 's':
                    freq_responses.append(self.network.s[:, i, j])
                elif parameter_type.lower() == 'z':
                    freq_responses.append(self.network.z[:, i, j])
                elif parameter_type.lower() == 'y':
                    freq_responses.append(self.network.y[:, i, j])
                else:
                    logging.warning('Invalid choice of matrix parameter type (S, Z, or Y); proceeding with S representation.')
                    freq_responses.append(self.network.s[:, i, j])
        freq_responses = np.array(freq_responses)

        # ITERATIVE FITTING OF POLES to the provided frequency responses
        # initial set of poles will be replaced with new poles after every iteration
        iterations = self.max_iterations
        self.d_res_history = []
        self.delta_max_history = []
        converged = False
        while iterations > 0:
            logging.info('Iteration {}'.format(self.max_iterations - iterations + 1))

            # generate coefficients of approximation function for each target frequency response
            # responses will be treated independently using QR decomposition
            # simplified coeff matrices of all responses will be stacked for least-squares solver
            A_matrix = []
            b_vector = []

            for freq_response in freq_responses:
                # calculate coefficients (row A_k in matrix) for each frequency sample s_k of the target response
                # row will be appended to submatrix A_sub of complete coeff matrix A_matrix
                # 2 rows per pole in result vector (1st for real part, 2nd for imaginary part)
                # --> 2 columns per pole in coeff matrix
                A_sub = []
                b_sub = []
                n_unused = 0

                for k, f_sample in enumerate(freqs_norm):
                    s_k = 2j * np.pi * f_sample
                    A_k = []
                    n_unused = 0

                    # add coefficients for a pair of complex conjugate poles
                    # part 1: first sum of rational functions (residue variable c)
                    for pole in poles:
                        # separate and stack real and imaginary part to preserve conjugacy of the pole pair
                        if np.imag(pole) == 0.0:
                            A_k.append(1 / (s_k - pole))
                            n_unused += 1
                        else:
                            # complex pole of a conjugated pair
                            A_k.append(1 / (s_k - pole) + 1 / (s_k - np.conjugate(pole)))       # real part of residue
                            A_k.append(1j / (s_k - pole) - 1j / (s_k - np.conjugate(pole)))     # imaginary part of residue
                            n_unused += 2

                    # part 2: constant (variable d) and proportional term (variable e)
                    if fit_constant:
                        A_k.append(1)
                        n_unused += 1
                    if fit_proportional:
                        A_k.append(s_k)
                        n_unused += 1

                    # part 3: second sum of rational functions (variable c_res)
                    for pole in poles:
                        # separate and stack real and imaginary part to preserve conjugacy of the pole pair
                        if np.imag(pole) == 0.0:
                            A_k.append(-1 * freq_response[k] / (s_k - pole))
                        else:
                            # complex pole of a conjugated pair
                            A_k.append(-1 * freq_response[k] / (s_k - pole)
                                       - freq_response[k] / (s_k - np.conjugate(pole)))
                            A_k.append(-1j * freq_response[k] / (s_k - pole)
                                       + 1j * freq_response[k] / (s_k - np.conjugate(pole)))

                    # part 4: constant (variable d_res)
                    A_k.append(-1 * freq_response[k])

                    A_sub.append(np.sqrt(weight_regular) * np.real(A_k))
                    A_sub.append(np.sqrt(weight_regular) * np.imag(A_k))
                    b_sub.append(np.sqrt(weight_regular) * 0.0)
                    b_sub.append(np.sqrt(weight_regular) * 0.0)

                # QR decomposition
                Q, R = np.linalg.qr(A_sub, 'reduced')

                # only R22 is required to solve for c_res and d_res
                R22 = R[n_unused:, n_unused:]
                # similarly, only right half of Q is required
                Q2 = Q[:, n_unused:]

                if len(A_matrix) == 0:
                    A_matrix = R22
                    b_vector = np.matmul(np.transpose(Q2), b_sub)
                else:
                    A_matrix = np.vstack((A_matrix, R22))
                    b_vector = np.append(b_vector, np.matmul(np.transpose(Q2), b_sub))

                # add extra equation to avoid trivial solution
                # use weight=1 for all equations, except for this extra equation
                weight_extra = np.linalg.norm(weight_regular * freq_response) / len(freq_response)
                A_k = np.zeros(np.shape(A_matrix)[1])
                for k, f_sample in enumerate(freqs_norm):
                    s_k = 2j * np.pi * f_sample
                    i = 0

                    # part 3: second sum of rational functions (variable c_res)
                    for pole in poles:
                        if np.imag(pole) == 0.0:
                            # real pole
                            A_k[i] += np.real(1 / (s_k - pole))
                            i += 1
                        else:
                            # complex pole of a conjugated pair
                            A_k[i] += np.real(1 / (s_k - pole) + 1 / (s_k - np.conjugate(pole)))
                            A_k[i+1] += np.real(1j / (s_k - pole) - 1j / (s_k - np.conjugate(pole)))
                            i += 2
                    # part 4: constant (d_res)
                    A_k[i] += 1

                A_matrix = np.vstack((A_matrix, np.sqrt(weight_extra) * A_k))
                b_vector = np.append(b_vector, np.sqrt(weight_extra) * len(freqs_norm))

            logging.info('A_matrix: condition number = {}'.format(np.linalg.cond(A_matrix)))

            # solve least squares for real parts
            x, residuals, rank, singular_vals = np.linalg.lstsq(A_matrix, b_vector, rcond=None)

            # assemble individual result vectors from single LS result x
            c_res = x[:-1]
            d_res = x[-1]

            # check if d_res is suited for zeros calculation
            tol_res = 1e-8
            if np.abs(d_res) < tol_res:
                # d_res is too small, discard solution and proceed the |d_res| = tol_res
                d_res = tol_res * (d_res / np.abs(d_res))
                logging.warning('Replacing d_res solution as it was too small')

            self.d_res_history.append(d_res)
            logging.info('d_res = {}'.format(d_res))

            # build test matrix H, which will hold the new poles as Eigenvalues
            A_matrix = np.zeros((len(c_res), len(c_res)))
            i = 0
            for pole in poles:
                # fill diagonal with previous poles
                if np.imag(pole) == 0.0:
                    # one row for a real pole
                    A_matrix[i, i] = np.real(pole)
                    A_matrix[i] -= c_res / d_res
                    i += 1
                else:
                    # two rows for a complex pole of a conjugated pair
                    A_matrix[i, i] = np.real(pole)
                    A_matrix[i, i + 1] = np.imag(pole)
                    A_matrix[i + 1, i] = -1 * np.imag(pole)
                    A_matrix[i + 1, i + 1] = np.real(pole)
                    A_matrix[i] -= 2 * c_res / d_res
                    i += 2
            poles_new = np.linalg.eigvals(A_matrix)

            # replace poles for next iteration
            poles = []
            for k, pole in enumerate(poles_new):
                # flip real part of unstable poles (real part needs to be negative for stability)
                if np.real(pole) > 0.0:
                    pole = -1 * np.real(pole) + 1j * np.imag(pole)
                if np.imag(pole) >= 0.0:
                    # complex poles need to come in complex conjugate pairs; append only the positive part
                    poles.append(pole)

            poles = np.sort_complex(poles)

            # calculate relative changes in the singular values; stop iteration loop once poles have converged
            new_max_singular = np.amax(singular_vals)
            delta_max = np.abs(1 - new_max_singular / max_singular)
            self.delta_max_history.append(delta_max)
            logging.info('Max. relative change in residues = {}\n'.format(delta_max))
            max_singular = new_max_singular

            stop = False
            if delta_max < self.max_tol:
                if converged:
                    # is really converged, finish
                    logging.info('Pole relocation process converged after {} iterations.'.format(
                        self.max_iterations - iterations + 1))
                    stop = True
                else:
                    # might be converged, but do one last run to be sure
                    converged = True
            else:
                if converged:
                    # is not really converged, continue
                    converged = False

            iterations -= 1

            if iterations == 0:
                if converged and stop is False:
                    logging.warning('Reached tolerance only after max. number of iterations (N_max = {}). '
                                    'Results might not have converged properly.'.format(self.max_iterations))
                else:
                    logging.warning('Reached maximum number of iterations (N_max = {}). '
                                    'Results did not converge.'.format(self.max_iterations))

            if stop:
                iterations = 0

        # ITERATIONS DONE
        poles = np.array(poles)

        logging.info('Initial poles before relocation:')
        logging.info(self.initial_poles)

        logging.info('Final poles:')
        logging.info(poles * norm)

        logging.info('\n### Starting zeros calculation process.\n')

        # finally, solve for the residues with the previously calculated poles
        zeros = []
        constant_coeff = []
        proportional_coeff = []

        for freq_response in freq_responses:
            # calculate coefficients (row A_k in matrix) for each frequency sample s_k of the target response
            # row will be appended to submatrix A_sub of complete coeff matrix A_matrix
            # 2 rows per pole in result vector (1st for real part, 2nd for imaginary part)
            # --> 2 columns per pole in coeff matrix
            A_matrix = []
            b_vector = []

            for k, f_sample in enumerate(freqs_norm):
                s_k = 2j * np.pi * f_sample
                A_k = []
                # add coefficients for a pair of complex conjugate poles
                # part 1: first sum of rational functions (residue variable c)
                for pole in poles:
                    # separate and stack real and imaginary part to preserve conjugacy of the pole pair
                    if np.imag(pole) == 0.0:
                        A_k.append(1 / (s_k - pole))
                    else:
                        A_k.append(1 / (s_k - pole) + 1 / (s_k - np.conjugate(pole)))       # real part of residue
                        A_k.append(1j / (s_k - pole) - 1j / (s_k - np.conjugate(pole)))     # imaginary part of residue

                # part 2: constant (variable d) and proportional term (variable e)
                if fit_constant:
                    A_k.append(1.0)
                if fit_proportional:
                    A_k.append(s_k)

                A_matrix.append(np.array(A_k))
                b_vector.append(np.array(freq_response[k]))

            A_matrix = np.vstack((np.real(A_matrix), np.imag(A_matrix)))
            b_vector = np.append(np.real(b_vector), np.imag(b_vector))

            logging.info('A_matrix: condition number = {}'.format(np.linalg.cond(A_matrix)))

            # solve least squares and obtain results as stack of real part vector and imaginary part vector
            x, residuals, rank, singular_vals = np.linalg.lstsq(A_matrix, b_vector, rcond=None)

            i = 0
            zeros_response = []
            for pole in poles:
                if np.imag(pole) == 0.0:
                    zeros_response.append(x[i] + 0j)
                    i += 1
                else:
                    zeros_response.append(x[i] + 1j * x[i + 1])
                    i += 2
            zeros.append(zeros_response)

            if fit_constant and fit_proportional:
                # both constant d and proportional e were fitted
                constant_coeff.append(x[-2])
                proportional_coeff.append(x[-1])
            elif fit_constant:
                # only constant d was fitted
                constant_coeff.append(x[-1])
                proportional_coeff.append(0.0)
            elif fit_proportional:
                # only proportional e was fitted
                constant_coeff.append(0.0)
                proportional_coeff.append(x[-1])
            else:
                # neither constant d nor proportional e was fitted
                constant_coeff.append(0.0)
                proportional_coeff.append(0.0)

        # save poles, zeros, d, e in actual frequencies (un-normalized)
        self.poles = poles * norm
        self.zeros = np.array(zeros) * norm
        self.constant_coeff = np.array(constant_coeff)
        self.proportional_coeff = np.array(proportional_coeff) / norm

        logging.info('\n### Vector fitting finished.\n')

    def _get_ABCDE(self):
        """
        Private method.
        Returns the real-valued system matrices of the state-space representation of the current rational model, as
        defined in [1]_.

        Returns
        -------
        A : ndarray
            State-space matrix A holding the poles on the diagonal as real values with imaginary parts on the sub-
            diagonal
        B : ndarray
            State-space matrix B holding coefficients (1, 2, or 0), depending on the respective type of pole in A
        C : ndarray
            State-space matrix C holding the residues (zeros)
        D : ndarray
            State-space matrix D holding the constants
        E : ndarray
            State-space matrix E holding the proportional coefficients (usually 0 in case of fitted S-parameters)

        References
        ----------
        .. [1] <NAME> and <NAME>, "Fast Passivity Assessment for S-Parameter Rational Models Via a Half-Size
            Test Matrix," in IEEE Transactions on Microwave Theory and Techniques, vol. 56, no. 12, pp. 2701-2708,
            Dec. 2008, DOI: 10.1109/TMTT.2008.2007319.
        """

        # initial checks
        if self.poles is None:
            logging.error('self.poles = None; nothing to do. You need to run vector_fit() first.')
            return
        if self.zeros is None:
            logging.error('self.zeros = None; nothing to do. You need to run vector_fit() first.')
            return
        if self.proportional_coeff is None:
            logging.error('self.proportional_coeff = None; nothing to do. You need to run vector_fit() first.')
            return
        if self.constant_coeff is None:
            logging.error('self.constant_coeff = None; nothing to do. You need to run vector_fit() first.')
            return

        # assemble real-valued state-space matrices A, B, C, D, E from fitted complex-valued pole-residue model

        # determine size of the matrix system
        n_ports = int(np.sqrt(len(self.constant_coeff)))
        n_poles_real = 0
        n_poles_cplx = 0
        for pole in self.poles:
            if np.imag(pole) == 0.0:
                n_poles_real += 1
            else:
                n_poles_cplx += 1
        n_matrix = (n_poles_real + 2 * n_poles_cplx) * n_ports

        # state-space matrix A holds the poles on the diagonal as real values with imaginary parts on the sub-diagonal
        # state-space matrix B holds coefficients (1, 2, or 0), depending on the respective type of pole in A
        # assemble A = [[poles_real,   0,                  0],
        #               [0,            real(poles_cplx),   imag(poles_cplx],
        #               [0,            -imag(poles_cplx),  real(poles_cplx]]
        A = np.identity(n_matrix)
        B = np.zeros(shape=(n_matrix, n_ports))
        i_A = 0  # index on diagonal of A
        for j in range(n_ports):
            for pole in self.poles:
                if np.imag(pole) == 0.0:
                    # adding a real pole
                    A[i_A, i_A] = np.real(pole)
                    B[i_A, j] = 1
                    i_A += 1
                else:
                    # adding a complex-conjugate pole
                    A[i_A, i_A] = np.real(pole)
                    A[i_A, i_A + 1] = np.imag(pole)
                    A[i_A + 1, i_A] = -1 * np.imag(pole)
                    A[i_A + 1, i_A + 1] = np.real(pole)
                    B[i_A, j] = 2
                    i_A += 2

        # state-space matrix C holds the residues (zeros)
        # assemble C = [[R1.11, R1.12, R1.13, ...], [R2.11, R2.12, R2.13, ...], ...]
        C = np.zeros(shape=(n_ports, n_matrix))
        for i in range(n_ports):
            for j in range(n_ports):
                # i: row index
                # j: column index
                i_response = i * n_ports + j

                j_zeros = 0
                for zero in self.zeros[i_response]:
                    if np.imag(zero) == 0.0:
                        C[i, j * (n_poles_real + 2 * n_poles_cplx) + j_zeros] = np.real(zero)
                        j_zeros += 1
                    else:
                        C[i, j * (n_poles_real + 2 * n_poles_cplx) + j_zeros] = np.real(zero)
                        C[i, j * (n_poles_real + 2 * n_poles_cplx) + j_zeros + 1] = np.imag(zero)
                        j_zeros += 2

        # state-space matrix D holds the constants
        # assemble D = [[d11, d12, ...], [d21, d22, ...], ...]
        D = np.zeros(shape=(n_ports, n_ports))
        for i in range(n_ports):
            for j in range(n_ports):
                # i: row index
                # j: column index
                i_response = i * n_ports + j
                D[i, j] = self.constant_coeff[i_response]

        # state-space matrix E holds the proportional coefficients (usually 0 in case of fitted S-parameters)
        # assemble E = [[e11, e12, ...], [e21, e22, ...], ...]
        E = np.zeros(shape=(n_ports, n_ports))
        for i in range(n_ports):
            for j in range(n_ports):
                # i: row index
                # j: column index
                i_response = i * n_ports + j
                E[i, j] = self.proportional_coeff[i_response]

        return A, B, C, D, E

    @staticmethod
    def _get_s_from_ABCDE(freq, A, B, C, D, E):
        """
        Private method.
        Returns the S-matrix of the vector fitted model calculated from the real-valued system matrices of the state-
        space representation, as provided by `_get_ABCDE()`.

        Parameters
        ----------
        freq : float
            Frequency (in Hz) at which to calculate the S-matrix.
        A : ndarray
        B : ndarray
        C : ndarray
        D : ndarray
        E : ndarray

        Returns
        -------
        ndarray
            Complex-valued S-matrix (NxN) calculated at frequency `freq`.
        """

        dim_A = np.shape(A)[0]
        stsp_poles = np.linalg.inv(2j * np.pi * freq * np.identity(dim_A) - A)
        stsp_S = np.matmul(np.matmul(C, stsp_poles), B)
        stsp_S += D + 2j * np.pi * freq * E
        return stsp_S

    def passivity_test(self, parameter_type='S'):
        """
        Evaluates the passivity of reciprocal vector fitted models by means of a half-size test matrix [1]_. Any
        existing frequency bands of passivity violations will be returned as a sorted list.

        Parameters
        ----------
        parameter_type: str, optional
            Representation type of the fitted frequency responses. Either *scattering* (:attr:`s` or :attr:`S`),
            *impedance* (:attr:`z` or :attr:`Z`) or *admittance* (:attr:`y` or :attr:`Y`). Currently, only scattering
            parameters are supported for passivity evaluation.

        Returns
        -------
        ndarray
            NumPy array with frequency bands of passivity violation:
            [[f_start_1, f_stop_1], [f_start_2, f_stop_2], ...].

        See Also
        --------
        is_passive : Query the model passivity as a boolean value.

        References
        ----------
        .. [1] <NAME> and <NAME>, "Fast Passivity Assessment for S-Parameter Rational Models Via a Half-Size
            Test Matrix," in IEEE Transactions on Microwave Theory and Techniques, vol. 56, no. 12, pp. 2701-2708,
            Dec. 2008, DOI: 10.1109/TMTT.2008.2007319.
        """

        if parameter_type.lower() != 's':
            logging.error('Passivity testing is currently only supported for scattering (S) parameters.')
            return
        if parameter_type.lower() == 's' and len(np.flatnonzero(self.proportional_coeff)) > 0:
            logging.error('Passivity testing of scattering parameters with nonzero proportional coefficients does not '
                          'make any sense; you need to run vector_fit() with option \'fit_proportional=False\' first.')
            return

        # # the network needs to be reciprocal for this passivity test method to work: S = transpose(S)
        # if not np.allclose(self.zeros, np.transpose(self.zeros)) or \
        #         not np.allclose(self.constant_coeff, np.transpose(self.constant_coeff)) or \
        #         not np.allclose(self.proportional_coeff, np.transpose(self.proportional_coeff)):
        #     logging.error('Passivity testing with unsymmetrical model parameters is not supported. '
        #                   'The model needs to be reciprocal.')
        #     return

        # get state-space matrices
        A, B, C, D, E = self._get_ABCDE()
        n_ports = np.shape(D)[0]

        # build half-size test matrix P from state-space matrices A, B, C, D
        inv_neg = np.linalg.inv(D - np.identity(n_ports))
        inv_pos = np.linalg.inv(D + np.identity(n_ports))
        prod_neg = np.matmul(np.matmul(B, inv_neg), C)
        prod_pos = np.matmul(np.matmul(B, inv_pos), C)
        P = np.matmul(A - prod_neg, A - prod_pos)

        # extract eigenvalues of P
        P_eigs = np.linalg.eigvals(P)

        # purely imaginary square roots of eigenvalues identify frequencies (2*pi*f) of borders of passivity violations
        freqs_violation = []
        for sqrt_eigenval in np.sqrt(P_eigs):
            if np.real(sqrt_eigenval) == 0.0:
                freqs_violation.append(np.imag(sqrt_eigenval) / 2 / np.pi)

        # sort the output from lower to higher frequencies
        freqs_violation = np.sort(freqs_violation)

        # identify frequency bands of passivity violations

        # sweep the bands between crossover frequencies and identify bands of passivity violations
        violation_bands = []
        for i, freq in enumerate(freqs_violation):
            if i == 0:
                f_start = 0
                f_stop = freq
            else:
                f_start = freqs_violation[i - 1]
                f_stop = freq

            # calculate singular values at the center frequency between crossover frequencies to identify violations
            f_center = 0.5 * (f_start + f_stop)
            s_center = self._get_s_from_ABCDE(f_center, A, B, C, D, E)
            u, sigma, vh = np.linalg.svd(s_center)
            passive = True
            for singval in sigma:
                if singval > 1:
                    # passivity violation in this band
                    passive = False
            if not passive:
                # add this band to the list of passivity violations
                if violation_bands is None:
                    violation_bands = [[f_start, f_stop]]
                else:
                    violation_bands.append([f_start, f_stop])

        return np.array(violation_bands)

    def is_passive(self, parameter_type='S'):
        """
        Provides the passivity status of the model as a boolean value.

        Parameters
        ----------
        parameter_type : str, optional
            Representation type of the fitted frequency responses. Either *scattering* (:attr:`s` or :attr:`S`),
            *impedance* (:attr:`z` or :attr:`Z`) or *admittance* (:attr:`y` or :attr:`Y`). Currently, only scattering
            parameters are supported for passivity evaluation.

        Returns
        -------
        bool
            :attr:`True` is model is passive, else :attr:`False`.
        """
        viol_bands = self.passivity_test(parameter_type)
        if len(viol_bands) == 0:
            return True
        else:
            return False

    def passivity_enforce(self, n_samples=100, parameter_type='S'):
        """
        Enforces the passivity of the vector fitted model, if required. This is an implementation of the method
        presented in [1]_.

        Parameters
        ----------
        n_samples : int, optional
            Number of linearly spaced frequency samples at which passivity will be evaluated and enforce. (Default: 100)

        parameter_type : str, optional
            Representation type of the fitted frequency responses. Either *scattering* (:attr:`s` or :attr:`S`),
            *impedance* (:attr:`z` or :attr:`Z`) or *admittance* (:attr:`y` or :attr:`Y`). Currently, only scattering
            parameters are supported for passivity evaluation.

        Returns
        -------
        None

        See Also
        --------
        passivity_test : More verbose passivity evaluation routine.
        plot_passivation : Convergence plot for passivity enforcement iterations.

        References
        ----------
        .. [1] <NAME>, <NAME> and <NAME>, "Efficient Algorithm for Passivity Enforcement of S-Parameter-
            Based Macromodels," in IEEE Transactions on Microwave Theory and Techniques, vol. 57, no. 2, pp. 415-420,
            Feb. 2009, DOI: 10.1109/TMTT.2008.2011201.
        """

        if parameter_type.lower() != 's':
            logging.error('Passivity testing is currently only supported for scattering (S) parameters.')
            return
        if parameter_type.lower() == 's' and len(np.flatnonzero(self.proportional_coeff)) > 0:
            logging.error('Passivity testing of scattering parameters with nonzero proportional coefficients does not '
                          'make any sense; you need to run vector_fit() with option \'fit_proportional=False\' first.')
            return

        # always run passivity test first; this will write 'self.violation_bands'
        violation_bands = self.passivity_test()

        if len(violation_bands) == 0:
            # model is already passive; do nothing and return
            logging.info('Passivity enforcement: The model is already passive. Nothing to do.')
            return

        freqs_eval = np.linspace(0, 1.2 * violation_bands[0, -1], n_samples)
        A, B, C, D, E = self._get_ABCDE()
        dim_A = np.shape(A)[0]
        C_t = C
        delta = 0.999   # predefined tolerance parameter (users should not need to change this)

        # iterative compensation of passivity violations
        t = 0
        self.history_max_sigma = []
        while t < self.max_iterations:
            logging.info('Passivity enforcement; Iteration {}'.format(t + 1))

            A_matrix = []
            b_vector = []
            sigma_max = 0

            # sweep through evaluation frequencies
            for i_eval, freq_eval in enumerate(freqs_eval):<commit_msg>Fixed frequency interval for passivity enforcement.<commit_after>"""
=========================================
VectorFitting (:mod:`skrf.vectorFitting`)
=========================================

.. autoclass:: VectorFitting
    :members:
    :member-order: groupwise
    :special-members: __init__

"""

import numpy as np
import os
from functools import wraps
try:
    from . import plotting    # will perform the correct setup for matplotlib before it is called below
    import matplotlib.pyplot as mplt
    from matplotlib.ticker import EngFormatter
except ImportError:
    mplt = None

import logging


def check_plotting(func):
    """
    This decorator checks if matplotlib.pyplot is available under the name mplt.
    If not, raise an RuntimeError.

    Raises
    ------
    RuntimeError
        When trying to run the decorated function without matplotlib
    """

    @wraps(func)
    def wrapper(*args, **kwargs):
        if mplt is None:
            raise RuntimeError('Plotting is not available')
        func(*args, **kwargs)

    return wrapper


class VectorFitting:
    """
    =========================================================
    VectorFitting (:class:`skrf.vectorFitting.VectorFitting`)
    =========================================================

    This class provides a Python implementation of the Vector Fitting algorithm and various functions for the fit
    analysis and export of SPICE equivalent circuits.

    Notes
    -----
    The fitting code is based on the original algorithm [1]_ and on two improvements for relaxed pole relocation [2]_
    and efficient (fast) solving [3]_. See also the Vector Fitting website [4]_ for further information and download of
    the papers listed below. A Matlab implementation is also available there for reference.

    References
    ----------
    .. [1] <NAME>, <NAME>, "Rational Approximation of Frequency Domain Responses by Vector Fitting", IEEE
        Transactions on Power Delivery, vol. 14, no. 3, pp. 1052-1061, July 1999, DOI: https://doi.org/10.1109/61.772353

    .. [2] <NAME>, "Improving the Pole Relocating Properties of Vector Fitting", IEEE Transactions on Power
        Delivery, vol. 21, no. 3, pp. 1587-1592, July 2006, DOI: https://doi.org/10.1109/TPWRD.2005.860281

    .. [3] <NAME>, <NAME>, <NAME>, <NAME>, "Marcomodeling of Multiport Systems Using a Fast
        Implementation of the Vector Fitting Method", IEEE Microwave and Wireless Components Letters, vol. 18, no. 6,
        pp. 383-385, June 2008, DOI: https://doi.org/10.1109/LMWC.2008.922585

    .. [4] Vector Fitting website: https://www.sintef.no/projectweb/vectorfitting/
    """

    def __init__(self, network):
        """
        Creates a VectorFitting instance based on a supplied :class:`skrf.network.Network` containing the frequency
        responses of the N-port.

        Parameters
        ----------
        network : :class:`skrf.network.Network`
            Network instance of the N-port holding the S-matrix to be fitted.
        """

        self.network = network
        self.initial_poles = None

        self.poles = None
        """ Instance variable holding the list of fitted poles. Will be initialized by :func:`vector_fit`. """

        self.zeros = None
        """ Instance variable holding the list of fitted zeros. Will be initialized by :func:`vector_fit`. """

        self.proportional_coeff = None
        """ Instance variable holding the list of fitted proportional coefficients. Will be initialized by 
        :func:`vector_fit`. """

        self.constant_coeff = None
        """ Instance variable holding the list of fitted constants. Will be initialized by :func:`vector_fit`. """

        self.max_iterations = 100
        """ Instance variable specifying the maximum number of iterations for the fitting process. To be changed by the
         user before calling :func:`vector_fit`. """

        self.max_tol = 1e-6
        """ Instance variable specifying the convergence criterion in terms of relative tolerance. To be changed by the
         user before calling :func:`vector_fit`. """

        self.d_res_history = []
        self.delta_max_history = []
        self.history_max_sigma = []

    def vector_fit(self, n_poles_real=2, n_poles_cmplx=2, init_pole_spacing='lin', parameter_type='S',
                   fit_constant=True, fit_proportional=False):
        """
        Main work routine performing the vector fit. The results will be stored in the class variables
        :attr:`poles`, :attr:`zeros`, :attr:`proportional_coeff` and :attr:`constant_coeff`.

        Parameters
        ----------
        n_poles_real : int, optional
            Number of initial real poles. See notes.

        n_poles_cmplx : int, optional
            Number of initial complex conjugate poles. See notes.

        init_pole_spacing : str, optional
            Type of initial pole spacing across the frequency interval of the S-matrix. Either linear (lin) or
            logarithmic (log).

        parameter_type : str, optional
            Representation type of the frequency responses to be fitted. Either *scattering* (:attr:`s` or :attr:`S`),
            *impedance* (:attr:`z` or :attr:`Z`) or *admittance* (:attr:`y` or :attr:`Y`). As scikit-rf can currently
            only read S parameters from a Touchstone file, the fit should also be performed on the original S
            parameters. Otherwise, scikit-rf will convert the responses from S to Z or Y, which might work for the fit
            but can cause other issues.

        fit_constant : bool, optional
            Include a constant term **d** in the fit.

        fit_proportional : bool, optional
            Include a proportional term **e** in the fit.

        Returns
        -------
        None

        Notes
        -----
        The required number of real or complex conjugate starting poles depends on the behaviour of the frequency
        responses. To fit a smooth response such as a low-pass characteristic, 1-3 real poles and no complex conjugate
        poles is usually sufficient. If resonances or other types of peaks are present in some or all of the responses,
        a similar number of complex conjugate poles is required. Be careful not to use too many poles, as excessive
        poles will not only increase the computation workload during the fitting and the subsequent use of the model,
        but they can also introduce unwanted resonances at frequencies well outside the fit interval.
        """

        # create initial poles and space them across the frequencies in the provided Touchstone file
        # use normalized frequencies during the iterations (seems to be more stable during least-squares fit)
        norm = np.average(self.network.f)
        freqs_norm = self.network.f / norm

        fmin = np.amin(freqs_norm)
        fmax = np.amax(freqs_norm)
        weight_regular = 1.0
        if init_pole_spacing == 'log':
            pole_freqs_real = np.geomspace(fmin, fmax, n_poles_real)
            pole_freqs_cmplx = np.geomspace(fmin, fmax, n_poles_cmplx)
        elif init_pole_spacing == 'lin':
            pole_freqs_real = np.linspace(fmin, fmax, n_poles_real)
            pole_freqs_cmplx = np.linspace(fmin, fmax, n_poles_cmplx)
        else:
            logging.warning('Invalid choice of initial pole spacing; proceeding with linear spacing')
            pole_freqs_real = np.linspace(fmin, fmax, n_poles_real)
            pole_freqs_cmplx = np.linspace(fmin, fmax, n_poles_cmplx)
        poles = []

        # add real poles
        for i, f in enumerate(pole_freqs_real):
            omega = 2 * np.pi * f
            poles.append((-1 / 100 + 0j) * omega)
        # add complex-conjugate poles (store only positive imaginary parts)
        for i, f in enumerate(pole_freqs_cmplx):
            omega = 2 * np.pi * f
            poles.append((-1 / 100 + 1j) * omega)
        poles = np.array(poles)

        # save initial poles (un-normalize first)
        self.initial_poles = poles * norm
        max_singular = 1

        logging.info('### Starting pole relocation process.\n')

        # stack frequency responses as a single vector
        # stacking order:
        # s11, s12, s13, ..., s21, s22, s23, ...
        freq_responses = []
        for i in range(self.network.nports):
            for j in range(self.network.nports):
                if parameter_type.lower() == 's':
                    freq_responses.append(self.network.s[:, i, j])
                elif parameter_type.lower() == 'z':
                    freq_responses.append(self.network.z[:, i, j])
                elif parameter_type.lower() == 'y':
                    freq_responses.append(self.network.y[:, i, j])
                else:
                    logging.warning('Invalid choice of matrix parameter type (S, Z, or Y); proceeding with S representation.')
                    freq_responses.append(self.network.s[:, i, j])
        freq_responses = np.array(freq_responses)

        # ITERATIVE FITTING OF POLES to the provided frequency responses
        # initial set of poles will be replaced with new poles after every iteration
        iterations = self.max_iterations
        self.d_res_history = []
        self.delta_max_history = []
        converged = False
        while iterations > 0:
            logging.info('Iteration {}'.format(self.max_iterations - iterations + 1))

            # generate coefficients of approximation function for each target frequency response
            # responses will be treated independently using QR decomposition
            # simplified coeff matrices of all responses will be stacked for least-squares solver
            A_matrix = []
            b_vector = []

            for freq_response in freq_responses:
                # calculate coefficients (row A_k in matrix) for each frequency sample s_k of the target response
                # row will be appended to submatrix A_sub of complete coeff matrix A_matrix
                # 2 rows per pole in result vector (1st for real part, 2nd for imaginary part)
                # --> 2 columns per pole in coeff matrix
                A_sub = []
                b_sub = []
                n_unused = 0

                for k, f_sample in enumerate(freqs_norm):
                    s_k = 2j * np.pi * f_sample
                    A_k = []
                    n_unused = 0

                    # add coefficients for a pair of complex conjugate poles
                    # part 1: first sum of rational functions (residue variable c)
                    for pole in poles:
                        # separate and stack real and imaginary part to preserve conjugacy of the pole pair
                        if np.imag(pole) == 0.0:
                            A_k.append(1 / (s_k - pole))
                            n_unused += 1
                        else:
                            # complex pole of a conjugated pair
                            A_k.append(1 / (s_k - pole) + 1 / (s_k - np.conjugate(pole)))       # real part of residue
                            A_k.append(1j / (s_k - pole) - 1j / (s_k - np.conjugate(pole)))     # imaginary part of residue
                            n_unused += 2

                    # part 2: constant (variable d) and proportional term (variable e)
                    if fit_constant:
                        A_k.append(1)
                        n_unused += 1
                    if fit_proportional:
                        A_k.append(s_k)
                        n_unused += 1

                    # part 3: second sum of rational functions (variable c_res)
                    for pole in poles:
                        # separate and stack real and imaginary part to preserve conjugacy of the pole pair
                        if np.imag(pole) == 0.0:
                            A_k.append(-1 * freq_response[k] / (s_k - pole))
                        else:
                            # complex pole of a conjugated pair
                            A_k.append(-1 * freq_response[k] / (s_k - pole)
                                       - freq_response[k] / (s_k - np.conjugate(pole)))
                            A_k.append(-1j * freq_response[k] / (s_k - pole)
                                       + 1j * freq_response[k] / (s_k - np.conjugate(pole)))

                    # part 4: constant (variable d_res)
                    A_k.append(-1 * freq_response[k])

                    A_sub.append(np.sqrt(weight_regular) * np.real(A_k))
                    A_sub.append(np.sqrt(weight_regular) * np.imag(A_k))
                    b_sub.append(np.sqrt(weight_regular) * 0.0)
                    b_sub.append(np.sqrt(weight_regular) * 0.0)

                # QR decomposition
                Q, R = np.linalg.qr(A_sub, 'reduced')

                # only R22 is required to solve for c_res and d_res
                R22 = R[n_unused:, n_unused:]
                # similarly, only right half of Q is required
                Q2 = Q[:, n_unused:]

                if len(A_matrix) == 0:
                    A_matrix = R22
                    b_vector = np.matmul(np.transpose(Q2), b_sub)
                else:
                    A_matrix = np.vstack((A_matrix, R22))
                    b_vector = np.append(b_vector, np.matmul(np.transpose(Q2), b_sub))

                # add extra equation to avoid trivial solution
                # use weight=1 for all equations, except for this extra equation
                weight_extra = np.linalg.norm(weight_regular * freq_response) / len(freq_response)
                A_k = np.zeros(np.shape(A_matrix)[1])
                for k, f_sample in enumerate(freqs_norm):
                    s_k = 2j * np.pi * f_sample
                    i = 0

                    # part 3: second sum of rational functions (variable c_res)
                    for pole in poles:
                        if np.imag(pole) == 0.0:
                            # real pole
                            A_k[i] += np.real(1 / (s_k - pole))
                            i += 1
                        else:
                            # complex pole of a conjugated pair
                            A_k[i] += np.real(1 / (s_k - pole) + 1 / (s_k - np.conjugate(pole)))
                            A_k[i+1] += np.real(1j / (s_k - pole) - 1j / (s_k - np.conjugate(pole)))
                            i += 2
                    # part 4: constant (d_res)
                    A_k[i] += 1

                A_matrix = np.vstack((A_matrix, np.sqrt(weight_extra) * A_k))
                b_vector = np.append(b_vector, np.sqrt(weight_extra) * len(freqs_norm))

            logging.info('A_matrix: condition number = {}'.format(np.linalg.cond(A_matrix)))

            # solve least squares for real parts
            x, residuals, rank, singular_vals = np.linalg.lstsq(A_matrix, b_vector, rcond=None)

            # assemble individual result vectors from single LS result x
            c_res = x[:-1]
            d_res = x[-1]

            # check if d_res is suited for zeros calculation
            tol_res = 1e-8
            if np.abs(d_res) < tol_res:
                # d_res is too small, discard solution and proceed the |d_res| = tol_res
                d_res = tol_res * (d_res / np.abs(d_res))
                logging.warning('Replacing d_res solution as it was too small')

            self.d_res_history.append(d_res)
            logging.info('d_res = {}'.format(d_res))

            # build test matrix H, which will hold the new poles as Eigenvalues
            A_matrix = np.zeros((len(c_res), len(c_res)))
            i = 0
            for pole in poles:
                # fill diagonal with previous poles
                if np.imag(pole) == 0.0:
                    # one row for a real pole
                    A_matrix[i, i] = np.real(pole)
                    A_matrix[i] -= c_res / d_res
                    i += 1
                else:
                    # two rows for a complex pole of a conjugated pair
                    A_matrix[i, i] = np.real(pole)
                    A_matrix[i, i + 1] = np.imag(pole)
                    A_matrix[i + 1, i] = -1 * np.imag(pole)
                    A_matrix[i + 1, i + 1] = np.real(pole)
                    A_matrix[i] -= 2 * c_res / d_res
                    i += 2
            poles_new = np.linalg.eigvals(A_matrix)

            # replace poles for next iteration
            poles = []
            for k, pole in enumerate(poles_new):
                # flip real part of unstable poles (real part needs to be negative for stability)
                if np.real(pole) > 0.0:
                    pole = -1 * np.real(pole) + 1j * np.imag(pole)
                if np.imag(pole) >= 0.0:
                    # complex poles need to come in complex conjugate pairs; append only the positive part
                    poles.append(pole)

            poles = np.sort_complex(poles)

            # calculate relative changes in the singular values; stop iteration loop once poles have converged
            new_max_singular = np.amax(singular_vals)
            delta_max = np.abs(1 - new_max_singular / max_singular)
            self.delta_max_history.append(delta_max)
            logging.info('Max. relative change in residues = {}\n'.format(delta_max))
            max_singular = new_max_singular

            stop = False
            if delta_max < self.max_tol:
                if converged:
                    # is really converged, finish
                    logging.info('Pole relocation process converged after {} iterations.'.format(
                        self.max_iterations - iterations + 1))
                    stop = True
                else:
                    # might be converged, but do one last run to be sure
                    converged = True
            else:
                if converged:
                    # is not really converged, continue
                    converged = False

            iterations -= 1

            if iterations == 0:
                if converged and stop is False:
                    logging.warning('Reached tolerance only after max. number of iterations (N_max = {}). '
                                    'Results might not have converged properly.'.format(self.max_iterations))
                else:
                    logging.warning('Reached maximum number of iterations (N_max = {}). '
                                    'Results did not converge.'.format(self.max_iterations))

            if stop:
                iterations = 0

        # ITERATIONS DONE
        poles = np.array(poles)

        logging.info('Initial poles before relocation:')
        logging.info(self.initial_poles)

        logging.info('Final poles:')
        logging.info(poles * norm)

        logging.info('\n### Starting zeros calculation process.\n')

        # finally, solve for the residues with the previously calculated poles
        zeros = []
        constant_coeff = []
        proportional_coeff = []

        for freq_response in freq_responses:
            # calculate coefficients (row A_k in matrix) for each frequency sample s_k of the target response
            # row will be appended to submatrix A_sub of complete coeff matrix A_matrix
            # 2 rows per pole in result vector (1st for real part, 2nd for imaginary part)
            # --> 2 columns per pole in coeff matrix
            A_matrix = []
            b_vector = []

            for k, f_sample in enumerate(freqs_norm):
                s_k = 2j * np.pi * f_sample
                A_k = []
                # add coefficients for a pair of complex conjugate poles
                # part 1: first sum of rational functions (residue variable c)
                for pole in poles:
                    # separate and stack real and imaginary part to preserve conjugacy of the pole pair
                    if np.imag(pole) == 0.0:
                        A_k.append(1 / (s_k - pole))
                    else:
                        A_k.append(1 / (s_k - pole) + 1 / (s_k - np.conjugate(pole)))       # real part of residue
                        A_k.append(1j / (s_k - pole) - 1j / (s_k - np.conjugate(pole)))     # imaginary part of residue

                # part 2: constant (variable d) and proportional term (variable e)
                if fit_constant:
                    A_k.append(1.0)
                if fit_proportional:
                    A_k.append(s_k)

                A_matrix.append(np.array(A_k))
                b_vector.append(np.array(freq_response[k]))

            A_matrix = np.vstack((np.real(A_matrix), np.imag(A_matrix)))
            b_vector = np.append(np.real(b_vector), np.imag(b_vector))

            logging.info('A_matrix: condition number = {}'.format(np.linalg.cond(A_matrix)))

            # solve least squares and obtain results as stack of real part vector and imaginary part vector
            x, residuals, rank, singular_vals = np.linalg.lstsq(A_matrix, b_vector, rcond=None)

            i = 0
            zeros_response = []
            for pole in poles:
                if np.imag(pole) == 0.0:
                    zeros_response.append(x[i] + 0j)
                    i += 1
                else:
                    zeros_response.append(x[i] + 1j * x[i + 1])
                    i += 2
            zeros.append(zeros_response)

            if fit_constant and fit_proportional:
                # both constant d and proportional e were fitted
                constant_coeff.append(x[-2])
                proportional_coeff.append(x[-1])
            elif fit_constant:
                # only constant d was fitted
                constant_coeff.append(x[-1])
                proportional_coeff.append(0.0)
            elif fit_proportional:
                # only proportional e was fitted
                constant_coeff.append(0.0)
                proportional_coeff.append(x[-1])
            else:
                # neither constant d nor proportional e was fitted
                constant_coeff.append(0.0)
                proportional_coeff.append(0.0)

        # save poles, zeros, d, e in actual frequencies (un-normalized)
        self.poles = poles * norm
        self.zeros = np.array(zeros) * norm
        self.constant_coeff = np.array(constant_coeff)
        self.proportional_coeff = np.array(proportional_coeff) / norm

        logging.info('\n### Vector fitting finished.\n')

    def _get_ABCDE(self):
        """
        Private method.
        Returns the real-valued system matrices of the state-space representation of the current rational model, as
        defined in [1]_.

        Returns
        -------
        A : ndarray
            State-space matrix A holding the poles on the diagonal as real values with imaginary parts on the sub-
            diagonal
        B : ndarray
            State-space matrix B holding coefficients (1, 2, or 0), depending on the respective type of pole in A
        C : ndarray
            State-space matrix C holding the residues (zeros)
        D : ndarray
            State-space matrix D holding the constants
        E : ndarray
            State-space matrix E holding the proportional coefficients (usually 0 in case of fitted S-parameters)

        References
        ----------
        .. [1] <NAME> and <NAME>, "Fast Passivity Assessment for S-Parameter Rational Models Via a Half-Size
            Test Matrix," in IEEE Transactions on Microwave Theory and Techniques, vol. 56, no. 12, pp. 2701-2708,
            Dec. 2008, DOI: 10.1109/TMTT.2008.2007319.
        """

        # initial checks
        if self.poles is None:
            logging.error('self.poles = None; nothing to do. You need to run vector_fit() first.')
            return
        if self.zeros is None:
            logging.error('self.zeros = None; nothing to do. You need to run vector_fit() first.')
            return
        if self.proportional_coeff is None:
            logging.error('self.proportional_coeff = None; nothing to do. You need to run vector_fit() first.')
            return
        if self.constant_coeff is None:
            logging.error('self.constant_coeff = None; nothing to do. You need to run vector_fit() first.')
            return

        # assemble real-valued state-space matrices A, B, C, D, E from fitted complex-valued pole-residue model

        # determine size of the matrix system
        n_ports = int(np.sqrt(len(self.constant_coeff)))
        n_poles_real = 0
        n_poles_cplx = 0
        for pole in self.poles:
            if np.imag(pole) == 0.0:
                n_poles_real += 1
            else:
                n_poles_cplx += 1
        n_matrix = (n_poles_real + 2 * n_poles_cplx) * n_ports

        # state-space matrix A holds the poles on the diagonal as real values with imaginary parts on the sub-diagonal
        # state-space matrix B holds coefficients (1, 2, or 0), depending on the respective type of pole in A
        # assemble A = [[poles_real,   0,                  0],
        #               [0,            real(poles_cplx),   imag(poles_cplx],
        #               [0,            -imag(poles_cplx),  real(poles_cplx]]
        A = np.identity(n_matrix)
        B = np.zeros(shape=(n_matrix, n_ports))
        i_A = 0  # index on diagonal of A
        for j in range(n_ports):
            for pole in self.poles:
                if np.imag(pole) == 0.0:
                    # adding a real pole
                    A[i_A, i_A] = np.real(pole)
                    B[i_A, j] = 1
                    i_A += 1
                else:
                    # adding a complex-conjugate pole
                    A[i_A, i_A] = np.real(pole)
                    A[i_A, i_A + 1] = np.imag(pole)
                    A[i_A + 1, i_A] = -1 * np.imag(pole)
                    A[i_A + 1, i_A + 1] = np.real(pole)
                    B[i_A, j] = 2
                    i_A += 2

        # state-space matrix C holds the residues (zeros)
        # assemble C = [[R1.11, R1.12, R1.13, ...], [R2.11, R2.12, R2.13, ...], ...]
        C = np.zeros(shape=(n_ports, n_matrix))
        for i in range(n_ports):
            for j in range(n_ports):
                # i: row index
                # j: column index
                i_response = i * n_ports + j

                j_zeros = 0
                for zero in self.zeros[i_response]:
                    if np.imag(zero) == 0.0:
                        C[i, j * (n_poles_real + 2 * n_poles_cplx) + j_zeros] = np.real(zero)
                        j_zeros += 1
                    else:
                        C[i, j * (n_poles_real + 2 * n_poles_cplx) + j_zeros] = np.real(zero)
                        C[i, j * (n_poles_real + 2 * n_poles_cplx) + j_zeros + 1] = np.imag(zero)
                        j_zeros += 2

        # state-space matrix D holds the constants
        # assemble D = [[d11, d12, ...], [d21, d22, ...], ...]
        D = np.zeros(shape=(n_ports, n_ports))
        for i in range(n_ports):
            for j in range(n_ports):
                # i: row index
                # j: column index
                i_response = i * n_ports + j
                D[i, j] = self.constant_coeff[i_response]

        # state-space matrix E holds the proportional coefficients (usually 0 in case of fitted S-parameters)
        # assemble E = [[e11, e12, ...], [e21, e22, ...], ...]
        E = np.zeros(shape=(n_ports, n_ports))
        for i in range(n_ports):
            for j in range(n_ports):
                # i: row index
                # j: column index
                i_response = i * n_ports + j
                E[i, j] = self.proportional_coeff[i_response]

        return A, B, C, D, E

    @staticmethod
    def _get_s_from_ABCDE(freq, A, B, C, D, E):
        """
        Private method.
        Returns the S-matrix of the vector fitted model calculated from the real-valued system matrices of the state-
        space representation, as provided by `_get_ABCDE()`.

        Parameters
        ----------
        freq : float
            Frequency (in Hz) at which to calculate the S-matrix.
        A : ndarray
        B : ndarray
        C : ndarray
        D : ndarray
        E : ndarray

        Returns
        -------
        ndarray
            Complex-valued S-matrix (NxN) calculated at frequency `freq`.
        """

        dim_A = np.shape(A)[0]
        stsp_poles = np.linalg.inv(2j * np.pi * freq * np.identity(dim_A) - A)
        stsp_S = np.matmul(np.matmul(C, stsp_poles), B)
        stsp_S += D + 2j * np.pi * freq * E
        return stsp_S

    def passivity_test(self, parameter_type='S'):
        """
        Evaluates the passivity of reciprocal vector fitted models by means of a half-size test matrix [1]_. Any
        existing frequency bands of passivity violations will be returned as a sorted list.

        Parameters
        ----------
        parameter_type: str, optional
            Representation type of the fitted frequency responses. Either *scattering* (:attr:`s` or :attr:`S`),
            *impedance* (:attr:`z` or :attr:`Z`) or *admittance* (:attr:`y` or :attr:`Y`). Currently, only scattering
            parameters are supported for passivity evaluation.

        Returns
        -------
        ndarray
            NumPy array with frequency bands of passivity violation:
            [[f_start_1, f_stop_1], [f_start_2, f_stop_2], ...].

        See Also
        --------
        is_passive : Query the model passivity as a boolean value.

        References
        ----------
        .. [1] <NAME> and <NAME>, "Fast Passivity Assessment for S-Parameter Rational Models Via a Half-Size
            Test Matrix," in IEEE Transactions on Microwave Theory and Techniques, vol. 56, no. 12, pp. 2701-2708,
            Dec. 2008, DOI: 10.1109/TMTT.2008.2007319.
        """

        if parameter_type.lower() != 's':
            logging.error('Passivity testing is currently only supported for scattering (S) parameters.')
            return
        if parameter_type.lower() == 's' and len(np.flatnonzero(self.proportional_coeff)) > 0:
            logging.error('Passivity testing of scattering parameters with nonzero proportional coefficients does not '
                          'make any sense; you need to run vector_fit() with option \'fit_proportional=False\' first.')
            return

        # # the network needs to be reciprocal for this passivity test method to work: S = transpose(S)
        # if not np.allclose(self.zeros, np.transpose(self.zeros)) or \
        #         not np.allclose(self.constant_coeff, np.transpose(self.constant_coeff)) or \
        #         not np.allclose(self.proportional_coeff, np.transpose(self.proportional_coeff)):
        #     logging.error('Passivity testing with unsymmetrical model parameters is not supported. '
        #                   'The model needs to be reciprocal.')
        #     return

        # get state-space matrices
        A, B, C, D, E = self._get_ABCDE()
        n_ports = np.shape(D)[0]

        # build half-size test matrix P from state-space matrices A, B, C, D
        inv_neg = np.linalg.inv(D - np.identity(n_ports))
        inv_pos = np.linalg.inv(D + np.identity(n_ports))
        prod_neg = np.matmul(np.matmul(B, inv_neg), C)
        prod_pos = np.matmul(np.matmul(B, inv_pos), C)
        P = np.matmul(A - prod_neg, A - prod_pos)

        # extract eigenvalues of P
        P_eigs = np.linalg.eigvals(P)

        # purely imaginary square roots of eigenvalues identify frequencies (2*pi*f) of borders of passivity violations
        freqs_violation = []
        for sqrt_eigenval in np.sqrt(P_eigs):
            if np.real(sqrt_eigenval) == 0.0:
                freqs_violation.append(np.imag(sqrt_eigenval) / 2 / np.pi)

        # sort the output from lower to higher frequencies
        freqs_violation = np.sort(freqs_violation)

        # identify frequency bands of passivity violations

        # sweep the bands between crossover frequencies and identify bands of passivity violations
        violation_bands = []
        for i, freq in enumerate(freqs_violation):
            if i == 0:
                f_start = 0
                f_stop = freq
            else:
                f_start = freqs_violation[i - 1]
                f_stop = freq

            # calculate singular values at the center frequency between crossover frequencies to identify violations
            f_center = 0.5 * (f_start + f_stop)
            s_center = self._get_s_from_ABCDE(f_center, A, B, C, D, E)
            u, sigma, vh = np.linalg.svd(s_center)
            passive = True
            for singval in sigma:
                if singval > 1:
                    # passivity violation in this band
                    passive = False
            if not passive:
                # add this band to the list of passivity violations
                if violation_bands is None:
                    violation_bands = [[f_start, f_stop]]
                else:
                    violation_bands.append([f_start, f_stop])

        return np.array(violation_bands)

    def is_passive(self, parameter_type='S'):
        """
        Provides the passivity status of the model as a boolean value.

        Parameters
        ----------
        parameter_type : str, optional
            Representation type of the fitted frequency responses. Either *scattering* (:attr:`s` or :attr:`S`),
            *impedance* (:attr:`z` or :attr:`Z`) or *admittance* (:attr:`y` or :attr:`Y`). Currently, only scattering
            parameters are supported for passivity evaluation.

        Returns
        -------
        bool
            :attr:`True` is model is passive, else :attr:`False`.
        """
        viol_bands = self.passivity_test(parameter_type)
        if len(viol_bands) == 0:
            return True
        else:
            return False

    def passivity_enforce(self, n_samples=100, parameter_type='S'):
        """
        Enforces the passivity of the vector fitted model, if required. This is an implementation of the method
        presented in [1]_.

        Parameters
        ----------
        n_samples : int, optional
            Number of linearly spaced frequency samples at which passivity will be evaluated and enforce. (Default: 100)

        parameter_type : str, optional
            Representation type of the fitted frequency responses. Either *scattering* (:attr:`s` or :attr:`S`),
            *impedance* (:attr:`z` or :attr:`Z`) or *admittance* (:attr:`y` or :attr:`Y`). Currently, only scattering
            parameters are supported for passivity evaluation.

        Returns
        -------
        None

        See Also
        --------
        passivity_test : More verbose passivity evaluation routine.
        plot_passivation : Convergence plot for passivity enforcement iterations.

        References
        ----------
        .. [1] <NAME>, <NAME> and <NAME>, "Efficient Algorithm for Passivity Enforcement of S-Parameter-
            Based Macromodels," in IEEE Transactions on Microwave Theory and Techniques, vol. 57, no. 2, pp. 415-420,
            Feb. 2009, DOI: 10.1109/TMTT.2008.2011201.
        """

        if parameter_type.lower() != 's':
            logging.error('Passivity testing is currently only supported for scattering (S) parameters.')
            return
        if parameter_type.lower() == 's' and len(np.flatnonzero(self.proportional_coeff)) > 0:
            logging.error('Passivity testing of scattering parameters with nonzero proportional coefficients does not '
                          'make any sense; you need to run vector_fit() with option \'fit_proportional=False\' first.')
            return

        # always run passivity test first; this will write 'self.violation_bands'
        violation_bands = self.passivity_test()

        if len(violation_bands) == 0:
            # model is already passive; do nothing and return
            logging.info('Passivity enforcement: The model is already passive. Nothing to do.')
            return

        freqs_eval = np.linspace(0, 1.2 * violation_bands[-1, -1], n_samples)
        A, B, C, D, E = self._get_ABCDE()
        dim_A = np.shape(A)[0]
        C_t = C
        delta = 0.999   # predefined tolerance parameter (users should not need to change this)

        # iterative compensation of passivity violations
        t = 0
        self.history_max_sigma = []
        while t < self.max_iterations:
            logging.info('Passivity enforcement; Iteration {}'.format(t + 1))

            A_matrix = []
            b_vector = []
            sigma_max = 0

            # sweep through evaluation frequencies
            for i_eval, freq_eval in enumerate(freqs_eval):
{'fraction_non_alphanumeric': 0.09493179247468622, 'fraction_numerical': 0.013000994835034547, 'mean_word_length': 1.7111505209487918, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 47, 'https://': 8, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1641071', 'n_tokens_mistral': 21669, 'n_tokens_neox': 19770, 'n_words': 6536}