412224039

Build Type

push

github

Committed by Jeffrey Hokanson @ Thor

Commit Message

Array-valued Vector Fitting

Run Details

306 of 375 branches covered (81.6%)

Branch coverage included in aggregate %.

10 of 10 new or added lines in 1 file covered. (100.0%)

6 existing lines in 1 file now uncovered.

1309 of 1418 relevant lines covered (92.31%)

2.73 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

93.28

/polyrat/vecfit.py

import numpy as np
from copy import deepcopy as copy
from .aaa import _build_cauchy
from .arnoldi import ArnoldiPolynomialBasis
from .linratfit import linearized_ratfit
from .polynomial import Polynomial
from .lagrange import LagrangePolynomialInterpolant
from .basis import MonomialPolynomialBasis, LegendrePolynomialBasis
from .rational import RationalFunction, RationalRatio
from iterprinter import IterationPrinter
import scipy.linalg

from .util import minimize_2norm_varpro, minimize_2norm_dense


def _solve_linearized_vecfit(num_basis, denom_basis, y):
        A = np.hstack([ num_basis, -(denom_basis[:,1:].T * y).T ])
        b = y
        x, res, rank, s = scipy.linalg.lstsq(A, b, overwrite_a = True, overwrite_b = False)
        a = x[0:num_basis.shape[1]]
        b = np.hstack([1, x[num_basis.shape[1]:]])
        return a, b, s[0]/s[-1]

# TODO: Implement a DGB15a variant of the linear solve that allows array-valued y

# TODO: Add matrix-weight option

def vecfit(X, y, num_degree, denom_degree, verbose = True, 
        Basis = ArnoldiPolynomialBasis, poles0 = 'linearized',
        maxiter = 50, ftol = 1e-10, btol = 1e-7):
        r"""Implements Vector Fitting 
        
        See: GS99

        Parameters
        ----------
        X: numpy array  (M, m)
                Input coordinates
        y: numpy array (M,*)
                Output values

        poles0: ['GS', 'linearized', array-like]
                Specifies how the initial poles are to be selected
                * 'GS': as recommeneded by Gustavsen & Semlyen, linearized spaced between largest imaginary value
                * 'linearized', perform a linearized rational fitting
                * array-like: specify an array of denom_degree initial poles
        """
        nout_dim = len(y.shape[1:])
        assert num_degree >= 0 and denom_degree >= 0, "numerator and denominator degrees must be nonnegative integers"
        assert num_degree + 1 >= denom_degree, "Vector fitting requires denominator degree to be at most one less than numerator degree"
        if isinstance(poles0, str):
                assert poles0 in ['linearized', 'GS']

        if verbose:
                printer = IterationPrinter(it = '4d', res = '20.10e', delta = '10.4e', bnorm = '10.4e', cond = '10.4e')
                printer.print_header(it = 'iter', res = 'residual norm', delta = 'Δ fit', bnorm = '‖b - e₀‖₂', cond = 'condion #') 

        if isinstance(poles0, str):
                if poles0 == 'linearized':
                        # Generate initial poles by one step of SK iteration (i.e., the linearized ratfit)
                        numerator, denominator = linearized_ratfit(X, y, num_degree, denom_degree, simultaneous = True)
                        poles = denominator.roots().flatten()
                elif poles0 == 'GS':
                        # Generate initial poles as recommened in GS99, Sec. 3.2 (eqns. 9-10)
                        im_max = np.max(np.abs(X.imag))
                        assert im_max > 0, "Must have a non-zero imaginary extent to use `poles0='GS'` initialization"
                        poles = -im_max/100 + 1j*np.linspace(-im_max, im_max, denom_degree)
        else:
                assert len(poles0) == denom_degree, "Number of poles must match the degree of the denominator"
                poles = np.array(poles0)        

        # Construct the Vandermonde matrix for the remaining terms
        if num_degree - denom_degree >= 0:
                bonus_basis = Basis(X, num_degree - denom_degree)
                V = bonus_basis.vandermonde_X
        else:
                bonus_basis = None
                V = np.zeros((len(y),0))

        r_old = np.zeros(y.shape)

        best_fit = {'residual_norm':np.inf}

        for it in range(maxiter):
                C = _build_cauchy(X, poles)

                num_basis = np.hstack([C, V])
                denom_basis = np.hstack([np.ones((len(X), 1)), C])

                #a, b, cond = _solve_linearized_vecfit(num_basis, denom_basis, y)
                a, b, cond = minimize_2norm_varpro(num_basis, denom_basis, y, P_orth = False, method = 'ls')
                b_norm = b[0] #np.linalg.norm(b)
                b /= b_norm
                a /= b_norm
        
                # Compute the rational approximation
                Pa = np.einsum('ij,j...->i...', num_basis, a)
                Qb = denom_basis @ b
                r = np.multiply(1./Qb.reshape(-1, *([1,]*nout_dim)), Pa)

                residual_norm = np.linalg.norm( (y - r).flatten(), 2)
                delta_norm = np.linalg.norm( (r_old - r).flatten(), 2)
                b_norm = np.linalg.norm(b[1:]/b[0], np.inf)

                if residual_norm < best_fit['residual_norm']:
                        best_fit['residual_norm'] = residual_norm
                        best_fit['a'] = a
                        best_fit['b'] = b
                        best_fit['poles'] = poles

                if verbose:
                        printer.print_iter(it = it, res = residual_norm, delta = delta_norm, bnorm = b_norm, cond = cond)        

                if it == maxiter - 1:
                        if verbose: print("maximum iteration limit reached")
                        break

                if delta_norm < ftol:
                        if verbose: print("terminated due to small change in approximation")
                        break

                if b_norm < btol:
                        if verbose: print("terminated due to denominator being approximately constant")
                        break
        
                # Update roots only if we are going to continue the iteration
                # (if we update these we change the polynomial values)
                # This is the root finding approach that Gustavsen takes in Vector Fitting
                # See Gus06: eq. 5
                poles = np.linalg.eigvals(np.diag(poles) - np.outer(np.ones(len(poles)), b[1:]))

                r_old = r                

        a = best_fit['a']
        b = best_fit['b']
        poles = best_fit['poles']

        if bonus_basis:
                bonus_poly = Polynomial(bonus_basis, a[len(poles):])
        else:
                bonus_poly = Polynomial(MonomialPolynomialBasis(X, 0), 0*a[0:1])

        return a[:len(poles)], b, poles, bonus_poly        


