Coverage for /builds/kinetik161/ase/ase/optimize/precon/precon.py: 85.19%
675 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""
2Implementation of the Precon abstract base class and subclasses
3"""
5import copy
6import sys
7import time
8import warnings
9from abc import ABC, abstractmethod
11import numpy as np
12from scipy import sparse
13from scipy.interpolate import CubicSpline
14from scipy.sparse.linalg import spsolve
16import ase.units as units
17import ase.utils.ff as ff
18from ase.constraints import FixAtoms
19from ase.filters import Filter
20from ase.geometry import find_mic
21from ase.neighborlist import neighbor_list
22from ase.optimize.precon.neighbors import (estimate_nearest_neighbour_distance,
23 get_neighbours)
24from ase.utils import longsum, tokenize_version
26try:
27 import pyamg
28except ImportError:
29 have_pyamg = False
30else:
31 have_pyamg = True
34def create_pyamg_solver(P, max_levels=15):
35 if not have_pyamg:
36 raise RuntimeError(
37 'pyamg not available: install with `pip install pyamg`')
38 filter_key = 'filter'
39 if tokenize_version(pyamg.__version__) >= tokenize_version('4.2.1'):
40 filter_key = 'filter_entries'
41 return pyamg.smoothed_aggregation_solver(
42 P, B=None,
43 strength=('symmetric', {'theta': 0.0}),
44 smooth=(
45 'jacobi', {filter_key: True, 'weighting': 'local'}),
46 improve_candidates=([('block_gauss_seidel',
47 {'sweep': 'symmetric', 'iterations': 4})] +
48 [None] * (max_levels - 1)),
49 aggregate='standard',
50 presmoother=('block_gauss_seidel',
51 {'sweep': 'symmetric', 'iterations': 1}),
52 postsmoother=('block_gauss_seidel',
53 {'sweep': 'symmetric', 'iterations': 1}),
54 max_levels=max_levels,
55 max_coarse=300,
56 coarse_solver='pinv')
59THz = 1e12 * 1. / units.s
62class Precon(ABC):
64 @abstractmethod
65 def make_precon(self, atoms, reinitialize=None):
66 """
67 Create a preconditioner matrix based on the passed set of atoms.
69 Creates a general-purpose preconditioner for use with optimization
70 algorithms, based on examining distances between pairs of atoms in the
71 lattice. The matrix will be stored in the attribute self.P and
72 returned.
74 Args:
75 atoms: the Atoms object used to create the preconditioner.
77 reinitialize: if True, parameters of the preconditioner
78 will be recalculated before the preconditioner matrix is
79 created. If False, they will be calculated only when they
80 do not currently have a value (ie, the first time this
81 function is called).
83 Returns:
84 P: A sparse scipy csr_matrix. BE AWARE that using
85 numpy.dot() with sparse matrices will result in
86 errors/incorrect results - use the .dot method directly
87 on the matrix instead.
88 """
89 ...
91 @abstractmethod
92 def Pdot(self, x):
93 """
94 Return the result of applying P to a vector x
95 """
96 ...
98 def dot(self, x, y):
99 """
100 Return the preconditioned dot product <P x, y>
102 Uses 128-bit floating point math for vector dot products
103 """
104 return longsum(self.Pdot(x) * y)
106 def norm(self, x):
107 """
108 Return the P-norm of x, where |x|_P = sqrt(<Px, x>)
109 """
110 return np.sqrt(self.dot(x, x))
112 def vdot(self, x, y):
113 return self.dot(x.reshape(-1),
114 y.reshape(-1))
116 @abstractmethod
117 def solve(self, x):
118 """
119 Solve the (sparse) linear system P x = y and return y
120 """
121 ...
123 def apply(self, forces, atoms):
124 """
125 Convenience wrapper that combines make_precon() and solve()
127 Parameters
128 ----------
129 forces: array
130 (len(atoms)*3) array of input forces
131 atoms: ase.atoms.Atoms
133 Returns
134 -------
135 precon_forces: array
136 (len(atoms), 3) array of preconditioned forces
137 residual: float
138 inf-norm of original forces, i.e. maximum absolute force
139 """
140 self.make_precon(atoms)
141 residual = np.linalg.norm(forces, np.inf)
142 precon_forces = self.solve(forces)
143 return precon_forces, residual
145 @abstractmethod
146 def copy(self):
147 ...
149 @abstractmethod
150 def asarray(self):
151 """
152 Array representation of preconditioner, as a dense matrix
153 """
154 ...
157class Logfile:
158 def __init__(self, logfile=None):
159 if isinstance(logfile, str):
160 if logfile == "-":
161 logfile = sys.stdout
162 else:
163 logfile = open(logfile, "a")
164 self.logfile = logfile
166 def write(self, *args):
167 if self.logfile is None:
168 return
169 self.logfile.write(*args)
172class SparsePrecon(Precon):
173 def __init__(self, r_cut=None, r_NN=None,
174 mu=None, mu_c=None,
175 dim=3, c_stab=0.1, force_stab=False,
176 reinitialize=False, array_convention='C',
177 solver="auto", solve_tol=1e-8,
178 apply_positions=True, apply_cell=True,
179 estimate_mu_eigmode=False, logfile=None, rng=None,
180 neighbour_list=neighbor_list):
181 """Initialise a preconditioner object based on passed parameters.
183 Parameters:
184 r_cut: float. This is a cut-off radius. The preconditioner matrix
185 will be created by considering pairs of atoms that are within a
186 distance r_cut of each other. For a regular lattice, this is
187 usually taken somewhere between the first- and second-nearest
188 neighbour distance. If r_cut is not provided, default is
189 2 * r_NN (see below)
190 r_NN: nearest neighbour distance. If not provided, this is
191 calculated
192 from input structure.
193 mu: float
194 energy scale for position degreees of freedom. If `None`, mu
195 is precomputed using finite difference derivatives.
196 mu_c: float
197 energy scale for cell degreees of freedom. Also precomputed
198 if None.
199 estimate_mu_eigmode:
200 If True, estimates mu based on the lowest eigenmodes of
201 unstabilised preconditioner. If False it uses the sine based
202 approach.
203 dim: int; dimensions of the problem
204 c_stab: float. The diagonal of the preconditioner matrix will have
205 a stabilisation constant added, which will be the value of
206 c_stab times mu.
207 force_stab:
208 If True, always add the stabilisation to diagnonal, regardless
209 of the presence of fixed atoms.
210 reinitialize: if True, the value of mu will be recalculated when
211 self.make_precon is called. This can be overridden in specific
212 cases with reinitialize argument in self.make_precon. If it
213 is set to True here, the value passed for mu will be
214 irrelevant unless reinitialize is set to False the first time
215 make_precon is called.
216 array_convention: Either 'C' or 'F' for Fortran; this will change
217 the preconditioner to reflect the ordering of the indices in
218 the vector it will operate on. The C convention assumes the
219 vector will be arranged atom-by-atom (ie [x1, y1, z1, x2, ...])
220 while the F convention assumes it will be arranged component
221 by component (ie [x1, x2, ..., y1, y2, ...]).
222 solver: One of "auto", "direct" or "pyamg", specifying whether to
223 use a direct sparse solver or PyAMG to solve P x = y.
224 Default is "auto" which uses PyAMG if available, falling
225 back to sparse solver if not. solve_tol: tolerance used for
226 PyAMG sparse linear solver, if available.
227 apply_positions: bool
228 if True, apply preconditioner to position DoF
229 apply_cell: bool
230 if True, apply preconditioner to cell DoF
231 logfile: file object or str
232 If *logfile* is a string, a file with that name will be opened.
233 Use '-' for stdout.
234 rng: None or np.random.RandomState instance
235 Random number generator to use for initialising pyamg solver
236 neighbor_list: function (optional). Optionally replace the built-in
237 ASE neighbour list with an alternative with the same call
238 signature, e.g. `matscipy.neighbours.neighbour_list`.
240 Raises:
241 ValueError for problem with arguments
243 """
244 self.r_NN = r_NN
245 self.r_cut = r_cut
246 self.mu = mu
247 self.mu_c = mu_c
248 self.estimate_mu_eigmode = estimate_mu_eigmode
249 self.c_stab = c_stab
250 self.force_stab = force_stab
251 self.array_convention = array_convention
252 self.reinitialize = reinitialize
253 self.P = None
254 self.old_positions = None
256 use_pyamg = False
257 if solver == "auto":
258 use_pyamg = have_pyamg
259 elif solver == "direct":
260 use_pyamg = False
261 elif solver == "pyamg":
262 if not have_pyamg:
263 raise RuntimeError("solver='pyamg', PyAMG can't be imported!")
264 use_pyamg = True
265 else:
266 raise ValueError('unknown solver - '
267 'should be "auto", "direct" or "pyamg"')
269 self.use_pyamg = use_pyamg
270 self.solve_tol = solve_tol
271 self.apply_positions = apply_positions
272 self.apply_cell = apply_cell
274 if dim < 1:
275 raise ValueError('Dimension must be at least 1')
276 self.dim = dim
277 self.logfile = Logfile(logfile)
279 if rng is None:
280 rng = np.random.RandomState()
281 self.rng = rng
283 self.neighbor_list = neighbor_list
285 def copy(self):
286 return copy.deepcopy(self)
288 def Pdot(self, x):
289 return self.P.dot(x)
291 def solve(self, x):
292 start_time = time.time()
293 if self.use_pyamg and have_pyamg:
294 y = self.ml.solve(x, x0=self.rng.random(self.P.shape[0]),
295 tol=self.solve_tol,
296 accel='cg',
297 maxiter=300,
298 cycle='W')
299 else:
300 y = spsolve(self.P, x)
301 self.logfile.write(
302 f'--- Precon applied in {(time.time() - start_time)} seconds ---\n')
303 return y
305 def estimate_mu(self, atoms, H=None):
306 r"""Estimate optimal preconditioner coefficient \mu
308 \mu is estimated from a numerical solution of
310 [dE(p+v) - dE(p)] \cdot v = \mu < P1 v, v >
312 with perturbation
314 v(x,y,z) = H P_lowest_nonzero_eigvec(x, y, z)
316 or
318 v(x,y,z) = H (sin(x / Lx), sin(y / Ly), sin(z / Lz))
320 After the optimal \mu is found, self.mu will be set to its value.
322 If `atoms` is an instance of Filter an additional \mu_c
323 will be computed for the cell degrees of freedom .
325 Args:
326 atoms: Atoms object for initial system
328 H: 3x3 array or None
329 Magnitude of deformation to apply.
330 Default is 1e-2*rNN*np.eye(3)
332 Returns:
333 mu : float
334 mu_c : float or None
335 """
336 logfile = self.logfile
338 if self.dim != 3:
339 raise ValueError('Automatic calculation of mu only possible for '
340 'three-dimensional preconditioners. Try setting '
341 'mu manually instead.')
343 if self.r_NN is None:
344 self.r_NN = estimate_nearest_neighbour_distance(atoms,
345 self.neighbor_list)
347 # deformation matrix, default is diagonal
348 if H is None:
349 H = 1e-2 * self.r_NN * np.eye(3)
351 # compute perturbation
352 p = atoms.get_positions()
354 if self.estimate_mu_eigmode:
355 self.mu = 1.0
356 self.mu_c = 1.0
357 c_stab = self.c_stab
358 self.c_stab = 0.0
360 if isinstance(atoms, Filter):
361 n = len(atoms.atoms)
362 else:
363 n = len(atoms)
364 self._make_sparse_precon(atoms, initial_assembly=True)
365 self.P = self.P[:3 * n, :3 * n]
366 eigvals, eigvecs = sparse.linalg.eigsh(self.P, k=4, which='SM')
368 logfile.write('estimate_mu(): lowest 4 eigvals = %f %f %f %f\n' %
369 (eigvals[0], eigvals[1], eigvals[2], eigvals[3]))
370 # check eigenvalues
371 if any(eigvals[0:3] > 1e-6):
372 raise ValueError('First 3 eigenvalues of preconditioner matrix'
373 'do not correspond to translational modes.')
374 elif eigvals[3] < 1e-6:
375 raise ValueError('Fourth smallest eigenvalue of '
376 'preconditioner matrix '
377 'is too small, increase r_cut.')
379 x = np.zeros(n)
380 for i in range(n):
381 x[i] = eigvecs[:, 3][3 * i]
382 x = x / np.linalg.norm(x)
383 if x[0] < 0:
384 x = -x
386 v = np.zeros(3 * len(atoms))
387 for i in range(n):
388 v[3 * i] = x[i]
389 v[3 * i + 1] = x[i]
390 v[3 * i + 2] = x[i]
391 v = v / np.linalg.norm(v)
392 v = v.reshape((-1, 3))
394 self.c_stab = c_stab
395 else:
396 Lx, Ly, Lz = (p[:, i].max() - p[:, i].min() for i in range(3))
397 logfile.write('estimate_mu(): Lx=%.1f Ly=%.1f Lz=%.1f\n' %
398 (Lx, Ly, Lz))
400 x, y, z = p.T
401 # sine_vr = [np.sin(x/Lx), np.sin(y/Ly), np.sin(z/Lz)], but we need
402 # to take into account the possibility that one of Lx/Ly/Lz is
403 # zero.
404 sine_vr = [x, y, z]
406 for i, L in enumerate([Lx, Ly, Lz]):
407 if L == 0:
408 warnings.warn(
409 'Cell length L[%d] == 0. Setting H[%d,%d] = 0.' %
410 (i, i, i))
411 H[i, i] = 0.0
412 else:
413 sine_vr[i] = np.sin(sine_vr[i] / L)
415 v = np.dot(H, sine_vr).T
417 natoms = len(atoms)
418 if isinstance(atoms, Filter):
419 natoms = len(atoms.atoms)
420 eps = H / self.r_NN
421 v[natoms:, :] = eps
423 v1 = v.reshape(-1)
425 # compute LHS
426 dE_p = -atoms.get_forces().reshape(-1)
427 atoms_v = atoms.copy()
428 atoms_v.calc = atoms.calc
429 if isinstance(atoms, Filter):
430 atoms_v = atoms.__class__(atoms_v)
431 if hasattr(atoms, 'constant_volume'):
432 atoms_v.constant_volume = atoms.constant_volume
433 atoms_v.set_positions(p + v)
434 dE_p_plus_v = -atoms_v.get_forces().reshape(-1)
436 # compute left hand side
437 LHS = (dE_p_plus_v - dE_p) * v1
439 # assemble P with \mu = 1
440 self.mu = 1.0
441 self.mu_c = 1.0
443 self._make_sparse_precon(atoms, initial_assembly=True)
445 # compute right hand side
446 RHS = self.P.dot(v1) * v1
448 # use partial sums to compute separate mu for positions and cell DoFs
449 self.mu = longsum(LHS[:3 * natoms]) / longsum(RHS[:3 * natoms])
450 if self.mu < 1.0:
451 logfile.write('estimate_mu(): mu (%.3f) < 1.0, '
452 'capping at mu=1.0' % self.mu)
453 self.mu = 1.0
455 if isinstance(atoms, Filter):
456 self.mu_c = longsum(LHS[3 * natoms:]) / longsum(RHS[3 * natoms:])
457 if self.mu_c < 1.0:
458 logfile.write('estimate_mu(): mu_c (%.3f) < 1.0, '
459 'capping at mu_c=1.0\n' % self.mu_c)
460 self.mu_c = 1.0
462 logfile.write('estimate_mu(): mu=%r, mu_c=%r\n' %
463 (self.mu, self.mu_c))
465 self.P = None # force a rebuild with new mu (there may be fixed atoms)
466 return (self.mu, self.mu_c)
468 def asarray(self):
469 return np.array(self.P.todense())
471 def one_dim_to_ndim(self, csc_P, N):
472 """
473 Expand an N x N precon matrix to self.dim*N x self.dim * N
475 Args:
476 csc_P (sparse matrix): N x N sparse matrix, in CSC format
477 """
478 start_time = time.time()
479 if self.dim == 1:
480 P = csc_P
481 elif self.array_convention == 'F':
482 csc_P = csc_P.tocsr()
483 P = csc_P
484 for i in range(self.dim - 1):
485 P = sparse.block_diag((P, csc_P)).tocsr()
486 else:
487 # convert back to triplet and read the arrays
488 csc_P = csc_P.tocoo()
489 i = csc_P.row * self.dim
490 j = csc_P.col * self.dim
491 z = csc_P.data
493 # N-dimensionalise, interlaced coordinates
494 I = np.hstack([i + d for d in range(self.dim)])
495 J = np.hstack([j + d for d in range(self.dim)])
496 Z = np.hstack([z for d in range(self.dim)])
497 P = sparse.csc_matrix((Z, (I, J)),
498 shape=(self.dim * N, self.dim * N))
499 P = P.tocsr()
500 self.logfile.write(
501 f'--- N-dim precon created in {(time.time() - start_time)} s ---\n')
502 return P
504 def create_solver(self):
505 if self.use_pyamg and have_pyamg:
506 start_time = time.time()
507 self.ml = create_pyamg_solver(self.P)
508 self.logfile.write(
509 f'--- multi grid solver created in {(time.time() - start_time)}'
510 ' ---\n')
513class SparseCoeffPrecon(SparsePrecon):
514 def _make_sparse_precon(self, atoms, initial_assembly=False,
515 force_stab=False):
516 """Create a sparse preconditioner matrix based on the passed atoms.
518 Creates a general-purpose preconditioner for use with optimization
519 algorithms, based on examining distances between pairs of atoms in the
520 lattice. The matrix will be stored in the attribute self.P and
521 returned. Note that this function will use self.mu, whatever it is.
523 Args:
524 atoms: the Atoms object used to create the preconditioner.
526 Returns:
527 A scipy.sparse.csr_matrix object, representing a d*N by d*N matrix
528 (where N is the number of atoms, and d is the value of self.dim).
529 BE AWARE that using numpy.dot() with this object will result in
530 errors/incorrect results - use the .dot method directly on the
531 sparse matrix instead.
533 """
534 logfile = self.logfile
535 logfile.write('creating sparse precon: initial_assembly=%r, '
536 'force_stab=%r, apply_positions=%r, apply_cell=%r\n' %
537 (initial_assembly, force_stab, self.apply_positions,
538 self.apply_cell))
540 N = len(atoms)
541 diag_i = np.arange(N, dtype=int)
542 start_time = time.time()
543 if self.apply_positions:
544 # compute neighbour list
545 i, j, rij, fixed_atoms = get_neighbours(
546 atoms, self.r_cut,
547 neighbor_list=self.neighbor_list)
548 logfile.write(
549 f'--- neighbour list created in {(time.time() - start_time)} s '
550 '--- \n')
552 # compute entries in triplet format: without the constraints
553 start_time = time.time()
554 coeff = self.get_coeff(rij)
555 diag_coeff = np.bincount(i, -coeff, minlength=N).astype(np.float64)
556 if force_stab or len(fixed_atoms) == 0:
557 logfile.write('adding stabilisation to precon')
558 diag_coeff += self.mu * self.c_stab
559 else:
560 diag_coeff = np.ones(N)
562 # precon is mu_c * identity for cell DoF
563 if isinstance(atoms, Filter):
564 if self.apply_cell:
565 diag_coeff[-3:] = self.mu_c
566 else:
567 diag_coeff[-3:] = 1.0
568 logfile.write(
569 f'--- computed triplet format in {(time.time() - start_time)} s '
570 '---\n')
572 if self.apply_positions and not initial_assembly:
573 # apply the constraints
574 start_time = time.time()
575 mask = np.ones(N)
576 mask[fixed_atoms] = 0.0
577 coeff *= mask[i] * mask[j]
578 diag_coeff[fixed_atoms] = 1.0
579 logfile.write(
580 f'--- applied fixed_atoms in {(time.time() - start_time)} s '
581 '---\n')
583 if self.apply_positions:
584 # remove zeros
585 start_time = time.time()
586 inz = np.nonzero(coeff)
587 i = np.hstack((i[inz], diag_i))
588 j = np.hstack((j[inz], diag_i))
589 coeff = np.hstack((coeff[inz], diag_coeff))
590 logfile.write(
591 f'--- remove zeros in {(time.time() - start_time)} s '
592 '---\n')
593 else:
594 i = diag_i
595 j = diag_i
596 coeff = diag_coeff
598 # create an N x N precon matrix in compressed sparse column (CSC) format
599 start_time = time.time()
600 csc_P = sparse.csc_matrix((coeff, (i, j)), shape=(N, N))
601 logfile.write(
602 f'--- created CSC matrix in {(time.time() - start_time)} s '
603 '---\n')
605 self.P = self.one_dim_to_ndim(csc_P, N)
606 self.create_solver()
608 def make_precon(self, atoms, reinitialize=None):
609 if self.r_NN is None:
610 self.r_NN = estimate_nearest_neighbour_distance(atoms,
611 self.neighbor_list)
613 if self.r_cut is None:
614 # This is the first time this function has been called, and no
615 # cutoff radius has been specified, so calculate it automatically.
616 self.r_cut = 2.0 * self.r_NN
617 elif self.r_cut < self.r_NN:
618 warning = ('WARNING: r_cut (%.2f) < r_NN (%.2f), '
619 'increasing to 1.1*r_NN = %.2f' % (self.r_cut,
620 self.r_NN,
621 1.1 * self.r_NN))
622 warnings.warn(warning)
623 self.r_cut = 1.1 * self.r_NN
625 if reinitialize is None:
626 # The caller has not specified whether or not to recalculate mu,
627 # so the Precon's setting is used.
628 reinitialize = self.reinitialize
630 if self.mu is None:
631 # Regardless of what the caller has specified, if we don't
632 # currently have a value of mu, then we need one.
633 reinitialize = True
635 if reinitialize:
636 self.estimate_mu(atoms)
638 if self.P is not None:
639 real_atoms = atoms
640 if isinstance(atoms, Filter):
641 real_atoms = atoms.atoms
642 if self.old_positions is None:
643 self.old_positions = real_atoms.positions
644 displacement, _ = find_mic(real_atoms.positions -
645 self.old_positions,
646 real_atoms.cell, real_atoms.pbc)
647 self.old_positions = real_atoms.get_positions()
648 max_abs_displacement = abs(displacement).max()
649 self.logfile.write('max(abs(displacements)) = %.2f A (%.2f r_NN)' %
650 (max_abs_displacement,
651 max_abs_displacement / self.r_NN))
652 if max_abs_displacement < 0.5 * self.r_NN:
653 return
655 start_time = time.time()
656 self._make_sparse_precon(atoms, force_stab=self.force_stab)
657 self.logfile.write(
658 f'--- Precon created in {(time.time() - start_time)} seconds '
659 '--- \n')
661 @abstractmethod
662 def get_coeff(self, r):
663 ...
666class Pfrommer(Precon):
667 """
668 Use initial guess for inverse Hessian from Pfrommer et al. as a
669 simple preconditioner
671 J. Comput. Phys. vol 131 p233-240 (1997)
672 """
674 def __init__(self, bulk_modulus=500 * units.GPa, phonon_frequency=50 * THz,
675 apply_positions=True, apply_cell=True):
676 """
677 Default bulk modulus is 500 GPa and default phonon frequency is 50 THz
678 """
679 self.bulk_modulus = bulk_modulus
680 self.phonon_frequency = phonon_frequency
681 self.apply_positions = apply_positions
682 self.apply_cell = apply_cell
683 self.H0 = None
685 def make_precon(self, atoms, reinitialize=None):
686 if self.H0 is not None:
687 # only build H0 on first call
688 return
690 variable_cell = False
691 if isinstance(atoms, Filter):
692 variable_cell = True
693 atoms = atoms.atoms
695 # position DoF
696 omega = self.phonon_frequency
697 mass = atoms.get_masses().mean()
698 block = np.eye(3) / (mass * omega**2)
699 blocks = [block] * len(atoms)
701 # cell DoF
702 if variable_cell:
703 coeff = 1.0
704 if self.apply_cell:
705 coeff = 1.0 / (3 * self.bulk_modulus)
706 blocks.append(np.diag([coeff] * 9))
708 self.H0 = sparse.block_diag(blocks, format='csr')
709 return
711 def Pdot(self, x):
712 return self.H0.solve(x)
714 def solve(self, x):
715 y = self.H0.dot(x)
716 return y
718 def copy(self):
719 return Pfrommer(self.bulk_modulus,
720 self.phonon_frequency,
721 self.apply_positions,
722 self.apply_cell)
724 def asarray(self):
725 return np.array(self.H0.todense())
728class IdentityPrecon(Precon):
729 """
730 Dummy preconditioner which does not modify forces
731 """
733 def make_precon(self, atoms, reinitialize=None):
734 self.atoms = atoms
736 def Pdot(self, x):
737 return x
739 def solve(self, x):
740 return x
742 def copy(self):
743 return IdentityPrecon()
745 def asarray(self):
746 return np.eye(3 * len(self.atoms))
749class C1(SparseCoeffPrecon):
750 """Creates matrix by inserting a constant whenever r_ij is less than r_cut.
751 """
753 def __init__(self, r_cut=None, mu=None, mu_c=None, dim=3, c_stab=0.1,
754 force_stab=False,
755 reinitialize=False, array_convention='C',
756 solver="auto", solve_tol=1e-9,
757 apply_positions=True, apply_cell=True, logfile=None):
758 super().__init__(r_cut=r_cut, mu=mu, mu_c=mu_c,
759 dim=dim, c_stab=c_stab,
760 force_stab=force_stab,
761 reinitialize=reinitialize,
762 array_convention=array_convention,
763 solver=solver, solve_tol=solve_tol,
764 apply_positions=apply_positions,
765 apply_cell=apply_cell,
766 logfile=logfile)
768 def get_coeff(self, r):
769 return -self.mu * np.ones_like(r)
772class Exp(SparseCoeffPrecon):
773 """Creates matrix with values decreasing exponentially with distance.
774 """
776 def __init__(self, A=3.0, r_cut=None, r_NN=None, mu=None, mu_c=None,
777 dim=3, c_stab=0.1,
778 force_stab=False, reinitialize=False, array_convention='C',
779 solver="auto", solve_tol=1e-9,
780 apply_positions=True, apply_cell=True,
781 estimate_mu_eigmode=False, logfile=None):
782 """
783 Initialise an Exp preconditioner with given parameters.
785 Args:
786 r_cut, mu, c_stab, dim, sparse, reinitialize, array_convention: see
787 precon.__init__()
788 A: coefficient in exp(-A*r/r_NN). Default is A=3.0.
789 """
790 super().__init__(r_cut=r_cut, r_NN=r_NN,
791 mu=mu, mu_c=mu_c, dim=dim, c_stab=c_stab,
792 force_stab=force_stab,
793 reinitialize=reinitialize,
794 array_convention=array_convention,
795 solver=solver,
796 solve_tol=solve_tol,
797 apply_positions=apply_positions,
798 apply_cell=apply_cell,
799 estimate_mu_eigmode=estimate_mu_eigmode,
800 logfile=logfile)
802 self.A = A
804 def get_coeff(self, r):
805 return -self.mu * np.exp(-self.A * (r / self.r_NN - 1))
808def ij_to_x(i, j):
809 x = [3 * i, 3 * i + 1, 3 * i + 2,
810 3 * j, 3 * j + 1, 3 * j + 2]
811 return x
814def ijk_to_x(i, j, k):
815 x = [3 * i, 3 * i + 1, 3 * i + 2,
816 3 * j, 3 * j + 1, 3 * j + 2,
817 3 * k, 3 * k + 1, 3 * k + 2]
818 return x
821def ijkl_to_x(i, j, k, l):
822 x = [3 * i, 3 * i + 1, 3 * i + 2,
823 3 * j, 3 * j + 1, 3 * j + 2,
824 3 * k, 3 * k + 1, 3 * k + 2,
825 3 * l, 3 * l + 1, 3 * l + 2]
826 return x
829def apply_fixed(atoms, P):
830 fixed_atoms = []
831 for constraint in atoms.constraints:
832 if isinstance(constraint, FixAtoms):
833 fixed_atoms.extend(list(constraint.index))
834 else:
835 raise TypeError(
836 'only FixAtoms constraints are supported by Precon class')
837 if len(fixed_atoms) != 0:
838 P = P.tolil()
839 for i in fixed_atoms:
840 P[i, :] = 0.0
841 P[:, i] = 0.0
842 P[i, i] = 1.0
843 return P
846class FF(SparsePrecon):
847 """Creates matrix using morse/bond/angle/dihedral force field parameters.
848 """
850 def __init__(self, dim=3, c_stab=0.1, force_stab=False,
851 array_convention='C', solver="auto", solve_tol=1e-9,
852 apply_positions=True, apply_cell=True,
853 hessian='spectral', morses=None, bonds=None, angles=None,
854 dihedrals=None, logfile=None):
855 """Initialise an FF preconditioner with given parameters.
857 Args:
858 dim, c_stab, force_stab, array_convention, use_pyamg, solve_tol:
859 see SparsePrecon.__init__()
860 morses: Morse instance
861 bonds: Bond instance
862 angles: Angle instance
863 dihedrals: Dihedral instance
864 """
866 if (morses is None and bonds is None and angles is None and
867 dihedrals is None):
868 raise ImportError(
869 'At least one of morses, bonds, angles or dihedrals must be '
870 'defined!')
872 super().__init__(dim=dim, c_stab=c_stab,
873 force_stab=force_stab,
874 array_convention=array_convention,
875 solver=solver,
876 solve_tol=solve_tol,
877 apply_positions=apply_positions,
878 apply_cell=apply_cell, logfile=logfile)
880 self.hessian = hessian
881 self.morses = morses
882 self.bonds = bonds
883 self.angles = angles
884 self.dihedrals = dihedrals
886 def make_precon(self, atoms, reinitialize=None):
887 start_time = time.time()
888 self._make_sparse_precon(atoms, force_stab=self.force_stab)
889 self.logfile.write(
890 f'--- Precon created in {(time.time() - start_time)} seconds '
891 '---\n')
893 def add_morse(self, morse, atoms, row, col, data, conn=None):
894 if self.hessian == 'reduced':
895 i, j, Hx = ff.get_morse_potential_reduced_hessian(
896 atoms, morse)
897 elif self.hessian == 'spectral':
898 i, j, Hx = ff.get_morse_potential_hessian(
899 atoms, morse, spectral=True)
900 else:
901 raise NotImplementedError('Not implemented hessian')
902 x = ij_to_x(i, j)
903 row.extend(np.repeat(x, 6))
904 col.extend(np.tile(x, 6))
905 data.extend(Hx.flatten())
906 if conn is not None:
907 conn[i, j] = True
908 conn[j, i] = True
910 def add_bond(self, bond, atoms, row, col, data, conn=None):
911 if self.hessian == 'reduced':
912 i, j, Hx = ff.get_bond_potential_reduced_hessian(
913 atoms, bond, self.morses)
914 elif self.hessian == 'spectral':
915 i, j, Hx = ff.get_bond_potential_hessian(
916 atoms, bond, self.morses, spectral=True)
917 else:
918 raise NotImplementedError('Not implemented hessian')
919 x = ij_to_x(i, j)
920 row.extend(np.repeat(x, 6))
921 col.extend(np.tile(x, 6))
922 data.extend(Hx.flatten())
923 if conn is not None:
924 conn[i, j] = True
925 conn[j, i] = True
927 def add_angle(self, angle, atoms, row, col, data, conn=None):
928 if self.hessian == 'reduced':
929 i, j, k, Hx = ff.get_angle_potential_reduced_hessian(
930 atoms, angle, self.morses)
931 elif self.hessian == 'spectral':
932 i, j, k, Hx = ff.get_angle_potential_hessian(
933 atoms, angle, self.morses, spectral=True)
934 else:
935 raise NotImplementedError('Not implemented hessian')
936 x = ijk_to_x(i, j, k)
937 row.extend(np.repeat(x, 9))
938 col.extend(np.tile(x, 9))
939 data.extend(Hx.flatten())
940 if conn is not None:
941 conn[i, j] = conn[i, k] = conn[j, k] = True
942 conn[j, i] = conn[k, i] = conn[k, j] = True
944 def add_dihedral(self, dihedral, atoms, row, col, data, conn=None):
945 if self.hessian == 'reduced':
946 i, j, k, l, Hx = \
947 ff.get_dihedral_potential_reduced_hessian(
948 atoms, dihedral, self.morses)
949 elif self.hessian == 'spectral':
950 i, j, k, l, Hx = ff.get_dihedral_potential_hessian(
951 atoms, dihedral, self.morses, spectral=True)
952 else:
953 raise NotImplementedError('Not implemented hessian')
954 x = ijkl_to_x(i, j, k, l)
955 row.extend(np.repeat(x, 12))
956 col.extend(np.tile(x, 12))
957 data.extend(Hx.flatten())
958 if conn is not None:
959 conn[i, j] = conn[i, k] = conn[i, l] = conn[
960 j, k] = conn[j, l] = conn[k, l] = True
961 conn[j, i] = conn[k, i] = conn[l, i] = conn[
962 k, j] = conn[l, j] = conn[l, k] = True
964 def _make_sparse_precon(self, atoms, initial_assembly=False,
965 force_stab=False):
966 N = len(atoms)
968 row = []
969 col = []
970 data = []
972 if self.morses is not None:
973 for morse in self.morses:
974 self.add_morse(morse, atoms, row, col, data)
976 if self.bonds is not None:
977 for bond in self.bonds:
978 self.add_bond(bond, atoms, row, col, data)
980 if self.angles is not None:
981 for angle in self.angles:
982 self.add_angle(angle, atoms, row, col, data)
984 if self.dihedrals is not None:
985 for dihedral in self.dihedrals:
986 self.add_dihedral(dihedral, atoms, row, col, data)
988 # add the diagonal
989 row.extend(range(self.dim * N))
990 col.extend(range(self.dim * N))
991 data.extend([self.c_stab] * self.dim * N)
993 # create the matrix
994 start_time = time.time()
995 self.P = sparse.csc_matrix(
996 (data, (row, col)), shape=(self.dim * N, self.dim * N))
997 self.logfile.write(
998 f'--- created CSC matrix in {(time.time() - start_time)} s ---\n')
1000 self.P = apply_fixed(atoms, self.P)
1001 self.P = self.P.tocsr()
1002 self.logfile.write(
1003 f'--- N-dim precon created in {(time.time() - start_time)} s ---\n')
1004 self.create_solver()
1007class Exp_FF(Exp, FF):
1008 """Creates matrix with values decreasing exponentially with distance.
1009 """
1011 def __init__(self, A=3.0, r_cut=None, r_NN=None, mu=None, mu_c=None,
1012 dim=3, c_stab=0.1,
1013 force_stab=False, reinitialize=False, array_convention='C',
1014 solver="auto", solve_tol=1e-9,
1015 apply_positions=True, apply_cell=True,
1016 estimate_mu_eigmode=False,
1017 hessian='spectral', morses=None, bonds=None, angles=None,
1018 dihedrals=None, logfile=None):
1019 """Initialise an Exp+FF preconditioner with given parameters.
1021 Args:
1022 r_cut, mu, c_stab, dim, reinitialize, array_convention: see
1023 precon.__init__()
1024 A: coefficient in exp(-A*r/r_NN). Default is A=3.0.
1025 """
1026 if (morses is None and bonds is None and angles is None and
1027 dihedrals is None):
1028 raise ImportError(
1029 'At least one of morses, bonds, angles or dihedrals must '
1030 'be defined!')
1032 SparsePrecon.__init__(self, r_cut=r_cut, r_NN=r_NN,
1033 mu=mu, mu_c=mu_c, dim=dim, c_stab=c_stab,
1034 force_stab=force_stab,
1035 reinitialize=reinitialize,
1036 array_convention=array_convention,
1037 solver=solver,
1038 solve_tol=solve_tol,
1039 apply_positions=apply_positions,
1040 apply_cell=apply_cell,
1041 estimate_mu_eigmode=estimate_mu_eigmode,
1042 logfile=logfile)
1044 self.A = A
1045 self.hessian = hessian
1046 self.morses = morses
1047 self.bonds = bonds
1048 self.angles = angles
1049 self.dihedrals = dihedrals
1051 def make_precon(self, atoms, reinitialize=None):
1052 if self.r_NN is None:
1053 self.r_NN = estimate_nearest_neighbour_distance(atoms,
1054 self.neighbor_list)
1056 if self.r_cut is None:
1057 # This is the first time this function has been called, and no
1058 # cutoff radius has been specified, so calculate it automatically.
1059 self.r_cut = 2.0 * self.r_NN
1060 elif self.r_cut < self.r_NN:
1061 warning = ('WARNING: r_cut (%.2f) < r_NN (%.2f), '
1062 'increasing to 1.1*r_NN = %.2f' % (self.r_cut,
1063 self.r_NN,
1064 1.1 * self.r_NN))
1065 warnings.warn(warning)
1066 self.r_cut = 1.1 * self.r_NN
1068 if reinitialize is None:
1069 # The caller has not specified whether or not to recalculate mu,
1070 # so the Precon's setting is used.
1071 reinitialize = self.reinitialize
1073 if self.mu is None:
1074 # Regardless of what the caller has specified, if we don't
1075 # currently have a value of mu, then we need one.
1076 reinitialize = True
1078 if reinitialize:
1079 self.estimate_mu(atoms)
1081 if self.P is not None:
1082 real_atoms = atoms
1083 if isinstance(atoms, Filter):
1084 real_atoms = atoms.atoms
1085 if self.old_positions is None:
1086 self.old_positions = real_atoms.positions
1087 displacement, _ = find_mic(real_atoms.positions -
1088 self.old_positions,
1089 real_atoms.cell, real_atoms.pbc)
1090 self.old_positions = real_atoms.get_positions()
1091 max_abs_displacement = abs(displacement).max()
1092 self.logfile.write('max(abs(displacements)) = %.2f A (%.2f r_NN)' %
1093 (max_abs_displacement,
1094 max_abs_displacement / self.r_NN))
1095 if max_abs_displacement < 0.5 * self.r_NN:
1096 return
1098 # Create the preconditioner:
1099 start_time = time.time()
1100 self._make_sparse_precon(atoms, force_stab=self.force_stab)
1101 self.logfile.write(
1102 f'--- Precon created in {(time.time() - start_time)} seconds ---\n')
1104 def _make_sparse_precon(self, atoms, initial_assembly=False,
1105 force_stab=False):
1106 """Create a sparse preconditioner matrix based on the passed atoms.
1108 Args:
1109 atoms: the Atoms object used to create the preconditioner.
1111 Returns:
1112 A scipy.sparse.csr_matrix object, representing a d*N by d*N matrix
1113 (where N is the number of atoms, and d is the value of self.dim).
1114 BE AWARE that using numpy.dot() with this object will result in
1115 errors/incorrect results - use the .dot method directly on the
1116 sparse matrix instead.
1118 """
1119 self.logfile.write('creating sparse precon: initial_assembly=%r, '
1120 'force_stab=%r, apply_positions=%r, '
1121 'apply_cell=%r\n' %
1122 (initial_assembly, force_stab,
1123 self.apply_positions, self.apply_cell))
1125 N = len(atoms)
1126 start_time = time.time()
1127 if self.apply_positions:
1128 # compute neighbour list
1129 i_list, j_list, rij_list, fixed_atoms = get_neighbours(
1130 atoms, self.r_cut, self.neighbor_list)
1131 self.logfile.write(
1132 f'--- neighbour list created in {(time.time() - start_time)} s '
1133 '---\n')
1135 row = []
1136 col = []
1137 data = []
1139 # precon is mu_c*identity for cell DoF
1140 start_time = time.time()
1141 if isinstance(atoms, Filter):
1142 i = N - 3
1143 j = N - 2
1144 k = N - 1
1145 x = ijk_to_x(i, j, k)
1146 row.extend(x)
1147 col.extend(x)
1148 if self.apply_cell:
1149 data.extend(np.repeat(self.mu_c, 9))
1150 else:
1151 data.extend(np.repeat(self.mu_c, 9))
1152 self.logfile.write(
1153 f'--- computed triplet format in {(time.time() - start_time)} s '
1154 '---\n')
1156 conn = sparse.lil_matrix((N, N), dtype=bool)
1158 if self.apply_positions and not initial_assembly:
1159 if self.morses is not None:
1160 for morse in self.morses:
1161 self.add_morse(morse, atoms, row, col, data, conn)
1163 if self.bonds is not None:
1164 for bond in self.bonds:
1165 self.add_bond(bond, atoms, row, col, data, conn)
1167 if self.angles is not None:
1168 for angle in self.angles:
1169 self.add_angle(angle, atoms, row, col, data, conn)
1171 if self.dihedrals is not None:
1172 for dihedral in self.dihedrals:
1173 self.add_dihedral(dihedral, atoms, row, col, data, conn)
1175 if self.apply_positions:
1176 for i, j, rij in zip(i_list, j_list, rij_list):
1177 if not conn[i, j]:
1178 coeff = self.get_coeff(rij)
1179 x = ij_to_x(i, j)
1180 row.extend(x)
1181 col.extend(x)
1182 data.extend(3 * [-coeff] + 3 * [coeff])
1184 row.extend(range(self.dim * N))
1185 col.extend(range(self.dim * N))
1186 if initial_assembly:
1187 data.extend([self.mu * self.c_stab] * self.dim * N)
1188 else:
1189 data.extend([self.c_stab] * self.dim * N)
1191 # create the matrix
1192 start_time = time.time()
1193 self.P = sparse.csc_matrix(
1194 (data, (row, col)), shape=(self.dim * N, self.dim * N))
1195 self.logfile.write(
1196 f'--- created CSC matrix in {(time.time() - start_time)} s ---\n')
1198 if not initial_assembly:
1199 self.P = apply_fixed(atoms, self.P)
1201 self.P = self.P.tocsr()
1202 self.create_solver()
1205def make_precon(precon, atoms=None, **kwargs):
1206 """
1207 Construct preconditioner from a string and optionally build for Atoms
1209 Parameters
1210 ----------
1211 precon - one of 'C1', 'Exp', 'Pfrommer', 'FF', 'Exp_FF', 'ID', None
1212 or an instance of a subclass of `ase.optimize.precon.Precon`
1214 atoms - ase.atoms.Atoms instance, optional
1215 If present, build apreconditioner for this Atoms object
1217 **kwargs - additional keyword arguments to pass to Precon constructor
1219 Returns
1220 -------
1221 precon - instance of relevant subclass of `ase.optimize.precon.Precon`
1222 """
1223 lookup = {
1224 'C1': C1,
1225 'Exp': Exp,
1226 'Pfrommer': Pfrommer,
1227 'FF': FF,
1228 'Exp_FF': Exp_FF,
1229 'ID': IdentityPrecon,
1230 'IdentityPrecon': IdentityPrecon,
1231 None: IdentityPrecon
1232 }
1233 if isinstance(precon, str) or precon is None:
1234 cls = lookup[precon]
1235 precon = cls(**kwargs)
1236 if atoms is not None:
1237 precon.make_precon(atoms)
1238 return precon
1241class SplineFit:
1242 """
1243 Fit a cubic spline fit to images
1244 """
1246 def __init__(self, s, x):
1247 self._s = s
1248 self._x_data = x
1249 self._x = CubicSpline(self._s, x, bc_type='not-a-knot')
1250 self._dx_ds = self._x.derivative()
1251 self._d2x_ds2 = self._x.derivative(2)
1253 @property
1254 def s(self):
1255 return self._s
1257 @property
1258 def x_data(self):
1259 return self._x_data
1261 @property
1262 def x(self):
1263 return self._x
1265 @property
1266 def dx_ds(self):
1267 return self._dx_ds
1269 @property
1270 def d2x_ds2(self):
1271 return self._d2x_ds2
1274class PreconImages:
1275 def __init__(self, precon, images, **kwargs):
1276 """
1277 Wrapper for a list of Precon objects and associated images
1279 This is used when preconditioning a NEB object. Equation references
1280 refer to Paper IV in the :class:`ase.mep.NEB` documentation, i.e.
1282 S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys.
1283 150, 094109 (2019)
1284 https://dx.doi.org/10.1063/1.5064465
1286 Args:
1287 precon (str or list): preconditioner(s) to use
1288 images (list of Atoms): Atoms objects that define the state
1290 """
1291 self.images = images
1292 self._spline = None
1294 if isinstance(precon, list):
1295 if len(precon) != len(images):
1296 raise ValueError(f'length mismatch: len(precon)={len(precon)} '
1297 f'!= len(images)={len(images)}')
1298 self.precon = precon
1299 return
1300 P0 = make_precon(precon, images[0], **kwargs)
1301 self.precon = [P0]
1302 for image in images[1:]:
1303 P = P0.copy()
1304 P.make_precon(image)
1305 self.precon.append(P)
1307 def __len__(self):
1308 return len(self.precon)
1310 def __iter__(self):
1311 return iter(self.precon)
1313 def __getitem__(self, index):
1314 return self.precon[index]
1316 def apply(self, all_forces, index=None):
1317 """Apply preconditioners to stored images
1319 Args:
1320 all_forces (array): forces on images, shape (nimages, natoms, 3)
1321 index (slice, optional): Which images to include. Defaults to all.
1323 Returns:
1324 precon_forces: array of preconditioned forces
1325 """
1326 if index is None:
1327 index = slice(None)
1328 precon_forces = []
1329 for precon, image, forces in zip(self.precon[index],
1330 self.images[index],
1331 all_forces):
1332 f_vec = forces.reshape(-1)
1333 pf_vec, _ = precon.apply(f_vec, image)
1334 precon_forces.append(pf_vec.reshape(-1, 3))
1336 return np.array(precon_forces)
1338 def average_norm(self, i, j, dx):
1339 """Average norm between images i and j
1341 Args:
1342 i (int): left image
1343 j (int): right image
1344 dx (array): vector
1346 Returns:
1347 norm: norm of vector wrt average of precons at i and j
1348 """
1349 return np.sqrt(0.5 * (self.precon[i].dot(dx, dx) +
1350 self.precon[j].dot(dx, dx)))
1352 def get_tangent(self, i):
1353 """Normalised tangent vector at image i
1355 Args:
1356 i (int): image of interest
1358 Returns:
1359 tangent: tangent vector, normalised with appropriate precon norm
1360 """
1361 tangent = self.spline.dx_ds(self.spline.s[i])
1362 tangent /= self.precon[i].norm(tangent)
1363 return tangent.reshape(-1, 3)
1365 def get_residual(self, i, imgforce):
1366 # residuals computed according to eq. 11 in the paper
1367 P_dot_imgforce = self.precon[i].Pdot(imgforce.reshape(-1))
1368 return np.linalg.norm(P_dot_imgforce, np.inf)
1370 def get_spring_force(self, i, k1, k2, tangent):
1371 """Spring force on image
1373 Args:
1374 i (int): image of interest
1375 k1 (float): spring constant for left spring
1376 k2 (float): spring constant for right spring
1377 tangent (array): tangent vector, shape (natoms, 3)
1379 Returns:
1380 eta: NEB spring forces, shape (natoms, 3)
1381 """
1382 # Definition following Eq. 9 in Paper IV
1383 nimages = len(self.images)
1384 k = 0.5 * (k1 + k2) / (nimages ** 2)
1385 curvature = self.spline.d2x_ds2(self.spline.s[i]).reshape(-1, 3)
1386 # complete Eq. 9 by including the spring force
1387 eta = k * self.precon[i].vdot(curvature, tangent) * tangent
1388 return eta
1390 def get_coordinates(self, positions=None):
1391 """Compute displacements wrt appropriate precon metric for each image
1393 Args:
1394 positions (list or array, optional) - images positions.
1395 Shape either (nimages * natoms, 3) or ((nimages-2)*natoms, 3)
1397 Returns:
1398 s : array shape (nimages,), reaction coordinates, in range [0, 1]
1399 x : array shape (nimages, 3 * natoms), flat displacement vectors
1400 """
1401 nimages = len(self.precon)
1402 natoms = len(self.images[0])
1403 d_P = np.zeros(nimages)
1404 x = np.zeros((nimages, 3 * natoms)) # flattened positions
1405 if positions is None:
1406 positions = [image.positions for image in self.images]
1407 elif isinstance(positions, np.ndarray) and len(positions.shape) == 2:
1408 positions = positions.reshape(-1, natoms, 3)
1409 positions = [positions[i, :, :] for i in range(len(positions))]
1410 if len(positions) == len(self.images) - 2:
1411 # prepend and append the non-moving images
1412 positions = ([self.images[0].positions] + positions +
1413 [self.images[-1].positions])
1414 assert len(positions) == len(self.images)
1416 x[0, :] = positions[0].reshape(-1)
1417 for i in range(1, nimages):
1418 x[i, :] = positions[i].reshape(-1)
1419 dx, _ = find_mic(positions[i] - positions[i - 1],
1420 self.images[i - 1].cell,
1421 self.images[i - 1].pbc)
1422 dx = dx.reshape(-1)
1423 d_P[i] = self.average_norm(i, i - 1, dx)
1425 s = d_P.cumsum() / d_P.sum() # Eq. A1 in paper IV
1426 return s, x
1428 def spline_fit(self, positions=None):
1429 """Fit 3 * natoms cubic splines as a function of reaction coordinate
1431 Returns:
1432 fit : :class:`ase.optimize.precon.SplineFit` object
1433 """
1434 s, x = self.get_coordinates(positions)
1435 return SplineFit(s, x)
1437 @property
1438 def spline(self):
1439 s, x = self.get_coordinates()
1440 if self._spline and (np.abs(s - self._old_s).max() < 1e-6 and
1441 np.abs(x - self._old_x).max() < 1e-6):
1442 return self._spline
1444 self._spline = self.spline_fit()
1445 self._old_s = s
1446 self._old_x = x
1447 return self._spline