<commit_before>"""Model order reduction for linear state space systems can be done with
Proper Orthogonal Decomposition (POD) methods. Some of them can be simply
applied when creating a system.

"""

from __future__ import division

import collections
import numpy as np
import scipy as sp
from numpy import array
from scipy.integrate import ode
# from scipy.signal.ltisys import abcd_normalize
# currently buggy. Will be fixed in scipy v0.15
# instead the following can be used temporarily
from futurescipy import abcd_normalize


def truncation_square_root(A, B, C,
                           k=0, tol=0.0,
                           balance=True, scale=True,
                           check_stability=True,
                           length_cache_array=None):
    """Perform truncation of a system. Scaling and balancing are optional

    This allows to reduce a linear state space system by either specifying it
    to a certain number of states `k` or by specifying a error tolerance `tol`
    relative to the input. In theory the most accurate results are achieved by
    using `balance` and `scale` but the size of the error strongly depends on
    the particular problem and scaling and balancing may in some cases cost
    too much.

    Parameters
    ----------
    A : array_like
    B : array_like
    C : array_like
        State-Space matrices of the system that should be reduced
    k : int, optional
        Order of the output system
    tol : float, optional
        Error of the output system, based on the Hankel Singular Values
    balance : Boolean, optional
        Balance the system before reducing it to make sure, that the error
        is kept small
    scale : Boolean, optional
        Scale the system
    check_stability : Boolean, optional
        Checks if all the real parts of the eigenvalues of A are in the left
        half of the complex plane.

    Returns
    -------
    Nr : int
        Actual size of the system, based on the error: If the Machine error
        would have been bigger than the error of the reduced system, it may
        happen that ``Nr < k`` and if the error would be inconsiderably bad,
        It might be the case that ``Nr > k``. In case `k` was never specified,
        this is purely based on `tol`
    Ar, Br, Cr : ndarray
        Reduced arrays
    hsv : ndarray
        Hankel singular values of the original system. The size of the error
        may be calculated based on this.

    Raises
    ------
    ValueError
        If the system that's provided is not stable (i.e. `A` has eigenvalues
        which have non-negative real parts)
    ImportError
        If the slycot subroutine `ab09ad` can't be found. Occurs if the
        slycot package is not installed.

    """
    if check_stability and not isStable(A):
        raise ValueError("This doesn't seem to be a stable system!")
    try:
        from slycot import ab09ad
    except ImportError:
        raise ImportError("can't find slycot subroutine ab09ad")

    if balance:
        job = 'B'
    else:
        job = 'N'

    if scale:
        equil = 'S'
    else:
        equil = 'N'

    dico = 'C'
    n = np.size(A,0)
    m = np.size(B,1)
    p = np.size(C,0)
    nr = k

    if not length_cache_array:
        length_cache_array = 5*n**2 + n*m + n*p + 10*n + m*p

    return ab09ad(dico,job,equil,n,m,p,A,B,C,nr,tol,ldwork=length_cache_array)

def controllabilityTruncation(A,B,C,k=None,check_stability=True):
    """Truncate the system based on the controllability Gramian

    Solves the Lyapunov Equation for ``AP + PA^H + B B^H`` and computes the
    eigenvalues and eigenvectors of `P` which are used to truncate the system.
    """
    if check_stability and not isStable(A):
        raise ValueError("This doesn't seem to be a stable system!")

    if not k:
        k = int(a.shape[0]/2)

    P = sp.linalg.solve_lyapunov(A, -np.dot(B, B.H))

    Lambdak, Uk = sp.linalg.eigh(P, eigvals=(0,k), check_finite=False,
                                 overwrite_a=True, overwrite_b=True)

    UkH = Uk.H

    A = np.dot(UkH, np.dot(A, Uk))
    B = np.dot(UkH, B)
    C = np.dot(C, Uk)

    return k, A, B, C, Lambdak

def isStable(A):
    """Check if all eigenvalues are in the left half of the complex plane"""
    D, V = np.linalg.eig(A)
    return (D.real < 0).all()

class lss(object):
    """linear time independent state space system.

    Default contstructor is called like ``lss(A,B,C,D)`` and a system can
    easily be copied by calling ``lss(sys)`` where `sys` is a<commit_msg>added docstring and fixed eigval order of gramian model reduction<commit_after>"""Model order reduction for linear state space systems can be done with
Proper Orthogonal Decomposition (POD) methods. Some of them can be simply
applied when creating a system.

