Page MenuHomec4science
Files F61662673

point_matching.py
No OneTemporary
Actions

File Metadata

Created
Wed, May 8, 04:20

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 .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 buildResiduesForDistance(res:Np2D, projMat:Np2D, supp:int,
projMapping : Np2D = None,
projMappingReal : bool = False) -> Np2D:
if projMapping is not None:
badidx = np.where(projMapping >= len(res))
if len(badidx[1]) > 0: projMapping = projMapping[:, : badidx[1][0]]
res = res[projMapping[0]] * res[projMapping[1]].conj()
if isinstance(projMat, (np.ndarray,)):
res = projMat[:, supp : supp + len(res)].dot(res)
if projMapping is not None and projMappingReal: res = np.real(res)
return res
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,
chordalRadius : Tuple[float, float] = [-1] * 2,
projMapping : Np2D = None,
projMappingReal : bool = False) \
-> 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.
chordalRadius(optional): Radius to be used in chordal metric. If <= 0,
Euclidean metric is used. Defaults to [-1, -1].
projMapping(optional): Mapping for projection based on projMap. Should
be assigned for nonlinear outputs. Defaults to None.
projMappingReal(optional): Whether projection based on projMap is
followed by collapse onto real part. Defaults to False.
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.logical_not(np.any(np.isinf(p)))
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:
resC = [buildResiduesForDistance(coeffs[j][: N].T, projMat, supps[j],
projMapping, projMappingReal)
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.logical_not(np.isinf(plsfree)))[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,
chordalRadius)
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