 r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# flake8: noqa

2import logging

3import math

4import numpy as np

6###CO <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

7try:

8 import scipy

9 import scipy.linalg

10 have_scipy = True

11except ImportError:

12 have_scipy = False

13#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

15from ase.utils import longsum

17logger = logging.getLogger(__name__)

19###CO <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

20class LinearPath:

21 """Describes a linear search path of the form t -> t g

22 """

24 def __init__(self, dirn):

25 """Initialise LinearPath object

27 Args:

28 dirn : search direction

29 """

30 self.dirn = dirn

32 def step(self, alpha):

33 return alpha * self.dirn

37def nullspace(A, myeps=1e-10):

38 """The RumPath class needs the ability to compute the null-space of

39 a small matrix. This is provided here. But we now also need scipy!

41 This routine was copy-pasted from

42 http://stackoverflow.com/questions/5889142/python-numpy-scipy-finding-the-null-space-of-a-matrix

43 How the h*** does numpy/scipy not have a null-space implemented?

44 """

45 u, s, vh = scipy.linalg.svd(A)

46 padding = max(0, np.shape(A) - np.shape(s))

47 null_mask = np.concatenate(((s <= myeps),

49 axis=0)

50 null_space = scipy.compress(null_mask, vh, axis=0)

51 return scipy.transpose(null_space)

55class RumPath:

56 """Describes a curved search path, taking into account information

57 about (near-) rigid unit motions (RUMs).

59 One can tag sub-molecules of the system, which are collections of

60 particles that form a (near-)rigid unit. Let x1, ... xn be the positions

61 of one such molecule, then we construct a path of the form

62 xi(t) = xi(0) + (exp(K t) - I) yi + t wi + t c

63 where yi = xi - <x>, c = <g> is a rigid translation, K is anti-symmetric

64 so that exp(tK) yi denotes a rotation about the centre of mass, and wi

65 is the remainind stretch of the molecule.

67 The following variables are stored:

68 * rotation_factors : array of acceleration factors

69 * rigid_units : array of molecule indices

70 * stretch : w

71 * K : list of K matrices

72 * y : list of y-vectors

73 """

75 def __init__(self, x_start, dirn, rigid_units, rotation_factors):

76 """Initialise a RumPath

78 Args:

79 x_start : vector containing the positions in d x nAt shape

80 dirn : search direction, same shape as x_start vector

81 rigid_units : array of arrays of molecule indices

82 rotation_factors : factor by which the rotation of each molecular

83 is accelerated; array of scalars, same length as

84 rigid_units

85 """

87 if not have_scipy:

88 raise RuntimeError("RumPath depends on scipy, which could not be imported")

90 # keep some stuff stored

91 self.rotation_factors = rotation_factors

92 self.rigid_units = rigid_units

93 # create storage for more stuff

94 self.K = []

95 self.y = []

96 # We need to reshape x_start and dirn since we want to apply

97 # rotations to individual position vectors!

98 # we will eventually store the stretch in w, X is just a reference

100 w = dirn.copy().reshape( [3, len(dirn)/3] )

101 X = x_start.reshape( [3, len(dirn)/3] )

103 for I in rigid_units: # I is a list of indices for one molecule

104 # get the positions of the i-th molecule, subtract mean

105 x = X[:, I]

106 y = x - x.mean(0).T # PBC?

107 # same for forces >>> translation component

108 g = w[:, I]

109 f = g - g.mean(0).T

110 # compute the system to solve for K (see accompanying note!)

111 # A = \sum_j Yj Yj'

112 # b = \sum_j Yj' fj

113 A = np.zeros((3,3))

114 b = np.zeros(3)

115 for j in range(len(I)):

116 Yj = np.array( [ [ y[1,j], 0.0, -y[2,j] ],

117 [ -y[0,j], y[2,j], 0.0 ],

118 [ 0.0, -y[1,j], y[0,j] ] ] )

119 A += np.dot(Yj.T, Yj)

120 b += np.dot(Yj.T, f[:, j])

121 # If the directions y[:,j] span all of R^3 (canonically this is true

122 # when there are at least three atoms in the molecule) but if

123 # not, then A is singular so we cannot solve A k = b. In this case

124 # we solve Ak = b in the space orthogonal to the null-space of A.

125 # TODO:

126 # this can get unstable if A is "near-singular"! We may

127 # need to revisit this idea at some point to get something

128 # more robust

129 N = nullspace(A)

130 b -= np.dot(np.dot(N, N.T), b)

131 A += np.dot(N, N.T)

132 k = scipy.linalg.solve(A, b, sym_pos=True)

133 K = np.array( [ [ 0.0, k, -k ],

134 [ -k, 0.0, k ],

135 [ k, -k, 0.0 ] ] )

136 # now remove the rotational component from the search direction

137 # ( we actually keep the translational component as part of w,

138 # but this could be changed as well! )

139 w[:, I] -= np.dot(K, y)

140 # store K and y

141 self.K.append(K)

142 self.y.append(y)

144 # store the stretch (no need to copy here, since w is already a copy)

145 self.stretch = w

148 def step(self, alpha):

149 """perform a step in the line-search, given a step-length alpha

151 Args:

152 alpha : step-length

154 Returns:

155 s : update for positions

156 """

157 # translation and stretch

158 s = alpha * self.stretch

159 # loop through rigid_units

160 for (I, K, y, rf) in zip(self.rigid_units, self.K, self.y,

161 self.rotation_factors):

162 # with matrix exponentials:

163 # s[:, I] += expm(K * alpha * rf) * p.y - p.y

164 # third-order taylor approximation:

165 # I + t K + 1/2 t^2 K^2 + 1/6 t^3 K^3 - I

166 # = t K (I + 1/2 t K (I + 1/3 t K))

167 aK = alpha * rf * K

168 s[:, I] += np.dot(aK, y + 0.5 * np.dot(aK, y + 1/3. * np.dot( aK, y )) )

170 return s.ravel()

171#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

175class LineSearchArmijo:

177 def __init__(self, func, c1=0.1, tol=1e-14):

178 """Initialise the linesearch with set parameters and functions.

180 Args:

181 func: the function we are trying to minimise (energy), which should

182 take an array of positions for its argument

183 c1: parameter for the sufficient decrease condition in (0.0 0.5)

184 tol: tolerance for evaluating equality

186 """

188 self.tol = tol

189 self.func = func

191 if not (0 < c1 < 0.5):

192 logger.error("c1 outside of allowed interval (0, 0.5). Replacing with "

193 "default value.")

194 print("Warning: C1 outside of allowed interval. Replacing with "

195 "default value.")

196 c1 = 0.1

198 self.c1 = c1

201 ###CO : added rigid_units and rotation_factors

202 def run(self, x_start, dirn, a_max=None, a_min=None, a1=None,

203 func_start=None, func_old=None, func_prime_start=None,

204 rigid_units=None, rotation_factors=None, maxstep=None):

206 """Perform a backtracking / quadratic-interpolation linesearch

207 to find an appropriate step length with Armijo condition.

208 NOTE THIS LINESEARCH DOES NOT IMPOSE WOLFE CONDITIONS!

210 The idea is to do backtracking via quadratic interpolation, stabilised

211 by putting a lower bound on the decrease at each linesearch step.

212 To ensure BFGS-behaviour, whenever "reasonable" we take 1.0 as the

213 starting step.

215 Since Armijo does not guarantee convergence of BFGS, the outer

216 BFGS algorithm must restart when the current search direction

217 ceases to be a descent direction.

219 Args:

220 x_start: vector containing the position to begin the linesearch

221 from (ie the current location of the optimisation)

222 dirn: vector pointing in the direction to search in (pk in [NW]).

223 Note that this does not have to be a unit vector, but the

224 function will return a value scaled with respect to dirn.

225 a_max: an upper bound on the maximum step length allowed. Default is 2.0.

226 a_min: a lower bound on the minimum step length allowed. Default is 1e-10.

227 A RuntimeError is raised if this bound is violated

228 during the line search.

229 a1: the initial guess for an acceptable step length. If no value is

230 given, this will be set automatically, using quadratic

231 interpolation using func_old, or "rounded" to 1.0 if the

232 initial guess lies near 1.0. (specifically for LBFGS)

233 func_start: the value of func at the start of the linesearch, ie

234 phi(0). Passing this information avoids potentially expensive

235 re-calculations

236 func_prime_start: the value of func_prime at the start of the

237 linesearch (this will be dotted with dirn to find phi_prime(0))

238 func_old: the value of func_start at the previous step taken in

239 the optimisation (this will be used to calculate the initial

240 guess for the step length if it is not provided)

241 rigid_units, rotationfactors : see documentation of RumPath, if it is

242 unclear what these parameters are, then leave them at None

243 maxstep: maximum allowed displacement in Angstrom. Default is 0.2.

245 Returns:

246 A tuple: (step, func_val, no_update)

248 step: the final chosen step length, representing the number of

249 multiples of the direction vector to move

250 func_val: the value of func after taking this step, ie phi(step)

251 no_update: true if the linesearch has not performed any updates of

252 phi or alpha, due to errors or immediate convergence

254 Raises:

255 ValueError for problems with arguments

256 RuntimeError for problems encountered during iteration

257 """

259 a1 = self.handle_args(x_start, dirn, a_max, a_min, a1, func_start,

260 func_old, func_prime_start, maxstep)

262 # DEBUG

263 logger.debug("a1(auto) = ", a1)

265 if abs(a1 - 1.0) <= 0.5:

266 a1 = 1.0

268 logger.debug("-----------NEW LINESEARCH STARTED---------")

270 a_final = None

271 phi_a_final = None

272 num_iter = 0

274 ###CO <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

275 # create a search-path

276 if rigid_units is None:

277 # standard linear search-path

278 logger.debug("-----using LinearPath-----")

279 path = LinearPath(dirn)

280 else:

281 logger.debug("-----using RumPath------")

282 # if rigid_units != None, but rotation_factors == None, then

283 # raise an error.

284 if rotation_factors == None:

285 raise RuntimeError('RumPath cannot be created since rotation_factors == None')

286 path = RumPath(x_start, dirn, rigid_units, rotation_factors)

287 #>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

290 while(True):

292 logger.debug("-----------NEW ITERATION OF LINESEARCH----------")

293 logger.debug("Number of linesearch iterations: %d", num_iter)

294 logger.debug("a1 = %e", a1)

296 ###CO replaced: func_a1 = self.func(x_start + a1 * self.dirn)

297 func_a1 = self.func(x_start + path.step(a1))

298 phi_a1 = func_a1

299 # compute sufficient decrease (Armijo) condition

300 suff_dec = (phi_a1 <= self.func_start+self.c1*a1*self.phi_prime_start)

302 # DEBUG

303 # print("c1*a1*phi_prime_start = ", self.c1*a1*self.phi_prime_start,

304 # " | phi_a1 - phi_0 = ", phi_a1 - self.func_start)

305 logger.info("a1 = %.3f, suff_dec = %r", a1, suff_dec)

306 if a1 < self.a_min:

307 raise RuntimeError('a1 < a_min, giving up')

308 if self.phi_prime_start > 0.0:

309 raise RuntimeError("self.phi_prime_start > 0.0")

311 # check sufficient decrease (Armijo condition)

312 if suff_dec:

313 a_final = a1

314 phi_a_final = phi_a1

315 logger.debug("Linesearch returned a = %e, phi_a = %e",

316 a_final, phi_a_final)

317 logger.debug("-----------LINESEARCH COMPLETE-----------")

318 return a_final, phi_a_final, num_iter==0

320 # we don't have sufficient decrease, so we need to compute a

321 # new trial step-length

322 at = - ((self.phi_prime_start * a1) /

323 (2*((phi_a1 - self.func_start)/a1 - self.phi_prime_start)))

324 logger.debug("quadratic_min: initial at = %e", at)

326 # because a1 does not satisfy Armijo it follows that at must

327 # lie between 0 and a1. In fact, more strongly,

328 # at \leq (2 (1-c1))^{-1} a1, which is a back-tracking condition

329 # therefore, we should now only check that at has not become too small,

330 # in which case it is likely that nonlinearity has played a big role

331 # here, so we take an ultra-conservative backtracking step

332 a1 = max( at, a1 / 10.0 )

333 if a1 > at:

334 logger.debug("at (%e) < a1/10: revert to backtracking a1/10", at)

336 # (end of while(True) line-search loop)

337 # (end of run())

341 def handle_args(self, x_start, dirn, a_max, a_min, a1, func_start, func_old,

342 func_prime_start, maxstep):

344 """Verify passed parameters and set appropriate attributes accordingly.

346 A suitable value for the initial step-length guess will be either

347 verified or calculated, stored in the attribute self.a_start, and

348 returned.

350 Args:

351 The args should be identical to those of self.run().

353 Returns:

354 The suitable initial step-length guess a_start

356 Raises:

357 ValueError for problems with arguments

359 """

361 self.a_max = a_max

362 self.a_min = a_min

363 self.x_start = x_start

364 self.dirn = dirn

365 self.func_old = func_old

366 self.func_start = func_start

367 self.func_prime_start = func_prime_start

369 if a_max is None:

370 a_max = 2.0

372 if a_max < self.tol:

373 logger.warning("a_max too small relative to tol. Reverting to "

374 "default value a_max = 2.0 (twice the <ideal> step).")

375 a_max = 2.0 # THIS ASSUMES NEWTON/BFGS TYPE BEHAVIOUR!

377 if self.a_min is None:

378 self.a_min = 1e-10

380 if func_start is None:

381 logger.debug("Setting func_start")

382 self.func_start = self.func(x_start)

384 self.phi_prime_start = longsum(self.func_prime_start * self.dirn)

385 if self.phi_prime_start >= 0:

386 logger.error("Passed direction which is not downhill. Aborting...")

387 raise ValueError("Direction is not downhill.")

388 elif math.isinf(self.phi_prime_start):

389 logger.error("Passed func_prime_start and dirn which are too big. "

390 "Aborting...")

391 raise ValueError("func_prime_start and dirn are too big.")

393 if a1 is None:

394 if func_old is not None:

395 # Interpolating a quadratic to func and func_old - see NW

396 # equation 3.60

397 a1 = 2*(self.func_start - self.func_old)/self.phi_prime_start

398 logger.debug("Interpolated quadratic, obtained a1 = %e", a1)

399 if a1 is None or a1 > a_max:

400 logger.debug("a1 greater than a_max. Reverting to default value "

401 "a1 = 1.0")

402 a1 = 1.0

403 if a1 is None or a1 < self.tol:

404 logger.debug("a1 is None or a1 < self.tol. Reverting to default value "

405 "a1 = 1.0")

406 a1 = 1.0

407 if a1 is None or a1 < self.a_min:

408 logger.debug("a1 is None or a1 < a_min. Reverting to default value "

409 "a1 = 1.0")

410 a1 = 1.0

412 if maxstep is None:

413 maxstep = 0.2

414 logger.debug("maxstep = %e", maxstep)

416 r = np.reshape(dirn, (-1, 3))

417 steplengths = ((a1*r)**2).sum(1)**0.5

418 maxsteplength = np.max(steplengths)

419 if maxsteplength >= maxstep:

420 a1 *= maxstep / maxsteplength

421 logger.debug("Rescaled a1 to fulfill maxstep criterion")

423 self.a_start = a1

425 logger.debug("phi_start = %e, phi_prime_start = %e", self.func_start,

426 self.phi_prime_start)

427 logger.debug("func_start = %s, self.func_old = %s", self.func_start,

428 self.func_old)

429 logger.debug("a1 = %e, a_max = %e, a_min = %e", a1, a_max, self.a_min)

431 return a1