4607244921

Build Type

push

github

Committed by GitHub

Commit Message

Adding examples (#77)

Run Details

11 of 91 new or added lines in 5 files covered. (12.09%)

1 existing line in 1 file now uncovered.

830 of 929 relevant lines covered (89.34%)

2.68 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

90.16

/src/SOAPify/analysis.py

"""Module that contains various analysis routines"""
import numpy
from numpy import ndarray
from .distances import simpleSOAPdistance
from MDAnalysis import Universe, AtomGroup
from MDAnalysis.lib.NeighborSearch import AtomNeighborSearch


def tempoSOAP(
    SOAPTrajectory: ndarray,
    window: int = 1,
    stride: int = None,
    backward: bool = False,
    distanceFunction: callable = simpleSOAPdistance,
) -> "tuple[ndarray, ndarray]":
    """performs the 'tempoSOAP' analysis on the given SOAP trajectory

        * Original author: Cristina Caruso
        * Mantainer: Daniele Rapetti
    Args:
        SOAPTrajectory (int):
            an array containg a soap trajectory, the shape is (frames, atom, SOAPdim)
        window (int):
            the dimension of the windows between each state confrontations.
            Defaults to 1.
        stride (int):
            the stride in frames between each state confrontation. **NOT IN USE**.
            Defaults to None.
        deltaT (int): number of frames to skip
        distanceFunction (callable, optional):
            the function that define the distance. Defaults to :func:`SOAPify.distances.simpleSOAPdistance`.

    Returns:
        tuple[numpy.ndarray,numpy.ndarray]: _description_
    """
    if stride is None:
        stride = window
    if stride > window:
        raise ValueError("the window must be bigger than the stride")
    if window >= SOAPTrajectory.shape[0] or stride >= SOAPTrajectory.shape[0]:
        raise ValueError("stride and window must be smaller than simulation lenght")
    timedSOAP = numpy.zeros((SOAPTrajectory.shape[0] - window, SOAPTrajectory.shape[1]))

    for frame in range(window, SOAPTrajectory.shape[0]):
        for molecule in range(0, SOAPTrajectory.shape[1]):
            x = SOAPTrajectory[frame, molecule, :]
            y = SOAPTrajectory[frame - window, molecule, :]
            # fill the matrix (each molecule for each frame)
            timedSOAP[frame - window, molecule] = distanceFunction(x, y)
    # vectorizedDistanceFunction = numpy.vectorize(
    #     distanceFunction, signature="(n),(n)->(1)"
    # )
    # print(SOAPTrajectory.shape)

    deltaTimedSOAP = numpy.diff(timedSOAP.T, axis=-1)

    return timedSOAP, deltaTimedSOAP


def tempoSOAPsimple(
    SOAPTrajectory: ndarray,
    window: int = 1,
    stride: int = None,
    backward: bool = False,
) -> "tuple[ndarray, ndarray]":
    """performs the 'tempoSOAP' analysis on the given SOAP trajectory

        this is optimized to use :func:`SOAPify.distances.simpleSOAPdistance`
        without calling it

        .. warning:: this function works **only** with normalized numpy.float64
          soap vectors!


        * Original author: Cristina Caruso
        * Mantainer: Daniele Rapetti
    Args:
        SOAPTrajectory (int):
            _description_
        window (int):
            the dimension of the windows between each state confrontations.
            Defaults to 1.
        stride (int):
            the stride in frames between each state confrontation. **NOT IN USE**.
            Defaults to None.
        deltaT (int): number of frames to skip

    Returns:
        tuple[numpy.ndarray,numpy.ndarray]:
            - **timedSOAP** the tempoSOAP values, shape(frames-1,natoms)
            - **deltaTimedSOAP** the derivatives of tempoSOAP, shape(natoms, frames-2)
    """
    if stride is None:
        stride = window
    if stride > window:
        raise ValueError("the window must be bigger than the stride")
    if window >= SOAPTrajectory.shape[0] or stride >= SOAPTrajectory.shape[0]:
        raise ValueError("stride and window must be smaller than simulation lenght")

    timedSOAP = numpy.zeros((SOAPTrajectory.shape[0] - window, SOAPTrajectory.shape[1]))
    prev = SOAPTrajectory[0]
    for frame in range(window, SOAPTrajectory.shape[0]):
        actual = SOAPTrajectory[frame]
        # this is equivalent to distance of two normalized SOAP vector
        timedSOAP[frame - window] = numpy.linalg.norm(actual - prev, axis=1)
        prev = actual

    deltaTimedSOAP = numpy.diff(timedSOAP.T, axis=-1)

    return timedSOAP, deltaTimedSOAP


def listNeighboursAlongTrajectory(
    inputUniverse: Universe, cutOff: float, trajSlice: slice = slice(None)
) -> "list[list[AtomGroup]]":
    """produce a per frame list of the neighbours, atom per atom

        * Original author: Martina Crippa
        * Mantainer: Daniele Rapetti
    Args:
        inputUniverse (Universe):
            the universe, or the atomgroup containing the trajectory
        cutOff (float):
            the maximum neighbour distance
        trajSlice (slice, optional):
            the slice of the trajectory to consider. Defaults to slice(None).

    Returns:
        list[list[AtomGroup]]:
            list of AtomGroup wint the neighbours of each atom for each frame
    """
    nnListPerFrame = []
    for ts in inputUniverse.universe.trajectory[trajSlice]:
        nnListPerAtom = []
        nnSearch = AtomNeighborSearch(inputUniverse.atoms, box=inputUniverse.dimensions)
        for atom in inputUniverse.atoms:
            nnListPerAtom.append(nnSearch.search(atom, cutOff))
        nnListPerFrame.append([at.ix for at in nnListPerAtom])
    return nnListPerFrame


def neighbourChangeInTime(
    nnListPerFrame: "list[list[AtomGroup]]",
) -> "tuple[ndarray,ndarray,ndarray,ndarray]":
    """return, listed per each atoms the parameters used in the LENS analysis

        * Original author: Martina Crippa
        * Mantainer: Daniele Rapetti
    Args:
        nnListPerFrame (list[list[AtomGroup]]):
             a frame by frame list of the neighbours of each atom: output of
             :func:`listNeighboursAlongTrajectory
    Returns:
        tuple[numpy.ndarray,numpy.ndarray,numpy.ndarray,numpy.ndarray]:
            - **lensArray** The calculated LENS parameter
            - **numberOfNeighs** the count of neighbours per frame
            - **lensNumerators** the numerators used for calculating LENS parameter
            - **lensDenominators** the denominators used for calculating LENS parameter
    """
    nAt = numpy.asarray(nnListPerFrame, dtype=object).shape[1]
    nFrames = numpy.asarray(nnListPerFrame, dtype=object).shape[0]
    # this is the number of common NN between frames
    lensArray = numpy.zeros((nAt, nFrames))
    # this is the number of NN at that frame
    numberOfNeighs = numpy.zeros((nAt, nFrames))
    # this is the numerator of LENS
    lensNumerators = numpy.zeros((nAt, nFrames))
    # this is the denominator of lens
    lensDenominators = numpy.zeros((nAt, nFrames))
    # each nnlist contains also the atom that generates them,
    # so 0 nn is a 1 element list
    for atomID in range(nAt):
        numberOfNeighs[atomID, 0] = nnListPerFrame[0][atomID].shape[0] - 1
        # let's calculate the numerators and the denominators
        for frame in range(1, nFrames):
            numberOfNeighs[atomID, frame] = nnListPerFrame[frame][atomID].shape[0] - 1
            lensDenominators[atomID, frame] = (
                nnListPerFrame[frame][atomID].shape[0]
                + nnListPerFrame[frame - 1][atomID].shape[0]
                - 2
            )
            lensNumerators[atomID, frame] = numpy.setxor1d(
                nnListPerFrame[frame][atomID], nnListPerFrame[frame - 1][atomID]
            ).shape[0]

    denIsNot0 = lensDenominators != 0
    # lens
    lensArray[denIsNot0] = lensNumerators[denIsNot0] / lensDenominators[denIsNot0]
    return lensArray, numberOfNeighs, lensNumerators, lensDenominators

1	"""Module that contains various analysis routines"""	3✔
2	import numpy	3✔
3	from numpy import ndarray	3✔
4	from .distances import simpleSOAPdistance	3✔
5	from MDAnalysis import Universe, AtomGroup	3✔
6	from MDAnalysis.lib.NeighborSearch import AtomNeighborSearch	3✔
7
8
9	def tempoSOAP(	3✔
10	SOAPTrajectory: ndarray,
11	window: int = 1,
12	stride: int = None,
13	backward: bool = False,
14	distanceFunction: callable = simpleSOAPdistance,
15	) -> "tuple[ndarray, ndarray]":
16	"""performs the 'tempoSOAP' analysis on the given SOAP trajectory
17
18	* Original author: Cristina Caruso
19	* Mantainer: Daniele Rapetti
20	Args:
21	SOAPTrajectory (int):
22	an array containg a soap trajectory, the shape is (frames, atom, SOAPdim)
23	window (int):
24	the dimension of the windows between each state confrontations.
25	Defaults to 1.
26	stride (int):
27	the stride in frames between each state confrontation. NOT IN USE.
28	Defaults to None.
29	deltaT (int): number of frames to skip
30	distanceFunction (callable, optional):
31	the function that define the distance. Defaults to :func:`SOAPify.distances.simpleSOAPdistance`.
32
33	Returns:
34	tuple[numpy.ndarray,numpy.ndarray]: _description_
35	"""
36	if stride is None:	3✔
37	stride = window	×
38	if stride > window:	3✔
39	raise ValueError("the window must be bigger than the stride")	×
40	if window >= SOAPTrajectory.shape[0] or stride >= SOAPTrajectory.shape[0]:	3✔
UNCOV 41	raise ValueError("stride and window must be smaller than simulation lenght")	×
42	timedSOAP = numpy.zeros((SOAPTrajectory.shape[0] - window, SOAPTrajectory.shape[1]))	3✔
43
44	for frame in range(window, SOAPTrajectory.shape[0]):	3✔
45	for molecule in range(0, SOAPTrajectory.shape[1]):	3✔
46	x = SOAPTrajectory[frame, molecule, :]	3✔
47	y = SOAPTrajectory[frame - window, molecule, :]	3✔
48	# fill the matrix (each molecule for each frame)
49	timedSOAP[frame - window, molecule] = distanceFunction(x, y)	3✔
50	# vectorizedDistanceFunction = numpy.vectorize(
51	# distanceFunction, signature="(n),(n)->(1)"
52	# )
53	# print(SOAPTrajectory.shape)
54
55	deltaTimedSOAP = numpy.diff(timedSOAP.T, axis=-1)	3✔
56
57	return timedSOAP, deltaTimedSOAP	3✔
58
59
60	def tempoSOAPsimple(	3✔
61	SOAPTrajectory: ndarray,
62	window: int = 1,
63	stride: int = None,
64	backward: bool = False,
65	) -> "tuple[ndarray, ndarray]":
66	"""performs the 'tempoSOAP' analysis on the given SOAP trajectory
67
68	this is optimized to use :func:`SOAPify.distances.simpleSOAPdistance`
69	without calling it
70
71	.. warning:: this function works only with normalized numpy.float64
72	soap vectors!
73
74
75	* Original author: Cristina Caruso
76	* Mantainer: Daniele Rapetti
77	Args:
78	SOAPTrajectory (int):
79	_description_
80	window (int):
81	the dimension of the windows between each state confrontations.
82	Defaults to 1.
83	stride (int):
84	the stride in frames between each state confrontation. NOT IN USE.
85	Defaults to None.
86	deltaT (int): number of frames to skip
87
88	Returns:
89	tuple[numpy.ndarray,numpy.ndarray]:
90	- timedSOAP the tempoSOAP values, shape(frames-1,natoms)
91	- deltaTimedSOAP the derivatives of tempoSOAP, shape(natoms, frames-2)
92	"""
93	if stride is None:	3✔
94	stride = window	×
95	if stride > window:	3✔
96	raise ValueError("the window must be bigger than the stride")	×
97	if window >= SOAPTrajectory.shape[0] or stride >= SOAPTrajectory.shape[0]:	3✔
98	raise ValueError("stride and window must be smaller than simulation lenght")	×
99
100	timedSOAP = numpy.zeros((SOAPTrajectory.shape[0] - window, SOAPTrajectory.shape[1]))	3✔
101	prev = SOAPTrajectory[0]	3✔
102	for frame in range(window, SOAPTrajectory.shape[0]):	3✔
103	actual = SOAPTrajectory[frame]	3✔
104	# this is equivalent to distance of two normalized SOAP vector
105	timedSOAP[frame - window] = numpy.linalg.norm(actual - prev, axis=1)	3✔
106	prev = actual	3✔
107
108	deltaTimedSOAP = numpy.diff(timedSOAP.T, axis=-1)	3✔
109
110	return timedSOAP, deltaTimedSOAP	3✔
111
112
113	def listNeighboursAlongTrajectory(	3✔
114	inputUniverse: Universe, cutOff: float, trajSlice: slice = slice(None)
115	) -> "list[list[AtomGroup]]":
116	"""produce a per frame list of the neighbours, atom per atom
117
118	* Original author: Martina Crippa
119	* Mantainer: Daniele Rapetti
120	Args:
121	inputUniverse (Universe):
122	the universe, or the atomgroup containing the trajectory
123	cutOff (float):
124	the maximum neighbour distance
125	trajSlice (slice, optional):
126	the slice of the trajectory to consider. Defaults to slice(None).
127
128	Returns:
129	list[list[AtomGroup]]:
130	list of AtomGroup wint the neighbours of each atom for each frame
131	"""
132	nnListPerFrame = []	3✔
133	for ts in inputUniverse.universe.trajectory[trajSlice]:	3✔
134	nnListPerAtom = []	3✔
135	nnSearch = AtomNeighborSearch(inputUniverse.atoms, box=inputUniverse.dimensions)	3✔
136	for atom in inputUniverse.atoms:	3✔
137	nnListPerAtom.append(nnSearch.search(atom, cutOff))	3✔
138	nnListPerFrame.append([at.ix for at in nnListPerAtom])	3✔
139	return nnListPerFrame	3✔
140
141
142	def neighbourChangeInTime(	3✔
143	nnListPerFrame: "list[list[AtomGroup]]",
144	) -> "tuple[ndarray,ndarray,ndarray,ndarray]":
145	"""return, listed per each atoms the parameters used in the LENS analysis
146
147	* Original author: Martina Crippa
148	* Mantainer: Daniele Rapetti
149	Args:
150	nnListPerFrame (list[list[AtomGroup]]):
151	a frame by frame list of the neighbours of each atom: output of
152	:func:`listNeighboursAlongTrajectory
153	Returns:
154	tuple[numpy.ndarray,numpy.ndarray,numpy.ndarray,numpy.ndarray]:
155	- lensArray The calculated LENS parameter
156	- numberOfNeighs the count of neighbours per frame
157	- lensNumerators the numerators used for calculating LENS parameter
158	- lensDenominators the denominators used for calculating LENS parameter
159	"""
160	nAt = numpy.asarray(nnListPerFrame, dtype=object).shape[1]	3✔
161	nFrames = numpy.asarray(nnListPerFrame, dtype=object).shape[0]	3✔
162	# this is the number of common NN between frames
163	lensArray = numpy.zeros((nAt, nFrames))	3✔
164	# this is the number of NN at that frame
165	numberOfNeighs = numpy.zeros((nAt, nFrames))	3✔
166	# this is the numerator of LENS
167	lensNumerators = numpy.zeros((nAt, nFrames))	3✔
168	# this is the denominator of lens
169	lensDenominators = numpy.zeros((nAt, nFrames))	3✔
170	# each nnlist contains also the atom that generates them,
171	# so 0 nn is a 1 element list
172	for atomID in range(nAt):	3✔
173	numberOfNeighs[atomID, 0] = nnListPerFrame[0][atomID].shape[0] - 1	3✔
174	# let's calculate the numerators and the denominators
175	for frame in range(1, nFrames):	3✔
176	numberOfNeighs[atomID, frame] = nnListPerFrame[frame][atomID].shape[0] - 1	3✔
177	lensDenominators[atomID, frame] = (	3✔
178	nnListPerFrame[frame][atomID].shape[0]
179	+ nnListPerFrame[frame - 1][atomID].shape[0]
180	- 2
181	)
182	lensNumerators[atomID, frame] = numpy.setxor1d(	3✔
183	nnListPerFrame[frame][atomID], nnListPerFrame[frame - 1][atomID]
184	).shape[0]
185
186	denIsNot0 = lensDenominators != 0	3✔
187	# lens
188	lensArray[denIsNot0] = lensNumerators[denIsNot0] / lensDenominators[denIsNot0]	3✔
189	return lensArray, numberOfNeighs, lensNumerators, lensDenominators	3✔

GMPavanLab / SOAPify / 4607244921

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous