412243578

Build Type

push

github

Committed by Jeffrey Hokanson @ Thor

Commit Message

Fixed bug in other norm approximation

Run Details

308 of 377 branches covered (81.7%)

Branch coverage included in aggregate %.

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

6 existing lines in 1 file now uncovered.

1313 of 1422 relevant lines covered (92.33%)

2.73 hits per line

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

90.1

/polyrat/skiter_stabilized.py

r""" A stabilized variant of the Sanathanan-Koerner iteration

"""

import numpy as np
from numpy.linalg import LinAlgError
from iterprinter import IterationPrinter
from .skiter import _minimize_2_norm, _minimize_inf_norm
from .arnoldi import *
from .rational import *

from .util import minimize_2norm_varpro 



def skfit_stabilized(X, y, num_degree, denom_degree, maxiter = 20, verbose = True, 
        xtol = 1e-7, history = False, denom0 = None, norm = 2):
        r""" The Stabilized Sanathanan-Koerner Iteration


        Parameters
        ----------
        X: np.array (M,dim)
                
        y: np.array (M,)


        Returns
        -------
                        
        """

        if history:
                hist = []

        if denom0 is None:
                denom = np.ones(X.shape[0], dtype = X.dtype)
        else:
                assert denom0.shape[0] == X.shape[0]
                denom = denom0


        if np.isclose(norm, 2):
                linearized_solution = _minimize_2_norm
        elif ~np.isfinite(norm):
                linearized_solution = _minimize_inf_norm
        else: 
                raise NotImplementedError
        
        if verbose:
                printer = IterationPrinter(it = '4d', res_norm = '21.15e', delta_fit = '8.2e', cond = '8.2e')
                printer.print_header(it = 'iter', res_norm = 'residual norm', delta_fit = 'delta fit', cond = 'cond')

        # As there are no guarntees about convergence,
        # we record the best iteration
        best_res_norm = np.inf
        best_sol = None        
        nout_dim = len(y.shape[1:])
        
        # For comparison with current iterate to determine termination
        fit_old = np.zeros(y.shape, dtype = X.dtype)

        for it in range(maxiter):
                try:
                        num_basis = ArnoldiPolynomialBasis(X, num_degree, weight = 1./denom)
                        denom_basis = ArnoldiPolynomialBasis(X, denom_degree, weight = 1./denom)
                        P = num_basis.vandermonde_X
                        Q = denom_basis.vandermonde_X
                        
                        if np.isclose(norm, 2):
                                a, b, cond = minimize_2norm_varpro(P, Q, y)
                        else:                
                                A = np.hstack([P, np.multiply(-y[:,None], Q) ])
                                x, cond = linearized_solution(A)
                                a = x[:P.shape[1]]
                                b = x[-Q.shape[1]:]
        
                        Pa = np.einsum('ij,j...->i...', P, a)
                        Qb = Q @ b

                        #fit = Pa/Qb
                        fit = np.multiply(1./Qb.reshape(-1, *([1,]*nout_dim)), Pa)

                        delta_fit = np.linalg.norm( (fit - fit_old).flatten(), norm)                
                        res_norm = np.linalg.norm( (fit - y).flatten(), norm)
                
                except (LinAlgError) as e:
                        if verbose: print(e)
                        break
        

                # If we have improved the fit, append this 
                if res_norm < best_res_norm:
                        numerator = Polynomial(num_basis, a)
                        denominator = Polynomial(denom_basis, b)
                        best_sol = [numerator, denominator]
                        best_res_norm = res_norm

                if history:
                        hist.append({'fit': fit, 'cond':cond})

                if verbose:
                        printer.print_iter(it = it, delta_fit = delta_fit, res_norm = res_norm, cond = cond) 

                if delta_fit < xtol:
                        break

                denom = np.abs(denom * Qb)
                denom[denom == 0.] = 1.
                fit_old = fit

        if history:        
                return best_sol + [hist]
        else:
                return best_sol



class StabilizedSKRationalApproximation(RationalApproximation, RationalRatio):
        r"""

        Parameters
        ----------
        
        """

        def __init__(self, num_degree, denom_degree, norm = 2, maxiter = 20, verbose = True, xtol = 1e-7):

                RationalApproximation.__init__(self, num_degree, denom_degree)
                self.norm = norm
                self.xtol = float(xtol)

                self.maxiter = int(maxiter)
                self.verbose = verbose
                
                self.numerator = None
                self.denominator = None

        def fit(self, X, y, denom0 = None):
                X = np.array(X)
                y = np.array(y)
                assert X.shape[0] == y.shape[0], "X and y do not have the same number of rows"

                self.numerator, self.denominator, self.hist = skfit_stabilized(
                        X, y, self.num_degree, self.denom_degree,
                        maxiter = self.maxiter, verbose = self.verbose, norm = self.norm,
                        history = True, xtol = self.xtol, denom0 = denom0,
                        )
                