class VectorFittingRationalFunction(RationalRatio):
        def __init__(self, a, b, poles, bonus_poly = None):
                self._a = np.copy(a)
                self._b = np.copy(b)
                self.poles = np.copy(poles)
        
                self.bonus_poly = copy(bonus_poly)

        def eval(self, X):
                C = _build_cauchy(X, self.poles)
                num = np.einsum('ij,j...->i...', C, self._a)
                denom = self._b[0] + C @ self._b[1:]
                
                if self.bonus_poly is not None:
                        num += self.bonus_poly(X)

                nout_dim = len(self._a.shape[1:])
                return np.multiply(1./denom.reshape(-1, *([1,]*nout_dim)), num)

                        

class VectorFittingRationalApproximation(VectorFittingRationalFunction):
        def __init__(self, num_degree, denom_degree, *args, **kwargs):
                self.num_degree = int(num_degree)
                self.denom_degree = int(denom_degree)
                self.args = args
                self.kwargs = kwargs
        
        def fit(self, X, y):
                self._a, self._b, self.poles, self.bonus_poly = vecfit(X, y, self.num_degree, self.denom_degree, *self.args, **self.kwargs)

        

1	import numpy as np	3✔
2	from copy import deepcopy as copy	3✔
3	from .aaa import _build_cauchy	3✔
4	from .arnoldi import ArnoldiPolynomialBasis	3✔
5	from .linratfit import linearized_ratfit	3✔
6	from .polynomial import Polynomial	3✔
7	from .lagrange import LagrangePolynomialInterpolant	3✔
8	from .basis import MonomialPolynomialBasis, LegendrePolynomialBasis	3✔
9	from .rational import RationalFunction, RationalRatio	3✔
10	from iterprinter import IterationPrinter	3✔
11	import scipy.linalg	3✔
12
13	from .util import minimize_2norm_varpro, minimize_2norm_dense	3✔
14
15
16	def _solve_linearized_vecfit(num_basis, denom_basis, y):	3✔
UNCOV 17	A = np.hstack([ num_basis, -(denom_basis[:,1:].T * y).T ])	×
UNCOV 18	b = y	×
UNCOV 19	x, res, rank, s = scipy.linalg.lstsq(A, b, overwrite_a = True, overwrite_b = False)	×
UNCOV 20	a = x[0:num_basis.shape[1]]	×
UNCOV 21	b = np.hstack([1, x[num_basis.shape[1]:]])	×
UNCOV 22	return a, b, s[0]/s[-1]	×
23
24	# TODO: Implement a DGB15a variant of the linear solve that allows array-valued y
25
26	# TODO: Add matrix-weight option
27
28	def vecfit(X, y, num_degree, denom_degree, verbose = True,	3✔
29	Basis = ArnoldiPolynomialBasis, poles0 = 'linearized',
30	maxiter = 50, ftol = 1e-10, btol = 1e-7):
31	r"""Implements Vector Fitting
32
33	See: GS99
34
35	Parameters
36	----------
37	X: numpy array (M, m)
38	Input coordinates
39	y: numpy array (M,*)
40	Output values
41
42	poles0: ['GS', 'linearized', array-like]
43	Specifies how the initial poles are to be selected
44	* 'GS': as recommeneded by Gustavsen & Semlyen, linearized spaced between largest imaginary value
45	* 'linearized', perform a linearized rational fitting
46	* array-like: specify an array of denom_degree initial poles
47	"""
48	nout_dim = len(y.shape[1:])	3✔
49	assert num_degree >= 0 and denom_degree >= 0, "numerator and denominator degrees must be nonnegative integers"	3✔
50	assert num_degree + 1 >= denom_degree, "Vector fitting requires denominator degree to be at most one less than numerator degree"	3✔
51	if isinstance(poles0, str):	3✔
52	assert poles0 in ['linearized', 'GS']	3✔
53
54	if verbose:	3✔
55	printer = IterationPrinter(it = '4d', res = '20.10e', delta = '10.4e', bnorm = '10.4e', cond = '10.4e')	3✔
56	printer.print_header(it = 'iter', res = 'residual norm', delta = 'Δ fit', bnorm = '‖b - e₀‖₂', cond = 'condion #')	3✔
57
58	if isinstance(poles0, str):	3✔
59	if poles0 == 'linearized':	3✔
60	# Generate initial poles by one step of SK iteration (i.e., the linearized ratfit)
61	numerator, denominator = linearized_ratfit(X, y, num_degree, denom_degree, simultaneous = True)	3✔
62	poles = denominator.roots().flatten()	3✔
63	elif poles0 == 'GS':	3!
64	# Generate initial poles as recommened in GS99, Sec. 3.2 (eqns. 9-10)
65	im_max = np.max(np.abs(X.imag))	3✔
66	assert im_max > 0, "Must have a non-zero imaginary extent to use `poles0='GS'` initialization"	3✔
67	poles = -im_max/100 + 1j*np.linspace(-im_max, im_max, denom_degree)	3✔
68	else:
69	assert len(poles0) == denom_degree, "Number of poles must match the degree of the denominator"	3✔
70	poles = np.array(poles0)	3✔
71
72	# Construct the Vandermonde matrix for the remaining terms
73	if num_degree - denom_degree >= 0:	3✔
74	bonus_basis = Basis(X, num_degree - denom_degree)	3✔
75	V = bonus_basis.vandermonde_X	3✔
76	else:
77	bonus_basis = None	3✔
78	V = np.zeros((len(y),0))	3✔
79
80	r_old = np.zeros(y.shape)	3✔
81
82	best_fit = {'residual_norm':np.inf}	3✔
83
84	for it in range(maxiter):	3!
85	C = _build_cauchy(X, poles)	3✔
86
87	num_basis = np.hstack([C, V])	3✔
88	denom_basis = np.hstack([np.ones((len(X), 1)), C])	3✔
89
90	#a, b, cond = _solve_linearized_vecfit(num_basis, denom_basis, y)
91	a, b, cond = minimize_2norm_varpro(num_basis, denom_basis, y, P_orth = False, method = 'ls')	3✔
92	b_norm = b[0] #np.linalg.norm(b)	3✔
93	b /= b_norm	3✔
94	a /= b_norm	3✔
95
96	# Compute the rational approximation
97	Pa = np.einsum('ij,j...->i...', num_basis, a)	3✔
98	Qb = denom_basis @ b	3✔
99	r = np.multiply(1./Qb.reshape(-1, ([1,]nout_dim)), Pa)	3✔
100
101	residual_norm = np.linalg.norm( (y - r).flatten(), 2)	3✔
102	delta_norm = np.linalg.norm( (r_old - r).flatten(), 2)	3✔
103	b_norm = np.linalg.norm(b[1:]/b[0], np.inf)	3✔
104
105	if residual_norm < best_fit['residual_norm']:	3✔
106	best_fit['residual_norm'] = residual_norm	3✔
107	best_fit['a'] = a	3✔
108	best_fit['b'] = b	3✔
109	best_fit['poles'] = poles	3✔
110
111	if verbose:	3✔
112	printer.print_iter(it = it, res = residual_norm, delta = delta_norm, bnorm = b_norm, cond = cond)	3✔
113
114	if it == maxiter - 1:	3✔
115	if verbose: print("maximum iteration limit reached")	3✔
116	break	3✔
117
118	if delta_norm < ftol:	3✔
119	if verbose: print("terminated due to small change in approximation")	3✔
120	break	3✔
121
122	if b_norm < btol:	3✔
123	if verbose: print("terminated due to denominator being approximately constant")	3✔
124	break	3✔
125
126	# Update roots only if we are going to continue the iteration
127	# (if we update these we change the polynomial values)
128	# This is the root finding approach that Gustavsen takes in Vector Fitting
129	# See Gus06: eq. 5
130	poles = np.linalg.eigvals(np.diag(poles) - np.outer(np.ones(len(poles)), b[1:]))	3✔
131
132	r_old = r	3✔
133
134	a = best_fit['a']	3✔
135	b = best_fit['b']	3✔
136	poles = best_fit['poles']	3✔
137
138	if bonus_basis:	3✔
139	bonus_poly = Polynomial(bonus_basis, a[len(poles):])	3✔
140	else:
141	bonus_poly = Polynomial(MonomialPolynomialBasis(X, 0), 0*a[0:1])	3✔
142
143	return a[:len(poles)], b, poles, bonus_poly	3✔
144
145
146	class VectorFittingRationalFunction(RationalRatio):	3✔
147	def __init__(self, a, b, poles, bonus_poly = None):	3✔
148	self._a = np.copy(a)	3✔
149	self._b = np.copy(b)	3✔
150	self.poles = np.copy(poles)	3✔
151
152	self.bonus_poly = copy(bonus_poly)	3✔
153
154	def eval(self, X):	3✔
155	C = _build_cauchy(X, self.poles)	3✔
156	num = np.einsum('ij,j...->i...', C, self._a)	3✔
157	denom = self._b[0] + C @ self._b[1:]	3✔
158
159	if self.bonus_poly is not None:	3!
160	num += self.bonus_poly(X)	3✔
161
162	nout_dim = len(self._a.shape[1:])	3✔
163	return np.multiply(1./denom.reshape(-1, ([1,]nout_dim)), num)	3✔
164
165
166
167	class VectorFittingRationalApproximation(VectorFittingRationalFunction):	3✔
168	def __init__(self, num_degree, denom_degree, args, *kwargs):	3✔
169	self.num_degree = int(num_degree)	3✔
170	self.denom_degree = int(denom_degree)	3✔
171	self.args = args	3✔
172	self.kwargs = kwargs	3✔
173
174	def fit(self, X, y):	3✔
175	self._a, self._b, self.poles, self.bonus_poly = vecfit(X, y, self.num_degree, self.denom_degree, self.args, *self.kwargs)	3✔
176
177

jeffrey-hokanson / polyrat / 412224039

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous