diff --git a/examples/qmlearn.py b/examples/qmlearn.py index 135ea4628..1907263bf 100644 --- a/examples/qmlearn.py +++ b/examples/qmlearn.py @@ -135,8 +135,7 @@ def models(): energies = np.loadtxt("../test/data/hof_qm7.txt", usecols=1) model = qmlearn.representations.CoulombMatrix(data) # Create 1000 random indices - indices = np.arange(1000) - np.random.shuffle(indices) + indices = np.random.choice(np.arange(len(energies)), size=1000, replace=False) representations = model.generate(indices) model = qmlearn.kernels.GaussianKernel(sigma='auto') @@ -184,8 +183,7 @@ def pipelines(): ) # Create 1000 random indices - indices = np.arange(1000) - np.random.shuffle(indices) + indices = np.random.choice(np.arange(len(energies)), size=1000, replace=False) model.fit(indices[:800]) scores = model.score(indices[800:]) @@ -203,8 +201,7 @@ def pipelines(): ) # Create 1000 random indices - indices = np.arange(1000) - np.random.shuffle(indices) + indices = np.random.choice(np.arange(len(energies)), size=1000, replace=False) model.fit(indices[:800]) scores = model.score(indices[800:]) @@ -293,8 +290,7 @@ def cross_validation(): ) # Create 1000 random indices - indices = np.arange(1000) - np.random.shuffle(indices) + indices = np.random.choice(np.arange(len(energies)), size=1000, replace=False) # 3-fold CV of a given model can easily be done scores = sklearn.model_selection.cross_validate(model, indices, cv=3) @@ -344,4 +340,4 @@ def cross_validation(): pipelines() cross_validation() pipelines_2() - pipelines_3() + pipelines_3() \ No newline at end of file diff --git a/qml/fchl/fchl_kernels.py b/qml/fchl/fchl_kernels.py deleted file mode 100644 index feb450ef0..000000000 --- a/qml/fchl/fchl_kernels.py +++ /dev/null @@ -1,1551 +0,0 @@ -# MIT License -# -# Copyright (c) 2017-2018 Felix Faber and Anders Steen Christensen -# -# 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. - -from __future__ import absolute_import - -import numpy as np -import copy - - -from .ffchl_module import fget_kernels_fchl -from .ffchl_module import fget_symmetric_kernels_fchl -from .ffchl_module import fget_global_kernels_fchl -from .ffchl_module import fget_global_symmetric_kernels_fchl -from .ffchl_module import fget_atomic_kernels_fchl -from .ffchl_module import fget_atomic_symmetric_kernels_fchl - -from .ffchl_module import fget_local_full_kernels_fchl -from .ffchl_module import fget_local_gradient_kernels_fchl -from .ffchl_module import fget_local_hessian_kernels_fchl -from .ffchl_module import fget_local_symmetric_hessian_kernels_fchl - -from .ffchl_module import fget_local_invariant_alphas_fchl -from .ffchl_module import fget_atomic_gradient_kernels_fchl -from .ffchl_module import fget_local_atomic_kernels_fchl - -from .ffchl_module import fget_kernels_fchl_ef -from .ffchl_module import fget_kernels_fchl_ef_2ndderiv -from .ffchl_module import fget_kernels_fchl_ef_deriv -from .ffchl_module import fget_gp_kernels_fchl_ef -from .ffchl_module import fget_smooth_atomic_gradient_kernels_fchl - -from .ffchl_module import fget_kernels_fchl_ef_field - -from .fchl_kernel_functions import get_kernel_parameters - -from qml.utils.alchemy import get_alchemy - - -def get_local_kernels(A, B, \ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - """ Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: - - :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\sigma^2} \\big)` - - Where :math:`A_{i}` and :math:`B_{j}` are FCHL representation vectors. - K is calculated analytically using an OpenMP parallel Fortran routine. - Note, that this kernel will ONLY work with FCHL representations as input. - - :param A: Array of FCHL representation - shape=(N, maxsize, 5, maxneighbors). - :type A: numpy array - :param B: Array of FCHL representation - shape=(M, maxsize, 5, maxneighbors). - :type B: numpy array - - :param two_body_scaling: Weight for 2-body terms. - :type two_body_scaling: float - :param three_body_scaling: Weight for 3-body terms. - :type three_body_scaling: float - - :param two_body_width: Gaussian width for 2-body terms - :type two_body_width: float - :param three_body_width: Gaussian width for 3-body terms. - :type three_body_width: float - - :param two_body_power: Powerlaw for :math:`r^{-n}` 2-body terms. - :type two_body_power: float - :param three_body_power: Powerlaw for Axilrod-Teller-Muto 3-body term - :type three_body_power: float - - :param cut_start: The fraction of the cut-off radius at which cut-off damping start. - :type cut_start: float - :param cut_distance: Cut-off radius. (default=5 angstrom) - :type cut_distance: float - - :param fourier_order: 3-body Fourier-expansion truncation order. - :type fourier_order: integer - :param alchemy: Type of alchemical interpolation ``"periodic-table"`` or ``"off"`` are possible options. Disabling alchemical interpolation can yield dramatic speedups. - :type alchemy: string - - :param alchemy_period_width: Gaussian width along periods (columns) in the periodic table. - :type alchemy_period_width: float - :param alchemy_group_width: Gaussian width along groups (rows) in the periodic table. - :type alchemy_group_width: float - - :return: Array of FCHL kernel matrices matrix - shape=(n_sigmas, N, M), - :rtype: numpy array - """ - - 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" - - nm1 = A.shape[0] - nm2 = B.shape[0] - - N1 = np.zeros((nm1),dtype=np.int32) - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - for a in range(nm2): - N2[a] = len(np.where(B[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - neighbors2 = np.zeros((nm2, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - for a, representation in enumerate(B): - ni = N2[a] - for i, x in enumerate(representation[:ni]): - neighbors2[a,i] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, n_kernels = get_kernel_parameters(kernel, kernel_args) - - return fget_kernels_fchl(A, B, N1, N2, neighbors1, neighbors2, nm1, nm2, n_kernels, \ - three_body_width, two_body_width, cut_start, cut_distance, fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, three_body_power, kernel_idx, kernel_parameters) - - -def get_local_symmetric_kernels(A, \ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - """ Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: - - :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - A_j\\|_2^2}{2\sigma^2} \\big)` - - Where :math:`A_{i}` and :math:`A_{j}` are FCHL representation vectors. - K is calculated analytically using an OpenMP parallel Fortran routine. - Note, that this kernel will ONLY work with FCHL representations as input. - - :param A: Array of FCHL representation - shape=(N, maxsize, 5, maxneighbors). - :type A: numpy array - - :param two_body_scaling: Weight for 2-body terms. - :type two_body_scaling: float - :param three_body_scaling: Weight for 3-body terms. - :type three_body_scaling: float - - :param two_body_width: Gaussian width for 2-body terms - :type two_body_width: float - :param three_body_width: Gaussian width for 3-body terms. - :type three_body_width: float - - :param two_body_power: Powerlaw for :math:`r^{-n}` 2-body terms. - :type two_body_power: float - :param three_body_power: Powerlaw for Axilrod-Teller-Muto 3-body term - :type three_body_power: float - - :param cut_start: The fraction of the cut-off radius at which cut-off damping start. - :type cut_start: float - :param cut_distance: Cut-off radius. (default=5 angstrom) - :type cut_distance: float - - :param fourier_order: 3-body Fourier-expansion truncation order. - :type fourier_order: integer - :param alchemy: Type of alchemical interpolation ``"periodic-table"`` or ``"off"`` are possible options. Disabling alchemical interpolation can yield dramatic speedups. - :type alchemy: string - - :param alchemy_period_width: Gaussian width along periods (columns) in the periodic table. - :type alchemy_period_width: float - :param alchemy_group_width: Gaussian width along groups (rows) in the periodic table. - :type alchemy_group_width: float - - :return: Array of FCHL kernel matrices matrix - shape=(n_sigmas, N, N), - :rtype: numpy array - """ - - atoms_max = A.shape[1] - neighbors_max = A.shape[3] - - nm1 = A.shape[0] - N1 = np.zeros((nm1),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - kernel_idx, kernel_parameters, n_kernels = get_kernel_parameters(kernel, kernel_args) - - return fget_symmetric_kernels_fchl(A, N1, neighbors1, nm1, n_kernels, \ - three_body_width, two_body_width, cut_start, cut_distance, fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, three_body_power, kernel_idx, kernel_parameters) - - -def get_global_symmetric_kernels(A, \ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - """ Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: - - :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - A_j\\|_2^2}{2\sigma^2} \\big)` - - Where :math:`A_{i}` and :math:`A_{j}` are FCHL representation vectors. - K is calculated analytically using an OpenMP parallel Fortran routine. - Note, that this kernel will ONLY work with FCHL representations as input. - - :param A: Array of FCHL representation - shape=(N, maxsize, 5, maxneighbors). - :type A: numpy array - - :param two_body_scaling: Weight for 2-body terms. - :type two_body_scaling: float - :param three_body_scaling: Weight for 3-body terms. - :type three_body_scaling: float - - :param two_body_width: Gaussian width for 2-body terms - :type two_body_width: float - :param three_body_width: Gaussian width for 3-body terms. - :type three_body_width: float - - :param two_body_power: Powerlaw for :math:`r^{-n}` 2-body terms. - :type two_body_power: float - :param three_body_power: Powerlaw for Axilrod-Teller-Muto 3-body term - :type three_body_power: float - - :param cut_start: The fraction of the cut-off radius at which cut-off damping start. - :type cut_start: float - :param cut_distance: Cut-off radius. (default=5 angstrom) - :type cut_distance: float - - :param fourier_order: 3-body Fourier-expansion truncation order. - :type fourier_order: integer - :param alchemy: Type of alchemical interpolation ``"periodic-table"`` or ``"off"`` are possible options. Disabling alchemical interpolation can yield dramatic speedups. - :type alchemy: string - - :param alchemy_period_width: Gaussian width along periods (columns) in the periodic table. - :type alchemy_period_width: float - :param alchemy_group_width: Gaussian width along groups (rows) in the periodic table. - :type alchemy_group_width: float - - :return: Array of FCHL kernel matrices matrix - shape=(n_sigmas, N, N), - :rtype: numpy array - """ - - atoms_max = A.shape[1] - neighbors_max = A.shape[3] - - nm1 = A.shape[0] - N1 = np.zeros((nm1),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - kernel_idx, kernel_parameters, nsigmas = get_kernel_parameters(kernel, kernel_args) - - return fget_global_symmetric_kernels_fchl(A, N1, neighbors1, nm1, nsigmas, \ - three_body_width, two_body_width, cut_start, cut_distance, fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, three_body_power, kernel_idx, kernel_parameters) - - -def get_global_kernels(A, B, \ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - """ Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: - - :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\sigma^2} \\big)` - - Where :math:`A_{i}` and :math:`B_{j}` are FCHL representation vectors. - K is calculated analytically using an OpenMP parallel Fortran routine. - Note, that this kernel will ONLY work with FCHL representations as input. - - :param A: Array of FCHL representation - shape=(N, maxsize, 5, maxneighbors). - :type A: numpy array - :param B: Array of FCHL representation - shape=(M, maxsize, 5, maxneighbors). - :type B: numpy array - - :param two_body_scaling: Weight for 2-body terms. - :type two_body_scaling: float - :param three_body_scaling: Weight for 3-body terms. - :type three_body_scaling: float - - :param two_body_width: Gaussian width for 2-body terms - :type two_body_width: float - :param three_body_width: Gaussian width for 3-body terms. - :type three_body_width: float - - :param two_body_power: Powerlaw for :math:`r^{-n}` 2-body terms. - :type two_body_power: float - :param three_body_power: Powerlaw for Axilrod-Teller-Muto 3-body term - :type three_body_power: float - - :param cut_start: The fraction of the cut-off radius at which cut-off damping start. - :type cut_start: float - :param cut_distance: Cut-off radius. (default=5 angstrom) - :type cut_distance: float - - :param fourier_order: 3-body Fourier-expansion truncation order. - :type fourier_order: integer - :param alchemy: Type of alchemical interpolation ``"periodic-table"`` or ``"off"`` are possible options. Disabling alchemical interpolation can yield dramatic speedups. - :type alchemy: string - - :param alchemy_period_width: Gaussian width along periods (columns) in the periodic table. - :type alchemy_period_width: float - :param alchemy_group_width: Gaussian width along groups (rows) in the periodic table. - :type alchemy_group_width: float - - :return: Array of FCHL kernel matrices matrix - shape=(n_sigmas, N, M), - :rtype: numpy array - """ - - 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!" - - nm1 = A.shape[0] - nm2 = B.shape[0] - - N1 = np.zeros((nm1),dtype=np.int32) - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - for a in range(nm2): - N2[a] = len(np.where(B[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - neighbors2 = np.zeros((nm2, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - for a, representation in enumerate(B): - ni = N2[a] - for i, x in enumerate(representation[:ni]): - neighbors2[a,i] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - kernel_idx, kernel_parameters, nsigmas = get_kernel_parameters(kernel, kernel_args) - - return fget_global_kernels_fchl(A, B, N1, N2, neighbors1, neighbors2, nm1, nm2, nsigmas, \ - three_body_width, two_body_width, cut_start, cut_distance, fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, three_body_power, kernel_idx, kernel_parameters) - - -def get_atomic_kernels(A, B, \ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - """ Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: - - :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\sigma^2} \\big)` - - Where :math:`A_{i}` and :math:`B_{j}` are FCHL representation vectors. - K is calculated analytically using an OpenMP parallel Fortran routine. - Note, that this kernel will ONLY work with FCHL representations as input. - - :param A: Array of FCHL representation - shape=(N, maxsize, 5, size). - :type A: numpy array - :param B: Array of FCHL representation - shape=(M, maxsize, 5, size). - :type B: numpy array - :param sigma: List of kernel-widths. - :type sigma: list - - :param two_body_scaling: Weight for 2-body terms. - :type two_body_scaling: float - :param three_body_scaling: Weight for 3-body terms. - :type three_body_scaling: float - - :param two_body_width: Gaussian width for 2-body terms - :type two_body_width: float - :param three_body_width: Gaussian width for 3-body terms. - :type three_body_width: float - - :param two_body_power: Powerlaw for :math:`r^{-n}` 2-body terms. - :type two_body_power: float - :param three_body_power: Powerlaw for Axilrod-Teller-Muto 3-body term - :type three_body_power: float - - :param cut_start: The fraction of the cut-off radius at which cut-off damping start. - :type cut_start: float - :param cut_distance: Cut-off radius. (default=5 angstrom) - :type cut_distance: float - - :param fourier_order: 3-body Fourier-expansion truncation order. - :type fourier_order: integer - :param alchemy: Type of alchemical interpolation ``"periodic-table"`` or ``"off"`` are possible options. Disabling alchemical interpolation can yield dramatic speedups. - :type alchemy: string - - :param alchemy_period_width: Gaussian width along periods (columns) in the periodic table. - :type alchemy_period_width: float - :param alchemy_group_width: Gaussian width along groups (rows) in the periodic table. - :type alchemy_group_width: float - - :return: Array of FCHL kernel matrices matrix - shape=(n_sigmas, N, M), - :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" - - na1 = A.shape[0] - na2 = B.shape[0] - - neighbors1 = np.zeros((na1), dtype=np.int32) - neighbors2 = np.zeros((na2), dtype=np.int32) - - for i, x in enumerate(A): - neighbors1[i] = len(np.where(x[0]< cut_distance)[0]) - - for i, x in enumerate(B): - neighbors2[i] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - kernel_idx, kernel_parameters, nsigmas = get_kernel_parameters(kernel, kernel_args) - - return fget_atomic_kernels_fchl(A, B, neighbors1, neighbors2, na1, na2, nsigmas, \ - three_body_width, two_body_width, cut_start, cut_distance, fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, three_body_power, kernel_idx, kernel_parameters) - - -def get_atomic_symmetric_kernels(A, \ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - """ Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: - - :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\sigma^2} \\big)` - - Where :math:`A_{i}` and :math:`B_{j}` are FCHL representation vectors. - K is calculated analytically using an OpenMP parallel Fortran routine. - Note, that this kernel will ONLY work with FCHL representations as input. - - :param A: Array of FCHL representation - shape=(N, maxsize, 5, size). - :type A: numpy array - :param sigma: List of kernel-widths. - :type sigma: list - - :param two_body_scaling: Weight for 2-body terms. - :type two_body_scaling: float - :param three_body_scaling: Weight for 3-body terms. - :type three_body_scaling: float - - :param two_body_width: Gaussian width for 2-body terms - :type two_body_width: float - :param three_body_width: Gaussian width for 3-body terms. - :type three_body_width: float - - :param two_body_power: Powerlaw for :math:`r^{-n}` 2-body terms. - :type two_body_power: float - :param three_body_power: Powerlaw for Axilrod-Teller-Muto 3-body term - :type three_body_power: float - - :param cut_start: The fraction of the cut-off radius at which cut-off damping start. - :type cut_start: float - :param cut_distance: Cut-off radius. (default=5 angstrom) - :type cut_distance: float - - :param fourier_order: 3-body Fourier-expansion truncation order. - :type fourier_order: integer - :param alchemy: Type of alchemical interpolation ``"periodic-table"`` or ``"off"`` are possible options. Disabling alchemical interpolation can yield dramatic speedups. - :type alchemy: string - - :param alchemy_period_width: Gaussian width along periods (columns) in the periodic table. - :type alchemy_period_width: float - :param alchemy_group_width: Gaussian width along groups (rows) in the periodic table. - :type alchemy_group_width: float - - :return: Array of FCHL kernel matrices matrix - shape=(n_sigmas, N, M), - :rtype: numpy array - """ - - assert len(A.shape) == 3 - - na1 = A.shape[0] - - neighbors1 = np.zeros((na1), dtype=np.int32) - - for i, x in enumerate(A): - neighbors1[i] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - kernel_idx, kernel_parameters, nsigmas = get_kernel_parameters(kernel, kernel_args) - - return fget_atomic_symmetric_kernels_fchl(A, neighbors1, na1, nsigmas, \ - three_body_width, two_body_width, cut_start, cut_distance, fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, three_body_power, kernel_idx, kernel_parameters) - - -def get_local_full_kernels(A, B, dx=0.005, \ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - - 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 - - - N1 = np.zeros((nm1),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm2): - N2[a] = len(np.where(B[a,0,0,0,:,1,0] > 0.0001)[0]) - - neighbors2 = np.zeros((nm2, 3, 2, atoms_max, atoms_max), dtype=np.int32) - - for m in range(nm2): - ni = N2[m] - for xyz in range(3): - for pm in range(2): - for i in range(ni): - for a, x in enumerate(B[m,xyz,pm,i,:ni]): - neighbors2[m,xyz,pm,i,a] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - naq2 = np.sum(N2) * 3 - - kernel_idx, kernel_parameters, nsigmas = get_kernel_parameters(kernel, kernel_args) - - return fget_local_full_kernels_fchl(A, B, N1, N2, neighbors1, neighbors2, \ - nm1, nm2, naq2, nsigmas, three_body_width, two_body_width, cut_start, cut_distance, \ - fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, three_body_power, dx, kernel_idx, kernel_parameters) - - -def get_local_gradient_kernels(A, B, dx=0.005, - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - - 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 - - - N1 = np.zeros((nm1),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm2): - N2[a] = len(np.where(B[a,0,0,0,:,1,0] > 0.0001)[0]) - - neighbors2 = np.zeros((nm2, 3, 2, atoms_max, atoms_max), dtype=np.int32) - - for m in range(nm2): - ni = N2[m] - for xyz in range(3): - for pm in range(2): - for i in range(ni): - for a, x in enumerate(B[m,xyz,pm,i,:ni]): - neighbors2[m,xyz,pm,i,a] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, nsigmas = get_kernel_parameters(kernel, kernel_args) - - naq2 = np.sum(N2) * 3 - - return fget_local_gradient_kernels_fchl(A, B, N1, N2, neighbors1, neighbors2, \ - nm1, nm2, naq2, nsigmas, three_body_width, two_body_width, cut_start, cut_distance, \ - fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, three_body_power, dx, kernel_idx, kernel_parameters) - - -def get_local_hessian_kernels(A, B, dx=0.005, - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - - 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 - - N1 = np.zeros((nm1),dtype=np.int32) - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,0,0,0,:,1,0] > 0.0001)[0]) - - for a in range(nm2): - N2[a] = len(np.where(B[a,0,0,0,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, 3, 2, atoms_max, atoms_max), dtype=np.int32) - neighbors2 = np.zeros((nm2, 3, 2, atoms_max, atoms_max), dtype=np.int32) - - for m in range(nm1): - ni = N1[m] - for xyz in range(3): - for pm in range(2): - for i in range(ni): - for a, x in enumerate(A[m,xyz,pm,i,:ni]): - neighbors1[m,xyz,pm,i,a] = len(np.where(x[0]< cut_distance)[0]) - - for m in range(nm2): - ni = N2[m] - for xyz in range(3): - for pm in range(2): - for i in range(ni): - for a, x in enumerate(B[m,xyz,pm,i,:ni]): - neighbors2[m,xyz,pm,i,a] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, nsigmas = get_kernel_parameters(kernel, kernel_args) - - - naq1 = np.sum(N1) * 3 - naq2 = np.sum(N2) * 3 - - # print naq1, naq2, nsigmas - return fget_local_hessian_kernels_fchl(A, B, N1, N2, neighbors1, neighbors2, \ - nm1, nm2, naq1, naq2, nsigmas, three_body_width, two_body_width, cut_start, cut_distance, \ - fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, three_body_power, dx, kernel_idx, kernel_parameters) - - -def get_local_symmetric_hessian_kernels(A, dx=0.005, - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - - 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 - - neighbors_max = A.shape[6] - - N1 = np.zeros((nm1),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,0,0,0,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, 3, 2, atoms_max, atoms_max), dtype=np.int32) - - for m in range(nm1): - ni = N1[m] - for xyz in range(3): - for pm in range(2): - for i in range(ni): - for a, x in enumerate(A[m,xyz,pm,i,:ni]): - neighbors1[m,xyz,pm,i,a] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, nsigmas = get_kernel_parameters(kernel, kernel_args) - - - naq1 = np.sum(N1) * 3 - - return fget_local_symmetric_hessian_kernels_fchl(A, N1, neighbors1, nm1, naq1, nsigmas, \ - three_body_width, two_body_width, cut_start, cut_distance, \ - fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, three_body_power, dx, kernel_idx, kernel_parameters) - - -def get_local_invariant_alphas(A, B, F, energy=None, dx=0.005, regularization=1e-7, - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - 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 - - - N1 = np.zeros((nm1),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm2): - N2[a] = len(np.where(B[a,0,0,0,:,1,0] > 0.0001)[0]) - - neighbors2 = np.zeros((nm2, 3, 2, atoms_max, atoms_max), dtype=np.int32) - - for m in range(nm2): - ni = N2[m] - for xyz in range(3): - for pm in range(2): - for i in range(ni): - for a, x in enumerate(B[m,xyz,pm,i,:ni]): - neighbors2[m,xyz,pm,i,a] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, nsigmas = get_kernel_parameters(kernel, kernel_args) - - na1 = np.sum(N1) - naq2 = np.sum(N2) * 3 - - E = np.zeros((nm1)) - if energy is not None: - E = energy - - return fget_local_invariant_alphas_fchl(A, B, F, E, N1, N2, neighbors1, neighbors2, \ - nm1, nm2, na1, nsigmas, three_body_width, two_body_width, cut_start, cut_distance, \ - fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, - three_body_power, dx, kernel_idx, kernel_parameters, regularization) - - -def get_atomic_gradient_kernels(A, B, dx = 0.005, - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - - 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 - - - N1 = np.zeros((nm1),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm2): - N2[a] = len(np.where(B[a,0,0,0,:,1,0] > 0.0001)[0]) - - neighbors2 = np.zeros((nm2, 3, 2, atoms_max, atoms_max), dtype=np.int32) - - for m in range(nm2): - ni = N2[m] - for xyz in range(3): - for pm in range(2): - for i in range(ni): - for a, x in enumerate(B[m,xyz,pm,i,:ni]): - neighbors2[m,xyz,pm,i,a] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, nsigmas = get_kernel_parameters(kernel, kernel_args) - - - na1 = np.sum(N1) - naq2 = np.sum(N2) * 3 - - return fget_atomic_gradient_kernels_fchl(A, B, N1, N2, neighbors1, neighbors2, \ - nm1, nm2, na1, naq2, nsigmas, three_body_width, two_body_width, cut_start, cut_distance, \ - fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, - three_body_power, dx, kernel_idx, kernel_parameters) - - -def get_local_atomic_kernels(A, B, \ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - """ Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: - - :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\sigma^2} \\big)` - - Where :math:`A_{i}` and :math:`B_{j}` are FCHL representation vectors. - K is calculated analytically using an OpenMP parallel Fortran routine. - Note, that this kernel will ONLY work with FCHL representations as input. - - :param A: Array of FCHL representation - shape=(N, maxsize, 5, maxneighbors). - :type A: numpy array - :param B: Array of FCHL representation - shape=(M, maxsize, 5, maxneighbors). - :type B: numpy array - :param sigma: List of kernel-widths. - :type sigma: list - :param t_width: Gaussian width for the angular (theta) terms. - :type t_width: float - :param d_width: Gaussian width for the distance terms. - :type d_width: float - :param cut_start: The fraction of the cut-off radius at which cut-off damping start - :type cut_start: float - :param cut_distance: Cut-off radius. - :type cut_distance: float - :param r_width: Gaussian width along rows in the periodic table. - :type r_width: float - :param c_width: Gaussian width along columns in the periodic table. - :type c_width: float - :param order: Fourier-expansion truncation order. - :type order: integer - :param scale_distance: Weight for distance-dependent terms. - :type scale_distance: float - :param scale_angular: Weight for angle-dependent terms. - :type scale_angular: float - - :return: Array of FCHL kernel matrices matrix - shape=(n_sigmas, N, M), - :rtype: numpy array - """ - - 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" - - nm1 = A.shape[0] - nm2 = B.shape[0] - - N1 = np.zeros((nm1),dtype=np.int32) - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - for a in range(nm2): - N2[a] = len(np.where(B[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - neighbors2 = np.zeros((nm2, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - for a, representation in enumerate(B): - ni = N2[a] - for i, x in enumerate(representation[:ni]): - neighbors2[a,i] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, nsigmas = get_kernel_parameters(kernel, kernel_args) - - na1 = np.sum(N1) - - return fget_local_atomic_kernels_fchl(A, B, N1, N2, neighbors1, neighbors2, \ - nm1, nm2, na1, nsigmas, three_body_width, two_body_width, cut_start, cut_distance, \ - fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, \ - three_body_power, kernel_idx, kernel_parameters) - - -def get_local_kernels_ef(A, B, df=1e-5, ef_scaling=1.0,\ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - """ Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: - - :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\sigma^2} \\big)` - - Where :math:`A_{i}` and :math:`B_{j}` are FCHL representation vectors. - K is calculated analytically using an OpenMP parallel Fortran routine. - Note, that this kernel will ONLY work with FCHL representations as input. - - :param A: Array of FCHL representation - shape=(N, maxsize, 5, maxneighbors). - :type A: numpy array - :param B: Array of FCHL representation - shape=(M, maxsize, 5, maxneighbors). - :type B: numpy array - - :param two_body_scaling: Weight for 2-body terms. - :type two_body_scaling: float - :param three_body_scaling: Weight for 3-body terms. - :type three_body_scaling: float - - :param two_body_width: Gaussian width for 2-body terms - :type two_body_width: float - :param three_body_width: Gaussian width for 3-body terms. - :type three_body_width: float - - :param two_body_power: Powerlaw for :math:`r^{-n}` 2-body terms. - :type two_body_power: float - :param three_body_power: Powerlaw for Axilrod-Teller-Muto 3-body term - :type three_body_power: float - - :param cut_start: The fraction of the cut-off radius at which cut-off damping start. - :type cut_start: float - :param cut_distance: Cut-off radius. (default=5 angstrom) - :type cut_distance: float - - :param fourier_order: 3-body Fourier-expansion truncation order. - :type fourier_order: integer - :param alchemy: Type of alchemical interpolation ``"periodic-table"`` or ``"off"`` are possible options. Disabling alchemical interpolation can yield dramatic speedups. - :type alchemy: string - - :param alchemy_period_width: Gaussian width along periods (columns) in the periodic table. - :type alchemy_period_width: float - :param alchemy_group_width: Gaussian width along groups (rows) in the periodic table. - :type alchemy_group_width: float - - :return: Array of FCHL kernel matrices matrix - shape=(n_sigmas, N, M), - :rtype: numpy array - """ - - 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" - - - nm1 = A.shape[0] - nm2 = B.shape[0] - - N1 = np.zeros((nm1),dtype=np.int32) - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - for a in range(nm2): - N2[a] = len(np.where(B[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - neighbors2 = np.zeros((nm2, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - for a, representation in enumerate(B): - ni = N2[a] - for i, x in enumerate(representation[:ni]): - neighbors2[a,i] = len(np.where(x[0]< cut_distance)[0]) - - # charge1 = np.zeros((nm1, atoms_max)) - # charge2 = np.zeros((nm2, atoms_max)) - # - # for a, representation in enumerate(A): - # ni = N1[a] - # for i, x in enumerate(representation[:ni]): - # charge1[a,i] = Q1[a][i] - # - # for a, representation in enumerate(B): - # ni = N2[a] - # for i, x in enumerate(representation[:ni]): - # charge2[a,i] = Q2[a][i] - - # print(charge1) - # print(charge2) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, n_kernels = get_kernel_parameters(kernel, kernel_args) - - return fget_kernels_fchl_ef(A, B, N1, N2, neighbors1, neighbors2, nm1, nm2, n_kernels, \ - three_body_width, two_body_width, cut_start, cut_distance, fourier_order, pd, two_body_scaling, - three_body_scaling, doalchemy, two_body_power, three_body_power, ef_scaling, df, kernel_idx, kernel_parameters) - -def get_local_kernels_ef_deriv(A, B, df=1e-5, ef_scaling=1.0,\ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - """ Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: - - :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\sigma^2} \\big)` - - Where :math:`A_{i}` and :math:`B_{j}` are FCHL representation vectors. - K is calculated analytically using an OpenMP parallel Fortran routine. - Note, that this kernel will ONLY work with FCHL representations as input. - - :param A: Array of FCHL representation - shape=(N, maxsize, 5, maxneighbors). - :type A: numpy array - :param B: Array of FCHL representation - shape=(M, maxsize, 5, maxneighbors). - :type B: numpy array - - :param two_body_scaling: Weight for 2-body terms. - :type two_body_scaling: float - :param three_body_scaling: Weight for 3-body terms. - :type three_body_scaling: float - - :param two_body_width: Gaussian width for 2-body terms - :type two_body_width: float - :param three_body_width: Gaussian width for 3-body terms. - :type three_body_width: float - - :param two_body_power: Powerlaw for :math:`r^{-n}` 2-body terms. - :type two_body_power: float - :param three_body_power: Powerlaw for Axilrod-Teller-Muto 3-body term - :type three_body_power: float - - :param cut_start: The fraction of the cut-off radius at which cut-off damping start. - :type cut_start: float - :param cut_distance: Cut-off radius. (default=5 angstrom) - :type cut_distance: float - - :param fourier_order: 3-body Fourier-expansion truncation order. - :type fourier_order: integer - :param alchemy: Type of alchemical interpolation ``"periodic-table"`` or ``"off"`` are possible options. Disabling alchemical interpolation can yield dramatic speedups. - :type alchemy: string - - :param alchemy_period_width: Gaussian width along periods (columns) in the periodic table. - :type alchemy_period_width: float - :param alchemy_group_width: Gaussian width along groups (rows) in the periodic table. - :type alchemy_group_width: float - - :return: Array of FCHL kernel matrices matrix - shape=(n_sigmas, N, M), - :rtype: numpy array - """ - - 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" - - - nm1 = A.shape[0] - nm2 = B.shape[0] - - N1 = np.zeros((nm1),dtype=np.int32) - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - for a in range(nm2): - N2[a] = len(np.where(B[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - neighbors2 = np.zeros((nm2, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - for a, representation in enumerate(B): - ni = N2[a] - for i, x in enumerate(representation[:ni]): - neighbors2[a,i] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, n_kernels = get_kernel_parameters(kernel, kernel_args) - - na1 = np.sum(N1) - - return fget_kernels_fchl_ef_deriv(A, B, N1, N2, neighbors1, neighbors2, nm1, nm2, na1, n_kernels, \ - three_body_width, two_body_width, cut_start, cut_distance, fourier_order, pd, two_body_scaling, - three_body_scaling, doalchemy, two_body_power, three_body_power, ef_scaling, df, kernel_idx, kernel_parameters) - - -def get_gp_kernels_ef(A, B, fields=None, df=1e-5, ef_scaling=1.0,\ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - """ Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: - - :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\sigma^2} \\big)` - - Where :math:`A_{i}` and :math:`B_{j}` are FCHL representation vectors. - K is calculated analytically using an OpenMP parallel Fortran routine. - Note, that this kernel will ONLY work with FCHL representations as input. - - :param A: Array of FCHL representation - shape=(N, maxsize, 5, maxneighbors). - :type A: numpy array - :param B: Array of FCHL representation - shape=(M, maxsize, 5, maxneighbors). - :type B: numpy array - - :param two_body_scaling: Weight for 2-body terms. - :type two_body_scaling: float - :param three_body_scaling: Weight for 3-body terms. - :type three_body_scaling: float - - :param two_body_width: Gaussian width for 2-body terms - :type two_body_width: float - :param three_body_width: Gaussian width for 3-body terms. - :type three_body_width: float - - :param two_body_power: Powerlaw for :math:`r^{-n}` 2-body terms. - :type two_body_power: float - :param three_body_power: Powerlaw for Axilrod-Teller-Muto 3-body term - :type three_body_power: float - - :param cut_start: The fraction of the cut-off radius at which cut-off damping start. - :type cut_start: float - :param cut_distance: Cut-off radius. (default=5 angstrom) - :type cut_distance: float - - :param fourier_order: 3-body Fourier-expansion truncation order. - :type fourier_order: integer - :param alchemy: Type of alchemical interpolation ``"periodic-table"`` or ``"off"`` are possible options. Disabling alchemical interpolation can yield dramatic speedups. - :type alchemy: string - - :param alchemy_period_width: Gaussian width along periods (columns) in the periodic table. - :type alchemy_period_width: float - :param alchemy_group_width: Gaussian width along groups (rows) in the periodic table. - :type alchemy_group_width: float - - :return: Array of FCHL kernel matrices matrix - shape=(n_sigmas, N, M), - :rtype: numpy array - """ - - 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" - - - nm1 = A.shape[0] - nm2 = B.shape[0] - - N1 = np.zeros((nm1),dtype=np.int32) - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - for a in range(nm2): - N2[a] = len(np.where(B[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - neighbors2 = np.zeros((nm2, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - for a, representation in enumerate(B): - ni = N2[a] - for i, x in enumerate(representation[:ni]): - neighbors2[a,i] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, n_kernels = get_kernel_parameters(kernel, kernel_args) - - F1 = np.zeros((nm1,3)) - F2 = np.zeros((nm2,3)) - - if (fields is not None): - - F1 = np.array(fields[0]) - F2 = np.array(fields[1]) - - return fget_gp_kernels_fchl_ef(A, B, F1, F2, N1, N2, neighbors1, neighbors2, nm1, nm2, n_kernels, \ - three_body_width, two_body_width, cut_start, cut_distance, fourier_order, pd, two_body_scaling, - three_body_scaling, doalchemy, two_body_power, three_body_power, ef_scaling, df, kernel_idx, kernel_parameters) - - -def get_local_kernels_pol(A, B, df=1e-5, ef_scaling=1.0,\ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - """ Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: - - :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\sigma^2} \\big)` - - Where :math:`A_{i}` and :math:`B_{j}` are FCHL representation vectors. - K is calculated analytically using an OpenMP parallel Fortran routine. - Note, that this kernel will ONLY work with FCHL representations as input. - - :param A: Array of FCHL representation - shape=(N, maxsize, 5, maxneighbors). - :type A: numpy array - :param B: Array of FCHL representation - shape=(M, maxsize, 5, maxneighbors). - :type B: numpy array - - :param two_body_scaling: Weight for 2-body terms. - :type two_body_scaling: float - :param three_body_scaling: Weight for 3-body terms. - :type three_body_scaling: float - - :param two_body_width: Gaussian width for 2-body terms - :type two_body_width: float - :param three_body_width: Gaussian width for 3-body terms. - :type three_body_width: float - - :param two_body_power: Powerlaw for :math:`r^{-n}` 2-body terms. - :type two_body_power: float - :param three_body_power: Powerlaw for Axilrod-Teller-Muto 3-body term - :type three_body_power: float - - :param cut_start: The fraction of the cut-off radius at which cut-off damping start. - :type cut_start: float - :param cut_distance: Cut-off radius. (default=5 angstrom) - :type cut_distance: float - - :param fourier_order: 3-body Fourier-expansion truncation order. - :type fourier_order: integer - :param alchemy: Type of alchemical interpolation ``"periodic-table"`` or ``"off"`` are possible options. Disabling alchemical interpolation can yield dramatic speedups. - :type alchemy: string - - :param alchemy_period_width: Gaussian width along periods (columns) in the periodic table. - :type alchemy_period_width: float - :param alchemy_group_width: Gaussian width along groups (rows) in the periodic table. - :type alchemy_group_width: float - - :return: Array of FCHL kernel matrices matrix - shape=(n_sigmas, N, M), - :rtype: numpy array - """ - - 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" - - - nm1 = A.shape[0] - nm2 = B.shape[0] - - N1 = np.zeros((nm1),dtype=np.int32) - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - for a in range(nm2): - N2[a] = len(np.where(B[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - neighbors2 = np.zeros((nm2, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - for a, representation in enumerate(B): - ni = N2[a] - for i, x in enumerate(representation[:ni]): - neighbors2[a,i] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, n_kernels = get_kernel_parameters(kernel, kernel_args) - - na1 = np.sum(N1) - - return fget_kernels_fchl_ef_2ndderiv(A, B, N1, N2, neighbors1, neighbors2, nm1, nm2, na1, n_kernels, \ - three_body_width, two_body_width, cut_start, cut_distance, fourier_order, pd, two_body_scaling, - three_body_scaling, doalchemy, two_body_power, three_body_power, ef_scaling, df, kernel_idx, kernel_parameters) - - -def get_kernels_ef_field(A, B, fields=None, df=1e-5, ef_scaling=1.0,\ - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - - 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" - - - nm1 = A.shape[0] - nm2 = B.shape[0] - - N1 = np.zeros((nm1),dtype=np.int32) - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - for a in range(nm2): - N2[a] = len(np.where(B[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - neighbors2 = np.zeros((nm2, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - for a, representation in enumerate(B): - ni = N2[a] - for i, x in enumerate(representation[:ni]): - neighbors2[a,i] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, n_kernels = get_kernel_parameters(kernel, kernel_args) - - F2 = np.zeros((nm2,3)) - - if (fields is not None): - - F2 = np.array(fields) - - na1 = np.sum(N1) - - return fget_kernels_fchl_ef_field(A, B, F2, N1, N2, neighbors1, neighbors2, nm1, nm2, n_kernels, na1, \ - three_body_width, two_body_width, cut_start, cut_distance, fourier_order, pd, two_body_scaling, - three_body_scaling, doalchemy, two_body_power, three_body_power, ef_scaling, kernel_idx, kernel_parameters) - - -def get_smooth_atomic_gradient_kernels(A, B, dx = 0.005, - two_body_scaling=np.sqrt(8), three_body_scaling=1.6, - two_body_width=0.2, three_body_width=np.pi, - two_body_power=4.0, three_body_power=2.0, - cut_start=1.0, cut_distance=5.0, - fourier_order=1, alchemy="periodic-table", - alchemy_period_width=1.6, alchemy_group_width=1.6, - kernel="gaussian", kernel_args=None): - - 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 - - - N1 = np.zeros((nm1),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(A[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - - for a, representation in enumerate(A): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0]< cut_distance)[0]) - - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm2): - N2[a] = len(np.where(B[a,0,0,0,:,1,0] > 0.0001)[0]) - - neighbors2 = np.zeros((nm2, 3, 5, atoms_max, atoms_max), dtype=np.int32) - - for m in range(nm2): - ni = N2[m] - for xyz in range(3): - for pm in range(5): - for i in range(ni): - for a, x in enumerate(B[m,xyz,pm,i,:ni]): - neighbors2[m,xyz,pm,i,a] = len(np.where(x[0]< cut_distance)[0]) - - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=alchemy_group_width, c_width=alchemy_period_width) - - kernel_idx, kernel_parameters, nsigmas = get_kernel_parameters(kernel, kernel_args) - - - na1 = np.sum(N1) - naq2 = np.sum(N2) * 3 - - return fget_smooth_atomic_gradient_kernels_fchl(A, B, N1, N2, neighbors1, neighbors2, \ - nm1, nm2, na1, naq2, nsigmas, three_body_width, two_body_width, cut_start, cut_distance, \ - fourier_order, pd, two_body_scaling, three_body_scaling, doalchemy, two_body_power, - three_body_power, dx, kernel_idx, kernel_parameters) diff --git a/qml/kernels/wrappers.py b/qml/kernels/wrappers.py index 35688567b..df057e033 100644 --- a/qml/kernels/wrappers.py +++ b/qml/kernels/wrappers.py @@ -31,7 +31,6 @@ from ..arad import get_local_kernels_arad from ..arad import get_local_symmetric_kernels_arad - def get_atomic_kernels_laplacian(mols1, mols2, sigmas): n1 = np.array([mol.natoms for mol in mols1], dtype=np.int32) @@ -87,31 +86,6 @@ def get_atomic_kernels_laplacian_symmetric(mols, sigmas): return fget_vector_kernels_laplacian(x1, n, sigmas, nm, nsigmas) - -def get_atomic_kernels_gaussian(mols, sigmas): - - n = np.array([mol.natoms for mol in mols], dtype=np.int32) - - max_atoms = np.max(n) - - nm = n.size - - cmat_size = mols[0].representation.shape[1] - - x1 = np.zeros((nm, max_atoms, cmat_size), dtype=np.float64, order="F") - - for imol in range(nm1): - x[imol,:n[imol],:cmat_size] = mols[imol].representation - - # Reorder for Fortran speed - x = np.swapaxes(x, 0, 2) - - sigmas = np.array(sigmas, dtype=np.float64) - nsigmas = sigmas.size - - return fget_vector_kernels_gaussian(x, n, sigmas, nm, nsigmas) - - def arad_local_kernels(mols1, mols2, sigmas, width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5): @@ -128,7 +102,6 @@ def arad_local_kernels(mols1, mols2, sigmas, return K - def arad_local_symmetric_kernels(mols1, sigmas, width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5): @@ -174,7 +147,6 @@ def get_atomic_kernels_laplacian(mols1, mols2, sigmas): return fget_vector_kernels_laplacian(x1, x2, n1, n2, sigmas, nm1, nm2, nsigmas) - def get_atomic_kernels_gaussian(mols1, mols2, sigmas): n1 = np.array([mol.natoms for mol in mols1], dtype=np.int32) @@ -206,3 +178,26 @@ def get_atomic_kernels_gaussian(mols1, mols2, sigmas): return fget_vector_kernels_gaussian(x1, x2, n1, n2, sigmas, nm1, nm2, nsigmas) + +def get_atomic_kernels_gaussian_symmetric(mols, sigmas): + + n = np.array([mol.natoms for mol in mols], dtype=np.int32) + + max_atoms = np.max(n) + + nm = n.size + + cmat_size = mols[0].representation.shape[1] + + x1 = np.zeros((nm, max_atoms, cmat_size), dtype=np.float64, order="F") + + for imol in range(nm1): + x[imol,:n[imol],:cmat_size] = mols[imol].representation + + # Reorder for Fortran speed + x = np.swapaxes(x, 0, 2) + + sigmas = np.array(sigmas, dtype=np.float64) + nsigmas = sigmas.size + + return fget_vector_kernels_gaussian_symmetric(x, n, sigmas, nm, nsigmas) diff --git a/qml/qmlearn/kernels.py b/qml/qmlearn/kernels.py index a1157af2d..bd13bb14a 100644 --- a/qml/qmlearn/kernels.py +++ b/qml/qmlearn/kernels.py @@ -28,19 +28,15 @@ from .data import Data -from ..kernels.fkernels import fgaussian_kernel -from ..kernels.fkernels import fgaussian_kernel_symmetric from ..kernels.fkernels import fget_vector_kernels_gaussian from ..kernels.fkernels import fget_vector_kernels_gaussian_symmetric -from ..kernels.fkernels import flaplacian_kernel -from ..kernels.fkernels import flaplacian_kernel_symmetric from ..kernels.fkernels import fget_vector_kernels_laplacian from ..kernels.fkernels import fget_vector_kernels_laplacian_symmetric -from ..fchl.ffchl_module import fget_kernels_fchl -from ..fchl.ffchl_module import fget_symmetric_kernels_fchl -from ..fchl.ffchl_module import fget_global_kernels_fchl -from ..fchl.ffchl_module import fget_global_symmetric_kernels_fchl +from ..fchl import get_global_symmetric_kernels +from ..fchl import get_global_kernels +from ..fchl import get_local_symmetric_kernels +from ..fchl import get_local_kernels from ..utils.alchemy import get_alchemy from ..utils import get_unique @@ -344,24 +340,11 @@ def generate(self, X, Y=None, representation_type='molecular'): return kernel def _generate_molecular(self, X, Y=None): - - X = np.asarray(X) - - # Note: Transposed for Fortran - n = X.shape[0] - if Y is None or X is Y: # Do symmetric matrix - K = np.empty((n, n), order='F') - fgaussian_kernel_symmetric(X.T, n, K, self.sigma) + return gaussian_kernel_symmetric(X, self.sigma) else: - Y = np.asarray(Y) - # Do asymmetric matrix - m = Y.shape[0] - K = np.empty((n, m), order='F') - fgaussian_kernel(X.T, n, Y.T, m, K, self.sigma) - - return K + return gaussian_kernel(X, Y, self.sigma) def _generate_atomic(self, X, Y=None): @@ -371,6 +354,8 @@ def _generate_atomic(self, X, Y=None): nm1 = n1.size + # This is a bit of a mess to get kernels without + # alchemy working for i in range(len(X)): rep_size = X[0].shape[1] if rep_size == 0: @@ -383,6 +368,8 @@ def _generate_atomic(self, X, Y=None): x1 = np.zeros((nm1, max1, rep_size), dtype=np.float64, order="F") + # This is a bit of a mess to get kernels without + # alchemy working for i in range(nm1): if X[i].shape[1] == 0: n1[i] = 0 @@ -406,6 +393,8 @@ def _generate_atomic(self, X, Y=None): max2 = np.max(n2) nm2 = n2.size x2 = np.zeros((nm2, max2, rep_size), dtype=np.float64, order="F") + # This is a bit of a mess to get kernels without + # alchemy working for i in range(nm2): if Y[i].shape[1] == 0: n2[i] = 0 @@ -534,24 +523,11 @@ def generate(self, X, Y=None, representation_type='molecular'): return kernel def _generate_molecular(self, X, Y=None): - - X = np.asarray(X) - - # Note: Transposed for Fortran - n = X.shape[0] - if Y is None or X is Y: # Do symmetric matrix - K = np.empty((n, n), order='F') - flaplacian_kernel_symmetric(X.T, n, K, self.sigma) + return laplacian_kernel_symmetric(X, self.sigma) else: - Y = np.asarray(Y) - # Do asymmetric matrix - m = Y.shape[0] - K = np.empty((n, m), order='F') - flaplacian_kernel(X.T, n, Y.T, m, K, self.sigma) - - return K + return laplacian_kernel(X, Y, self.sigma) def _generate_atomic(self, X, Y=None): @@ -561,6 +537,8 @@ def _generate_atomic(self, X, Y=None): nm1 = n1.size + # This is a bit of a mess to get kernels without + # alchemy working for i in range(len(X)): rep_size = X[0].shape[1] if rep_size == 0: @@ -573,6 +551,8 @@ def _generate_atomic(self, X, Y=None): x1 = np.zeros((nm1, max1, rep_size), dtype=np.float64, order="F") + # This is a bit of a mess to get kernels without + # alchemy working for i in range(nm1): if X[i].shape[1] == 0: n1[i] = 0 @@ -599,6 +579,8 @@ def _generate_atomic(self, X, Y=None): max2 = np.max(n2) nm2 = n2.size x2 = np.zeros((nm2, max2, rep_size), dtype=np.float64, order="F") + # This is a bit of a mess to get kernels without + # alchemy working for i in range(nm2): if Y[i].shape[1] == 0: n2[i] = 0 @@ -686,7 +668,7 @@ def _quick_estimate_sigma(self, X, sigma_init=1, count=1): # Generate kernel for a random subset indices = np.random.choice(np.arange(len(X)), size=n, replace=False) - kernel = self.generate(X)#[indices]) + kernel = self.generate([X[i] for i in indices]) if not self.local: # min smallest kernel element @@ -736,102 +718,51 @@ def generate(self, X, Y=None): def _generate_molecular(self, X, Y=None): - atoms_max = X.shape[1] - neighbors_max = X.shape[3] - nm1 = X.shape[0] - N1 = np.zeros((nm1),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(X[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - - for a, representation in enumerate(X): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0] < self.cutoff)[0]) - if self.alchemy: alchemy = 'periodic-table' else: alchemy = 'off' - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=self.alchemy_group_width, c_width=self.alchemy_period_width) - if Y is None or X is Y: - # Do symmetric kernel - return fget_global_symmetric_kernels_fchl(X, N1, neighbors1, [self.sigma], - nm1, 1, self.three_body_width, self.two_body_width, self.damping_start, - self.cutoff, self.fourier_order, pd, self.two_body_scaling, self.three_body_scaling, - doalchemy, self.two_body_power, self.three_body_power)[0] + return get_global_symmetric_kernels(X, two_body_scaling=self.two_body_scaling, + three_body_scaling=self.three_body_scaling, two_body_width=self.two_body_width, + three_body_width=self.three_body_width, two_body_power=self.two_body_power, + three_body_power=self.three_body_power, cut_start=self.damping_start, + cut_distance=self.cutoff, fourier_order=self.fourier_order, alchemy=alchemy, + alchemy_period_width=self.alchemy_period_width, alchemy_group_width=self.alchemy_group_width, + kernel_args={'sigma': [self.sigma]})[0] else: - # Do asymmetric kernel - nm2 = Y.shape[0] - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm2): - N2[a] = len(np.where(Y[a,:,1,0] > 0.0001)[0]) - - neighbors2 = np.zeros((nm2, atoms_max), dtype=np.int32) - - for a, representation in enumerate(Y): - ni = N2[a] - for i, x in enumerate(representation[:ni]): - neighbors2[a,i] = len(np.where(x[0] < self.cutoff)[0]) - - return fget_global_kernels_fchl(X, Y, N1, N2, neighbors1, neighbors2, [self.sigma], - nm1, nm2, 1, self.three_body_width, self.two_body_width, self.damping_start, - self.cutoff, self.fourier_order, pd, self.two_body_scaling, self.three_body_scaling, - doalchemy, self.two_body_power, self.three_body_power)[0] + return get_global_kernels(X, Y, two_body_scaling=self.two_body_scaling, + three_body_scaling=self.three_body_scaling, two_body_width=self.two_body_width, + three_body_width=self.three_body_width, two_body_power=self.two_body_power, + three_body_power=self.three_body_power, cut_start=self.damping_start, + cut_distance=self.cutoff, fourier_order=self.fourier_order, alchemy=alchemy, + alchemy_period_width=self.alchemy_period_width, alchemy_group_width=self.alchemy_group_width, + kernel_args={'sigma': [self.sigma]})[0] def _generate_atomic(self, X, Y=None): - - atoms_max = X.shape[1] - neighbors_max = X.shape[3] - - nm1 = X.shape[0] - N1 = np.zeros((nm1),dtype=np.int32) - - for a in range(nm1): - N1[a] = len(np.where(X[a,:,1,0] > 0.0001)[0]) - - neighbors1 = np.zeros((nm1, atoms_max), dtype=np.int32) - - for a, representation in enumerate(X): - ni = N1[a] - for i, x in enumerate(representation[:ni]): - neighbors1[a,i] = len(np.where(x[0] < self.cutoff)[0]) - if self.alchemy: alchemy = 'periodic-table' else: alchemy = 'off' - doalchemy, pd = get_alchemy(alchemy, emax=100, r_width=self.alchemy_group_width, c_width=self.alchemy_period_width) if Y is None or X is Y: - return fget_symmetric_kernels_fchl(X, N1, neighbors1, [self.sigma], - nm1, 1, self.three_body_width, self.two_body_width, self.damping_start, - self.cutoff, self.fourier_order, pd, self.two_body_scaling, self.three_body_scaling, - doalchemy, self.two_body_power, self.three_body_power)[0] + return get_local_symmetric_kernels(X, two_body_scaling=self.two_body_scaling, + three_body_scaling=self.three_body_scaling, two_body_width=self.two_body_width, + three_body_width=self.three_body_width, two_body_power=self.two_body_power, + three_body_power=self.three_body_power, cut_start=self.damping_start, + cut_distance=self.cutoff, fourier_order=self.fourier_order, alchemy=alchemy, + alchemy_period_width=self.alchemy_period_width, alchemy_group_width=self.alchemy_group_width, + kernel_args={'sigma': [self.sigma]})[0] else: - nm2 = Y.shape[0] - N2 = np.zeros((nm2),dtype=np.int32) - - for a in range(nm2): - N2[a] = len(np.where(Y[a,:,1,0] > 0.0001)[0]) - - neighbors2 = np.zeros((nm2, atoms_max), dtype=np.int32) - - for a, representation in enumerate(Y): - ni = N2[a] - for i, x in enumerate(representation[:ni]): - neighbors2[a,i] = len(np.where(x[0] < self.cutoff)[0]) - - return fget_kernels_fchl(X, Y, N1, N2, neighbors1, neighbors2, [self.sigma], - nm1, nm2, 1, self.three_body_width, self.two_body_width, self.damping_start, - self.cutoff, self.fourier_order, pd, self.two_body_scaling, self.three_body_scaling, - doalchemy, self.two_body_power, self.three_body_power)[0] + return get_local_kernels(X, Y, two_body_scaling=self.two_body_scaling, + three_body_scaling=self.three_body_scaling, two_body_width=self.two_body_width, + three_body_width=self.three_body_width, two_body_power=self.two_body_power, + three_body_power=self.three_body_power, cut_start=self.damping_start, + cut_distance=self.cutoff, fourier_order=self.fourier_order, alchemy=alchemy, + alchemy_period_width=self.alchemy_period_width, alchemy_group_width=self.alchemy_group_width, + kernel_args={'sigma': [self.sigma]})[0] def fit_transform(self, X, y=None): diff --git a/qml/qmlearn/representations.py b/qml/qmlearn/representations.py index fb9c13bde..510311892 100644 --- a/qml/qmlearn/representations.py +++ b/qml/qmlearn/representations.py @@ -468,7 +468,7 @@ def _transform(self, X, local): generate_slatm(xyz, charge, self.element_pairs, local=local, sigmas=[self.sigma2, self.sigma3], dgrids=[self.dgrid2, self.dgrid3], rcut=self.rcut, - alchemy=self.alchemy, rpower=-self.rpower))) + alchemy=self.alchemy, rpower=self.rpower))) data._representations = np.asarray(representations) @@ -483,7 +483,7 @@ class GlobalSLATM(_SLATM, _MolecularRepresentation): _representation_short_name = "slatm" def __init__(self, data=None, sigma2=0.05, sigma3=0.05, dgrid2=0.03, - dgrid3=0.03, rcut=4.8, alchemy=False, rpower=-6, elements='auto', + dgrid3=0.03, rcut=4.8, alchemy=False, rpower=6, elements='auto', element_pairs='auto'): """ Generate Spectrum of London and Axillrod-Teller-Muto potential (SLATM) representation. @@ -559,7 +559,7 @@ class AtomicSLATM(_SLATM, _AtomicRepresentation): _representation_short_name = "aslatm" def __init__(self, data=None, sigma2=0.05, sigma3=0.05, dgrid2=0.03, - dgrid3=0.03, rcut=4.8, alchemy=False, rpower=-6, elements='auto', + dgrid3=0.03, rcut=4.8, alchemy=False, rpower=6, elements='auto', element_pairs='auto'): """ Generate Spectrum of London and Axillrod-Teller-Muto potential (SLATM) representation. diff --git a/test/test_fchl_electric_field.py b/test/test_fchl_electric_field.py index 504fa0f97..b920385ad 100644 --- a/test/test_fchl_electric_field.py +++ b/test/test_fchl_electric_field.py @@ -470,5 +470,4 @@ def test_gaussian_process_field_dependent(): test_generate_representation() test_gaussian_process() test_gaussian_process_field_dependent() - test_explicit_electric_field()