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