"""

from __future__ import division

import collections
import numpy as np
import scipy as sp
from numpy import array
from scipy.integrate import ode
# from scipy.signal.ltisys import abcd_normalize
# currently buggy. Will be fixed in scipy v0.15
# instead the following can be used temporarily
from futurescipy import abcd_normalize


def truncation_square_root(A, B, C,
                           k=0, tol=0.0,
                           balance=True, scale=True,
                           check_stability=True,
                           length_cache_array=None):
    """Perform truncation of a system. Scaling and balancing are optional

    This allows to reduce a linear state space system by either specifying it
    to a certain number of states `k` or by specifying a error tolerance `tol`
    relative to the input. In theory the most accurate results are achieved by
    using `balance` and `scale` but the size of the error strongly depends on
    the particular problem and scaling and balancing may in some cases cost
    too much.

    Parameters
    ----------
    A : array_like
    B : array_like
    C : array_like
        State-Space matrices of the system that should be reduced
    k : int, optional
        Order of the output system
    tol : float, optional
        Error of the output system, based on the Hankel Singular Values
    balance : Boolean, optional
        Balance the system before reducing it to make sure, that the error
        is kept small
    scale : Boolean, optional
        Scale the system
    check_stability : Boolean, optional
        Checks if all the real parts of the eigenvalues of A are in the left
        half of the complex plane.

    Returns
    -------
    Nr : int
        Actual size of the system, based on the error: If the Machine error
        would have been bigger than the error of the reduced system, it may
        happen that ``Nr < k`` and if the error would be inconsiderably bad,
        It might be the case that ``Nr > k``. In case `k` was never specified,
        this is purely based on `tol`
    Ar, Br, Cr : ndarray
        Reduced arrays
    hsv : ndarray
        Hankel singular values of the original system. The size of the error
        may be calculated based on this.

    Raises
    ------
    ValueError
        If the system that's provided is not stable (i.e. `A` has eigenvalues
        which have non-negative real parts)
    ImportError
        If the slycot subroutine `ab09ad` can't be found. Occurs if the
        slycot package is not installed.

    """
    if check_stability and not isStable(A):
        raise ValueError("This doesn't seem to be a stable system!")
    try:
        from slycot import ab09ad
    except ImportError:
        raise ImportError("can't find slycot subroutine ab09ad")

    if balance:
        job = 'B'
    else:
        job = 'N'

    if scale:
        equil = 'S'
    else:
        equil = 'N'

    dico = 'C'
    n = np.size(A,0)
    m = np.size(B,1)
    p = np.size(C,0)
    nr = k

    if not length_cache_array:
        length_cache_array = 5*n**2 + n*m + n*p + 10*n + m*p

    return ab09ad(dico,job,equil,n,m,p,A,B,C,nr,tol,ldwork=length_cache_array)

def controllabilityTruncation(A,B,C,k,check_stability=True):
    """Truncate the system based on the controllability Gramian

    Solves the Lyapunov Equation for ``AP + PA^H + B B^H`` and computes the
    eigenvalues and eigenvectors of `P` which are used to truncate the system.

    Parameters
    ----------
    A, B, C : array_like
        State-Space matrices of the system that should be reduced
    k : int
        Order of the reduced system
    check_stability : boolean, optional
        Should be checked if A is stable and hence the controllability Gramian
        exists

    Returns
    -------
    Nr : int
        Actual size of the system, based on the error: If the Machine error
        would have been bigger than the error of the reduced system, it may
        happen that ``Nr < k`` and if the error would be inconsiderably bad,
        It might be the case that ``Nr > k``. In case `k` was never specified,
        this is purely based on `tol`
    Ar, Br, Cr : ndarray
        Reduced arrays
    Lambdak : ndarray
        largest k eigenvalues of the Gramian.
    """
    if check_stability and not isStable(A):
        raise ValueError("This doesn't seem to be a stable system!")

    N = A.shape[0]

    P = sp.linalg.solve_lyapunov(A, -np.dot(B, B.H))

    Lambdak, Uk = sp.linalg.eigh(P, eigvals=(N-k,N-1), check_finite=False,
                                 overwrite_a=True, overwrite_b=True)

    UkH = Uk.H

    A = np.dot(UkH, np.dot(A, Uk))
    B = np.dot(UkH, B)
    C = np.dot(C, Uk)

    return k, A, B, C, Lambdak

def isStable(A):
    """Check if all eigenvalues are in the left half of the complex plane"""
    D, V = np.linalg.eig(A)
    return (D.real < 0).all()

class lss(object):
    """linear time independent state space system.

    Default contstructor is called like ``lss(A,B,C,D)`` and a system can
    easily be copied by calling ``lss(sys)`` where `sys` is a
{'fraction_non_alphanumeric': 0.07680662091239403, 'fraction_numerical': 0.004945498587000404, 'mean_word_length': 2.4004804392587507, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1003953', 'n_tokens_mistral': 3038, 'n_tokens_neox': 2912, 'n_words': 1239}