Coverage for /builds/debichem-team/python-ase/ase/optimize/precon/precon.py: 85.19%
675 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +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 (
23 estimate_nearest_neighbour_distance,
24 get_neighbours,
25)
26from ase.utils import longsum, tokenize_version
28try:
29 import pyamg
30except ImportError:
31 have_pyamg = False
32else:
33 have_pyamg = True
36def create_pyamg_solver(P, max_levels=15):
37 if not have_pyamg:
38 raise RuntimeError(
39 'pyamg not available: install with `pip install pyamg`')
40 filter_key = 'filter'
41 if tokenize_version(pyamg.__version__) >= tokenize_version('4.2.1'):
42 filter_key = 'filter_entries'
43 return pyamg.smoothed_aggregation_solver(
44 P, B=None,
45 strength=('symmetric', {'theta': 0.0}),
46 smooth=(
47 'jacobi', {filter_key: True, 'weighting': 'local'}),
48 improve_candidates=([('block_gauss_seidel',
49 {'sweep': 'symmetric', 'iterations': 4})] +
50 [None] * (max_levels - 1)),
51 aggregate='standard',
52 presmoother=('block_gauss_seidel',
53 {'sweep': 'symmetric', 'iterations': 1}),
54 postsmoother=('block_gauss_seidel',
55 {'sweep': 'symmetric', 'iterations': 1}),
56 max_levels=max_levels,
57 max_coarse=300,
58 coarse_solver='pinv')
61THz = 1e12 * 1. / units.s
64class Precon(ABC):
66 @abstractmethod
67 def make_precon(self, atoms, reinitialize=None):
68 """
69 Create a preconditioner matrix based on the passed set of atoms.
71 Creates a general-purpose preconditioner for use with optimization
72 algorithms, based on examining distances between pairs of atoms in the
73 lattice. The matrix will be stored in the attribute self.P and
74 returned.
76 Args:
77 atoms: the Atoms object used to create the preconditioner.
79 reinitialize: if True, parameters of the preconditioner
80 will be recalculated before the preconditioner matrix is
81 created. If False, they will be calculated only when they
82 do not currently have a value (ie, the first time this
83 function is called).
85 Returns:
86 P: A sparse scipy csr_matrix. BE AWARE that using
87 numpy.dot() with sparse matrices will result in
88 errors/incorrect results - use the .dot method directly
89 on the matrix instead.
90 """
91 ...
93 @abstractmethod
94 def Pdot(self, x):
95 """
96 Return the result of applying P to a vector x
97 """
98 ...
100 def dot(self, x, y):
101 """
102 Return the preconditioned dot product <P x, y>
104 Uses 128-bit floating point math for vector dot products
105 """
106 return longsum(self.Pdot(x) * y)
108 def norm(self, x):
109 """
110 Return the P-norm of x, where |x|_P = sqrt(<Px, x>)
111 """
112 return np.sqrt(self.dot(x, x))
114 def vdot(self, x, y):
115 return self.dot(x.reshape(-1),
116 y.reshape(-1))
118 @abstractmethod
119 def solve(self, x):
120 """
121 Solve the (sparse) linear system P x = y and return y
122 """
123 ...
125 def apply(self, forces, atoms):
126 """
127 Convenience wrapper that combines make_precon() and solve()
129 Parameters
130 ----------
131 forces: array
132 (len(atoms)*3) array of input forces
133 atoms: ase.atoms.Atoms
135 Returns
136 -------
137 precon_forces: array
138 (len(atoms), 3) array of preconditioned forces
139 residual: float
140 inf-norm of original forces, i.e. maximum absolute force
141 """
142 self.make_precon(atoms)
143 residual = np.linalg.norm(forces, np.inf)
144 precon_forces = self.solve(forces)
145 return precon_forces, residual
147 @abstractmethod
148 def copy(self):
149 ...
151 @abstractmethod
152 def asarray(self):
153 """
154 Array representation of preconditioner, as a dense matrix
155 """
156 ...
159class Logfile:
160 def __init__(self, logfile=None):
161 if isinstance(logfile, str):
162 if logfile == "-":
163 logfile = sys.stdout
164 else:
165 logfile = open(logfile, "a")
166 self.logfile = logfile
168 def write(self, *args):
169 if self.logfile is None:
170 return
171 self.logfile.write(*args)
174class SparsePrecon(Precon):
175 def __init__(self, r_cut=None, r_NN=None,
176 mu=None, mu_c=None,
177 dim=3, c_stab=0.1, force_stab=False,
178 reinitialize=False, array_convention='C',
179 solver="auto", solve_tol=1e-8,
180 apply_positions=True, apply_cell=True,
181 estimate_mu_eigmode=False, logfile=None, rng=None,
182 neighbour_list=neighbor_list):
183 """Initialise a preconditioner object based on passed parameters.
185 Parameters:
186 r_cut: float. This is a cut-off radius. The preconditioner matrix
187 will be created by considering pairs of atoms that are within a
188 distance r_cut of each other. For a regular lattice, this is
189 usually taken somewhere between the first- and second-nearest
190 neighbour distance. If r_cut is not provided, default is
191 2 * r_NN (see below)
192 r_NN: nearest neighbour distance. If not provided, this is
193 calculated
194 from input structure.
195 mu: float
196 energy scale for position degreees of freedom. If `None`, mu
197 is precomputed using finite difference derivatives.
198 mu_c: float
199 energy scale for cell degreees of freedom. Also precomputed
200 if None.
201 estimate_mu_eigmode:
202 If True, estimates mu based on the lowest eigenmodes of
203 unstabilised preconditioner. If False it uses the sine based
204 approach.
205 dim: int; dimensions of the problem
206 c_stab: float. The diagonal of the preconditioner matrix will have
207 a stabilisation constant added, which will be the value of
208 c_stab times mu.
209 force_stab:
210 If True, always add the stabilisation to diagnonal, regardless
211 of the presence of fixed atoms.
212 reinitialize: if True, the value of mu will be recalculated when
213 self.make_precon is called. This can be overridden in specific
214 cases with reinitialize argument in self.make_precon. If it
215 is set to True here, the value passed for mu will be
216 irrelevant unless reinitialize is set to False the first time
217 make_precon is called.
218 array_convention: Either 'C' or 'F' for Fortran; this will change
219 the preconditioner to reflect the ordering of the indices in
220 the vector it will operate on. The C convention assumes the
221 vector will be arranged atom-by-atom (ie [x1, y1, z1, x2, ...])
222 while the F convention assumes it will be arranged component
223 by component (ie [x1, x2, ..., y1, y2, ...]).
224 solver: One of "auto", "direct" or "pyamg", specifying whether to
225 use a direct sparse solver or PyAMG to solve P x = y.
226 Default is "auto" which uses PyAMG if available, falling
227 back to sparse solver if not. solve_tol: tolerance used for
228 PyAMG sparse linear solver, if available.
229 apply_positions: bool
230 if True, apply preconditioner to position DoF
231 apply_cell: bool
232 if True, apply preconditioner to cell DoF
233 logfile: file object or str
234 If *logfile* is a string, a file with that name will be opened.
235 Use '-' for stdout.
236 rng: None or np.random.RandomState instance
237 Random number generator to use for initialising pyamg solver
238 neighbor_list: function (optional). Optionally replace the built-in
239 ASE neighbour list with an alternative with the same call
240 signature, e.g. `matscipy.neighbours.neighbour_list`.
242 Raises:
243 ValueError for problem with arguments
245 """
246 self.r_NN = r_NN
247 self.r_cut = r_cut
248 self.mu = mu
249 self.mu_c = mu_c
250 self.estimate_mu_eigmode = estimate_mu_eigmode
251 self.c_stab = c_stab
252 self.force_stab = force_stab
253 self.array_convention = array_convention
254 self.reinitialize = reinitialize
255 self.P = None
256 self.old_positions = None
258 use_pyamg = False
259 if solver == "auto":
260 use_pyamg = have_pyamg
261 elif solver == "direct":
262 use_pyamg = False
263 elif solver == "pyamg":
264 if not have_pyamg:
265 raise RuntimeError("solver='pyamg', PyAMG can't be imported!")
266 use_pyamg = True
267 else:
268 raise ValueError('unknown solver - '
269 'should be "auto", "direct" or "pyamg"')
271 self.use_pyamg = use_pyamg
272 self.solve_tol = solve_tol
273 self.apply_positions = apply_positions
274 self.apply_cell = apply_cell
276 if dim < 1:
277 raise ValueError('Dimension must be at least 1')
278 self.dim = dim
279 self.logfile = Logfile(logfile)
281 if rng is None:
282 rng = np.random.RandomState()
283 self.rng = rng
285 self.neighbor_list = neighbor_list
287 def copy(self):
288 return copy.deepcopy(self)
290 def Pdot(self, x):
291 return self.P.dot(x)
293 def solve(self, x):
294 start_time = time.time()
295 if self.use_pyamg and have_pyamg:
296 y = self.ml.solve(x, x0=self.rng.random(self.P.shape[0]),
297 tol=self.solve_tol,
298 accel='cg',
299 maxiter=300,
300 cycle='W')
301 else:
302 y = spsolve(self.P, x)
303 self.logfile.write(
304 f'--- Precon applied in {(time.time() - start_time)} seconds ---\n')
305 return y
307 def estimate_mu(self, atoms, H=None):
308 r"""Estimate optimal preconditioner coefficient \mu
310 \mu is estimated from a numerical solution of
312 [dE(p+v) - dE(p)] \cdot v = \mu < P1 v, v >
314 with perturbation
316 v(x,y,z) = H P_lowest_nonzero_eigvec(x, y, z)
318 or
320 v(x,y,z) = H (sin(x / Lx), sin(y / Ly), sin(z / Lz))
322 After the optimal \mu is found, self.mu will be set to its value.
324 If `atoms` is an instance of Filter an additional \mu_c
325 will be computed for the cell degrees of freedom .
327 Args:
328 atoms: Atoms object for initial system
330 H: 3x3 array or None
331 Magnitude of deformation to apply.
332 Default is 1e-2*rNN*np.eye(3)
334 Returns:
335 mu : float
336 mu_c : float or None
337 """
338 logfile = self.logfile
340 if self.dim != 3:
341 raise ValueError('Automatic calculation of mu only possible for '
342 'three-dimensional preconditioners. Try setting '
343 'mu manually instead.')
345 if self.r_NN is None:
346 self.r_NN = estimate_nearest_neighbour_distance(atoms,
347 self.neighbor_list)
349 # deformation matrix, default is diagonal
350 if H is None:
351 H = 1e-2 * self.r_NN * np.eye(3)
353 # compute perturbation
354 p = atoms.get_positions()
356 if self.estimate_mu_eigmode:
357 self.mu = 1.0
358 self.mu_c = 1.0
359 c_stab = self.c_stab
360 self.c_stab = 0.0
362 if isinstance(atoms, Filter):
363 n = len(atoms.atoms)
364 else:
365 n = len(atoms)
366 self._make_sparse_precon(atoms, initial_assembly=True)
367 self.P = self.P[:3 * n, :3 * n]
368 eigvals, eigvecs = sparse.linalg.eigsh(self.P, k=4, which='SM')
370 logfile.write('estimate_mu(): lowest 4 eigvals = %f %f %f %f\n' %
371 (eigvals[0], eigvals[1], eigvals[2], eigvals[3]))
372 # check eigenvalues
373 if any(eigvals[0:3] > 1e-6):
374 raise ValueError('First 3 eigenvalues of preconditioner matrix'
375 'do not correspond to translational modes.')
376 elif eigvals[3] < 1e-6:
377 raise ValueError('Fourth smallest eigenvalue of '
378 'preconditioner matrix '
379 'is too small, increase r_cut.')
381 x = np.zeros(n)
382 for i in range(n):
383 x[i] = eigvecs[:, 3][3 * i]
384 x = x / np.linalg.norm(x)
385 if x[0] < 0:
386 x = -x
388 v = np.zeros(3 * len(atoms))
389 for i in range(n):
390 v[3 * i] = x[i]
391 v[3 * i + 1] = x[i]
392 v[3 * i + 2] = x[i]
393 v = v / np.linalg.norm(v)
394 v = v.reshape((-1, 3))
396 self.c_stab = c_stab
397 else:
398 Lx, Ly, Lz = (p[:, i].max() - p[:, i].min() for i in range(3))
399 logfile.write('estimate_mu(): Lx=%.1f Ly=%.1f Lz=%.1f\n' %
400 (Lx, Ly, Lz))
402 x, y, z = p.T
403 # sine_vr = [np.sin(x/Lx), np.sin(y/Ly), np.sin(z/Lz)], but we need
404 # to take into account the possibility that one of Lx/Ly/Lz is
405 # zero.
406 sine_vr = [x, y, z]
408 for i, L in enumerate([Lx, Ly, Lz]):
409 if L == 0:
410 warnings.warn(
411 'Cell length L[%d] == 0. Setting H[%d,%d] = 0.' %
412 (i, i, i))
413 H[i, i] = 0.0
414 else:
415 sine_vr[i] = np.sin(sine_vr[i] / L)
417 v = np.dot(H, sine_vr).T
419 natoms = len(atoms)
420 if isinstance(atoms, Filter):
421 natoms = len(atoms.atoms)
422 eps = H / self.r_NN
423 v[natoms:, :] = eps
425 v1 = v.reshape(-1)
427 # compute LHS
428 dE_p = -atoms.get_forces().reshape(-1)
429 atoms_v = atoms.copy()
430 atoms_v.calc = atoms.calc
431 if isinstance(atoms, Filter):
432 atoms_v = atoms.__class__(atoms_v)
433 if hasattr(atoms, 'constant_volume'):
434 atoms_v.constant_volume = atoms.constant_volume
435 atoms_v.set_positions(p + v)
436 dE_p_plus_v = -atoms_v.get_forces().reshape(-1)
438 # compute left hand side
439 LHS = (dE_p_plus_v - dE_p) * v1
441 # assemble P with \mu = 1
442 self.mu = 1.0
443 self.mu_c = 1.0
445 self._make_sparse_precon(atoms, initial_assembly=True)
447 # compute right hand side
448 RHS = self.P.dot(v1) * v1
450 # use partial sums to compute separate mu for positions and cell DoFs
451 self.mu = longsum(LHS[:3 * natoms]) / longsum(RHS[:3 * natoms])
452 if self.mu < 1.0:
453 logfile.write('estimate_mu(): mu (%.3f) < 1.0, '
454 'capping at mu=1.0' % self.mu)
455 self.mu = 1.0
457 if isinstance(atoms, Filter):
458 self.mu_c = longsum(LHS[3 * natoms:]) / longsum(RHS[3 * natoms:])
459 if self.mu_c < 1.0:
460 logfile.write('estimate_mu(): mu_c (%.3f) < 1.0, '
461 'capping at mu_c=1.0\n' % self.mu_c)
462 self.mu_c = 1.0
464 logfile.write('estimate_mu(): mu=%r, mu_c=%r\n' %
465 (self.mu, self.mu_c))
467 self.P = None # force a rebuild with new mu (there may be fixed atoms)
468 return (self.mu, self.mu_c)
470 def asarray(self):
471 return np.array(self.P.todense())
473 def one_dim_to_ndim(self, csc_P, N):
474 """
475 Expand an N x N precon matrix to self.dim*N x self.dim * N
477 Args:
478 csc_P (sparse matrix): N x N sparse matrix, in CSC format
479 """
480 start_time = time.time()
481 if self.dim == 1:
482 P = csc_P
483 elif self.array_convention == 'F':
484 csc_P = csc_P.tocsr()
485 P = csc_P
486 for _ in range(self.dim - 1):
487 P = sparse.block_diag((P, csc_P)).tocsr()
488 else:
489 # convert back to triplet and read the arrays
490 csc_P = csc_P.tocoo()
491 i = csc_P.row * self.dim
492 j = csc_P.col * self.dim
493 z = csc_P.data
495 # N-dimensionalise, interlaced coordinates
496 I = np.hstack([i + d for d in range(self.dim)])
497 J = np.hstack([j + d for d in range(self.dim)])
498 Z = np.hstack([z for _ in range(self.dim)])
499 P = sparse.csc_matrix((Z, (I, J)),
500 shape=(self.dim * N, self.dim * N))
501 P = P.tocsr()
502 self.logfile.write(
503 f'--- N-dim precon created in {(time.time() - start_time)} s ---\n')
504 return P
506 def create_solver(self):
507 if self.use_pyamg and have_pyamg:
508 start_time = time.time()
509 self.ml = create_pyamg_solver(self.P)
510 self.logfile.write(
511 f'--- multi grid solver created in {(time.time() - start_time)}'
512 ' ---\n')
515class SparseCoeffPrecon(SparsePrecon):
516 def _make_sparse_precon(self, atoms, initial_assembly=False,
517 force_stab=False):
518 """Create a sparse preconditioner matrix based on the passed atoms.
520 Creates a general-purpose preconditioner for use with optimization
521 algorithms, based on examining distances between pairs of atoms in the
522 lattice. The matrix will be stored in the attribute self.P and
523 returned. Note that this function will use self.mu, whatever it is.
525 Args:
526 atoms: the Atoms object used to create the preconditioner.
528 Returns:
529 A scipy.sparse.csr_matrix object, representing a d*N by d*N matrix
530 (where N is the number of atoms, and d is the value of self.dim).
531 BE AWARE that using numpy.dot() with this object will result in
532 errors/incorrect results - use the .dot method directly on the
533 sparse matrix instead.
535 """
536 logfile = self.logfile
537 logfile.write('creating sparse precon: initial_assembly=%r, '
538 'force_stab=%r, apply_positions=%r, apply_cell=%r\n' %
539 (initial_assembly, force_stab, self.apply_positions,
540 self.apply_cell))
542 N = len(atoms)
543 diag_i = np.arange(N, dtype=int)
544 start_time = time.time()
545 if self.apply_positions:
546 # compute neighbour list
547 i, j, rij, fixed_atoms = get_neighbours(
548 atoms, self.r_cut,
549 neighbor_list=self.neighbor_list)
550 logfile.write(
551 f'--- neighbour list created in {(time.time() - start_time)} s '
552 '--- \n')
554 # compute entries in triplet format: without the constraints
555 start_time = time.time()
556 coeff = self.get_coeff(rij)
557 diag_coeff = np.bincount(i, -coeff, minlength=N).astype(np.float64)
558 if force_stab or len(fixed_atoms) == 0:
559 logfile.write('adding stabilisation to precon')
560 diag_coeff += self.mu * self.c_stab
561 else:
562 diag_coeff = np.ones(N)
564 # precon is mu_c * identity for cell DoF
565 if isinstance(atoms, Filter):
566 if self.apply_cell:
567 diag_coeff[-3:] = self.mu_c
568 else:
569 diag_coeff[-3:] = 1.0
570 logfile.write(
571 f'--- computed triplet format in {(time.time() - start_time)} s '
572 '---\n')
574 if self.apply_positions and not initial_assembly:
575 # apply the constraints
576 start_time = time.time()
577 mask = np.ones(N)
578 mask[fixed_atoms] = 0.0
579 coeff *= mask[i] * mask[j]
580 diag_coeff[fixed_atoms] = 1.0
581 logfile.write(
582 f'--- applied fixed_atoms in {(time.time() - start_time)} s '
583 '---\n')
585 if self.apply_positions:
586 # remove zeros
587 start_time = time.time()
588 inz = np.nonzero(coeff)
589 i = np.hstack((i[inz], diag_i))
590 j = np.hstack((j[inz], diag_i))
591 coeff = np.hstack((coeff[inz], diag_coeff))
592 logfile.write(
593 f'--- remove zeros in {(time.time() - start_time)} s '
594 '---\n')
595 else:
596 i = diag_i
597 j = diag_i
598 coeff = diag_coeff
600 # create an N x N precon matrix in compressed sparse column (CSC) format
601 start_time = time.time()
602 csc_P = sparse.csc_matrix((coeff, (i, j)), shape=(N, N))
603 logfile.write(
604 f'--- created CSC matrix in {(time.time() - start_time)} s '
605 '---\n')
607 self.P = self.one_dim_to_ndim(csc_P, N)
608 self.create_solver()
610 def make_precon(self, atoms, reinitialize=None):
611 if self.r_NN is None:
612 self.r_NN = estimate_nearest_neighbour_distance(atoms,
613 self.neighbor_list)
615 if self.r_cut is None:
616 # This is the first time this function has been called, and no
617 # cutoff radius has been specified, so calculate it automatically.
618 self.r_cut = 2.0 * self.r_NN
619 elif self.r_cut < self.r_NN:
620 warning = ('WARNING: r_cut (%.2f) < r_NN (%.2f), '
621 'increasing to 1.1*r_NN = %.2f' % (self.r_cut,
622 self.r_NN,
623 1.1 * self.r_NN))
624 warnings.warn(warning)
625 self.r_cut = 1.1 * self.r_NN
627 if reinitialize is None:
628 # The caller has not specified whether or not to recalculate mu,
629 # so the Precon's setting is used.
630 reinitialize = self.reinitialize
632 if self.mu is None:
633 # Regardless of what the caller has specified, if we don't
634 # currently have a value of mu, then we need one.
635 reinitialize = True
637 if reinitialize:
638 self.estimate_mu(atoms)
640 if self.P is not None:
641 real_atoms = atoms
642 if isinstance(atoms, Filter):
643 real_atoms = atoms.atoms
644 if self.old_positions is None:
645 self.old_positions = real_atoms.positions
646 displacement, _ = find_mic(real_atoms.positions -
647 self.old_positions,
648 real_atoms.cell, real_atoms.pbc)
649 self.old_positions = real_atoms.get_positions()
650 max_abs_displacement = abs(displacement).max()
651 self.logfile.write('max(abs(displacements)) = %.2f A (%.2f r_NN)' %
652 (max_abs_displacement,
653 max_abs_displacement / self.r_NN))
654 if max_abs_displacement < 0.5 * self.r_NN:
655 return
657 start_time = time.time()
658 self._make_sparse_precon(atoms, force_stab=self.force_stab)
659 self.logfile.write(
660 f'--- Precon created in {(time.time() - start_time)} seconds '
661 '--- \n')
663 @abstractmethod
664 def get_coeff(self, r):
665 ...
668class Pfrommer(Precon):
669 """
670 Use initial guess for inverse Hessian from Pfrommer et al. as a
671 simple preconditioner
673 J. Comput. Phys. vol 131 p233-240 (1997)
674 """
676 def __init__(self, bulk_modulus=500 * units.GPa, phonon_frequency=50 * THz,
677 apply_positions=True, apply_cell=True):
678 """
679 Default bulk modulus is 500 GPa and default phonon frequency is 50 THz
680 """
681 self.bulk_modulus = bulk_modulus
682 self.phonon_frequency = phonon_frequency
683 self.apply_positions = apply_positions
684 self.apply_cell = apply_cell
685 self.H0 = None
687 def make_precon(self, atoms, reinitialize=None):
688 if self.H0 is not None:
689 # only build H0 on first call
690 return
692 variable_cell = False
693 if isinstance(atoms, Filter):
694 variable_cell = True
695 atoms = atoms.atoms
697 # position DoF
698 omega = self.phonon_frequency
699 mass = atoms.get_masses().mean()
700 block = np.eye(3) / (mass * omega**2)
701 blocks = [block] * len(atoms)
703 # cell DoF
704 if variable_cell:
705 coeff = 1.0
706 if self.apply_cell:
707 coeff = 1.0 / (3 * self.bulk_modulus)
708 blocks.append(np.diag([coeff] * 9))
710 self.H0 = sparse.block_diag(blocks, format='csr')
711 return
713 def Pdot(self, x):
714 return self.H0.solve(x)
716 def solve(self, x):
717 y = self.H0.dot(x)
718 return y
720 def copy(self):
721 return Pfrommer(self.bulk_modulus,
722 self.phonon_frequency,
723 self.apply_positions,
724 self.apply_cell)
726 def asarray(self):
727 return np.array(self.H0.todense())
730class IdentityPrecon(Precon):
731 """
732 Dummy preconditioner which does not modify forces
733 """
735 def make_precon(self, atoms, reinitialize=None):
736 self.atoms = atoms
738 def Pdot(self, x):
739 return x
741 def solve(self, x):
742 return x
744 def copy(self):
745 return IdentityPrecon()
747 def asarray(self):
748 return np.eye(3 * len(self.atoms))
751class C1(SparseCoeffPrecon):
752 """Creates matrix by inserting a constant whenever r_ij is less than r_cut.
753 """
755 def __init__(self, r_cut=None, mu=None, mu_c=None, dim=3, c_stab=0.1,
756 force_stab=False,
757 reinitialize=False, array_convention='C',
758 solver="auto", solve_tol=1e-9,
759 apply_positions=True, apply_cell=True, logfile=None):
760 super().__init__(r_cut=r_cut, mu=mu, mu_c=mu_c,
761 dim=dim, c_stab=c_stab,
762 force_stab=force_stab,
763 reinitialize=reinitialize,
764 array_convention=array_convention,
765 solver=solver, solve_tol=solve_tol,
766 apply_positions=apply_positions,
767 apply_cell=apply_cell,
768 logfile=logfile)
770 def get_coeff(self, r):
771 return -self.mu * np.ones_like(r)
774class Exp(SparseCoeffPrecon):
775 """Creates matrix with values decreasing exponentially with distance.
776 """
778 def __init__(self, A=3.0, r_cut=None, r_NN=None, mu=None, mu_c=None,
779 dim=3, c_stab=0.1,
780 force_stab=False, reinitialize=False, array_convention='C',
781 solver="auto", solve_tol=1e-9,
782 apply_positions=True, apply_cell=True,
783 estimate_mu_eigmode=False, logfile=None):
784 """
785 Initialise an Exp preconditioner with given parameters.
787 Args:
788 r_cut, mu, c_stab, dim, sparse, reinitialize, array_convention: see
789 precon.__init__()
790 A: coefficient in exp(-A*r/r_NN). Default is A=3.0.
791 """
792 super().__init__(r_cut=r_cut, r_NN=r_NN,
793 mu=mu, mu_c=mu_c, dim=dim, c_stab=c_stab,
794 force_stab=force_stab,
795 reinitialize=reinitialize,
796 array_convention=array_convention,
797 solver=solver,
798 solve_tol=solve_tol,
799 apply_positions=apply_positions,
800 apply_cell=apply_cell,
801 estimate_mu_eigmode=estimate_mu_eigmode,
802 logfile=logfile)
804 self.A = A
806 def get_coeff(self, r):
807 return -self.mu * np.exp(-self.A * (r / self.r_NN - 1))
810def ij_to_x(i, j):
811 x = [3 * i, 3 * i + 1, 3 * i + 2,
812 3 * j, 3 * j + 1, 3 * j + 2]
813 return x
816def ijk_to_x(i, j, k):
817 x = [3 * i, 3 * i + 1, 3 * i + 2,
818 3 * j, 3 * j + 1, 3 * j + 2,
819 3 * k, 3 * k + 1, 3 * k + 2]
820 return x
823def ijkl_to_x(i, j, k, l):
824 x = [3 * i, 3 * i + 1, 3 * i + 2,
825 3 * j, 3 * j + 1, 3 * j + 2,
826 3 * k, 3 * k + 1, 3 * k + 2,
827 3 * l, 3 * l + 1, 3 * l + 2]
828 return x
831def apply_fixed(atoms, P):
832 fixed_atoms = []
833 for constraint in atoms.constraints:
834 if isinstance(constraint, FixAtoms):
835 fixed_atoms.extend(list(constraint.index))
836 else:
837 raise TypeError(
838 'only FixAtoms constraints are supported by Precon class')
839 if len(fixed_atoms) != 0:
840 P = P.tolil()
841 for i in fixed_atoms:
842 P[i, :] = 0.0
843 P[:, i] = 0.0
844 P[i, i] = 1.0
845 return P
848class FF(SparsePrecon):
849 """Creates matrix using morse/bond/angle/dihedral force field parameters.
850 """
852 def __init__(self, dim=3, c_stab=0.1, force_stab=False,
853 array_convention='C', solver="auto", solve_tol=1e-9,
854 apply_positions=True, apply_cell=True,
855 hessian='spectral', morses=None, bonds=None, angles=None,
856 dihedrals=None, logfile=None):
857 """Initialise an FF preconditioner with given parameters.
859 Args:
860 dim, c_stab, force_stab, array_convention, use_pyamg, solve_tol:
861 see SparsePrecon.__init__()
862 morses: Morse instance
863 bonds: Bond instance
864 angles: Angle instance
865 dihedrals: Dihedral instance
866 """
868 if (morses is None and bonds is None and angles is None and
869 dihedrals is None):
870 raise ImportError(
871 'At least one of morses, bonds, angles or dihedrals must be '
872 'defined!')
874 super().__init__(dim=dim, c_stab=c_stab,
875 force_stab=force_stab,
876 array_convention=array_convention,
877 solver=solver,
878 solve_tol=solve_tol,
879 apply_positions=apply_positions,
880 apply_cell=apply_cell, logfile=logfile)
882 self.hessian = hessian
883 self.morses = morses
884 self.bonds = bonds
885 self.angles = angles
886 self.dihedrals = dihedrals
888 def make_precon(self, atoms, reinitialize=None):
889 start_time = time.time()
890 self._make_sparse_precon(atoms, force_stab=self.force_stab)
891 self.logfile.write(
892 f'--- Precon created in {(time.time() - start_time)} seconds '
893 '---\n')
895 def add_morse(self, morse, atoms, row, col, data, conn=None):
896 if self.hessian == 'reduced':
897 i, j, Hx = ff.get_morse_potential_reduced_hessian(
898 atoms, morse)
899 elif self.hessian == 'spectral':
900 i, j, Hx = ff.get_morse_potential_hessian(
901 atoms, morse, spectral=True)
902 else:
903 raise NotImplementedError('Not implemented hessian')
904 x = ij_to_x(i, j)
905 row.extend(np.repeat(x, 6))
906 col.extend(np.tile(x, 6))
907 data.extend(Hx.flatten())
908 if conn is not None:
909 conn[i, j] = True
910 conn[j, i] = True
912 def add_bond(self, bond, atoms, row, col, data, conn=None):
913 if self.hessian == 'reduced':
914 i, j, Hx = ff.get_bond_potential_reduced_hessian(
915 atoms, bond, self.morses)
916 elif self.hessian == 'spectral':
917 i, j, Hx = ff.get_bond_potential_hessian(
918 atoms, bond, self.morses, spectral=True)
919 else:
920 raise NotImplementedError('Not implemented hessian')
921 x = ij_to_x(i, j)
922 row.extend(np.repeat(x, 6))
923 col.extend(np.tile(x, 6))
924 data.extend(Hx.flatten())
925 if conn is not None:
926 conn[i, j] = True
927 conn[j, i] = True
929 def add_angle(self, angle, atoms, row, col, data, conn=None):
930 if self.hessian == 'reduced':
931 i, j, k, Hx = ff.get_angle_potential_reduced_hessian(
932 atoms, angle, self.morses)
933 elif self.hessian == 'spectral':
934 i, j, k, Hx = ff.get_angle_potential_hessian(
935 atoms, angle, self.morses, spectral=True)
936 else:
937 raise NotImplementedError('Not implemented hessian')
938 x = ijk_to_x(i, j, k)
939 row.extend(np.repeat(x, 9))
940 col.extend(np.tile(x, 9))
941 data.extend(Hx.flatten())
942 if conn is not None:
943 conn[i, j] = conn[i, k] = conn[j, k] = True
944 conn[j, i] = conn[k, i] = conn[k, j] = True
946 def add_dihedral(self, dihedral, atoms, row, col, data, conn=None):
947 if self.hessian == 'reduced':
948 i, j, k, l, Hx = \
949 ff.get_dihedral_potential_reduced_hessian(
950 atoms, dihedral, self.morses)
951 elif self.hessian == 'spectral':
952 i, j, k, l, Hx = ff.get_dihedral_potential_hessian(
953 atoms, dihedral, self.morses, spectral=True)
954 else:
955 raise NotImplementedError('Not implemented hessian')
956 x = ijkl_to_x(i, j, k, l)
957 row.extend(np.repeat(x, 12))
958 col.extend(np.tile(x, 12))
959 data.extend(Hx.flatten())
960 if conn is not None:
961 conn[i, j] = conn[i, k] = conn[i, l] = conn[
962 j, k] = conn[j, l] = conn[k, l] = True
963 conn[j, i] = conn[k, i] = conn[l, i] = conn[
964 k, j] = conn[l, j] = conn[l, k] = True
966 def _make_sparse_precon(self, atoms, initial_assembly=False,
967 force_stab=False):
968 N = len(atoms)
970 row = []
971 col = []
972 data = []
974 if self.morses is not None:
975 for morse in self.morses:
976 self.add_morse(morse, atoms, row, col, data)
978 if self.bonds is not None:
979 for bond in self.bonds:
980 self.add_bond(bond, atoms, row, col, data)
982 if self.angles is not None:
983 for angle in self.angles:
984 self.add_angle(angle, atoms, row, col, data)
986 if self.dihedrals is not None:
987 for dihedral in self.dihedrals:
988 self.add_dihedral(dihedral, atoms, row, col, data)
990 # add the diagonal
991 row.extend(range(self.dim * N))
992 col.extend(range(self.dim * N))
993 data.extend([self.c_stab] * self.dim * N)
995 # create the matrix
996 start_time = time.time()
997 self.P = sparse.csc_matrix(
998 (data, (row, col)), shape=(self.dim * N, self.dim * N))
999 self.logfile.write(
1000 f'--- created CSC matrix in {(time.time() - start_time)} s ---\n')
1002 self.P = apply_fixed(atoms, self.P)
1003 self.P = self.P.tocsr()
1004 self.logfile.write(
1005 f'--- N-dim precon created in {(time.time() - start_time)} s ---\n')
1006 self.create_solver()
1009class Exp_FF(Exp, FF):
1010 """Creates matrix with values decreasing exponentially with distance.
1011 """
1013 def __init__(self, A=3.0, r_cut=None, r_NN=None, mu=None, mu_c=None,
1014 dim=3, c_stab=0.1,
1015 force_stab=False, reinitialize=False, array_convention='C',
1016 solver="auto", solve_tol=1e-9,
1017 apply_positions=True, apply_cell=True,
1018 estimate_mu_eigmode=False,
1019 hessian='spectral', morses=None, bonds=None, angles=None,
1020 dihedrals=None, logfile=None):
1021 """Initialise an Exp+FF preconditioner with given parameters.
1023 Args:
1024 r_cut, mu, c_stab, dim, reinitialize, array_convention: see
1025 precon.__init__()
1026 A: coefficient in exp(-A*r/r_NN). Default is A=3.0.
1027 """
1028 if (morses is None and bonds is None and angles is None and
1029 dihedrals is None):
1030 raise ImportError(
1031 'At least one of morses, bonds, angles or dihedrals must '
1032 'be defined!')
1034 SparsePrecon.__init__(self, r_cut=r_cut, r_NN=r_NN,
1035 mu=mu, mu_c=mu_c, dim=dim, c_stab=c_stab,
1036 force_stab=force_stab,
1037 reinitialize=reinitialize,
1038 array_convention=array_convention,
1039 solver=solver,
1040 solve_tol=solve_tol,
1041 apply_positions=apply_positions,
1042 apply_cell=apply_cell,
1043 estimate_mu_eigmode=estimate_mu_eigmode,
1044 logfile=logfile)
1046 self.A = A
1047 self.hessian = hessian
1048 self.morses = morses
1049 self.bonds = bonds
1050 self.angles = angles
1051 self.dihedrals = dihedrals
1053 def make_precon(self, atoms, reinitialize=None):
1054 if self.r_NN is None:
1055 self.r_NN = estimate_nearest_neighbour_distance(atoms,
1056 self.neighbor_list)
1058 if self.r_cut is None:
1059 # This is the first time this function has been called, and no
1060 # cutoff radius has been specified, so calculate it automatically.
1061 self.r_cut = 2.0 * self.r_NN
1062 elif self.r_cut < self.r_NN:
1063 warning = ('WARNING: r_cut (%.2f) < r_NN (%.2f), '
1064 'increasing to 1.1*r_NN = %.2f' % (self.r_cut,
1065 self.r_NN,
1066 1.1 * self.r_NN))
1067 warnings.warn(warning)
1068 self.r_cut = 1.1 * self.r_NN
1070 if reinitialize is None:
1071 # The caller has not specified whether or not to recalculate mu,
1072 # so the Precon's setting is used.
1073 reinitialize = self.reinitialize
1075 if self.mu is None:
1076 # Regardless of what the caller has specified, if we don't
1077 # currently have a value of mu, then we need one.
1078 reinitialize = True
1080 if reinitialize:
1081 self.estimate_mu(atoms)
1083 if self.P is not None:
1084 real_atoms = atoms
1085 if isinstance(atoms, Filter):
1086 real_atoms = atoms.atoms
1087 if self.old_positions is None:
1088 self.old_positions = real_atoms.positions
1089 displacement, _ = find_mic(real_atoms.positions -
1090 self.old_positions,
1091 real_atoms.cell, real_atoms.pbc)
1092 self.old_positions = real_atoms.get_positions()
1093 max_abs_displacement = abs(displacement).max()
1094 self.logfile.write('max(abs(displacements)) = %.2f A (%.2f r_NN)' %
1095 (max_abs_displacement,
1096 max_abs_displacement / self.r_NN))
1097 if max_abs_displacement < 0.5 * self.r_NN:
1098 return
1100 # Create the preconditioner:
1101 start_time = time.time()
1102 self._make_sparse_precon(atoms, force_stab=self.force_stab)
1103 self.logfile.write(
1104 f'--- Precon created in {(time.time() - start_time)} seconds ---\n')
1106 def _make_sparse_precon(self, atoms, initial_assembly=False,
1107 force_stab=False):
1108 """Create a sparse preconditioner matrix based on the passed atoms.
1110 Args:
1111 atoms: the Atoms object used to create the preconditioner.
1113 Returns:
1114 A scipy.sparse.csr_matrix object, representing a d*N by d*N matrix
1115 (where N is the number of atoms, and d is the value of self.dim).
1116 BE AWARE that using numpy.dot() with this object will result in
1117 errors/incorrect results - use the .dot method directly on the
1118 sparse matrix instead.
1120 """
1121 self.logfile.write('creating sparse precon: initial_assembly=%r, '
1122 'force_stab=%r, apply_positions=%r, '
1123 'apply_cell=%r\n' %
1124 (initial_assembly, force_stab,
1125 self.apply_positions, self.apply_cell))
1127 N = len(atoms)
1128 start_time = time.time()
1129 if self.apply_positions:
1130 # compute neighbour list
1131 i_list, j_list, rij_list, _fixed_atoms = get_neighbours(
1132 atoms, self.r_cut, self.neighbor_list)
1133 self.logfile.write(
1134 f'--- neighbour list created in {(time.time() - start_time)} s '
1135 '---\n')
1137 row = []
1138 col = []
1139 data = []
1141 # precon is mu_c*identity for cell DoF
1142 start_time = time.time()
1143 if isinstance(atoms, Filter):
1144 i = N - 3
1145 j = N - 2
1146 k = N - 1
1147 x = ijk_to_x(i, j, k)
1148 row.extend(x)
1149 col.extend(x)
1150 if self.apply_cell:
1151 data.extend(np.repeat(self.mu_c, 9))
1152 else:
1153 data.extend(np.repeat(self.mu_c, 9))
1154 self.logfile.write(
1155 f'--- computed triplet format in {(time.time() - start_time)} s '
1156 '---\n')
1158 conn = sparse.lil_matrix((N, N), dtype=bool)
1160 if self.apply_positions and not initial_assembly:
1161 if self.morses is not None:
1162 for morse in self.morses:
1163 self.add_morse(morse, atoms, row, col, data, conn)
1165 if self.bonds is not None:
1166 for bond in self.bonds:
1167 self.add_bond(bond, atoms, row, col, data, conn)
1169 if self.angles is not None:
1170 for angle in self.angles:
1171 self.add_angle(angle, atoms, row, col, data, conn)
1173 if self.dihedrals is not None:
1174 for dihedral in self.dihedrals:
1175 self.add_dihedral(dihedral, atoms, row, col, data, conn)
1177 if self.apply_positions:
1178 for i, j, rij in zip(i_list, j_list, rij_list):
1179 if not conn[i, j]:
1180 coeff = self.get_coeff(rij)
1181 x = ij_to_x(i, j)
1182 row.extend(x)
1183 col.extend(x)
1184 data.extend(3 * [-coeff] + 3 * [coeff])
1186 row.extend(range(self.dim * N))
1187 col.extend(range(self.dim * N))
1188 if initial_assembly:
1189 data.extend([self.mu * self.c_stab] * self.dim * N)
1190 else:
1191 data.extend([self.c_stab] * self.dim * N)
1193 # create the matrix
1194 start_time = time.time()
1195 self.P = sparse.csc_matrix(
1196 (data, (row, col)), shape=(self.dim * N, self.dim * N))
1197 self.logfile.write(
1198 f'--- created CSC matrix in {(time.time() - start_time)} s ---\n')
1200 if not initial_assembly:
1201 self.P = apply_fixed(atoms, self.P)
1203 self.P = self.P.tocsr()
1204 self.create_solver()
1207def make_precon(precon, atoms=None, **kwargs):
1208 """
1209 Construct preconditioner from a string and optionally build for Atoms
1211 Parameters
1212 ----------
1213 precon - one of 'C1', 'Exp', 'Pfrommer', 'FF', 'Exp_FF', 'ID', None
1214 or an instance of a subclass of `ase.optimize.precon.Precon`
1216 atoms - ase.atoms.Atoms instance, optional
1217 If present, build apreconditioner for this Atoms object
1219 **kwargs - additional keyword arguments to pass to Precon constructor
1221 Returns
1222 -------
1223 precon - instance of relevant subclass of `ase.optimize.precon.Precon`
1224 """
1225 lookup = {
1226 'C1': C1,
1227 'Exp': Exp,
1228 'Pfrommer': Pfrommer,
1229 'FF': FF,
1230 'Exp_FF': Exp_FF,
1231 'ID': IdentityPrecon,
1232 'IdentityPrecon': IdentityPrecon,
1233 None: IdentityPrecon
1234 }
1235 if isinstance(precon, str) or precon is None:
1236 cls = lookup[precon]
1237 precon = cls(**kwargs)
1238 if atoms is not None:
1239 precon.make_precon(atoms)
1240 return precon
1243class SplineFit:
1244 """
1245 Fit a cubic spline fit to images
1246 """
1248 def __init__(self, s, x):
1249 self._s = s
1250 self._x_data = x
1251 self._x = CubicSpline(self._s, x, bc_type='not-a-knot')
1252 self._dx_ds = self._x.derivative()
1253 self._d2x_ds2 = self._x.derivative(2)
1255 @property
1256 def s(self):
1257 return self._s
1259 @property
1260 def x_data(self):
1261 return self._x_data
1263 @property
1264 def x(self):
1265 return self._x
1267 @property
1268 def dx_ds(self):
1269 return self._dx_ds
1271 @property
1272 def d2x_ds2(self):
1273 return self._d2x_ds2
1276class PreconImages:
1277 def __init__(self, precon, images, **kwargs):
1278 """
1279 Wrapper for a list of Precon objects and associated images
1281 This is used when preconditioning a NEB object. Equation references
1282 refer to Paper IV in the :class:`ase.mep.NEB` documentation, i.e.
1284 S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys.
1285 150, 094109 (2019)
1286 https://dx.doi.org/10.1063/1.5064465
1288 Args:
1289 precon (str or list): preconditioner(s) to use
1290 images (list of Atoms): Atoms objects that define the state
1292 """
1293 self.images = images
1294 self._spline = None
1296 if isinstance(precon, list):
1297 if len(precon) != len(images):
1298 raise ValueError(f'length mismatch: len(precon)={len(precon)} '
1299 f'!= len(images)={len(images)}')
1300 self.precon = precon
1301 return
1302 P0 = make_precon(precon, images[0], **kwargs)
1303 self.precon = [P0]
1304 for image in images[1:]:
1305 P = P0.copy()
1306 P.make_precon(image)
1307 self.precon.append(P)
1309 def __len__(self):
1310 return len(self.precon)
1312 def __iter__(self):
1313 return iter(self.precon)
1315 def __getitem__(self, index):
1316 return self.precon[index]
1318 def apply(self, all_forces, index=None):
1319 """Apply preconditioners to stored images
1321 Args:
1322 all_forces (array): forces on images, shape (nimages, natoms, 3)
1323 index (slice, optional): Which images to include. Defaults to all.
1325 Returns:
1326 precon_forces: array of preconditioned forces
1327 """
1328 if index is None:
1329 index = slice(None)
1330 precon_forces = []
1331 for precon, image, forces in zip(self.precon[index],
1332 self.images[index],
1333 all_forces):
1334 f_vec = forces.reshape(-1)
1335 pf_vec, _ = precon.apply(f_vec, image)
1336 precon_forces.append(pf_vec.reshape(-1, 3))
1338 return np.array(precon_forces)
1340 def average_norm(self, i, j, dx):
1341 """Average norm between images i and j
1343 Args:
1344 i (int): left image
1345 j (int): right image
1346 dx (array): vector
1348 Returns:
1349 norm: norm of vector wrt average of precons at i and j
1350 """
1351 return np.sqrt(0.5 * (self.precon[i].dot(dx, dx) +
1352 self.precon[j].dot(dx, dx)))
1354 def get_tangent(self, i):
1355 """Normalised tangent vector at image i
1357 Args:
1358 i (int): image of interest
1360 Returns:
1361 tangent: tangent vector, normalised with appropriate precon norm
1362 """
1363 tangent = self.spline.dx_ds(self.spline.s[i])
1364 tangent /= self.precon[i].norm(tangent)
1365 return tangent.reshape(-1, 3)
1367 def get_residual(self, i, imgforce):
1368 # residuals computed according to eq. 11 in the paper
1369 P_dot_imgforce = self.precon[i].Pdot(imgforce.reshape(-1))
1370 return np.linalg.norm(P_dot_imgforce, np.inf)
1372 def get_spring_force(self, i, k1, k2, tangent):
1373 """Spring force on image
1375 Args:
1376 i (int): image of interest
1377 k1 (float): spring constant for left spring
1378 k2 (float): spring constant for right spring
1379 tangent (array): tangent vector, shape (natoms, 3)
1381 Returns:
1382 eta: NEB spring forces, shape (natoms, 3)
1383 """
1384 # Definition following Eq. 9 in Paper IV
1385 nimages = len(self.images)
1386 k = 0.5 * (k1 + k2) / (nimages ** 2)
1387 curvature = self.spline.d2x_ds2(self.spline.s[i]).reshape(-1, 3)
1388 # complete Eq. 9 by including the spring force
1389 eta = k * self.precon[i].vdot(curvature, tangent) * tangent
1390 return eta
1392 def get_coordinates(self, positions=None):
1393 """Compute displacements wrt appropriate precon metric for each image
1395 Args:
1396 positions (list or array, optional) - images positions.
1397 Shape either (nimages * natoms, 3) or ((nimages-2)*natoms, 3)
1399 Returns:
1400 s : array shape (nimages,), reaction coordinates, in range [0, 1]
1401 x : array shape (nimages, 3 * natoms), flat displacement vectors
1402 """
1403 nimages = len(self.precon)
1404 natoms = len(self.images[0])
1405 d_P = np.zeros(nimages)
1406 x = np.zeros((nimages, 3 * natoms)) # flattened positions
1407 if positions is None:
1408 positions = [image.positions for image in self.images]
1409 elif isinstance(positions, np.ndarray) and len(positions.shape) == 2:
1410 positions = positions.reshape(-1, natoms, 3)
1411 positions = [positions[i, :, :] for i in range(len(positions))]
1412 if len(positions) == len(self.images) - 2:
1413 # prepend and append the non-moving images
1414 positions = ([self.images[0].positions] + positions +
1415 [self.images[-1].positions])
1416 assert len(positions) == len(self.images)
1418 x[0, :] = positions[0].reshape(-1)
1419 for i in range(1, nimages):
1420 x[i, :] = positions[i].reshape(-1)
1421 dx, _ = find_mic(positions[i] - positions[i - 1],
1422 self.images[i - 1].cell,
1423 self.images[i - 1].pbc)
1424 dx = dx.reshape(-1)
1425 d_P[i] = self.average_norm(i, i - 1, dx)
1427 s = d_P.cumsum() / d_P.sum() # Eq. A1 in paper IV
1428 return s, x
1430 def spline_fit(self, positions=None):
1431 """Fit 3 * natoms cubic splines as a function of reaction coordinate
1433 Returns:
1434 fit : :class:`ase.optimize.precon.SplineFit` object
1435 """
1436 s, x = self.get_coordinates(positions)
1437 return SplineFit(s, x)
1439 @property
1440 def spline(self):
1441 s, x = self.get_coordinates()
1442 if self._spline and (np.abs(s - self._old_s).max() < 1e-6 and
1443 np.abs(x - self._old_x).max() < 1e-6):
1444 return self._spline
1446 self._spline = self.spline_fit()
1447 self._old_s = s
1448 self._old_x = x
1449 return self._spline