Page MenuHomec4science
Files F61026869

point_matching.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 4, 01:55

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 numpy as np
from scipy.optimize import linear_sum_assignment as LSA
from rrompy.utilities.base.types import Tuple, Np1D, Np2D, HFEng
from rrompy.utilities.exception_manager import RROMPyWarning
__all__ = ['pointMatching', 'potential', 'angleTable', 'chordalMetricTable',
'chordalMetricAdjusted']
def pointMatching(distanceMatrix:Np2D) -> Tuple[Np1D, Np1D]:
return LSA(distanceMatrix)
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) -> Np2D:
if HFEngine is None:
innerT = Y.dot(X.T.conj())
normX = np.linalg.norm(X, axis = 1)
normY = np.linalg.norm(Y, axis = 1)
else:
innerT = HFEngine.innerProduct(X, Y, is_state = is_state)
normX = HFEngine.norm(X, is_state = is_state)
normY = HFEngine.norm(Y, is_state = is_state)
xInf = np.where(np.isclose(normX, 0.))[0]
normX[xInf] = 1.
return normY - np.abs(innerT / normX).T
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