diff --git a/.coveragerc b/.coveragerc new file mode 100644 index 00000000..fe9847e6 --- /dev/null +++ b/.coveragerc @@ -0,0 +1,5 @@ +[report] +exclude_also = + def __repr__ + raise ValueError + raise NotImplementedError diff --git a/Makefile b/Makefile index 0e0dd17a..0e63be59 100644 --- a/Makefile +++ b/Makefile @@ -34,7 +34,7 @@ types: ${python} -m monkeytype list-modules | grep ${pkg} | parallel -j${j} "${python} -m monkeytype apply {} > /dev/null && echo {}" cov: - ${python} -m pytest -vrs --cov=${pkg} --cov-report html tests + ${python} -m pytest --cov=${pkg} --cov-config .coveragerc --cov-report html tests compile: ${python} _compile.py diff --git a/environment_dev.yaml b/environment_dev.yaml index 0ea84576..c6a6d25c 100644 --- a/environment_dev.yaml +++ b/environment_dev.yaml @@ -14,6 +14,7 @@ dependencies: - pre-commit - pyarrow - pytest + - pytest-cov - scikit-learn - scipy # build diff --git a/src/qmllib/kernels/fdistance.f90 b/src/qmllib/kernels/fdistance.f90 index c096f6af..0bf7b531 100644 --- a/src/qmllib/kernels/fdistance.f90 +++ b/src/qmllib/kernels/fdistance.f90 @@ -1,4 +1,3 @@ - subroutine fmanhattan_distance(A, B, D) implicit none diff --git a/src/qmllib/kernels/fgradient_kernels.f90 b/src/qmllib/kernels/fgradient_kernels.f90 index 1d6f05e9..f27a1124 100644 --- a/src/qmllib/kernels/fgradient_kernels.f90 +++ b/src/qmllib/kernels/fgradient_kernels.f90 @@ -1,25 +1,3 @@ - - - - - - - - - - - - - - - - - - - - - - subroutine fglobal_kernel(x1, x2, q1, q2, n1, n2, nm1, nm2, sigma, kernel) implicit none diff --git a/src/qmllib/kernels/fkernels.f90 b/src/qmllib/kernels/fkernels.f90 index 44e8dc39..5335ba2a 100644 --- a/src/qmllib/kernels/fkernels.f90 +++ b/src/qmllib/kernels/fkernels.f90 @@ -1,4 +1,3 @@ - subroutine fget_local_kernels_gaussian(q1, q2, n1, n2, sigmas, & & nm1, nm2, nsigmas, kernels) diff --git a/src/qmllib/kernels/fkpca.f90 b/src/qmllib/kernels/fkpca.f90 index 2c2e5159..9a6b52f8 100644 --- a/src/qmllib/kernels/fkpca.f90 +++ b/src/qmllib/kernels/fkpca.f90 @@ -1,4 +1,3 @@ - subroutine fkpca(k, n, centering, kpca) implicit none diff --git a/src/qmllib/kernels/fkwasserstein.f90 b/src/qmllib/kernels/fkwasserstein.f90 index 73820f54..90f4ed93 100644 --- a/src/qmllib/kernels/fkwasserstein.f90 +++ b/src/qmllib/kernels/fkwasserstein.f90 @@ -1,4 +1,3 @@ - module searchtools implicit none diff --git a/src/qmllib/kernels/gradient_kernels.py b/src/qmllib/kernels/gradient_kernels.py index fa105d48..07780b4f 100644 --- a/src/qmllib/kernels/gradient_kernels.py +++ b/src/qmllib/kernels/gradient_kernels.py @@ -59,12 +59,10 @@ def get_global_kernel( N1 = np.array([len(Q) for Q in Q1], dtype=np.int32) N2 = np.array([len(Q) for Q in Q2], dtype=np.int32) - assert ( - N1.shape[0] == X1.shape[0] - ), "Error: List of charges does not match shape of representations" - assert ( - N2.shape[0] == X2.shape[0] - ), "Error: List of charges does not match shape of representations" + if not (N1.shape[0] == X1.shape[0]): + raise ValueError("List of charges does not match shape of representations") + if not (N2.shape[0] == X2.shape[0]): + raise ValueError("Error: List of charges does not match shape of representations") Q1_input = np.zeros((max(N1), X1.shape[0]), dtype=np.int32) Q2_input = np.zeros((max(N2), X2.shape[0]), dtype=np.int32) @@ -114,12 +112,10 @@ def get_local_kernels( N1 = np.array([len(Q) for Q in Q1], dtype=np.int32) N2 = np.array([len(Q) for Q in Q2], dtype=np.int32) - assert ( - N1.shape[0] == X1.shape[0] - ), "Error: List of charges does not match shape of representations" - assert ( - N2.shape[0] == X2.shape[0] - ), "Error: List of charges does not match shape of representations" + if not (N1.shape[0] == X1.shape[0]): + raise ValueError("Error: List of charges does not match shape of representations") + if not (N2.shape[0] == X2.shape[0]): + raise ValueError("Error: List of charges does not match shape of representations") Q1_input = np.zeros((max(N1), X1.shape[0]), dtype=np.int32) Q2_input = np.zeros((max(N2), X2.shape[0]), dtype=np.int32) @@ -176,12 +172,10 @@ def get_local_kernel( N1 = np.array([len(Q) for Q in Q1], dtype=np.int32) N2 = np.array([len(Q) for Q in Q2], dtype=np.int32) - assert ( - N1.shape[0] == X1.shape[0] - ), "Error: List of charges does not match shape of representations" - assert ( - N2.shape[0] == X2.shape[0] - ), "Error: List of charges does not match shape of representations" + if not (N1.shape[0] == X1.shape[0]): + raise ValueError("List of charges does not match shape of representations") + if not (N2.shape[0] == X2.shape[0]): + raise ValueError("List of charges does not match shape of representations") Q1_input = np.zeros((max(N1), X1.shape[0]), dtype=np.int32) Q2_input = np.zeros((max(N2), X2.shape[0]), dtype=np.int32) @@ -228,9 +222,8 @@ def get_local_symmetric_kernels(X1: ndarray, Q1: List[List[int]], SIGMAS: List[f N1 = np.array([len(Q) for Q in Q1], dtype=np.int32) - assert ( - N1.shape[0] == X1.shape[0] - ), "Error: List of charges does not match shape of representations" + if not (N1.shape[0] == X1.shape[0]): + raise ValueError("Error: List of charges does not match shape of representations") Q1_input = np.zeros((max(N1), X1.shape[0]), dtype=np.int32) for i, q in enumerate(Q1): @@ -275,9 +268,8 @@ def get_local_symmetric_kernel( N1 = np.array([len(Q) for Q in Q1], dtype=np.int32) - assert ( - N1.shape[0] == X1.shape[0] - ), "Error: List of charges does not match shape of representations" + if not (N1.shape[0] == X1.shape[0]): + raise ValueError("Error: List of charges does not match shape of representations") Q1_input = np.zeros((max(N1), X1.shape[0]), dtype=np.int32) for i, q in enumerate(Q1): @@ -329,12 +321,10 @@ def get_atomic_local_kernel( N1 = np.array([len(Q) for Q in Q1], dtype=np.int32) N2 = np.array([len(Q) for Q in Q2], dtype=np.int32) - assert ( - N1.shape[0] == X1.shape[0] - ), "Error: List of charges does not match shape of representations" - assert ( - N2.shape[0] == X2.shape[0] - ), "Error: List of charges does not match shape of representations" + if not (N1.shape[0] == X1.shape[0]): + raise ValueError("List of charges does not match shape of representations") + if not (N2.shape[0] == X2.shape[0]): + raise ValueError("List of charges does not match shape of representations") Q1_input = np.zeros((max(N1), X1.shape[0]), dtype=np.int32) Q2_input = np.zeros((max(N2), X2.shape[0]), dtype=np.int32) @@ -394,12 +384,10 @@ def get_atomic_local_gradient_kernel( N1 = np.array([len(Q) for Q in Q1], dtype=np.int32) N2 = np.array([len(Q) for Q in Q2], dtype=np.int32) - assert ( - N1.shape[0] == X1.shape[0] - ), "Error: List of charges does not match shape of representations" - assert ( - N2.shape[0] == X2.shape[0] - ), "Error: List of charges does not match shape of representations" + if not (N1.shape[0] == X1.shape[0]): + raise ValueError("List of charges does not match shape of representations") + if not (N2.shape[0] == X2.shape[0]): + raise ValueError("List of charges does not match shape of representations") Q1_input = np.zeros((max(N1), X1.shape[0]), dtype=np.int32) Q2_input = np.zeros((max(N2), X2.shape[0]), dtype=np.int32) @@ -475,12 +463,10 @@ def get_local_gradient_kernel( N1 = np.array([len(Q) for Q in Q1], dtype=np.int32) N2 = np.array([len(Q) for Q in Q2], dtype=np.int32) - assert ( - N1.shape[0] == X1.shape[0] - ), "Error: List of charges does not match shape of representations" - assert ( - N2.shape[0] == X2.shape[0] - ), "Error: List of charges does not match shape of representations" + if not (N1.shape[0] == X1.shape[0]): + raise ValueError("List of charges does not match shape of representations") + if not (N2.shape[0] == X2.shape[0]): + raise ValueError("List of charges does not match shape of representations") Q1_input = np.zeros((max(N1), X1.shape[0]), dtype=np.int32) Q2_input = np.zeros((max(N2), X2.shape[0]), dtype=np.int32) @@ -552,12 +538,10 @@ def get_gdml_kernel( N1 = np.array([len(Q) for Q in Q1], dtype=np.int32) N2 = np.array([len(Q) for Q in Q2], dtype=np.int32) - assert ( - N1.shape[0] == X1.shape[0] - ), "Error: List of charges does not match shape of representations" - assert ( - N2.shape[0] == X2.shape[0] - ), "Error: List of charges does not match shape of representations" + if not (N1.shape[0] == X1.shape[0]): + raise ValueError("List of charges does not match shape of representations") + if not (N2.shape[0] == X2.shape[0]): + raise ValueError("List of charges does not match shape of representations") Q1_input = np.zeros((max(N1), X1.shape[0]), dtype=np.int32) Q2_input = np.zeros((max(N2), X2.shape[0]), dtype=np.int32) @@ -627,9 +611,8 @@ def get_symmetric_gdml_kernel( N1 = np.array([len(Q) for Q in Q1], dtype=np.int32) - assert ( - N1.shape[0] == X1.shape[0] - ), "Error: List of charges does not match shape of representations" + if not (N1.shape[0] == X1.shape[0]): + raise ValueError("List of charges does not match shape of representations") Q1_input = np.zeros((max(N1), X1.shape[0]), dtype=np.int32) @@ -692,12 +675,10 @@ def get_gp_kernel( N1 = np.array([len(Q) for Q in Q1], dtype=np.int32) N2 = np.array([len(Q) for Q in Q2], dtype=np.int32) - assert ( - N1.shape[0] == X1.shape[0] - ), "Error: List of charges does not match shape of representations" - assert ( - N2.shape[0] == X2.shape[0] - ), "Error: List of charges does not match shape of representations" + if not (N1.shape[0] == X1.shape[0]): + raise ValueError("List of charges does not match shape of representations") + if not (N2.shape[0] == X2.shape[0]): + raise ValueError("List of charges does not match shape of representations") Q1_input = np.zeros((max(N1), X1.shape[0]), dtype=np.int32) Q2_input = np.zeros((max(N2), X2.shape[0]), dtype=np.int32) @@ -765,9 +746,8 @@ def get_symmetric_gp_kernel( N1 = np.array([len(Q) for Q in Q1], dtype=np.int32) - assert ( - N1.shape[0] == X1.shape[0] - ), "Error: List of charges does not match shape of representations" + if not (N1.shape[0] == X1.shape[0]): + raise ValueError("List of charges does not match shape of representations") Q1_input = np.zeros((max(N1), X1.shape[0]), dtype=np.int32) diff --git a/src/qmllib/kernels/kernels.py b/src/qmllib/kernels/kernels.py index 6064bb00..0f41d206 100644 --- a/src/qmllib/kernels/kernels.py +++ b/src/qmllib/kernels/kernels.py @@ -265,25 +265,28 @@ def matern_kernel( """ if metric == "l1": + if order == 0: gammas = [] + elif order == 1: gammas = [1] sigma /= np.sqrt(3) + elif order == 2: gammas = [1, 1 / 3.0] sigma /= np.sqrt(5) + else: - print("Order:%d not implemented in Matern Kernel" % order) - raise SystemExit + raise ValueError(f"Order '{order}' not implemented in Matern Kernel") return sargan_kernel(A, B, sigma, gammas) elif metric == "l2": pass + else: - print("Error: Unknown distance metric %s in Matern kernel" % str(metric)) - raise SystemExit + raise ValueError(f"Unknown distance metric {metric} in Matern kernel") na = A.shape[0] nb = B.shape[0] @@ -326,10 +329,13 @@ def get_local_kernels_gaussian( :rtype: numpy array """ - assert np.sum(na) == A.shape[0], "Error in A input" - assert np.sum(nb) == B.shape[0], "Error in B input" + if np.sum(na) != A.shape[0]: + raise ValueError("Error in A input") + if np.sum(nb) != B.shape[0]: + raise ValueError("Error in B input") - assert A.shape[1] == B.shape[1], "Error in representation sizes" + if A.shape[1] != B.shape[1]: + raise ValueError("Error in representation sizes") nma = len(na) nmb = len(nb) @@ -370,10 +376,13 @@ def get_local_kernels_laplacian( :rtype: numpy array """ - assert np.sum(na) == A.shape[0], "Error in A input" - assert np.sum(nb) == B.shape[0], "Error in B input" + if np.sum(na) != A.shape[0]: + raise ValueError("Error in A input") + if np.sum(nb) != B.shape[0]: + raise ValueError("Error in B input") - assert A.shape[1] == B.shape[1], "Error in representation sizes" + if A.shape[1] != B.shape[1]: + raise ValueError("Error in representation sizes") nma = len(na) nmb = len(nb) @@ -403,9 +412,12 @@ def kpca(K: ndarray, n: int = 2, centering: bool = True) -> ndarray: :rtype: numpy array """ - assert K.shape[0] == K.shape[1], "ERROR: Square matrix required for Kernel PCA." - assert np.allclose(K, K.T, atol=1e-8), "ERROR: Symmetric matrix required for Kernel PCA." - assert n <= K.shape[0], "ERROR: Requested more principal components than matrix size." + if K.shape[0] != K.shape[1]: + raise ValueError("Square matrix required for Kernel PCA.") + if not np.allclose(K, K.T, atol=1e-8): + raise ValueError("Symmetric matrix required for Kernel PCA.") + if not n <= K.shape[0]: + raise ValueError("Requested more principal components than matrix size.") size = K.shape[0] pca = fkpca(K, size, centering) diff --git a/src/qmllib/representations/arad/arad.py b/src/qmllib/representations/arad/arad.py index bebc6837..731cb7f4 100644 --- a/src/qmllib/representations/arad/arad.py +++ b/src/qmllib/representations/arad/arad.py @@ -23,8 +23,9 @@ def getAngle(sp: ndarray, norms: ndarray) -> ndarray: return angles + def generate_arad( - coordinates: ndarray, nuclear_charges: ndarray, size: int = 23, cut_distance: float = 5.0 + nuclear_charges: ndarray, coordinates: ndarray, size: int = 23, cut_distance: float = 5.0 ) -> ndarray: """Generates a representation for the ARAD kernel module. @@ -149,9 +150,12 @@ def get_global_kernels_arad( amax = X1.shape[1] - assert X1.shape[3] == amax, "ERROR: Check ARAD decriptor sizes! code = 1" - assert X2.shape[1] == amax, "ERROR: Check ARAD decriptor sizes! code = 2" - assert X2.shape[3] == amax, "ERROR: Check ARAD decriptor sizes! code = 3" + if not X1.shape[3] == amax: + raise ValueError("Check ARAD decriptor sizes") + if not X2.shape[1] == amax: + raise ValueError("Check ARAD decriptor sizes") + if not X2.shape[3] == amax: + raise ValueError("Check ARAD decriptor sizes") nm1 = X1.shape[0] nm2 = X2.shape[0] @@ -257,9 +261,12 @@ def get_local_kernels_arad( amax = X1.shape[1] - assert X1.shape[3] == amax, "ERROR: Check ARAD decriptor sizes! code = 1" - assert X2.shape[1] == amax, "ERROR: Check ARAD decriptor sizes! code = 2" - assert X2.shape[3] == amax, "ERROR: Check ARAD decriptor sizes! code = 3" + if not X1.shape[3] == amax: + raise ValueError("Check ARAD decriptor sizes") + if not X2.shape[1] == amax: + raise ValueError("Check ARAD decriptor sizes") + if not X2.shape[3] == amax: + raise ValueError("Check ARAD decriptor sizes") nm1 = X1.shape[0] nm2 = X2.shape[0] @@ -363,8 +370,10 @@ def get_atomic_kernels_arad( :rtype: numpy array """ - assert len(X1.shape) == 3 - assert len(X2.shape) == 3 + if not len(X1.shape) == 3: + raise ValueError("Expected different shape") + if not len(X2.shape) == 3: + raise ValueError("Expected different shape") na1 = X1.shape[0] na2 = X2.shape[0] @@ -427,7 +436,8 @@ def get_atomic_symmetric_kernels_arad( :rtype: numpy array """ - assert len(X1.shape) == 3 + if not len(X1.shape) == 3: + raise ValueError("Expected different shape") na1 = X1.shape[0] N1 = np.empty(na1, dtype=np.int32) diff --git a/src/qmllib/representations/arad/farad_kernels.f90 b/src/qmllib/representations/arad/farad_kernels.f90 index 9ec5da89..72e082f2 100644 --- a/src/qmllib/representations/arad/farad_kernels.f90 +++ b/src/qmllib/representations/arad/farad_kernels.f90 @@ -1,25 +1,3 @@ -! MIT License -! -! Copyright (c) 2016 Anders Steen Christensen, Felix Faber -! -! Permission is hereby granted, free of charge, to any person obtaining a copy -! of this software and associated documentation files (the "Software"), to deal -! in the Software without restriction, including without limitation the rights -! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -! copies of the Software, and to permit persons to whom the Software is -! furnished to do so, subject to the following conditions: -! -! The above copyright notice and this permission notice shall be included in all -! copies or substantial portions of the Software. -! -! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -! SOFTWARE. - module arad implicit none diff --git a/src/qmllib/representations/facsf.f90 b/src/qmllib/representations/facsf.f90 index 76163791..5f8d85e6 100644 --- a/src/qmllib/representations/facsf.f90 +++ b/src/qmllib/representations/facsf.f90 @@ -1,25 +1,3 @@ - - - - - - - - - - - - - - - - - - - - - - module acsf_utils implicit none diff --git a/src/qmllib/representations/fchl/fchl_electric_field_kernels.py b/src/qmllib/representations/fchl/fchl_electric_field_kernels.py index eae81649..b3cd9e83 100644 --- a/src/qmllib/representations/fchl/fchl_electric_field_kernels.py +++ b/src/qmllib/representations/fchl/fchl_electric_field_kernels.py @@ -78,10 +78,12 @@ def get_atomic_local_electric_field_gradient_kernels( """ atoms_max = A.shape[1] - # TODO Unused neighbors_max = A.shape[3] + neighbors_max = A.shape[3] - # assert B.shape[1] == atoms_max, "ERROR: Check FCHL representation sizes! code = 2" - # assert B.shape[3] == neighbors_max, "ERROR: Check FCHL representation sizes! code = 3" + if not B.shape[1] == atoms_max: + raise ValueError("Check FCHL representation sizes") + if not B.shape[3] == neighbors_max: + raise ValueError("Check FCHL representation sizes") nm1 = A.shape[0] nm2 = B.shape[0] @@ -218,8 +220,10 @@ def get_gaussian_process_electric_field_kernels( atoms_max = A.shape[1] neighbors_max = A.shape[3] - assert B.shape[1] == atoms_max, "ERROR: Check FCHL representation sizes! code = 2" - assert B.shape[3] == neighbors_max, "ERROR: Check FCHL representation sizes! code = 3" + if not B.shape[1] == atoms_max: + raise ValueError("Check FCHL representation sizes") + if not B.shape[3] == neighbors_max: + raise ValueError("Check FCHL representation sizes") nm1 = A.shape[0] nm2 = B.shape[0] @@ -318,8 +322,8 @@ def get_kernels_ef_field( # atoms_max = A.shape[1] # neighbors_max = A.shape[3] - # assert B.shape[1] == atoms_max, "ERROR: Check FCHL representation sizes! code = 2" - # assert B.shape[3] == neighbors_max, "ERROR: Check FCHL representation sizes! code = 3" + # assert B.shape[1] == atoms_max, "Check FCHL representation sizes! code = 2" + # assert B.shape[3] == neighbors_max, "Check FCHL representation sizes! code = 3" # nm1 = A.shape[0] # nm2 = B.shape[0] @@ -449,8 +453,8 @@ def get_kernels_ef_field( # atoms_max = A.shape[1] # neighbors_max = A.shape[3] -# assert B.shape[1] == atoms_max, "ERROR: Check FCHL representation sizes! code = 2" -# assert B.shape[3] == neighbors_max, "ERROR: Check FCHL representation sizes! code = 3" +# assert B.shape[1] == atoms_max, "Check FCHL representation sizes! code = 2" +# assert B.shape[3] == neighbors_max, "Check FCHL representation sizes! code = 3" # nm1 = A.shape[0] diff --git a/src/qmllib/representations/fchl/fchl_force_kernels.py b/src/qmllib/representations/fchl/fchl_force_kernels.py index dea6e139..ac7817d7 100644 --- a/src/qmllib/representations/fchl/fchl_force_kernels.py +++ b/src/qmllib/representations/fchl/fchl_force_kernels.py @@ -1,5 +1,3 @@ -from __future__ import absolute_import - import numpy as np from qmllib.utils.alchemy import get_alchemy @@ -40,19 +38,19 @@ def get_gaussian_process_kernels( nm1 = A.shape[0] nm2 = B.shape[0] - assert B.shape[1] == 3 - assert B.shape[2] == 2 - assert B.shape[5] == 5 - assert A.shape[2] == 5 - # assert B.shape[2] == 2 - # assert B.shape[5] == 5 - atoms_max = B.shape[4] - assert A.shape[1] == atoms_max - assert B.shape[3] == atoms_max - neighbors_max = B.shape[6] - assert A.shape[3] == neighbors_max + + if ( + not B.shape[1] == 3 + or not B.shape[2] == 2 + or not B.shape[5] == 5 + or not A.shape[2] == 5 + or not A.shape[1] == atoms_max + or not B.shape[3] == atoms_max + or not A.shape[3] == neighbors_max + ): + raise ValueError("Representations have wrong shape") N1 = np.zeros((nm1), dtype=np.int32) @@ -142,19 +140,20 @@ def get_local_gradient_kernels( nm1 = A.shape[0] nm2 = B.shape[0] - assert B.shape[1] == 3 - assert B.shape[2] == 2 - assert B.shape[5] == 5 - assert A.shape[2] == 5 - # assert B.shape[2] == 2 - # assert B.shape[5] == 5 - atoms_max = B.shape[4] - assert A.shape[1] == atoms_max - assert B.shape[3] == atoms_max neighbors_max = B.shape[6] - assert A.shape[3] == neighbors_max + + if ( + not B.shape[1] == 3 + or not B.shape[2] == 2 + or not B.shape[5] == 5 + or not A.shape[2] == 5 + or not A.shape[1] == atoms_max + or not B.shape[3] == atoms_max + or not A.shape[3] == neighbors_max + ): + raise ValueError("Representations have wrong shape") N1 = np.zeros((nm1), dtype=np.int32) @@ -244,19 +243,22 @@ def get_local_hessian_kernels( nm1 = A.shape[0] nm2 = B.shape[0] - assert A.shape[1] == 3 - assert A.shape[2] == 2 - assert A.shape[5] == 5 - assert B.shape[1] == 3 - assert B.shape[2] == 2 - assert B.shape[5] == 5 - atoms_max = A.shape[4] - assert A.shape[3] == atoms_max - assert B.shape[3] == atoms_max neighbors_max = A.shape[6] - assert B.shape[6] == neighbors_max + + if ( + not A.shape[1] == 3 + or not A.shape[2] == 2 + or not A.shape[5] == 5 + or not B.shape[1] == 3 + or not B.shape[2] == 2 + or not B.shape[5] == 5 + or not A.shape[3] == atoms_max + or not B.shape[3] == atoms_max + or not B.shape[6] == neighbors_max + ): + raise ValueError("Representations have wrong shape") N1 = np.zeros((nm1), dtype=np.int32) N2 = np.zeros((nm2), dtype=np.int32) @@ -348,15 +350,18 @@ def get_local_symmetric_hessian_kernels( nm1 = A.shape[0] - assert A.shape[1] == 3 - assert A.shape[2] == 2 - assert A.shape[5] == 5 - atoms_max = A.shape[4] - assert A.shape[3] == atoms_max # Unused neighbors_max = A.shape[6] + if ( + not A.shape[1] == 3 + or not A.shape[2] == 2 + or not A.shape[5] == 5 + or not A.shape[3] == atoms_max + ): + raise ValueError("Representations have wrong shape") + N1 = np.zeros((nm1), dtype=np.int32) for a in range(nm1): @@ -431,19 +436,20 @@ def get_force_alphas( nm1 = A.shape[0] nm2 = B.shape[0] - assert B.shape[1] == 3 - assert B.shape[2] == 2 - assert B.shape[5] == 5 - assert A.shape[2] == 5 - # assert B.shape[2] == 2 - # assert B.shape[5] == 5 - atoms_max = B.shape[4] - assert A.shape[1] == atoms_max - assert B.shape[3] == atoms_max neighbors_max = B.shape[6] - assert A.shape[3] == neighbors_max + + if ( + not B.shape[1] == 3 + or not B.shape[2] == 2 + or not B.shape[5] == 5 + or not A.shape[2] == 5 + or not A.shape[1] == atoms_max + or not B.shape[3] == atoms_max + or not A.shape[3] == neighbors_max + ): + raise ValueError("Representations have wrong shape") N1 = np.zeros((nm1), dtype=np.int32) @@ -541,17 +547,20 @@ def get_atomic_local_gradient_kernels( nm1 = A.shape[0] nm2 = B.shape[0] - assert B.shape[1] == 3 - assert B.shape[2] == 2 - assert B.shape[5] == 5 - assert A.shape[2] == 5 - atoms_max = B.shape[4] - assert A.shape[1] == atoms_max - assert B.shape[3] == atoms_max neighbors_max = B.shape[6] - assert A.shape[3] == neighbors_max + + if ( + not B.shape[1] == 3 + or not B.shape[2] == 2 + or not B.shape[5] == 5 + or not A.shape[2] == 5 + or not A.shape[1] == atoms_max + or not B.shape[3] == atoms_max + or not A.shape[3] == neighbors_max + ): + raise ValueError("Representations have wrong shape") N1 = np.zeros((nm1), dtype=np.int32) @@ -643,17 +652,20 @@ def get_atomic_local_gradient_5point_kernels( nm1 = A.shape[0] nm2 = B.shape[0] - assert B.shape[1] == 3 - assert B.shape[2] == 5 - assert B.shape[5] == 5 - assert A.shape[2] == 5 - atoms_max = B.shape[4] - assert A.shape[1] == atoms_max - assert B.shape[3] == atoms_max neighbors_max = B.shape[6] - assert A.shape[3] == neighbors_max + + if ( + not B.shape[1] == 3 + or not B.shape[2] == 5 + or not B.shape[5] == 5 + or not A.shape[2] == 5 + or not A.shape[1] == atoms_max + or not B.shape[3] == atoms_max + or not A.shape[3] == neighbors_max + ): + raise ValueError("Representations have wrong shape") N1 = np.zeros((nm1), dtype=np.int32) diff --git a/src/qmllib/representations/fchl/fchl_kernel_functions.py b/src/qmllib/representations/fchl/fchl_kernel_functions.py index 1af3cc7b..6f605915 100644 --- a/src/qmllib/representations/fchl/fchl_kernel_functions.py +++ b/src/qmllib/representations/fchl/fchl_kernel_functions.py @@ -1,5 +1,3 @@ -from __future__ import division, print_function - from typing import Dict, List, Optional, Tuple, Union import numpy as np @@ -52,8 +50,8 @@ def get_polynomial_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, nd tags = {"alpha": [1.0], "c": [0.0], "d": [1.0]} parameters = np.array([tags["alpha"], tags["c"], tags["d"]]).T - assert len(tags["alpha"]) == len(tags["c"]) - assert len(tags["alpha"]) == len(tags["d"]) + if not (len(tags["alpha"]) == len(tags["c"]) or len(tags["alpha"]) == len(tags["d"])): + raise ValueError("Unexpected parameter dimensions") n_kernels = len(tags["alpha"]) return kt.polynomial, parameters, n_kernels @@ -73,7 +71,9 @@ def get_sigmoid_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, ndarr tags["c"], ] ).T - assert len(tags["alpha"]) == len(tags["c"]) + if not len(tags["alpha"]) == len(tags["c"]): + raise ValueError("Unexpected parameter dimensions") + n_kernels = len(tags["alpha"]) return kt.sigmoid, parameters, n_kernels @@ -125,8 +125,8 @@ def get_bessel_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, ndarra tags = {"sigma": [1.0], "v": [1.0], "n": [1.0]} parameters = np.array([tags["sigma"], tags["v"], tags["n"]]).T - assert len(tags["sigma"]) == len(tags["v"]) - assert len(tags["sigma"]) == len(tags["n"]) + if not (len(tags["sigma"]) == len(tags["v"]) or len(tags["sigma"]) == len(tags["n"])): + raise ValueError("Unexpected parameter dimensions") n_kernels = len(tags["sigma"]) @@ -147,7 +147,8 @@ def get_l2_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, ndarray, i tags["c"], ] ).T - assert len(tags["alpha"]) == len(tags["c"]) + if not len(tags["alpha"]) == len(tags["c"]): + raise ValueError("Unexpected parameter dimensions") n_kernels = len(tags["alpha"]) return kt.l2, parameters, n_kernels @@ -161,7 +162,9 @@ def get_matern_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, ndarra "n": [2.0], } - assert len(tags["sigma"]) == len(tags["n"]) + if not len(tags["sigma"]) == len(tags["n"]): + raise ValueError("Unexpected parameter dimensions") + n_kernels = len(tags["sigma"]) n_max = int(max(tags["n"])) + 1 @@ -260,8 +263,6 @@ def get_kernel_parameters( idx, parameters, n_kernels = get_polynomial2_parameters(tags) else: - - print("QML ERROR: Unsupported kernel specification,", name) - exit() + raise ValueError(f"Unsupported kernel specification '{name}") return idx, parameters, n_kernels diff --git a/src/qmllib/representations/fchl/fchl_representations.py b/src/qmllib/representations/fchl/fchl_representations.py index 39541047..e47c2973 100644 --- a/src/qmllib/representations/fchl/fchl_representations.py +++ b/src/qmllib/representations/fchl/fchl_representations.py @@ -221,13 +221,6 @@ def generate_fchl18_electric_field( # If a list is given, assume these are the fictitious charges - print(fictitious_charges) - print(type(fictitious_charges)) - print(nuclear_charges) - - print(len(fictitious_charges)) - print(len(nuclear_charges)) - if isinstance(fictitious_charges, list) or isinstance(fictitious_charges, np.ndarray): if len(fictitious_charges) != len(nuclear_charges): @@ -267,9 +260,7 @@ def generate_fchl18_electric_field( # partial_charges = [atom.partialcharge for atom in mol] else: - # print("QML ERROR: Unable to parse argument for fictitious charges", fictitious_charges) - # exit() - raise ValueError("Missing charges") + raise ValueError("Missing fictitious charges") size = max_size neighbors = size diff --git a/src/qmllib/representations/fchl/fchl_scalar_kernels.py b/src/qmllib/representations/fchl/fchl_scalar_kernels.py index 9ae05112..1cc07dcf 100644 --- a/src/qmllib/representations/fchl/fchl_scalar_kernels.py +++ b/src/qmllib/representations/fchl/fchl_scalar_kernels.py @@ -1,5 +1,3 @@ -from __future__ import absolute_import - from typing import Dict, List, Optional, Union import numpy as np @@ -88,8 +86,10 @@ def get_local_kernels( atoms_max = A.shape[1] neighbors_max = A.shape[3] - assert B.shape[1] == atoms_max, "ERROR: Check FCHL representation sizes! code = 2" - assert B.shape[3] == neighbors_max, "ERROR: Check FCHL representation sizes! code = 3" + if not B.shape[1] == atoms_max: + raise ValueError("Check FCHL representation sizes") + if not B.shape[3] == neighbors_max: + raise ValueError("Check FCHL representation sizes") nm1 = A.shape[0] nm2 = B.shape[0] @@ -432,8 +432,10 @@ def get_global_kernels( atoms_max = A.shape[1] neighbors_max = A.shape[3] - assert B.shape[1] == atoms_max, "ERROR: Check FCHL representation sizes!" - assert B.shape[3] == neighbors_max, "ERROR: Check FCHL representation sizes!" + if not B.shape[1] == atoms_max: + raise ValueError("Check FCHL representation sizes!") + if not B.shape[3] == neighbors_max: + raise ValueError("Check FCHL representation sizes!") nm1 = A.shape[0] nm2 = B.shape[0] @@ -560,11 +562,8 @@ def get_atomic_kernels( :rtype: numpy array """ - assert len(A.shape) == 3 - assert len(B.shape) == 3 - - # assert B.shape[1] == atoms_max, "ERROR: Check FCHL representation sizes! code = 2" - # assert B.shape[3] == neighbors_max, "ERROR: Check FCHL representation sizes! code = 3" + if not len(A.shape) == 3 or not len(B.shape) == 3: + raise ValueError("Unexpected shapes") na1 = A.shape[0] na2 = B.shape[0] @@ -673,7 +672,8 @@ def get_atomic_symmetric_kernels( :rtype: numpy array """ - assert len(A.shape) == 3 + if not len(A.shape) == 3: + raise ValueError("Unexpected shape") na1 = A.shape[0] @@ -768,8 +768,10 @@ def get_atomic_local_kernels( atoms_max = A.shape[1] neighbors_max = A.shape[3] - assert B.shape[1] == atoms_max, "ERROR: Check FCHL representation sizes! code = 2" - assert B.shape[3] == neighbors_max, "ERROR: Check FCHL representation sizes! code = 3" + if not B.shape[1] == atoms_max: + raise ValueError("Check FCHL representation sizes") + if not B.shape[3] == neighbors_max: + raise ValueError("Check FCHL representation sizes") nm1 = A.shape[0] nm2 = B.shape[0] diff --git a/src/qmllib/representations/frepresentations.f90 b/src/qmllib/representations/frepresentations.f90 index eea7eaca..a0f06158 100644 --- a/src/qmllib/representations/frepresentations.f90 +++ b/src/qmllib/representations/frepresentations.f90 @@ -1,4 +1,3 @@ - module representations implicit none diff --git a/src/qmllib/representations/fslatm.f90 b/src/qmllib/representations/fslatm.f90 index 846614fd..824dd88f 100644 --- a/src/qmllib/representations/fslatm.f90 +++ b/src/qmllib/representations/fslatm.f90 @@ -1,4 +1,3 @@ - module slatm_utils implicit none diff --git a/src/qmllib/representations/representations.py b/src/qmllib/representations/representations.py index c65a3e4b..a1c486d5 100644 --- a/src/qmllib/representations/representations.py +++ b/src/qmllib/representations/representations.py @@ -32,8 +32,7 @@ def vector_to_matrix(v): """ if not (np.sqrt(8 * v.shape[0] + 1) == int(np.sqrt(8 * v.shape[0] + 1))): - print("ERROR: Can not make a square matrix.") - exit(1) + raise ValueError("Can not make a square matrix.") n = v.shape[0] l = (-1 + int(np.sqrt(8 * n + 1))) // 2 @@ -101,8 +100,7 @@ def generate_coulomb_matrix( return fgenerate_unsorted_coulomb_matrix(nuclear_charges, coordinates, size) else: - print("ERROR: Unknown sorting scheme requested") - raise SystemExit + raise ValueError("Unknown sorting scheme requested") def generate_coulomb_matrix_atomic( @@ -205,8 +203,7 @@ def generate_coulomb_matrix_atomic( return np.zeros((0, 0)) else: - print("ERROR: Unknown value %s given for 'indices' variable" % indices) - raise SystemExit + raise ValueError("Unknown value %s given for 'indices' variable" % indices) else: indices = np.asarray(indices, dtype=int) + 1 nindices = indices.size @@ -240,8 +237,7 @@ def generate_coulomb_matrix_atomic( ) else: - print("ERROR: Unknown sorting scheme requested") - raise SystemExit + raise ValueError("Unknown sorting scheme requested") def generate_coulomb_matrix_eigenvalue( @@ -401,8 +397,8 @@ def get_slatm_mbtypes(nuclear_charges: List[ndarray], pbc: str = "000") -> List[ def generate_slatm( - coordinates: ndarray, nuclear_charges: ndarray, + coordinates: ndarray, mbtypes: List[List[int64]], unit_cell: None = None, local: bool = False, @@ -421,10 +417,10 @@ def generate_slatm( NOTE: You will need to run the ``get_slatm_mbtypes()`` function to get the ``mbtypes`` input (or generate it manually). - :param coordinates: Input coordinates - :type coordinates: numpy array :param nuclear_charges: List of nuclear charges. :type nuclear_charges: numpy array + :param coordinates: Input coordinates + :type coordinates: numpy array :param mbtypes: Many-body types for the whole dataset, including 1-, 2- and 3-body types. Could be obtained by calling ``get_slatm_mbtypes()``. :type mbtypes: list :param local: Generate a local representation. Defaulted to False (i.e., global representation); otherwise, atomic version. @@ -451,8 +447,9 @@ def generate_slatm( c = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) if pbc != "000": - # print(' -- handling systems with periodic boundary condition') - assert c is not None, "ERROR: Please specify unit cell for SLATM" + if c is None: + raise ValueError("Please specify unit cell for SLATM") + # ======================================================================= # PBC may introduce new many-body terms, so at the stage of get statistics # info from db, we've already considered this point by letting maximal number @@ -470,7 +467,6 @@ def generate_slatm( mbs = [] X2Ns = [] for ia in range(na): - # if iprt: print ' -- ia = ', ia + 1 n1 = 0 n2 = 0 n3 = 0 @@ -486,7 +482,6 @@ def generate_slatm( ] ), ) - # print ' -- mbsi = ', mbsi if alchemy: n1 = 1 n1_0 = mbs_ia.shape[0] @@ -500,7 +495,6 @@ def generate_slatm( n1 += len(mbsi) mbs_ia = np.concatenate((mbs_ia, mbsi), axis=0) elif len(mbtype) == 2: - # print ' 001, pbc = ', pbc mbsi = get_sbop( mbtype, obj, @@ -513,7 +507,6 @@ def generate_slatm( rpower=rpower, ) mbsi *= 0.5 # only for the two-body parts, local rpst - # print ' 002' if alchemy: n2 = len(mbsi) n2_0 = mbs_ia.shape[0] @@ -557,7 +550,9 @@ def generate_slatm( X2N = [n1, n2, n3] if X2N not in X2Ns: X2Ns.append(X2N) - assert len(X2Ns) == 1, "#ERROR: multiple `X2N ???" + + if len(X2Ns) != 1: + raise ValueError("multiple `X2N ???") else: n1 = 0 @@ -837,8 +832,7 @@ def generate_fchl19( else: if nFourier > 1: - print("Error: FCHL-ACSF only supports nFourier=1, requested", nFourier) - exit() + raise ValueError(f"FCHL-ACSF only supports nFourier=1, requested {nFourier}") (rep, grad) = fgenerate_fchl_acsf_and_gradients( coordinates, diff --git a/src/qmllib/representations/slatm.py b/src/qmllib/representations/slatm.py index 98422660..ddcfdb8f 100644 --- a/src/qmllib/representations/slatm.py +++ b/src/qmllib/representations/slatm.py @@ -131,8 +131,8 @@ def get_sbop( z1, z2 = mbtype zs, coords, c = obj - if iloc: - assert ia is not None, "#ERROR: plz specify `za and `ia " + if iloc and ia is None: + raise ValueError("Please specify za and ia") if pbc != "000": raise NotImplementedError("Periodic boundary conditions not implemented") @@ -179,8 +179,8 @@ def get_sbot( z1, z2, z3 = mbtype zs, coords, c = obj - if iloc: - assert ia is not None, "#ERROR: plz specify `za and `ia " + if iloc and ia is None: + raise ValueError("Please specify za and ia") if pbc != "000": raise NotImplementedError("Periodic boundary conditions not implemented") diff --git a/src/qmllib/solvers/__init__.py b/src/qmllib/solvers/__init__.py index eab39b52..a9ba6c61 100644 --- a/src/qmllib/solvers/__init__.py +++ b/src/qmllib/solvers/__init__.py @@ -231,9 +231,9 @@ def condition_number(A, method="cholesky"): raise ValueError("expected square matrix") if method.lower() == "cholesky": - assert np.allclose( - A, A.T - ), "ERROR: Can't use a Cholesky-decomposition for a non-symmetric matrix." + + if not np.allclose(A, A.T): + raise ValueError("Can't use a Cholesky-decomposition for a non-symmetric matrix.") cond = fcond(A) diff --git a/src/qmllib/solvers/fsolvers.f90 b/src/qmllib/solvers/fsolvers.f90 index 4a3e5aa2..4131f3da 100644 --- a/src/qmllib/solvers/fsolvers.f90 +++ b/src/qmllib/solvers/fsolvers.f90 @@ -1,4 +1,3 @@ - subroutine fcho_solve(A, y, x) implicit none diff --git a/src/qmllib/utils/alchemy.py b/src/qmllib/utils/alchemy.py index fad96104..b8f1e416 100644 --- a/src/qmllib/utils/alchemy.py +++ b/src/qmllib/utils/alchemy.py @@ -307,7 +307,7 @@ def get_alchemy( return doalchemy, pd - raise NotImplementedError("QML ERROR: Unknown alchemical method specified:", alchemy) + raise NotImplementedError(f"Unknown alchemical method specified: {alchemy}") def QNum_distance(a, b, n_width, m_width, l_width, s_width): @@ -404,10 +404,12 @@ def check_if_unique(iterator): num_dims = [] for k, v in e_vec.items(): - assert isinstance(k, int), "Error! Keys need to be int" + if not isinstance(k, int): + raise ValueError("Keys need to be type int") num_dims.append(len(v)) - assert check_if_unique(num_dims), "Error! Unequal number of dimensions" + if not check_if_unique(num_dims): + raise ValueError("Unequal number of dimensions") tmp = np.zeros((emax, num_dims[0])) diff --git a/src/qmllib/utils/fib.f90 b/src/qmllib/utils/fib.f90 deleted file mode 100644 index ced53fd4..00000000 --- a/src/qmllib/utils/fib.f90 +++ /dev/null @@ -1,902 +0,0 @@ -! fib.f90 -subroutine fib(a, n) - USE OMP_LIB - use iso_c_binding - integer(c_int), intent(in) :: n - integer(c_int), intent(out) :: a(n) - - !$OMP PARALLEL - PRINT *, "Hello from process: ", OMP_GET_THREAD_NUM() -!$OMP END PARALLEL - - do i = 1, n - if (i .eq. 1) then - a(i) = 0.0d0 - elseif (i .eq. 2) then - a(i) = 1.0d0 - else - a(i) = a(i - 1) + a(i - 2) - end if - end do -end - -function searchsorted(all_values, sorted) result(cdf_idx) - - implicit none - - double precision, dimension(:), intent(in) :: all_values - double precision, dimension(:), intent(in) :: sorted - - integer, allocatable, dimension(:) :: cdf_idx - - double precision :: val - - integer :: i, j, n, m - - n = size(all_values) - 1 - m = size(sorted) - - allocate (cdf_idx(n)) - - cdf_idx(:) = 0 - - do i = 1, n - - val = all_values(i) - - do j = 1, m - - !write (*,*) i, j, sorted(j), val - - ! if ((sorted(j) <= val) .and. (val < sorted(j+1))) then - if (sorted(j) > val) then - - cdf_idx(i) = j - 1 - !write(*,*) "found" - exit - - ! endif - else !iif (val > maxval(sorted)) then - cdf_idx(i) = m - end if - - end do - - end do - -end function searchsorted - -recursive subroutine quicksort(a, first, last) - implicit none - double precision :: a(*), x, t - integer first, last - integer i, j - - x = a((first + last)/2) - i = first - j = last - do - do while (a(i) < x) - i = i + 1 - end do - do while (x < a(j)) - j = j - 1 - end do - if (i >= j) exit - t = a(i); a(i) = a(j); a(j) = t - i = i + 1 - j = j - 1 - end do - if (first < i - 1) call quicksort(a, first, i - 1) - if (j + 1 < last) call quicksort(a, j + 1, last) -end subroutine quicksort - -subroutine fget_local_kernels_gaussian(q1, q2, n1, n2, sigmas, & - & nm1, nm2, nsigmas, kernels) - - implicit none - - double precision, dimension(:, :), intent(in) :: q1 - double precision, dimension(:, :), intent(in) :: q2 - - ! List of numbers of atoms in each molecule - integer, dimension(:), intent(in) :: n1 - integer, dimension(:), intent(in) :: n2 - - ! Sigma in the Gaussian kernel - double precision, dimension(:), intent(in) :: sigmas - - ! Number of molecules - integer, intent(in) :: nm1 - integer, intent(in) :: nm2 - - ! Number of sigmas - integer, intent(in) :: nsigmas - - ! -1.0 / sigma^2 for use in the kernel - double precision, dimension(nsigmas) :: inv_sigma2 - - ! Resulting alpha vector - double precision, dimension(nsigmas, nm1, nm2), intent(out) :: kernels - - ! Internal counters - integer :: a, b, i, j, k, ni, nj - - ! Temporary variables necessary for parallelization - double precision, allocatable, dimension(:, :) :: atomic_distance - - integer, allocatable, dimension(:) :: i_starts - integer, allocatable, dimension(:) :: j_starts - - allocate (i_starts(nm1)) - allocate (j_starts(nm2)) - - !$OMP PARALLEL DO - do i = 1, nm1 - i_starts(i) = sum(n1(:i)) - n1(i) - end do - !$OMP END PARALLEL DO - - !$OMP PARALLEL DO - do j = 1, nm2 - j_starts(j) = sum(n2(:j)) - n2(j) - end do - !$OMP END PARALLEL DO - - inv_sigma2(:) = -0.5d0/(sigmas(:))**2 - kernels(:, :, :) = 0.0d0 - - allocate (atomic_distance(maxval(n1), maxval(n2))) - atomic_distance(:, :) = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj) SCHEDULE(dynamic) COLLAPSE(2) - do a = 1, nm1 - do b = 1, nm2 - nj = n2(b) - ni = n1(a) - - atomic_distance(:, :) = 0.0d0 - do i = 1, ni - do j = 1, nj - - atomic_distance(i, j) = sum((q1(:, i + i_starts(a)) - q2(:, j + j_starts(b)))**2) - - end do - end do - - do k = 1, nsigmas - kernels(k, a, b) = sum(exp(atomic_distance(:ni, :nj)*inv_sigma2(k))) - end do - - end do - end do - !$OMP END PARALLEL DO - - deallocate (atomic_distance) - deallocate (i_starts) - deallocate (j_starts) - -end subroutine fget_local_kernels_gaussian - -subroutine fget_local_kernels_laplacian(q1, q2, n1, n2, sigmas, & - & nm1, nm2, nsigmas, kernels) - - implicit none - - double precision, dimension(:, :), intent(in) :: q1 - double precision, dimension(:, :), intent(in) :: q2 - - ! List of numbers of atoms in each molecule - integer, dimension(:), intent(in) :: n1 - integer, dimension(:), intent(in) :: n2 - - ! Sigma in the Gaussian kernel - double precision, dimension(:), intent(in) :: sigmas - - ! Number of molecules - integer, intent(in) :: nm1 - integer, intent(in) :: nm2 - - ! Number of sigmas - integer, intent(in) :: nsigmas - - ! -1.0 / sigma^2 for use in the kernel - double precision, dimension(nsigmas) :: inv_sigma2 - - ! Resulting alpha vector - double precision, dimension(nsigmas, nm1, nm2), intent(out) :: kernels - - ! Internal counters - integer :: a, b, i, j, k, ni, nj - - ! Temporary variables necessary for parallelization - double precision, allocatable, dimension(:, :) :: atomic_distance - - integer, allocatable, dimension(:) :: i_starts - integer, allocatable, dimension(:) :: j_starts - - allocate (i_starts(nm1)) - allocate (j_starts(nm2)) - - !$OMP PARALLEL DO - do i = 1, nm1 - i_starts(i) = sum(n1(:i)) - n1(i) - end do - !$OMP END PARALLEL DO - - !$OMP PARALLEL DO - do j = 1, nm2 - j_starts(j) = sum(n2(:j)) - n2(j) - end do - !$OMP END PARALLEL DO - - inv_sigma2(:) = -1.0d0/sigmas(:) - kernels(:, :, :) = 0.0d0 - - allocate (atomic_distance(maxval(n1), maxval(n2))) - atomic_distance(:, :) = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj) SCHEDULE(dynamic) COLLAPSE(2) - do a = 1, nm1 - do b = 1, nm2 - nj = n2(b) - ni = n1(a) - - atomic_distance(:, :) = 0.0d0 - do i = 1, ni - do j = 1, nj - - atomic_distance(i, j) = sum(abs(q1(:, i + i_starts(a)) - q2(:, j + j_starts(b)))) - - end do - end do - - do k = 1, nsigmas - kernels(k, a, b) = sum(exp(atomic_distance(:ni, :nj)*inv_sigma2(k))) - end do - - end do - end do - !$OMP END PARALLEL DO - - deallocate (atomic_distance) - deallocate (i_starts) - deallocate (j_starts) - -end subroutine fget_local_kernels_laplacian - -subroutine fget_vector_kernels_laplacian(q1, q2, n1, n2, sigmas, & - & nm1, nm2, nsigmas, kernels) - - implicit none - - ! Descriptors for the training set - double precision, dimension(:, :, :), intent(in) :: q1 - double precision, dimension(:, :, :), intent(in) :: q2 - - ! List of numbers of atoms in each molecule - integer, dimension(:), intent(in) :: n1 - integer, dimension(:), intent(in) :: n2 - - ! Sigma in the Gaussian kernel - double precision, dimension(:), intent(in) :: sigmas - - ! Number of molecules - integer, intent(in) :: nm1 - integer, intent(in) :: nm2 - - ! Number of sigmas - integer, intent(in) :: nsigmas - - ! -1.0 / sigma^2 for use in the kernel - double precision, dimension(nsigmas) :: inv_sigma - - ! Resulting alpha vector - double precision, dimension(nsigmas, nm1, nm2), intent(out) :: kernels - - ! Internal counters - integer :: i, j, k, ni, nj, ia, ja - - ! Temporary variables necessary for parallelization - double precision, allocatable, dimension(:, :) :: atomic_distance - - inv_sigma(:) = -1.0d0/sigmas(:) - - kernels(:, :, :) = 0.0d0 - - allocate (atomic_distance(maxval(n1), maxval(n2))) - atomic_distance(:, :) = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj) SCHEDULE(dynamic) COLLAPSE(2) - do j = 1, nm2 - do i = 1, nm1 - ni = n1(i) - nj = n2(j) - - atomic_distance(:, :) = 0.0d0 - - do ja = 1, nj - do ia = 1, ni - - atomic_distance(ia, ja) = sum(abs(q1(:, ia, i) - q2(:, ja, j))) - - end do - end do - - do k = 1, nsigmas - kernels(k, i, j) = sum(exp(atomic_distance(:ni, :nj)*inv_sigma(k))) - end do - - end do - end do - !$OMP END PARALLEL DO - - deallocate (atomic_distance) - -end subroutine fget_vector_kernels_laplacian - -subroutine fget_vector_kernels_gaussian(q1, q2, n1, n2, sigmas, & - & nm1, nm2, nsigmas, kernels) - - implicit none - - ! Representations (n_samples, n_max_atoms, rep_size) - double precision, dimension(:, :, :), intent(in) :: q1 - double precision, dimension(:, :, :), intent(in) :: q2 - - ! List of numbers of atoms in each molecule - integer, dimension(:), intent(in) :: n1 - integer, dimension(:), intent(in) :: n2 - - ! Sigma in the Gaussian kernel - double precision, dimension(:), intent(in) :: sigmas - - ! Number of molecules - integer, intent(in) :: nm1 - integer, intent(in) :: nm2 - - ! Number of sigmas - integer, intent(in) :: nsigmas - - ! -1.0 / sigma^2 for use in the kernel - double precision, dimension(nsigmas) :: inv_sigma2 - - ! Resulting alpha vector - double precision, dimension(nsigmas, nm1, nm2), intent(out) :: kernels - - ! Internal counters - integer :: i, j, k, ni, nj, ia, ja - - ! Temporary variables necessary for parallelization - double precision, allocatable, dimension(:, :) :: atomic_distance - - inv_sigma2(:) = -0.5d0/(sigmas(:))**2 - - kernels(:, :, :) = 0.0d0 - - allocate (atomic_distance(maxval(n1), maxval(n2))) - atomic_distance(:, :) = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj,ja,ia) SCHEDULE(dynamic) COLLAPSE(2) - do j = 1, nm2 - do i = 1, nm1 - ni = n1(i) - nj = n2(j) - - atomic_distance(:, :) = 0.0d0 - - do ja = 1, nj - do ia = 1, ni - - atomic_distance(ia, ja) = sum((q1(:, ia, i) - q2(:, ja, j))**2) - - end do - end do - - do k = 1, nsigmas - kernels(k, i, j) = sum(exp(atomic_distance(:ni, :nj)*inv_sigma2(k))) - end do - - end do - end do - !$OMP END PARALLEL DO - - deallocate (atomic_distance) - -end subroutine fget_vector_kernels_gaussian - -subroutine fget_vector_kernels_gaussian_symmetric(q, n, sigmas, & - & nm, nsigmas, kernels) - - implicit none - - ! Representations (rep_size, n_samples, n_max_atoms) - double precision, dimension(:, :, :), intent(in) :: q - - ! List of numbers of atoms in each molecule - integer, dimension(:), intent(in) :: n - - ! Sigma in the Gaussian kernel - double precision, dimension(:), intent(in) :: sigmas - - ! Number of molecules - integer, intent(in) :: nm - - ! Number of sigmas - integer, intent(in) :: nsigmas - - ! Resulting kernels - double precision, dimension(nsigmas, nm, nm), intent(out) :: kernels - - ! Temporary variables necessary for parallelization - double precision, allocatable, dimension(:, :) :: atomic_distance - double precision, allocatable, dimension(:) :: inv_sigma2 - - ! Internal counters - integer :: i, j, k, ni, nj, ia, ja - double precision :: val - - allocate (inv_sigma2(nsigmas)) - - inv_sigma2 = -0.5d0/(sigmas)**2 - - kernels = 1.0d0 - - i = size(q, dim=3) - allocate (atomic_distance(i, i)) - atomic_distance(:, :) = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj,ja,ia,val) SCHEDULE(dynamic) COLLAPSE(2) - do j = 1, nm - do i = 1, nm - if (i .lt. j) cycle - ni = n(i) - nj = n(j) - - atomic_distance(:, :) = 0.0d0 - - do ja = 1, nj - do ia = 1, ni - - atomic_distance(ia, ja) = sum((q(:, ia, i) - q(:, ja, j))**2) - - end do - end do - - do k = 1, nsigmas - val = sum(exp(atomic_distance(:ni, :nj)*inv_sigma2(k))) - kernels(k, i, j) = val - kernels(k, j, i) = val - end do - - end do - end do - !$OMP END PARALLEL DO - - deallocate (atomic_distance) - deallocate (inv_sigma2) - -end subroutine fget_vector_kernels_gaussian_symmetric - -subroutine fget_vector_kernels_laplacian_symmetric(q, n, sigmas, & - & nm, nsigmas, kernels) - - implicit none - - ! Representations (rep_size, n_samples, n_max_atoms) - double precision, dimension(:, :, :), intent(in) :: q - - ! List of numbers of atoms in each molecule - integer, dimension(:), intent(in) :: n - - ! Sigma in the Laplacian kernel - double precision, dimension(:), intent(in) :: sigmas - - ! Number of molecules - integer, intent(in) :: nm - - ! Number of sigmas - integer, intent(in) :: nsigmas - - ! Resulting kernels - double precision, dimension(nsigmas, nm, nm), intent(out) :: kernels - - ! Temporary variables necessary for parallelization - double precision, allocatable, dimension(:, :) :: atomic_distance - double precision, allocatable, dimension(:) :: inv_sigma2 - - ! Internal counters - integer :: i, j, k, ni, nj, ia, ja - double precision :: val - - allocate (inv_sigma2(nsigmas)) - - inv_sigma2 = -1.0d0/sigmas - - kernels = 1.0d0 - - i = size(q, dim=3) - allocate (atomic_distance(i, i)) - atomic_distance(:, :) = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj,ja,ia,val) SCHEDULE(dynamic) COLLAPSE(2) - do j = 1, nm - do i = 1, nm - if (i .lt. j) cycle - ni = n(i) - nj = n(j) - - atomic_distance(:, :) = 0.0d0 - - do ja = 1, nj - do ia = 1, ni - - atomic_distance(ia, ja) = sum(abs(q(:, ia, i) - q(:, ja, j))) - - end do - end do - - do k = 1, nsigmas - val = sum(exp(atomic_distance(:ni, :nj)*inv_sigma2(k))) - kernels(k, i, j) = val - kernels(k, j, i) = val - end do - - end do - end do - !$OMP END PARALLEL DO - - deallocate (atomic_distance) - deallocate (inv_sigma2) - -end subroutine fget_vector_kernels_laplacian_symmetric - -subroutine fgaussian_kernel(a, na, b, nb, k, sigma) - - implicit none - - double precision, dimension(:, :), intent(in) :: a - double precision, dimension(:, :), intent(in) :: b - - integer, intent(in) :: na, nb - - double precision, dimension(:, :), intent(inout) :: k - double precision, intent(in) :: sigma - - double precision, allocatable, dimension(:) :: temp - - double precision :: inv_sigma - integer :: i, j - - inv_sigma = -0.5d0/(sigma*sigma) - - allocate (temp(size(a, dim=1))) - - !$OMP PARALLEL DO PRIVATE(temp) COLLAPSE(2) - do i = 1, nb - do j = 1, na - temp(:) = a(:, j) - b(:, i) - k(j, i) = exp(inv_sigma*dot_product(temp, temp)) - end do - end do - !$OMP END PARALLEL DO - - deallocate (temp) - -end subroutine fgaussian_kernel - -subroutine fgaussian_kernel_symmetric(x, n, k, sigma) - - implicit none - - double precision, dimension(:, :), intent(in) :: x - - integer, intent(in) :: n - - double precision, dimension(:, :), intent(inout) :: k - double precision, intent(in) :: sigma - - double precision, allocatable, dimension(:) :: temp - double precision :: val - - double precision :: inv_sigma - integer :: i, j - - inv_sigma = -0.5d0/(sigma*sigma) - - k = 1.0d0 - - allocate (temp(size(x, dim=1))) - - !$OMP PARALLEL DO PRIVATE(temp, val) SCHEDULE(dynamic) - do i = 1, n - do j = i, n - temp = x(:, j) - x(:, i) - val = exp(inv_sigma*dot_product(temp, temp)) - k(j, i) = val - k(i, j) = val - end do - end do - !$OMP END PARALLEL DO - - deallocate (temp) - -end subroutine fgaussian_kernel_symmetric - -subroutine flaplacian_kernel(a, na, b, nb, k, sigma) - - implicit none - - double precision, dimension(:, :), intent(in) :: a - double precision, dimension(:, :), intent(in) :: b - - integer, intent(in) :: na, nb - - double precision, dimension(:, :), intent(inout) :: k - double precision, intent(in) :: sigma - - double precision :: inv_sigma - - integer :: i, j - - inv_sigma = -1.0d0/sigma - - !$OMP PARALLEL DO COLLAPSE(2) - do i = 1, nb - do j = 1, na - k(j, i) = exp(inv_sigma*sum(abs(a(:, j) - b(:, i)))) - end do - end do - !$OMP END PARALLEL DO - -end subroutine flaplacian_kernel - -subroutine flaplacian_kernel_symmetric(x, n, k, sigma) - - implicit none - - double precision, dimension(:, :), intent(in) :: x - - integer, intent(in) :: n - - double precision, dimension(:, :), intent(inout) :: k - double precision, intent(in) :: sigma - - double precision :: val - - double precision :: inv_sigma - integer :: i, j - - inv_sigma = -1.0d0/sigma - - k = 1.0d0 - - !$OMP PARALLEL DO PRIVATE(val) SCHEDULE(dynamic) - do i = 1, n - do j = i, n - val = exp(inv_sigma*sum(abs(x(:, j) - x(:, i)))) - k(j, i) = val - k(i, j) = val - end do - end do - !$OMP END PARALLEL DO - -end subroutine flaplacian_kernel_symmetric - -subroutine flinear_kernel(a, na, b, nb, k) - - implicit none - - double precision, dimension(:, :), intent(in) :: a - double precision, dimension(:, :), intent(in) :: b - - integer, intent(in) :: na, nb - - double precision, dimension(:, :), intent(inout) :: k - - integer :: i, j - -!$OMP PARALLEL DO COLLAPSE(2) - do i = 1, nb - do j = 1, na - k(j, i) = dot_product(a(:, j), b(:, i)) - end do - end do -!$OMP END PARALLEL DO - -end subroutine flinear_kernel - -subroutine fmatern_kernel_l2(a, na, b, nb, k, sigma, order) - - implicit none - - double precision, dimension(:, :), intent(in) :: a - double precision, dimension(:, :), intent(in) :: b - - integer, intent(in) :: na, nb - - double precision, dimension(:, :), intent(inout) :: k - double precision, intent(in) :: sigma - integer, intent(in) :: order - - double precision, allocatable, dimension(:) :: temp - - double precision :: inv_sigma, inv_sigma2, d, d2 - integer :: i, j - - allocate (temp(size(a, dim=1))) - - if (order == 0) then - inv_sigma = -1.0d0/sigma - - !$OMP PARALLEL DO PRIVATE(temp) COLLAPSE(2) - do i = 1, nb - do j = 1, na - temp(:) = a(:, j) - b(:, i) - k(j, i) = exp(inv_sigma*sqrt(sum(temp*temp))) - end do - end do - !$OMP END PARALLEL DO - else if (order == 1) then - inv_sigma = -sqrt(3.0d0)/sigma - - !$OMP PARALLEL DO PRIVATE(temp, d) COLLAPSE(2) - do i = 1, nb - do j = 1, na - temp(:) = a(:, j) - b(:, i) - d = sqrt(sum(temp*temp)) - k(j, i) = exp(inv_sigma*d)*(1.0d0 - inv_sigma*d) - end do - end do - !$OMP END PARALLEL DO - else - inv_sigma = -sqrt(5.0d0)/sigma - inv_sigma2 = 5.0d0/(3.0d0*sigma*sigma) - - !$OMP PARALLEL DO PRIVATE(temp, d, d2) COLLAPSE(2) - do i = 1, nb - do j = 1, na - temp(:) = a(:, j) - b(:, i) - d2 = sum(temp*temp) - d = sqrt(d2) - k(j, i) = exp(inv_sigma*d)*(1.0d0 - inv_sigma*d + inv_sigma2*d2) - end do - end do - !$OMP END PARALLEL DO - end if - - deallocate (temp) - -end subroutine fmatern_kernel_l2 - -subroutine fsargan_kernel(a, na, b, nb, k, sigma, gammas, ng) - - implicit none - - double precision, dimension(:, :), intent(in) :: a - double precision, dimension(:, :), intent(in) :: b - double precision, dimension(:), intent(in) :: gammas - - integer, intent(in) :: na, nb, ng - - double precision, dimension(:, :), intent(inout) :: k - double precision, intent(in) :: sigma - - double precision, allocatable, dimension(:) :: prefactor - double precision :: inv_sigma - double precision :: d - - integer :: i, j, m - - inv_sigma = -1.0d0/sigma - - ! Allocate temporary - allocate (prefactor(ng)) - - !$OMP PARALLEL DO PRIVATE(d, prefactor) SCHEDULE(dynamic) COLLAPSE(2) - do i = 1, nb - do j = 1, na - d = sum(abs(a(:, j) - b(:, i))) - do m = 1, ng - prefactor(m) = gammas(m)*(-inv_sigma*d)**m - end do - k(j, i) = exp(inv_sigma*d)*(1 + sum(prefactor(:))) - end do - end do - !$OMP END PARALLEL DO - - ! Clean up - deallocate (prefactor) - -end subroutine fsargan_kernel - -subroutine fwasserstein_kernel(a, na, b, nb, k, sigma, p, q) - - implicit none - - double precision, dimension(:, :), intent(in) :: a - double precision, dimension(:, :), intent(in) :: b - - double precision, allocatable, dimension(:, :) :: asorted - double precision, allocatable, dimension(:, :) :: bsorted - - double precision, allocatable, dimension(:) :: rep - - integer, intent(in) :: na, nb - - double precision, dimension(:, :), intent(inout) :: k - double precision, intent(in) :: sigma - - integer, intent(in) :: p - integer, intent(in) :: q - - double precision :: inv_sigma - - integer :: i, j, l - integer :: rep_size - - double precision, allocatable, dimension(:) :: deltas - double precision, allocatable, dimension(:) :: all_values - - double precision, allocatable, dimension(:) :: a_cdf - double precision, allocatable, dimension(:) :: b_cdf - integer, allocatable, dimension(:) :: a_cdf_idx - integer, allocatable, dimension(:) :: b_cdf_idx - - rep_size = size(a, dim=1) - allocate (asorted(rep_size, na)) - allocate (bsorted(rep_size, nb)) - allocate (rep(rep_size)) - - allocate (all_values(rep_size*2)) - allocate (deltas(rep_size*2 - 1)) - - allocate (a_cdf(rep_size*2 - 1)) - allocate (b_cdf(rep_size*2 - 1)) - - allocate (a_cdf_idx(rep_size*2 - 1)) - allocate (b_cdf_idx(rep_size*2 - 1)) - - asorted(:, :) = a(:, :) - bsorted(:, :) = b(:, :) - - do i = 1, na - rep(:) = asorted(:, i) - call quicksort(rep, 1, rep_size) - asorted(:, i) = rep(:) - end do - - do i = 1, nb - rep(:) = bsorted(:, i) - call quicksort(rep, 1, rep_size) - bsorted(:, i) = rep(:) - end do - - !$OMP PARALLEL DO PRIVATE(all_values,a_cdf_idx,b_cdf_idx,a_cdf,b_cdf,deltas) - do j = 1, nb - do i = 1, na - - all_values(:rep_size) = asorted(:, i) - all_values(rep_size + 1:) = bsorted(:, j) - - call quicksort(all_values, 1, 2*rep_size) - - do l = 1, 2*rep_size - 1 - deltas(l) = all_values(l + 1) - all_values(l) - end do - - a_cdf_idx = searchsorted(all_values, asorted(:, i)) - b_cdf_idx = searchsorted(all_values, bsorted(:, j)) - - a_cdf(:) = a_cdf_idx(:) - b_cdf(:) = b_cdf_idx(:) - a_cdf(:) = a_cdf(:)/rep_size - b_cdf(:) = b_cdf(:)/rep_size - - ! k(i,j) = exp(-sum(abs(a_cdf-b_cdf)*deltas)/sigma) - k(i, j) = exp(-(sum((abs(a_cdf - b_cdf)**p)*deltas)**(1.0d0/p))**q/sigma) - - end do - end do - !$OMP END PARALLEL DO -end subroutine fwasserstein_kernel diff --git a/src/qmllib/utils/fsettings.f90 b/src/qmllib/utils/fsettings.f90 index 2a43a145..455604ec 100644 --- a/src/qmllib/utils/fsettings.f90 +++ b/src/qmllib/utils/fsettings.f90 @@ -1,4 +1,3 @@ - subroutine check_openmp(compiled_with_openmp) implicit none diff --git a/tests/test_arad.py b/tests/test_arad.py index de48b207..f92f1060 100644 --- a/tests/test_arad.py +++ b/tests/test_arad.py @@ -30,7 +30,7 @@ def test_arad(): properties.append(data[filename]) for coord, atoms in molecules: - rep = generate_arad(coord, atoms) + rep = generate_arad(atoms, coord) representations.append(rep) representations = np.array(representations) @@ -73,7 +73,7 @@ def test_arad(): molid = 5 coordinates, atoms = molecules[molid] natoms = len(atoms) - X1 = generate_arad(coordinates, atoms, size=natoms) + X1 = generate_arad(atoms, coordinates, size=natoms) XA = X1[:natoms] K_atomic_asymm = get_atomic_kernels_arad(XA, XA, sigmas) diff --git a/tests/test_slatm.py b/tests/test_slatm.py index 4bb5e67b..451d25fa 100644 --- a/tests/test_slatm.py +++ b/tests/test_slatm.py @@ -32,7 +32,7 @@ def test_slatm_global_representation(): representations = [] for coord, atoms in mols: - slatm_vector = generate_slatm(coord, atoms, mbtypes) + slatm_vector = generate_slatm(atoms, coord, mbtypes) representations.append(slatm_vector) X_qml = np.array([rep for rep in representations]) @@ -68,7 +68,7 @@ def test_slatm_local_representation(): for _, mol in enumerate(mols): coord, atoms = mol - slatm_vector = generate_slatm(coord, atoms, mbtypes, local=True) + slatm_vector = generate_slatm(atoms, coord, mbtypes, local=True) local_representations.append(slatm_vector)