test_gpr.py 9.18 KB
Newer Older
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
1
2
3
4
import numpy as np
import unittest
from gpr import gpr

5
6
7
from ase.io import read

# old gpr
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
8
from kernels import RBF, ConstantKernel as C, WhiteKernel
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
9
from descriptor.fingerprint import Fingerprint
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
from delta_functions_multi.delta import delta as deltaFunc
from GPR import GPR

def initialize_old_gpr(atoms):
    ### Set up feature ###

    # Radial part
    Rc1 = 6
    binwidth1 = 0.2
    sigma1 = 0.2
    
    # Angular part
    Rc2 = 4
    Nbins2 = 30
    sigma2 = 0.2
    gamma = 2
    
    # Radial/angular weighting
    eta = 20
    use_angular = True
    
    # Initialize feature
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
32
    featureCalculator = Fingerprint(Rc1=Rc1, Rc2=Rc2, binwidth1=binwidth1, Nbins2=Nbins2, sigma1=sigma1, sigma2=sigma2, gamma=gamma, eta=eta, use_angular=use_angular)
33

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
34
35
    kernel = C(10, (1e1, 1e6)) * RBF(10, (1,1000))
    #kernel = C(10, (1e1, 1e6)) * (RBF(10, (1,1000)) + C(0.01, (0.01, 0.01)) * RBF(10, (1,1000)) + WhiteKernel(1e-5, (1e-6,1e-2)))
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
36
    delta = deltaFunc(atoms=atoms, rcut=6)
37
38
39
40
41
42
43
44
45
    gpr = GPR(kernel=kernel,
          featureCalculator=featureCalculator,
          delta_function=delta,
          bias_func=None,
          optimize=False,
          n_restarts_optimizer=1)

    return gpr

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
46
47
48
49
50
51
52
def get_E_with_std(traj, gpr):
    E = []
    F = []
    for a in traj:
        e = gpr.predict_energy(a, )

class test_gpr(unittest.TestCase):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
53
54
55
56
57
58
59
60
61
62
63
64

    @classmethod
    def setUpClass(cls):
        print('setupClass')

    @classmethod
    def tearDownClass(cls):
        print('teardownClass')

    def setUp(self):
        print('setUp')
        #self.kernel = gauss_kernel()
65
66
        a = read('structures.traj', index='0')
        self.gpr_old = initialize_old_gpr(a)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
67
68
        #self.gpr = gpr()
        self.gpr = gpr(kernel='single')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
69
70
71
72

    def tearDown(self):
        print('tearDown\n')
    
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
73
74
75
76
77
78
79
80
81
82
83
84
85
86
    def test_compare_training_with_old(self):
        traj = read('structures.traj', index=':50')
        traj_train = traj[:40]
        traj_predict = traj[40:]

        self.gpr_old.train(traj_train, optimize=False)
        self.gpr.train(traj_train)

        np.testing.assert_almost_equal(self.gpr.alpha, self.gpr_old.alpha)

        E_old = np.array([self.gpr_old.predict_energy(a, return_error=True)[:2] for a in traj_predict])
        E = np.array([self.gpr.predict_energy(a, eval_std=True) for a in traj_predict])
        np.testing.assert_almost_equal(E, E_old)

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
87
        F_old = np.array([self.gpr_old.predict_force(a).reshape((-1,3)) for a in traj_predict])
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
88
89
90
        F = np.array([self.gpr.predict_forces(a) for a in traj_predict])
        np.testing.assert_almost_equal(F, F_old)
        
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
91
        Fstd_old = np.array([self.gpr_old.predict_force(a, return_error=True)[1].reshape((-1,3)) for a in traj_predict])
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
92
93
94
95
96
        Fstd = np.array([self.gpr.predict_forces(a, eval_std=True)[1] for a in traj_predict])
        np.testing.assert_almost_equal(Fstd, Fstd_old)


    def test_compare_lml_with_old(self):
97
98
99
100
        traj = read('structures.traj', index=':50')
        traj_train = traj[:40]
        traj_predict = traj[40:]

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
101
102
103
104
105
106
        self.gpr_old.train(traj_train)
        self.gpr.train(traj_train)

        lml_old = self.gpr_old.log_marginal_likelihood_value_
        lml_new = -self.gpr.neg_log_marginal_likelihood(eval_gradient=False)
        np.testing.assert_almost_equal(lml_new, lml_old)
107

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
108
109
110
111
112
113
114
115
116
117
    def test_compare_lml_gradient_with_old(self):
        traj = read('structures.traj', index=':50')
        traj_train = traj[:40]
        traj_predict = traj[40:]

        self.gpr_old.train(traj_train)
        self.gpr.train(traj_train)

        _, lml_ddTheta_old = self.gpr_old.log_marginal_likelihood(self.gpr_old.kernel.theta, eval_gradient=True)
        _, lml_ddTheta = self.gpr.neg_log_marginal_likelihood(eval_gradient=True)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
118
        np.testing.assert_almost_equal(-lml_ddTheta, lml_ddTheta_old)
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
119
120
121
122
123
124
125
126
127
128
129
130
131

    def test_lml_gradient(self):
        traj = read('structures.traj', index=':50')
        traj_train = traj[:40]
        traj_predict = traj[40:]

        self.gpr_old.train(traj_train)
        self.gpr.train(traj_train)

        _, lml_ddTheta = self.gpr.neg_log_marginal_likelihood(eval_gradient=True)
        lml_ddTheta_numeric = self.gpr.numerical_neg_lml()
        np.testing.assert_almost_equal(lml_ddTheta, lml_ddTheta_numeric)
    