1	r""" A stabilized variant of the Sanathanan-Koerner iteration
2
3	"""
4
5	import numpy as np	3✔
6	from numpy.linalg import LinAlgError	3✔
7	from iterprinter import IterationPrinter	3✔
8	from .skiter import _minimize_2_norm, _minimize_inf_norm	3✔
9	from .arnoldi import *	3✔
10	from .rational import *	3✔
11
12	from .util import minimize_2norm_varpro	3✔
13
14
15
16	def skfit_stabilized(X, y, num_degree, denom_degree, maxiter = 20, verbose = True,	3✔
17	xtol = 1e-7, history = False, denom0 = None, norm = 2):
18	r""" The Stabilized Sanathanan-Koerner Iteration
19
20
21	Parameters
22	----------
23	X: np.array (M,dim)
24
25	y: np.array (M,)
26
27
28	Returns
29	-------
30
31	"""
32
33	if history:	3✔
34	hist = []	3✔
35
36	if denom0 is None:	3!
37	denom = np.ones(X.shape[0], dtype = X.dtype)	3✔
38	else:
UNCOV 39	assert denom0.shape[0] == X.shape[0]	×
UNCOV 40	denom = denom0	×
41
42
43	if np.isclose(norm, 2):	3✔
44	linearized_solution = _minimize_2_norm	3✔
45	elif ~np.isfinite(norm):	3!
46	linearized_solution = _minimize_inf_norm	3✔
47	else:
UNCOV 48	raise NotImplementedError	×
49
50	if verbose:	3!
51	printer = IterationPrinter(it = '4d', res_norm = '21.15e', delta_fit = '8.2e', cond = '8.2e')	3✔
52	printer.print_header(it = 'iter', res_norm = 'residual norm', delta_fit = 'delta fit', cond = 'cond')	3✔
53
54	# As there are no guarntees about convergence,
55	# we record the best iteration
56	best_res_norm = np.inf	3✔
57	best_sol = None	3✔
58	nout_dim = len(y.shape[1:])	3✔
59
60	# For comparison with current iterate to determine termination
61	fit_old = np.zeros(y.shape, dtype = X.dtype)	3✔
62
63	for it in range(maxiter):	3✔
64	try:	3✔
65	num_basis = ArnoldiPolynomialBasis(X, num_degree, weight = 1./denom)	3✔
66	denom_basis = ArnoldiPolynomialBasis(X, denom_degree, weight = 1./denom)	3✔
67	P = num_basis.vandermonde_X	3✔
68	Q = denom_basis.vandermonde_X	3✔
69
70	if np.isclose(norm, 2):	3✔
71	a, b, cond = minimize_2norm_varpro(P, Q, y)	3✔
72	else:
73	A = np.hstack([P, np.multiply(-y[:,None], Q) ])	3✔
74	x, cond = linearized_solution(A)	3✔
75	a = x[:P.shape[1]]	3✔
76	b = x[-Q.shape[1]:]	3✔
77
78	Pa = np.einsum('ij,j...->i...', P, a)	3✔
79	Qb = Q @ b	3✔
80
81	#fit = Pa/Qb
82	fit = np.multiply(1./Qb.reshape(-1, ([1,]nout_dim)), Pa)	3✔
83
84	delta_fit = np.linalg.norm( (fit - fit_old).flatten(), norm)	3✔
85	res_norm = np.linalg.norm( (fit - y).flatten(), norm)	3✔
86
UNCOV 87	except (LinAlgError) as e:	×
UNCOV 88	if verbose: print(e)	×
UNCOV 89	break	×
90
91
92	# If we have improved the fit, append this
93	if res_norm < best_res_norm:	3✔
94	numerator = Polynomial(num_basis, a)	3✔
95	denominator = Polynomial(denom_basis, b)	3✔
96	best_sol = [numerator, denominator]	3✔
97	best_res_norm = res_norm	3✔
98
99	if history:	3✔
100	hist.append({'fit': fit, 'cond':cond})	3✔
101
102	if verbose:	3!
103	printer.print_iter(it = it, delta_fit = delta_fit, res_norm = res_norm, cond = cond)	3✔
104
105	if delta_fit < xtol:	3✔
106	break	3✔
107
108	denom = np.abs(denom * Qb)	3✔
109	denom[denom == 0.] = 1.	3✔
110	fit_old = fit	3✔
111
112	if history:	3✔
113	return best_sol + [hist]	3✔
114	else:
115	return best_sol	3✔
116
117
118
119	class StabilizedSKRationalApproximation(RationalApproximation, RationalRatio):	3✔
120	r"""
121
122	Parameters
123	----------
124
125	"""
126
127	def __init__(self, num_degree, denom_degree, norm = 2, maxiter = 20, verbose = True, xtol = 1e-7):	3✔
128
129	RationalApproximation.__init__(self, num_degree, denom_degree)	3✔
130	self.norm = norm	3✔
131	self.xtol = float(xtol)	3✔
132
133	self.maxiter = int(maxiter)	3✔
134	self.verbose = verbose	3✔
135
136	self.numerator = None	3✔
137	self.denominator = None	3✔
138
139	def fit(self, X, y, denom0 = None):	3✔
140	X = np.array(X)	3✔
141	y = np.array(y)	3✔
142	assert X.shape[0] == y.shape[0], "X and y do not have the same number of rows"	3✔
143
144	self.numerator, self.denominator, self.hist = skfit_stabilized(	3✔
145	X, y, self.num_degree, self.denom_degree,
146	maxiter = self.maxiter, verbose = self.verbose, norm = self.norm,
147	history = True, xtol = self.xtol, denom0 = denom0,
148	)
149

jeffrey-hokanson / polyrat / 412243578

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