(git:0de0cc2)
pao_ml_gaussprocess.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Gaussian Process implementation
10 !> \author Ole Schuett
11 ! **************************************************************************************************
13  USE kinds, ONLY: dp
14  USE pao_types, ONLY: pao_env_type,&
15  training_matrix_type
16 #include "./base/base_uses.f90"
17 
18  IMPLICIT NONE
19 
20  PRIVATE
21 
22  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_ml_gaussprocess'
23 
25 
26 CONTAINS
27 
28 ! **************************************************************************************************
29 !> \brief Builds the covariance matrix
30 !> \param pao ...
31 ! **************************************************************************************************
32  SUBROUTINE pao_ml_gp_train(pao)
33  TYPE(pao_env_type), POINTER :: pao
34 
35  INTEGER :: i, ikind, info, j, npoints
36  REAL(dp), DIMENSION(:), POINTER :: idescr, jdescr
37  TYPE(training_matrix_type), POINTER :: training_matrix
38 
39  ! TODO this could be parallelized over ranks
40  DO ikind = 1, SIZE(pao%ml_training_matrices)
41  training_matrix => pao%ml_training_matrices(ikind)
42  npoints = SIZE(training_matrix%inputs, 2) ! number of points
43  cpassert(SIZE(training_matrix%outputs, 2) == npoints)
44  IF (npoints == 0) cycle ! have no training data
45 
46  IF (pao%iw > 0) WRITE (pao%iw, *) "PAO|ML| Building covariance matrix for kind: ", &
47  trim(training_matrix%kindname), " from ", npoints, "training points."
48 
49  ! build covariance matrix
50  ALLOCATE (training_matrix%GP(npoints, npoints))
51  DO i = 1, npoints
52  DO j = i, npoints
53  idescr => training_matrix%inputs(:, i)
54  jdescr => training_matrix%inputs(:, j)
55  training_matrix%GP(i, j) = kernel(pao%gp_scale, idescr, jdescr)
56  training_matrix%GP(j, i) = training_matrix%GP(i, j)
57  END DO
58  END DO
59 
60  ! add noise of training data
61  DO i = 1, npoints
62  training_matrix%GP(i, i) = training_matrix%GP(i, i) + pao%gp_noise_var**2
63  END DO
64 
65  ! compute cholesky decomposition of covariance matrix
66  CALL dpotrf("U", npoints, training_matrix%GP, npoints, info)
67  cpassert(info == 0)
68  END DO
69 
70  END SUBROUTINE pao_ml_gp_train
71 
72 ! **************************************************************************************************
73 !> \brief Uses covariance matrix to make prediction
74 !> \param pao ...
75 !> \param ikind ...
76 !> \param descriptor ...
77 !> \param output ...
78 !> \param variance ...
79 ! **************************************************************************************************
80  SUBROUTINE pao_ml_gp_predict(pao, ikind, descriptor, output, variance)
81  TYPE(pao_env_type), POINTER :: pao
82  INTEGER, INTENT(IN) :: ikind
83  REAL(dp), DIMENSION(:), INTENT(IN) :: descriptor
84  REAL(dp), DIMENSION(:), INTENT(OUT) :: output
85  REAL(dp), INTENT(OUT) :: variance
86 
87  INTEGER :: i, info, npoints
88  REAL(dp), ALLOCATABLE, DIMENSION(:) :: cov, weights
89  TYPE(training_matrix_type), POINTER :: training_matrix
90 
91  training_matrix => pao%ml_training_matrices(ikind)
92  npoints = SIZE(training_matrix%outputs, 2)
93 
94  ! calculate covariances between descriptor and training-points
95  ALLOCATE (cov(npoints))
96  DO i = 1, npoints
97  cov(i) = kernel(pao%gp_scale, descriptor, training_matrix%inputs(:, i))
98  END DO
99 
100  ! calculate weights
101  ALLOCATE (weights(npoints))
102  weights(:) = cov(:)
103  CALL dpotrs("U", npoints, 1, training_matrix%GP, npoints, weights, npoints, info)
104  cpassert(info == 0)
105 
106  ! calculate predicted output
107  output = 0.0_dp
108  DO i = 1, npoints
109  output(:) = output + weights(i)*training_matrix%outputs(:, i)
110  END DO
111 
112  ! calculate prediction's variance
113  variance = kernel(pao%gp_scale, descriptor, descriptor) - dot_product(weights, cov)
114 
115  IF (variance < 0.0_dp) &
116  cpabort("PAO gaussian process found negative variance")
117 
118  DEALLOCATE (cov, weights)
119  END SUBROUTINE pao_ml_gp_predict
120 
121 ! **************************************************************************************************
122 !> \brief Calculate gradient of Gaussian process
123 !> \param pao ...
124 !> \param ikind ...
125 !> \param descriptor ...
126 !> \param outer_deriv ...
127 !> \param gradient ...
128 ! **************************************************************************************************
129  SUBROUTINE pao_ml_gp_gradient(pao, ikind, descriptor, outer_deriv, gradient)
130  TYPE(pao_env_type), POINTER :: pao
131  INTEGER, INTENT(IN) :: ikind
132  REAL(dp), DIMENSION(:), INTENT(IN), TARGET :: descriptor
133  REAL(dp), DIMENSION(:), INTENT(IN) :: outer_deriv
134  REAL(dp), DIMENSION(:), INTENT(OUT) :: gradient
135 
136  INTEGER :: i, info, npoints
137  REAL(dp), ALLOCATABLE, DIMENSION(:) :: cov_deriv, weights_deriv
138  REAL(dp), DIMENSION(SIZE(descriptor)) :: kg
139  TYPE(training_matrix_type), POINTER :: training_matrix
140 
141  training_matrix => pao%ml_training_matrices(ikind)
142  npoints = SIZE(training_matrix%outputs, 2)
143 
144  ! calculate derivative of weights
145  ALLOCATE (weights_deriv(npoints))
146  DO i = 1, npoints
147  weights_deriv(i) = sum(outer_deriv*training_matrix%outputs(:, i))
148  END DO
149 
150  ! calculate derivative of covariances
151  ALLOCATE (cov_deriv(npoints))
152  cov_deriv(:) = weights_deriv(:)
153  CALL dpotrs("U", npoints, 1, training_matrix%GP, npoints, cov_deriv, npoints, info)
154  cpassert(info == 0)
155 
156  ! calculate derivative of kernel
157  gradient(:) = 0.0_dp
158  DO i = 1, npoints
159  kg = kernel_grad(pao%gp_scale, descriptor, training_matrix%inputs(:, i))
160  gradient(:) = gradient(:) + kg(:)*cov_deriv(i)
161  END DO
162 
163  DEALLOCATE (cov_deriv, weights_deriv)
164  END SUBROUTINE pao_ml_gp_gradient
165 
166 ! **************************************************************************************************
167 !> \brief Gaussian kernel used to measure covariance between two descriptors.
168 !> \param scale ...
169 !> \param descr1 ...
170 !> \param descr2 ...
171 !> \return ...
172 ! **************************************************************************************************
173  PURE FUNCTION kernel(scale, descr1, descr2) RESULT(cov)
174  REAL(dp), INTENT(IN) :: scale
175  REAL(dp), DIMENSION(:), INTENT(IN) :: descr1, descr2
176  REAL(dp) :: cov
177 
178  REAL(dp) :: fdist2
179  REAL(dp), DIMENSION(SIZE(descr1)) :: diff
180 
181  diff = descr1 - descr2
182  fdist2 = sum((diff/scale)**2)
183  cov = exp(-fdist2/2.0_dp)
184  END FUNCTION kernel
185 
186 ! **************************************************************************************************
187 !> \brief Gradient of Gaussian kernel wrt descr1
188 !> \param scale ...
189 !> \param descr1 ...
190 !> \param descr2 ...
191 !> \return ...
192 ! **************************************************************************************************
193  PURE FUNCTION kernel_grad(scale, descr1, descr2) RESULT(grad)
194  REAL(dp), INTENT(IN) :: scale
195  REAL(dp), DIMENSION(:), INTENT(IN) :: descr1, descr2
196  REAL(dp), DIMENSION(SIZE(descr1)) :: grad
197 
198  REAL(dp) :: cov, fdist2
199  REAL(dp), DIMENSION(SIZE(descr1)) :: diff
200 
201  diff = descr1 - descr2
202  fdist2 = sum((diff/scale)**2)
203  cov = exp(-fdist2/2.0_dp)
204  grad(:) = cov*(-diff/scale**2)
205 
206  END FUNCTION kernel_grad
207 
208 END MODULE pao_ml_gaussprocess
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Gaussian Process implementation.
subroutine, public pao_ml_gp_gradient(pao, ikind, descriptor, outer_deriv, gradient)
Calculate gradient of Gaussian process.
subroutine, public pao_ml_gp_train(pao)
Builds the covariance matrix.
subroutine, public pao_ml_gp_predict(pao, ikind, descriptor, output, variance)
Uses covariance matrix to make prediction.
Types used by the PAO machinery.
Definition: pao_types.F:12