Page MenuHomec4science
Files F60928420

point_matching.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 3, 11:00

point_matching.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import warnings
import numpy as np
from scipy.optimize import linear_sum_assignment as LSA
from rrompy.utilities.base.types import Tuple, List, ListAny, Np1D, Np2D, HFEng
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyWarning,
RROMPyAssert)
__all__ = ['pointMatching', 'rationalFunctionMatching', 'potential',
'angleTable', 'chordalMetricTable', 'chordalMetricAdjusted']
def pointMatching(distanceMatrix:Np2D) -> Tuple[Np1D, Np1D]:
return LSA(distanceMatrix)
def rationalFunctionMatching(poles:List[Np1D], coeffs:List[Np2D],
featPts:Np2D, matchingWeight:float,
matchingMode:str, supps:ListAny, projMat:Np2D,
HFEngine : HFEng = None, is_state : bool = True) \
-> Tuple[List[Np1D], List[Np2D]]:
M, N = len(featPts), len(poles[0])
RROMPyAssert(len(poles), M, "Number of rational functions to be matched")
RROMPyAssert(len(coeffs), M, "Number of rational functions to be matched")
if M <= 1: return poles, coeffs
matchingMode = matchingMode.upper().strip().replace(" ", "")
if matchingMode != "NONE":
if matchingMode[: 5] != "SHIFT":
raise RROMPyException("Prescribed matching mode not recognized.")
if "-" in matchingMode:
shiftdeg = int(matchingMode.split("-")[-1])
else:
shiftdeg = 1
if matchingMode == "SHIFT":
avg = [np.mean(pls[np.logical_not(np.isinf(pls))]) for pls in poles]
with warnings.catch_warnings():
from rrompy.utilities.poly_fitting.polynomial import (
PolynomialInterpolator as PI)
for deg in range(shiftdeg, 0, -1):
try:
shift = PI()
shift.setupByInterpolation(featPts, np.array(avg), deg,
verbose = False)
break
except: pass
else:
shift = lambda _: np.mean(avg)
else: #if matchingMode == "NONE":
shift = lambda _: 0.
featDist = np.sum(np.abs(np.repeat(featPts, M, 0)
- np.tile(featPts, [M, 1])), axis = 1)
free = list(range(M))
fixed = [free.pop(np.argpartition(featDist, M)[M] % M)]
featDist = featDist.reshape(M, M)
for j in range(M - 1, 0, -1):
idx = np.argmin(featDist[np.ix_(fixed, free)].flatten())
Ifix = fixed[idx // j]
fixed += [free.pop(idx % j)]
Ifree = fixed[-1]
plsfix = poles[Ifix]
plsfree = (poles[Ifree] + shift([featPts[Ifix]])
- shift([featPts[Ifree]]))
resfix, resfree = None, None
if matchingWeight != 0:
resfix, resfree = coeffs[Ifix][: N].T, coeffs[Ifree][: N].T
if isinstance(projMat, (np.ndarray,)):
suppfix, suppfree = supps[Ifix], supps[Ifree]
resfix = projMat[:, suppfix : suppfix + len(resfix)].dot(
resfix)
resfree = projMat[:, suppfree : suppfree + len(resfree)].dot(
resfree)
distj = chordalMetricAdjusted(plsfix, plsfree, matchingWeight, resfix,
resfree, HFEngine, is_state)
reordering = pointMatching(distj)[1]
poles[Ifree] = poles[Ifree][reordering]
coeffs[Ifree][: N] = coeffs[Ifree][reordering]
return poles, coeffs
def potential(x:Np1D, foci : Tuple[float, float] = [- 1., 1.]) -> Np1D:
mu0 = np.mean(foci)
musig = foci[0] - mu0
isInf = np.isinf(x)
dist = np.empty(len(x))
dist[isInf] = np.inf
xEffR = x[np.logical_not(isInf)] - mu0
if np.isclose(musig, 0.):
if foci[0] != foci[1]:
RROMPyWarning("Collapsing different but numerically equal foci.")
dist[np.logical_not(isInf)] = np.abs(xEffR)
else:
xEffR /= musig
bernEff = (xEffR ** 2. - 1) ** .5
dist[np.logical_not(isInf)] = np.max(np.vstack((
np.abs(xEffR + bernEff), np.abs(xEffR - bernEff)
)), axis = 0)
return dist
def angleTable(X:Np2D, Y:Np2D, HFEngine : HFEng = None,
is_state : bool = True, radius : float = None) -> Np2D:
if HFEngine is None:
innerT = np.real(Y.T.conj().dot(X))
norm2X = np.sum(np.abs(X) ** 2., axis = 0)
norm2Y = np.sum(np.abs(Y) ** 2., axis = 0)
else:
innerT = np.real(HFEngine.innerProduct(X, Y, is_state = is_state))
norm2X = HFEngine.norm(X, is_state = is_state) ** 2.
norm2Y = HFEngine.norm(Y, is_state = is_state) ** 2.
xInf = np.where(np.isclose(norm2X, 0.))[0]
yInf = np.where(np.isclose(norm2Y, 0.))[0]
if radius is None: radius = np.mean(norm2Y) ** .5
dist2T = (np.tile(norm2Y.reshape(-1, 1), len(norm2X))
+ norm2X.reshape(1, -1) - 2 * innerT)
dist2T[:, xInf], dist2T[yInf, :] = 1., 1.
dist2T[np.ix_(yInf, xInf)] = 0.
dist2T[dist2T < 0.] = 0.
return radius * ((dist2T / (norm2X + radius ** 2.)).T
/ (norm2Y + radius ** 2.)) ** .5
def chordalMetricTable(x:Np1D, y:Np1D, radius : float = 1.) -> Np2D:
x, y = np.array(x), np.array(y)
xInf, yInf = np.where(np.isinf(x))[0], np.where(np.isinf(y))[0]
x[xInf], y[yInf] = 0., 0.
distT = np.abs(np.tile(y.reshape(-1, 1), len(x)) - x.reshape(1, -1))
distT[:, xInf], distT[yInf, :] = 1., 1.
distT[np.ix_(yInf, xInf)] = 0.
return radius * ((distT / (np.abs(x) ** 2. + radius ** 2.) ** .5).T
/ (np.abs(y) ** 2. + radius ** 2.) ** .5)
def chordalMetricAdjusted(x:Np1D, y:Np1D, w : float = 0, X : Np2D = None,
Y : Np2D = None, HFEngine : HFEng = None,
is_state : bool = True) -> Np2D:
dist = chordalMetricTable(x, y)
if w == 0: return dist
distAdj = angleTable(X, Y, HFEngine, is_state)
return (dist + w * distAdj) / (1. + w)

Event Timeline

c4science · Help