132
    def test_forces(self):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
        traj = read('structures.traj', index=':50')
        traj_train = traj[:40]
        traj_predict = traj[40:]

        self.gpr_old.train(traj_train)
        self.gpr.train(traj_train)

        a = traj_predict[0]
        F = self.gpr.predict_forces(a)
        F_numeric = self.gpr.numerical_forces(a)
        np.testing.assert_almost_equal(F, F_numeric)

    def test_forces_std(self):
        traj = read('structures.traj', index=':50')
        traj_train = traj[:40]
        traj_predict = traj[40:]

        self.gpr_old.train(traj_train)
        self.gpr.train(traj_train)

        a = traj_predict[0]
        _, Fstd = self.gpr.predict_forces(a, eval_std=True)
        _, Fstd_numeric = self.gpr.numerical_forces(a, eval_std=True)
        np.testing.assert_almost_equal(Fstd, Fstd_numeric)

    def test_optimize_hyperparameters(self):
        traj = read('structures.traj', index=':50')
        traj_train = traj[:40]
        traj_predict = traj[40:]

        self.gpr_old.train(traj_train)
        self.gpr.train(traj_train)

        self.gpr.optimize_hyperparameters()
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
167
168

if __name__ == '__main__':
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
    #unittest.main()

    import matplotlib.pyplot as plt
    from ase import Atoms
    from ase.visualize import view
    from descriptor.fingerprint import Fingerprint
    from custom_calculators import doubleLJ_calculator
    from gpr import gpr as GPR
    
    def finite_diff(krr, a, dx=1e-5, eval_std=False):
        Natoms, dim = a.positions.shape
        F = np.zeros((Natoms, dim))
        Fstd = np.zeros((Natoms, dim))
        for ia in range(Natoms):
            for idim in range(dim):
                a_up = a.copy()
                a_down = a.copy()
                a_up.positions[ia,idim] += dx/2
                a_down.positions[ia,idim] -= dx/2


                if not eval_std:
                    E_up = krr.predict_energy(a_up, eval_std=False)
                    E_down = krr.predict_energy(a_down, eval_std=False)
                    F[ia,idim] = -(E_up - E_down)/dx
                else:
                    E_up, err_up = krr.predict_energy(a_up, eval_std=True)
                    E_down, err_down = krr.predict_energy(a_down, eval_std=True)
                    
                    F[ia,idim] = -(E_up - E_down)/dx
                    Fstd[ia,idim] = -(err_up - err_down)/dx
        if eval_std:
            return F[1,0], Fstd[1,0]
        else:
            return F
    
    def createData(r):
        positions = np.array([[0,0,0],[r,0,0]])
        a = Atoms('2H', positions, cell=[3,3,1], pbc=[0,0,0])
        calc = doubleLJ_calculator()
        a.set_calculator(calc)
        return a

    def test1():
        a_train = [createData(r) for r in [0.9,1,1.3,2,3]]

        view(a_train[3])
        
        E_train = np.array([a.get_potential_energy() for a in a_train])
        Natoms = a_train[0].get_number_of_atoms()
        
        Rc1 = 5
        binwidth1 = 0.2
        sigma1 = 0.2
        
        Rc2 = 4
        Nbins2 = 30
        sigma2 = 0.2
        
        gamma = 1
        eta = 30
        use_angular = False
        
        descriptor = Fingerprint(Rc1=Rc1, Rc2=Rc2, binwidth1=binwidth1, Nbins2=Nbins2, sigma1=sigma1, sigma2=sigma2, gamma=gamma, eta=eta, use_angular=use_angular)
        
        # Set up KRR-model
        gpr = GPR(kernel='single', descriptor=descriptor)
        
        gpr.train(atoms_list=a_train)
        
        Ntest = 500
        r_test = np.linspace(0.87, 3.5, Ntest)
        E_test = np.zeros(Ntest)
        err_test = np.zeros(Ntest)
        F_test = np.zeros(Ntest)
        Fstd_test = np.zeros(Ntest)
        
        E_true = np.zeros(Ntest)
        F_true = np.zeros(Ntest)
        F_num = np.zeros(Ntest)
        Fstd_num = np.zeros(Ntest)
        for i, r in enumerate(r_test):
            ai = createData(r)
            E, err = gpr.predict_energy(ai, eval_std=True)
            E_test[i] = E
            err_test[i] = err
            
            result_test = gpr.predict_forces(ai, eval_std=True)
            F_test[i] = result_test[0][1,0]
            Fstd_test[i] = result_test[1][1,0]
            F_num[i], Fstd_num[i] = finite_diff(gpr, ai, eval_std=True)
            
            E_true[i] = ai.get_potential_energy()
            F_true[i] = ai.get_forces()[1,0]
            
            
        plt.figure()
        plt.title('Energy')
        plt.xlabel('r')
        plt.ylabel('E')
        plt.plot(r_test, E_true, color='darkgreen', label='true')
        plt.plot(r_test, E_test, color='steelblue', label='model')
        plt.fill_between(r_test, E_test-err_test, E_test+err_test, color='steelblue', alpha=0.3, label='model')
        plt.legend()
        
        plt.figure()
        plt.title('Force')
        plt.xlabel('r')
        plt.ylabel('F')
        plt.plot(r_test, F_true, color='darkgreen', label='true')
        plt.plot(r_test, F_test, color='steelblue', label='model')
        plt.plot(r_test, F_num, 'k:', label='num')
        plt.plot(r_test, F_test-Fstd_test, color='crimson', label='model')
        plt.plot(r_test, F_num-Fstd_num, 'k:', label='num')
        plt.legend()

    test1()
    plt.show()