Page MenuHomec4science
Files F93367755

point_distances.py
No OneTemporary
Actions

File Metadata

Created
Thu, Nov 28, 06:11

point_distances.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/>.
#
import numpy as np
from scipy.sparse import spmatrix
from rrompy.utilities.base.types import List, Np1D, Np2D, HFEng
__all__ = ['baseDistanceMatrix', 'vectorDistanceMatrix',
'doubleDistanceMatrix']
def baseDistanceMatrix(x:Np2D, y : Np2D = None, npar : int = None,
magnitude : bool = True, weights : Np1D = None) -> Np2D:
if npar is None: npar = x.shape[1] if x.ndim > 1 else 1
if y is None: y = x
if x.ndim != 3 or x.shape[1] != npar: x = x.reshape(-1, 1, npar)
if y.ndim != 2 or y.shape[1] != npar: y = y.reshape(-1, npar)
dist = np.repeat(x, len(y), axis = 1) - y
if weights is not None: dist *= np.array(weights).flatten()
if magnitude:
if dist.shape[2] == 1:
dist = np.abs(dist)[..., 0]
else:
dist = np.sum(np.abs(dist) ** 2., axis = 2) ** .5
return dist
def vectorDistanceMatrix(X:Np2D, Y:Np2D, HFEngine : HFEng = None,
is_state : bool = True, Xbad : List[bool] = None,
Ybad : List[bool] = None) -> Np2D:
if HFEngine is None:
innerT = np.real(Y.T.conj().dot(X))
if isinstance(X, (spmatrix,)):
if hasattr(X, "data"):
norm2X = np.sum(np.abs(X.data) ** 2., axis = 0)
else:
norm2X = np.sum(np.abs(X.todense()) ** 2., axis = 0)
else:
norm2X = np.sum(np.abs(X) ** 2., axis = 0)
if isinstance(Y, (spmatrix,)):
if hasattr(Y, "data"):
norm2Y = np.sum(np.abs(Y.data) ** 2., axis = 0)
else:
norm2Y = np.sum(np.abs(Y.todense()) ** 2., axis = 0)
else:
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.
if Xbad is None: Xbad = np.where(np.isinf(norm2X))[0]
if Ybad is None: Ybad = np.where(np.isinf(norm2Y))[0]
dist2T = (np.tile(norm2Y.reshape(-1, 1), len(norm2X))
+ norm2X.reshape(1, -1) - 2 * innerT)
dist2T[:, Xbad], dist2T[Ybad, :] = np.inf, np.inf
dist2T[np.ix_(Ybad, Xbad)] = 0.
dist2T[dist2T < 0.] = 0.
return dist2T.T ** .5
def doubleDistanceMatrix(x:Np1D, y:Np1D, w : float = 0, X : Np2D = None,
Y : Np2D = None, HFEngine : HFEng = None,
is_state : bool = True) -> Np2D:
Xbad, Ybad = np.where(np.isinf(x))[0], np.where(np.isinf(y))[0]
dist = vectorDistanceMatrix(np.reshape(x, [1, -1]), np.reshape(y, [1, -1]),
Xbad = Xbad, Ybad = Ybad)
if w == 0: return dist
distAdj = vectorDistanceMatrix(X, Y, HFEngine, is_state, Xbad = Xbad,
Ybad = Ybad)
return ((dist ** 2. + w * distAdj ** 2.) / (1 + w)) ** .5

Event Timeline

c4science · Help