Page MenuHomec4science
Files F61254146

point_matching.py
No OneTemporary
Actions

File Metadata

Created
Sun, May 5, 12:28

point_matching.py
View Options

# Copyright (C) 2018-2020 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/>.
#
from copy import deepcopy as copy
import numpy as np
from scipy.optimize import linear_sum_assignment as LSA
from .tensor_la import dot
from .point_distances import baseDistanceMatrix, doubleDistanceMatrix
from rrompy.utilities.base.types import Tuple, List, ListAny, Np1D, Np2D, HFEng
from rrompy.utilities.exception_manager import RROMPyAssert
__all__ = ['pointMatching', 'rationalFunctionMatching']
def pointMatching(distMatrix:Np2D) -> Tuple[Np1D, Np1D]:
return LSA(distMatrix)
def rationalFunctionMatching(poles:List[Np1D], coeffs:List[Np2D],
featPts:Np2D, matchingWeight:float, supps:ListAny,
projMat:Np2D, HFEngine : HFEng = None,
is_state : bool = True, root : int = None) \
-> Tuple[List[Np1D], List[Np2D]]:
"""
Match poles and residues of a set of rational functions.
Args:
poles: List of (lists of) poles.
coeffs: List of (lists of) residues.
featPts: Marginal parameters corresponding to rational models.
matchingWeight: Matching weight in distance computation.
supps: Support indices for projection matrix.
projMat: Projection matrix for residues.
HFEngine(optional): Engine for distance evaluation. Defaults to None,
i.e. Euclidean metric.
is_state(optional): Whether residues are of system state. Defaults to
True.
root(optional): Root of search tree. Defaults to None, i.e.
automatically chosen.
Returns:
Matched list of (lists of) poles and list of (lists of) residues.
"""
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
featDist = baseDistanceMatrix(featPts)
free = list(range(M))
if root is None:
#start from sample point with closest neighbor,
#among those with no inf pole
notInfPls = np.where([np.any(np.isinf(p)) == False for p in poles])[0]
MEff = len(notInfPls)
if MEff == 1:
root = notInfPls[0]
else:
featDistEff = featDist[notInfPls][:, notInfPls]
root = notInfPls[np.argpartition(featDistEff.flatten(),
MEff)[MEff] % MEff]
polesC = copy(poles)
if matchingWeight != 0:
if hasattr(projMat, "shape"):
resC = [dot(projMat[:, supps[j] : supps[j] + coeffs[j].shape[1]],
coeffs[j][: N].T) for j in range(M)]
else:
resC = [dot(projMat, coeffs[j][: N].T) for j in range(M)]
fixed = [free.pop(root)]
for j in range(M - 1, 0, -1):
#find closest point
idx = np.argmin(featDist[np.ix_(fixed, free)].flatten())
Ifix = fixed[idx // j]
fixed += [free.pop(idx % j)]
Ifree = fixed[-1]
plsfix, plsfree = polesC[Ifix], polesC[Ifree]
freeInf = np.where(np.isinf(plsfree))[0]
freeNotInf = np.where(np.isinf(plsfree) == False)[0]
plsfree = plsfree[freeNotInf]
if matchingWeight == 0:
resfix, resfree = None, None
else:
resfix, resfree = resC[Ifix], resC[Ifree][:, freeNotInf]
#build assignment distance matrix
distj = doubleDistanceMatrix(plsfree, plsfix, matchingWeight, resfree,
resfix, HFEngine, is_state)
reordering = pointMatching(distj)[1]
reorderingInf = [x for x in range(N) if x not in reordering]
#reorder good poles
poles[Ifree][reordering], poles[Ifree][reorderingInf] = (
poles[Ifree][freeNotInf], poles[Ifree][freeInf])
coeffs[Ifree][reordering], coeffs[Ifree][reorderingInf] = (
coeffs[Ifree][freeNotInf], coeffs[Ifree][freeInf])
#transfer missing poles over
polesC[Ifree][reordering], polesC[Ifree][reorderingInf] = (
polesC[Ifree][freeNotInf], polesC[Ifix][reorderingInf])
if matchingWeight != 0:
resC[Ifree][:, reordering], resC[Ifree][:, reorderingInf] = (
resC[Ifree][:, freeNotInf], resC[Ifix][:, reorderingInf])
return poles, coeffs

Event Timeline

c4science · Help