gmisclib.opt

1 """This module is a Levenberg-Marquardt optimizer with 2 stem-size control for numeric differentiation. 3 It just needs a function that produces residuals. 4 It is multi-threaded, so calculations of residuals can 5 be farmed out to many processors. 6 7 UPDATED 12/2009 GPK: NOT TESTED! 8 """ 9 10 11 12 import os 13 import sys 14 import math 15 import random 16 import traceback 17 import threading 18 from gmisclib import die 19 from gmisclib import avio 20 from gmisclib import gpkmisc 21 from gmisclib import Num 22 RA = Num.RA 23 LA = Num.LA 24

25 -class OptError(RuntimeError):

26 - def __init__(self, s=None):

27 RuntimeError.__init__(self, s)

28

29 -class NoDownhill(OptError):

30 - def __init__(self, s=None):

31 OptError.__init__(self, s)

32

33 -class NoDerivative(OptError):

34 - def __init__(self, s=None):

35 OptError.__init__(self, s)

36

37 -class BadParamError(ValueError):

38 __doc__ = """Raised when you give and opt instance some invalid control parameter.""" 39

40 - def __init__(self, s=None):

41 ValueError.__init__(self, s)

42

43 -class BadResult(OptError):

44 __doc__ = """The function to be optimized returns some illegal result. 45 For instance, an array of the wrong size or type, 46 or it returns None for the initial evaluation.""" 47

48 - def __init__(self, s=None):

49 OptError.__init__(self, s)

50 51 52

53 -def vec_equal(a, b):

54 """Check that two Numeric vectors are exactly equal.""" 55 return Num.alltrue(Num.equal(a, b))

56 57

58 -def sumsq(a):

59 """This is the overall error measure.""" 60 # print "sumsq: a=", a, a.shape 61 # return Num.sum(a**2) 62 return Num.sum(Num.sort(a**2))

63 64

65 -def maxabs(a):

66 """This is used as a measure of how much the function changed as a 67 result of a parameter change.""" 68 absa = Num.absolute(a) 69 return absa[Num.argmax(absa)]

70 71

72 -def wtf(r_dz, delta, quantum):

73 """This shows how important different measurements are to 74 the final derivitive estimate. The weight must be small 75 when the change of a parameter is smaller than the quantum 76 (i.e. when r_delta < quantum). It must also be small when 77 the change in z is much larger than T (r_dz>>1). 78 It is largest when r_dz is about 1. 79 """ 80 z_term = abs(r_dz)/(1.0 + r_dz**2) 81 if quantum > 0: 82 r_delta = delta/quantum 83 return (r_delta**2/(1.0 + r_delta**2)) * z_term 84 return z_term

85 86

87 -def near_duplicate(guess, tmp, quantum):

88 for (delt, dz) in tmp: 89 if abs(delt-guess) < 0.5*quantum: 90 return True 91 return False

92 93

94 -def symmetry_point(tmp, T, quantum):

95 """Given a list of (delta, z) in a differentiation, 96 pick one more point that makes the differentiation more symmetric. 97 This is done by picking a new delta that brings 98 sum{delta_i * wtf()} closer to zero, 99 where wtf() is the weight given to 100 each point in the differentiation. 101 """ 102 sumwd, sumd, sumdzsize, sumwt = (0.0, 0.0, 0.0, 0.0) 103 for (delt, dz) in tmp: 104 if dz is not None and delt != 0.0: 105 dzsize = maxabs(dz) 106 wt = wtf(dzsize/T, delt, quantum) 107 sumwt += wt 108 sumwd += delt*wt 109 # print " ( delt=", delt, "wt=", wt, "dzsize=", dzsize, "sumwd=", sumwd, "dw=", delt*wt, ")" 110 sumdzsize += dzsize 111 sumd += abs(delt) 112 if sumwt <= 0.0: 113 return None 114 115 asymmetry = sumwd/sumwt 116 dzsize_vs_delt = sumdzsize/sumd 117 # print "symmetry_point: asymmetry=", asymmetry, "dzsize_vs_delt=", dzsize_vs_delt 118 best = None 119 best_asym_after_guess = abs(asymmetry) 120 for q in (-0.2, -0.4, -0.5, -0.6, -0.7, -0.8, -0.9, -1.0, -1.1, -1.2, -1.35, -1.5, -1.75, -2, -2.5, -3): 121 guess = q * asymmetry 122 if near_duplicate(guess, tmp, quantum): 123 continue 124 125 wt_guess = wtf(dzsize_vs_delt*guess/T, guess, quantum) 126 asym_after_guess = abs(sumwd + guess * wt_guess) / (sumwt + wt_guess) 127 # print " guess=", guess, "sumwd=", sumwd, "wt_guess=", wt_guess, "asym_after_guess=", asym_after_guess 128 if asym_after_guess < best_asym_after_guess: 129 best_asym_after_guess = asym_after_guess 130 best = q*asymmetry 131 return best

132 133 134

135 -def scale_est(tmp, T):

136 if len(tmp) < 2: 137 die.warn("len(tmp)<2 in scale_est") 138 return None 139 sumd, sumsts = (0.0, 0.0) 140 for (delt, dz) in tmp: 141 if dz is not None: 142 sumsts += maxabs(dz) 143 sumd += abs(delt) 144 if abs(sumsts) < 1e-10: 145 die.warn("sumsts is small (%g/%g) in scale_est" % (sumsts, sumd)) 146 if sumsts == 0.0: 147 return "HUGE" 148 s_vs_delt = sumsts/sumd 149 return T/s_vs_delt

150 151

152 -def explored_region(tmp):

153 sumad = 0.0 154 for (delt, dz) in tmp: 155 sumad += abs(delt) 156 return sumad/len(tmp)

157 158

159 -def clz_point(qq, T, quantum):

160 """Given a list of (delta, z) tuples in a differentiation attempt, 161 find another delta that (a) fills a gap in the sequence, 162 (b) is preferably has a sign opposite most of the {delta} values, 163 and (c) has a large wtf(). 164 Points outside the sequence are also considered. 165 """ 166 if len(qq) == 1: 167 return None 168 169 FUDGEFACTOR = 2.0 170 tmp = qq[:] 171 tmp.sort() 172 sumd = 0.0 173 for (delt, dz) in tmp: 174 if dz is not None: 175 wt = wtf(maxabs(dz)/T, delt, quantum) 176 sumd += delt*wt # Find out which side to prefer. 177 178 top = None 179 bot = None 180 possibilities = [] 181 for i in range(1, len(tmp)): 182 delti, dzi = tmp[i] 183 deltim, dzim = tmp[i-1] 184 if dzi is None and dzim is None: 185 continue 186 187 if dzi is None: 188 top = True 189 dzc = dzim 190 elif dzim is None: 191 bot = True 192 dzc = dzi 193 else: 194 dzc = 0.5*(dzi + dzim) 195 196 ctr = 0.5*(delti + deltim) 197 if not near_duplicate(ctr, tmp, quantum): 198 gap = (delti - deltim)/(abs(delti) + abs(deltim) + quantum) 199 wtc = wtf(maxabs(dzc)/T, ctr, quantum) 200 gap *= math.sqrt(wtc) 201 if ctr * sumd < 0: 202 gap *= FUDGEFACTOR 203 possibilities.append( (gap, ctr) ) 204 205 if not top and len(tmp)>1: 206 delti, dzi = tmp[-1] 207 assert dzi is not None 208 deltim, dzim = tmp[-2] 209 ctr = delti + 0.5*(delti-deltim) 210 if not near_duplicate(ctr, tmp, quantum): 211 gap = (delti - deltim)/(abs(delti) + abs(deltim) + quantum) 212 if abs(delti) > abs(deltim): 213 dzc = dzi * (ctr/delti) 214 else: 215 dzc = dzim * (ctr/deltim) 216 wtc = wtf(maxabs(dzc)/T, ctr, quantum) 217 gap *= math.sqrt(wtc) 218 if ctr * sumd < 0: 219 gap *= FUDGEFACTOR 220 possibilities.append( (gap, ctr) ) 221 222 if not bot and len(tmp)>1: 223 delti, dzi = tmp[0] 224 assert dzi is not None 225 deltip, dzip = tmp[1] 226 ctr = delti - 0.5*(deltip-delti) 227 if not near_duplicate(ctr, tmp, quantum): 228 gap = (deltip - delti)/(abs(delti) + abs(deltip) + quantum) 229 if abs(delti) > abs(deltip): 230 dzc = dzi * (ctr/delti) 231 else: 232 dzc = dzip * (ctr/deltip) 233 wtc = wtf(maxabs(dzc)/T, ctr, quantum) 234 gap *= math.sqrt(wtc) 235 if ctr * sumd < 0: 236 gap *= FUDGEFACTOR 237 possibilities.append( (gap, ctr) ) 238 239 if len(possibilities)==0: 240 return None 241 242 possibilities.sort() 243 return possibilities[-1][1]

244 245 246 247

248 -def deriv_estimate(tmp, T, quantum):

249 # This is pretty crude, and could be improved. 250 if len(tmp) < 2: 251 raise NoDerivative 252 sumd = Num.zeros(tmp[0][1].shape, Num.Float) 253 sumwt = 0.0 254 for (delt, dz) in tmp: 255 wt = wtf(maxabs(dz)/T, delt, quantum) 256 # print "#WT=", wt, "stepsize=", maxabs(dz), "T=", T, "delt=", delt 257 increment = (dz*(delt*wt)).astype(Num.Float) 258 Num.add(sumd, increment, sumd) 259 sumwt += delt*delt*wt 260 assert sumwt > 0, "Whoops!" 261 #sumd.savespace(1) 262 return sumd/sumwt

263 264

265 -def need_more_diff_pts(tmp, T, quantum):

266 n = 0 267 nT = 0 268 nq = 0 269 nTl = 0 270 for (delt, dz) in tmp: 271 if dz is not None: 272 n += 1 273 if 0.5*T < maxabs(dz) < 2*T: 274 nT += 1 275 if 0.2*T < maxabs(dz) < 5*T: 276 nTl += 1 277 if abs(delt) > 0.9*quantum: 278 nq += 1 279 return n<3 or nT<1 or nq<2 or nTl<2

280 281 282

283 -def diff_one(x, i, sem):

284 """Differentiate z with respect to parameter i. 285 We enter this function with the semaphore 'sem' already acquired. 286 """ 287 die.note("differentiating", i) 288 # print "# DIFF", i 289 z = x.z() 290 zraw = x.p.zraw() 291 assert z is not None 292 # print "at p=", x.p, "z=", z 293 delta = x.scale[i] 294 # print "delta=", delta 295 v = Num.zeros((x.np,), Num.Float) 296 v[i] = 1.0 297 # print "v=", v 298 tmp = [ (0.0, Num.zeros(z.shape, Num.Float)) ] 299 # print "tmp=", tmp 300 quantum = x.deriv_quantum[i] 301 302 non_None = 0 303 for sign in (-1, 1, -0.5, 0.5, -0.25, 0.25, -0.125, 0.125, -0.0625, 0.0625): 304 delta = sign*max(x.scale[i], quantum) 305 if near_duplicate(delta, tmp, quantum): 306 continue 307 # print "# scale=", x.scale[i], "quantum=", quantum, "delta=", delta, "sign=", sign 308 q = x.constrain(x.p.p, x.p.p + delta*v, x.args) 309 if q is None or vec_equal(q, x.p.p): 310 die.warn("Constraint complains: delta=%g" % delta) 311 continue 312 z1 = prms(q, x.fn, x.args, (sem.id, x.steps), diffparam=i, diffctr=zraw).z() 313 if z1 is not None: 314 if z1.shape != z.shape: 315 sem.release() 316 raise BadResult, 'Size mismatch: %s vs. %s' % (str(z1.shape), 317 str(z.shape)) 318 try: 319 delta_z = z1 - z 320 non_None += 1 321 except OverflowError: 322 die.warn("Overflow in differentiation attempt") 323 delta_z = None 324 325 tmp.append((delta, delta_z)) 326 327 if non_None > 0: 328 break 329 330 if len(tmp) < 2: 331 sem.release() 332 raise NoDerivative, 'len(tmp)<2, initial' 333 334 parity = True 335 while need_more_diff_pts(tmp, x.Td(), quantum): 336 if parity: 337 delta = clz_point(tmp, x.Td(), quantum) 338 if delta is None: 339 delta = symmetry_point(tmp, x.Td(), quantum) 340 else: 341 parity = not parity 342 else: 343 delta = symmetry_point(tmp, x.Td(), quantum) 344 if delta is None: 345 delta = clz_point(tmp, x.Td(), quantum) 346 else: 347 parity = not parity 348 if delta is None: 349 die.warn("Neither symmetry_point nor clz_point") 350 break 351 q = x.constrain(x.p.p, x.p.p + delta*v, x.args) 352 if q is not None and not vec_equal(q, x.p.p): 353 z1 = prms(q, x.fn, x.args, (sem.id, x.steps), diffparam=i, diffctr=zraw).z() 354 tmp.append((delta, z1-z)) 355 else: 356 die.warn("Constrain complains 2") 357 358 # We calculate a derivative from the tmp array. 359 x.diff[i,:] = deriv_estimate(tmp, x.Td(), quantum) 360 361 # Re-estimate scale: 362 newscale = scale_est(tmp, x.Td()) 363 if newscale == "HUGE": 364 die.info("scale_est is huge on %d %g %s" % (i, x.scale[i], str(tmp))) 365 x.newscale[i] = 3.0*x.scale[i] 366 elif newscale is not None: 367 x.newscale[i] = gpkmisc.limit(x.scale[i]*0.1, gpkmisc.geo_mean(x.scale[i], newscale), 368 x.scale[i]*10.0) 369 else: 370 die.info("scale_est is None on %d %g %s" % (i, x.scale[i], str(tmp))) 371 x.newscale[i] = x.scale[i] 372 x.newscale[i] = max(x.newscale[i], quantum) 373 374 # Used for termination conditions. 375 x.explored[i] = explored_region(tmp) 376 # print "# newscale=", x.newscale[i], "quantum=", quantum, "explored=", x.explored[i] 377 sys.stderr.flush() 378 sys.stdout.flush() 379 sem.release()

380 381 382

383 -def eval_lambda(processor, a, b, startp, lamb, out, opt):

384 rmat = RA.standard_normal(a.shape) 385 # The factor of self.np sets the geometric mean of the eigenvalues to 1, 386 # so that a typical eigenvalue is close to the eigenvalues of 387 # the identity matrix. 388 rsq = Num.matrixmultiply(rmat, Num.transpose(rmat)) / a.shape[0] 389 aa = a + rsq * lamb 390 391 try: 392 # s = LA.solve_linear_equations(aa, b) * opt.scale 393 s = LA.solve(aa, b) * opt.scale 394 except LA.LinAlgError: # Singular matrix 395 processor.release() 396 return 397 # print "s=", s 398 # print "startp.p=", startp.p 399 newp = opt.constrain(startp.p, startp.p - s, opt.args) # Trim the step. 400 if newp is None: 401 die.warn("eval_lambda: Step outside of valid region" ) 402 pnew = None 403 else: 404 pnew = prms(newp, opt.fn, opt.args, (processor.id, opt.steps)) 405 pnew.sumsq() # Force evaluation to fill cache. 406 out.append((lamb, pnew, rsq)) 407 processor.release()

408 409 410 411

412 -def lamb_correct(newscale, scale):

413 assert len(newscale.shape) == 1 414 assert newscale.shape == scale.shape 415 return math.exp(2*Num.sum(Num.log(newscale/scale))/newscale.shape[0])

416 417 # threading._VERBOSE = 1 418

419 -class mysem:

420 __doc__ = """This is a semaphore that is automatically released when 421 it is de-allocated. It also contains a processor ID.""" 422

423 - def __init__(self, id, sem):

424 self.id = id 425 self.sem = sem

426

427 - def release(self):

428 self.sem._release(self.id) 429 self.sem = None

430

431 - def __del__(self):

432 if self.sem: 433 self.release()

434 435

436 -class semclass:

437 __doc__ = """This is a semaphore to control the number of 438 simultaneous computations. It hands out a 439 mysem class which is used to release the semaphore.""" 440

441 - def __init__(self, n):

442 self.s = threading.Semaphore(n) 443 self.p_list = range(n) 444 self.nthreads = n

445 # print "SEMCLASS:", self.p_list 446

447 - def acquire(self):

448 # print "SEMCLASS: acquire_wait" 449 self.s.acquire() 450 tmp = self.p_list.pop(-1) 451 # print "SEMCLASS: acquired", tmp, "remaining:", self.p_list 452 return mysem(tmp, self)

453

454 - def _release(self, id):

455 self.p_list.append(id) 456 self.s.release()

457 458 459 460 _NumType = type(Num.array([1], Num.Float)) 461 _tupletype = type(()) 462 _sequenceType = type([]) 463 _floatType = type(1.0) 464

465 -def _zarray(z):

466 if type(z) == _NumType: 467 return Num.ravel(z) 468 assert type(z) == _floatType, "Opt needs an array of mixed floats and Num arrays" 469 return Num.array([z], Num.Float)

470 471

472 -class prms:

473 __doc__ = """This class records parameters and caches function evaluations. 474 All function evaluations happen here.""" 475

476 - def __init__(self, p, fn, args, processor, diffparam=None, diffctr=None):

477 assert len(p.shape) == 1 478 self.p = p 479 self.lock = threading.RLock() 480 self.fn = fn 481 self.args = args 482 self.processor = processor 483 self.diffparam = diffparam 484 self.diffctr = diffctr

485 486

487 - def zraw(self):

488 self.lock.acquire() 489 if not hasattr(self, 'zr_cache'): 490 try: 491 self.zr_cache = (self.fn)(self.p, self.args, self.diffparam, 492 self.diffctr, self.processor) 493 except: 494 traceback.print_exc(file=sys.stdout) 495 die.die("Uncaught exception in self.fn") 496 self.lock.release() 497 return self.zr_cache

498 499

500 - def z(self):

501 self.lock.acquire() 502 if not hasattr(self, 'z_cache'): 503 z = self.zraw() 504 ztype = type(z) 505 if z is None: 506 self.z_cache = None 507 elif ztype == _NumType: 508 self.z_cache = Num.ravel(z) 509 elif ztype == _sequenceType or ztype == _tupletype: # A list of floats mixed with Num arrays 510 self.z_cache = Num.concatenate(map(_zarray, z)) 511 else: 512 raise BadResult, "Unknown type(%s) returned from function" % repr(type(z)) 513 self.lock.release() 514 return self.z_cache

515 516

517 - def sumsq(self):

518 self.lock.acquire() 519 if not hasattr(self, 'sumsq_cache'): 520 z = self.z() 521 if z is None: 522 self.sumsq_cache = None 523 else: 524 self.sumsq_cache = sumsq(z) 525 self.lock.release() 526 return self.sumsq_cache

527

528 - def __len__(self):

529 return self.p.shape[0]

530

531 - def __str__(self):

532 return str(self.p)

533

534 - def __repr__(self):

535 return repr(self.p)

536 537

538 -class LockedList:

539 - def __init__(self, start=None):

540 if start is None: 541 self.x = [] 542 else: 543 self.x = start 544 self.l = threading.Lock()

545

546 - def append(self, x):

547 self.l.acquire() 548 self.x.append(x) 549 self.l.release()

550

551 - def __len__(self):

552 return len(self.x)

553

554 - def pop(self, i = -1):

555 self.l.acquire() 556 try: 557 tmp = self.x.pop(i) 558 except: 559 self.l.release() 560 raise 561 self.l.release() 562 return tmp

563

564 - def items(self):

565 self.l.acquire() 566 tmp = self.x 567 self.l.release() 568 return tmp

569 570

571 -def _dE(T):

572 return random.expovariate(1.0) * T

573 574

575 -def anneal_guts(p, rv, T, lamb, sem, opt):

576 ossq = p.sumsq() 577 if ossq is None: 578 die.warn("anneal_guts: ossq=None") 579 sem.release() 580 return 581 582 np = len(p) 583 # eval, evec = lapack_dsyevd.dsyevd(opt.curv) 584 # eval, evec = LA.Heigenvectors(opt.curv) 585 eval, evec = LA.eigh(opt.curv) 586 r = RA.standard_normal((np,)) * math.sqrt(2*T/np) 587 e = 1.0/Num.sqrt(eval + lamb) 588 s = Num.matrixmultiply(r*e, evec) * opt.curv_scale * opt.active 589 q = opt.constrain(p.p, p.p + s, opt.args) 590 if q is None: 591 die.warn("anneal_guts: started outside constraint region") 592 sem.release() 593 return 594 p = prms(q, opt.fn, opt.args, (sem.id, opt.steps)) 595 ssq = p.sumsq() 596 if ssq is None: 597 die.warn("anneal_guts: ssq=None") 598 sem.release() 599 return 600 prediction = sumsq(r) 601 if opt.verbose: 602 print "AGUTS: ssq-ossq=", ssq-ossq, "quad_predict=", prediction, "lambda=", lamb 603 if ssq < ossq + 2 * _dE(T): 604 rv.append( (p, lamb) ) # Return the result here. 605 sys.stdout.flush() 606 sem.release()

607 608

609 -def _evec_write(log_fd, a_vec, a_val, note):

610 Nlim = 500 611 o = ['note=%s' % note, 'eigenvalue=%g' % a_val] 612 # XXX Below error checks shouldn't be needed... 613 try: 614 sq = Num.absolute(a_vec)**2 615 except OverflowError: 616 die.info("note=" + note) 617 die.info("a=" + str(a_vec)) 618 raise 619 idx = Num.argsort(sq) 620 total = Num.sum(sq) 621 ssq = 0.0 622 n = 0 623 for i in range(idx.shape[0]-1, -1, -1): 624 j = idx[i] 625 o.append( 'c%d=%.3f' % (j, a_vec[j])) 626 ssq += sq[j] 627 n += 1 628 if ssq > total*(1-n/Nlim): 629 break 630 log_fd.write('#Evec: ' + '; '.join(o) + '\n')

631 632 633

634 -def _opt_eval(maybe_fcn, arg):

635 if callable(maybe_fcn): 636 return maybe_fcn(arg) 637 return float(maybe_fcn)

638 639 640

641 -class opt:

642 __doc__ = """A class that implements a optimizer. 643 @ivar scale: estimate of the step size to use in differentiation. 644 @ivar deriv_quantum: minimum step size used in differentiation. 645 @ivar T: simulated annealing temperature. 646 @ivar active: vector of which parameters to optimize. Only the active 647 parameters change their value. 648 @ivar sem: a semaphore to control the number of threads to use 649 (if you want to set the # of threads do self.sem = semclass(nnn)). 650 @ivar constrain: a function that constrains the step. 651 You can use it to do a constrained fit by having it project 652 the proposed step back into the legal volume. 653 Called as self.constrain(old_prms, proposed_prms, self.args). 654 It returns the constrained step, or None if there isn't any. 655 """ 656

657 - def __init__(self, p, T, fcn, args=None):

658 """Create an optimizer object. 659 You can also set the following variables between calls to self.step(): 660 C{scale}, C{deriv_quantum}, C{T}, C{active}, C{sem}, C{constrain} 661 662 @param p: initial parameter vector. 663 @param T: temperature for simulated annealing and numeric differentiation. 664 @param fcn: fcn(p, args, diffparam, diffctr, processor), where 665 C{p} = parameters, 666 C{args} = arbitrary reference passed to function, 667 C{diffparam} = L{None} (unless differentiating), 668 C{diffctr} = L{None} (unless differentiating), 669 C{processor} = (proc_id, i_steps), 670 where proc_id is in C{range(self.nthreads)} 671 and specifies which processor ought 672 to be assigned to the task, and 673 C{i_steps} is the number of steps the 674 optimization has gone through. 675 The function is assumed to return either a Num Python array 676 or a sequence of arrays and floats, mixed. 677 """ 678 self.fn = fcn 679 self.args = args 680 self.steps = 0 681 if isinstance(p, prms): 682 self.p = p 683 else: 684 self.p = prms(Num.array(p, Num.Float, copy=True), self.fn, self.args, (0, self.steps)) 685 self.last_p = self.p # Maintained by terminate(). 686 self.np = len(self.p.p) # Number of parameters 687 self.scale = Num.ones((self.np,), Num.Float) 688 self.deriv_quantum = None 689 self.T = T 690 self.active = Num.ones((self.np,), Num.Int) 691 try: 692 nthreads = os.sysconf('SC_NPROCESSORS_ONLN') 693 except (ValueError, AttributeError): 694 nthreads = 1 695 self.lamb = 1.0 696 self.sem = semclass(nthreads) 697 self.verbose = 1 698 self.constrain = lambda x, y, args: y

699 700 701 702

703 - def set_nthreads(self, n):

704 """Set the number of simultaneous threads to be allowed. 705 Calculations of the residuals are farmed out to threads. 706 """ 707 self.sem = semclass(min(self.np, max(0, n)))

708 709 710 711

712 - def z(self):

713 """What is the current residual?""" 714 return self.p.z()

715 716

717 - def sumsq(self):

718 """What is the current error?""" 719 return self.p.sumsq()

720 721

722 - def _differentiate(self, prms_to_measure = None):

723 """Sets self.diff to the Jacobian matrix. Jacobian is in real units, with 724 no parameter scaling. Also sets self.scale. 725 726 The differentials with respect to unmeasured parameters are left unchanged 727 from whatever they were. This is good under some circumstances, where 728 you might be optimizing just a single (of many) local parameter which only 729 affects the solution locally, and you have some global parameters which 730 you don't expect to change much as the local parameter changes. 731 However, sometimes you might want to zero out self.diff to get rid of 732 old information. 733 """ 734 735 sys.stderr.flush() 736 sys.stdout.flush() 737 if self.deriv_quantum is None: 738 self.deriv_quantum = Num.zeros((self.np,), Num.Float) 739 if prms_to_measure is None: 740 prms_to_measure = self.active 741 zctr = self.z() 742 if zctr is None: 743 raise BadResult, "Function returns None at starting point." 744 no = zctr.shape[0] 745 self.explored = Num.ones((self.np,), Num.Float) 746 self.newscale = Num.ones((self.np,), Num.Float) 747 if not hasattr(self, 'diff') or self.diff is None or self.diff.shape != (self.np, no): 748 self.diff = None # Free the memory. 749 self.diff = Num.zeros((self.np, no), Num.Float) 750 # self.diff.savespace(1) 751 threads = [] 752 # All the work is done in the threads. 753 # This obscure little bit below starts up at most self.nthreads 754 # at a time to do self.np differentiations, starting new ones as old 755 # ones terminate. 756 for i in range(self.np): 757 # We start threads to calculate derivatives, 758 # making sure not to have 759 # more then self.nthreads operating at one time. 760 if prms_to_measure[i]: 761 processor = self.sem.acquire() 762 thr = threading.Thread(target=diff_one, 763 args=(self, i, processor)) 764 threads.append(thr) 765 thr.start() 766 767 for thr in threads: 768 thr.join()

769 770

771 - def _diff_update(self, oldp, newp):

772 """Adjust the derivitive matrix to match a newly recorded newp.z() 773 value. This is a secant update.""" 774 if self.verbose: 775 print "# diff_update" 776 step = newp.p - oldp.p 777 step2 = sumsq(step/self.curv_scale) 778 if step2 < 1e-10: 779 return 0 # Parameters changed by trivial amount. 780 # We expect the residuals to change by step*diff 781 782 chg = newp.z() - oldp.z() 783 chg_sz = sumsq(chg) 784 785 ex = Num.matrixmultiply(step, self.diff) 786 ex_sz = sumsq(ex) 787 cor_sz = sumsq(chg - ex) 788 # If the residuals don't change much, then reduce the update, 789 # because the secant isn't long enough to get a good measurement. 790 f = chg_sz/(chg_sz + 0.5*self.Td()) 791 # We are cautious, and make only a partial update when the secant 792 # disagrees badly with our initial estimate. 793 f *= (ex_sz / (ex_sz + cor_sz)) * ( chg_sz / (chg_sz + cor_sz)) 794 if self.verbose: 795 print " f=", f 796 797 if self.verbose: 798 print " change in z=", chg_sz, "expected=", ex_sz, "correction=", cor_sz 799 # Correction is a BIG matrix... 800 correction = Num.outerproduct(step, chg-ex) * (f/step2) 801 Num.add(self.diff, correction, self.diff) 802 del correction # Free some memory... 803 self.calc_curv() 804 return 1

805 806 807

808 - def _curv_update(self, oldp, newp):

809 """Adjust the curvature matrix to match a newly recorded newp.sumsq() 810 value. This is a secant update.""" 811 # expect sumsq to change by 2*resid*diff*s + s*diff*diff*s 812 if self.verbose: 813 print "# curv_update" 814 step = (newp.p - oldp.p)/self.curv_scale 815 step2 = sumsq(step) 816 EPS = 1e-10 817 if step2 < EPS: 818 return 0.0 819 820 step_chg = sumsq(newp.z() - oldp.z()) 821 # If the residuals don't change much, then reduce the update, 822 # because the secant isn't long enough to get a good measurement. 823 f = step_chg/(step_chg + 10*self.Td()) 824 825 b = Num.matrixmultiply( self.diff, oldp.z()) 826 parab = Num.dot(step, Num.matrixmultiply(self.curv, step)) 827 if self.verbose: 828 print "# parabolic part:", parab 829 print "# downhill part:", 2*Num.matrixmultiply(step*self.curv_scale, b) 830 ex = 2*Num.matrixmultiply(step*self.curv_scale, b) + parab 831 if self.verbose: 832 print "# expect a change of", ex, "get", newp.sumsq() - oldp.sumsq() 833 delta = (newp.sumsq() - oldp.sumsq()) - ex 834 if delta <= 0: 835 return 0.0 # Don't let eigenvalues get negative. 836 print "# delta=", delta 837 correction = Num.outerproduct(step, step) * (f*delta/step2) 838 Num.add(self.curv, correction, self.curv) 839 return delta

840 841 842 LOG_NAME = 'a.tmp' 843

844 - def _log(self):

845 if not hasattr(self, 'log_fd'): 846 self.log_fd = open(self.LOG_NAME, "w") 847 return self.log_fd

848 849

850 - def print_ev_diag(self, misc={}, vectors=0):

851 log_fd = self._log() 852 a = self.curv 853 854 if vectors: 855 # eval, evec = lapack_dsyevd.dsyevd(a) 856 # eval, evec = LA.Heigenvectors(a) 857 eval, evec = LA.eigh(a) 858 else: 859 # eval, evec = lapack_dsyevd.dsy(a, jobz='N') 860 # eval = LA.Heigenvalues(a) 861 eval = LA.eigh(a) 862 evec = None 863 864 mx_e = gpkmisc.N_maximum(eval) 865 o = misc.copy() 866 o['Max_EV'] = mx_e 867 o['Min_EV'] = gpkmisc.N_minimum(eval) 868 lg_thr = mx_e/10.0 869 sm_thr = mx_e/1e6 870 o['N_large_ev'] = Num.sum(Num.greater(eval, lg_thr)) 871 o['large_thr'] = lg_thr 872 o['N_small_ev'] = Num.sum(Num.less(eval, sm_thr)) 873 o['small_thr'] = sm_thr 874 o['Median_EV'] = gpkmisc.N_median(eval) 875 o['N_sm_lambda'] = Num.sum(Num.less(eval, self.lamb)) 876 o['lambda'] = self.lamb 877 log_fd.write('#EV: ' + avio.concoct(o) + '\n') 878 if vectors: 879 n = eval.shape[0] 880 for i in range(n): 881 if eval[i] > lg_thr: 882 _evec_write(log_fd, evec[i,:], eval[i], "L") 883 if eval[i] < sm_thr: 884 _evec_write(log_fd, evec[i,:], eval[i], "S") 885 log_fd.flush()

886 887 888

889 - def Td(self):

890 """This allows the differentiation temperature to be dynamically set.""" 891 if callable(self.T): 892 return self.T(self) 893 return float(self.T)

894 895

896 - def dE(self):

897 return _dE(self.Td())

898

899 - def quick_lambda(self, a, b, startp, lamb, ntries=10):

900 if self.verbose: 901 print "# quick lambda: lamb=", lamb 902 attempt = 0 903 bestp, bestl, bestrsq = (None, 0.0, None) 904 NTRIES = 3 + int(round(math.sqrt(self.np))) * self.sem.nthreads 905 EPS = 1e-13 906 ev_sum = Num.trace(a) 907 lamb = min(lamb/5.0, ev_sum*5.0) * math.exp(-self.sem.nthreads * 0.2) 908 rv = LockedList() 909 while attempt < NTRIES: 910 processor = self.sem.acquire() 911 tmp = threading.Thread(target=eval_lambda, 912 args=(processor, a, b, startp, lamb, rv, self)) 913 tmp.start() 914 915 processor = self.sem.acquire() # A thread may have finished. 916 while len(rv) > 0: 917 bestl1, bestp1, bestrsq1 = rv.pop(0) 918 if bestp1 is None: 919 continue 920 if bestp1.z() is None: 921 continue 922 if bestp is None or bestp1.sumsq()<bestp.sumsq(): 923 bestp, bestl, bestrsq = (bestp1, bestl1, bestrsq1) 924 if bestp1.sumsq() < startp.sumsq() + self.dE(): 925 # Good enough. Don't start any more. 926 attempt = NTRIES 927 elif sumsq(bestp1.z() - startp.z()) < 0.2 * self.Td(): 928 # We're too close to the starting point. Go no closer. 929 attempt = NTRIES 930 processor.release() 931 attempt += 1 932 lamb = lamb*2.0*math.pow(10.0, 1.0/self.sem.nthreads) + EPS*ev_sum 933 # THIS FUNCTION CAN EXIT WITH CHILD THREADS STILL RUNNING. 934 # That's OK; they will eventually terminate and clean themselves up. 935 return (bestl, bestp, bestrsq)

936 937

938 - def search_with_update(self, lamb, newp):

939 if self.verbose: 940 print "search with update: lamb=", lamb, "ssq=", newp.sumsq() 941 # self._curv_update(self.p, newp) 942 self._diff_update(self.p, newp) 943 # self._curv_update(self.p, newp) 944 prm_per_proc = float(self.np)/float(self.sem.nthreads) 945 tthresh = 2*int(math.ceil(math.sqrt(prm_per_proc/7.0))) 946 impthresh = newp.sumsq()/prm_per_proc 947 tries = 0 948 improvement = 0.0 949 bestp = newp 950 while tries < tthresh: 951 if self.verbose: 952 print " swu: tries=", tries 953 # ds = self.scale[:,Num.NewAxis] * self.diff 954 # b = Num.matrixmultiply(ds, bestp.z()) 955 b = self.scale * Num.matrixmultiply(self.diff, bestp.z()) 956 # Make sure this is the same as above commented out... 957 958 bestl, bestp1, bestrsq1 = self.quick_lambda(self.curv, b, bestp, lamb, 3) 959 lamb = max(1e-10, bestl) 960 if bestp1 is None: 961 if self.verbose: 962 print " None" 963 break 964 if self.verbose: 965 print " bestp1=", bestp1.sumsq(), "bestp=", bestp.sumsq(), "T=", self.Td() 966 if bestp1.sumsq() <= bestp.sumsq(): 967 self._diff_update(self.p, newp) 968 # self._curv_update(bestp, bestp1) 969 improvement = bestp.sumsq() - bestp1.sumsq() 970 if self.verbose: 971 print " improvement=", improvement, "impthresh=", impthresh 972 bestp = bestp1 973 else: 974 improvement = 0.0 975 976 if improvement < self.Td(): 977 if self.verbose: 978 print " Not much better" 979 break 980 elif improvement <= impthresh: 981 tries += 1 982 return (bestp, bestrsq1)

983 984

985 - def calc_curv(self):

986 if self.verbose: 987 print "scale=", self.scale 988 # DS is a BIG matrix... 989 # ds = self.scale[:,Num.NewAxis] * self.diff 990 # a = Num.matrixmultiply(ds, Num.transpose(ds)) 991 a = self.scale[:,Num.NewAxis] \ 992 * Num.matrixmultiply(self.diff, Num.transpose(self.diff)) \ 993 * self.scale[Num.NewAxis, :] 994 assert a.shape == (self.np, self.np) 995 # Next is kluge to keep inactive parameters from generating numeric problems: 996 ev_avg = Num.trace(a)/Num.sum(self.active) 997 for i in Num.nonzero(1 - self.active): 998 a[i,i] += ev_avg 999 1000 # print "a=", a 1001 self.curv = a # Use curv_scale with a. 1002 self.curv_scale = self.scale

1003

1004 - def measure(self):

1005 """If you decrease the number of active parameters, it might behoove 1006 you to set self.diff=None, before you call this function, 1007 so that old derivitives get flushed. 1008 """ 1009 self._differentiate() 1010 self.calc_curv() 1011 self.last = 'measure'

1012 1013

1014 - def step(self):

1015 """A minimization step.""" 1016 if not self.Td() > 0: 1017 raise BadParamError, 'Temperature must be positive.' 1018 if self.verbose: 1019 print "BEGIN STEP" 1020 # ds = self.measure() 1021 # b = Num.matrixmultiply(ds, self.z()) 1022 self.measure() 1023 b = self.scale * Num.matrixmultiply(self.diff, self.z()) 1024 # print "b=", b 1025 1026 self.lamb, bestp, bestrsq = self.quick_lambda(self.curv, b, self.p, self.lamb) 1027 1028 if bestp is None or bestp.sumsq() is None: 1029 die.warn("No Downhill, quick_lambda. Will anneal.") 1030 # We use Td instead of Ta below, because this is essentially 1031 # an attempt to recover from an error, not a simulated 1032 # anealing step to get error bars. We want to take a step 1033 # comparable to the size of the region we've explored during 1034 # the differentiation, and that region is set by Td. 1035 self.one_anneal(self.Td()) 1036 return 1037 if not ( bestp.sumsq() < self.sumsq() ): 1038 die.warn("No downhill after quick_lambda: sumsq %g -> %g" 1039 % (self.sumsq(), bestp.sumsq())) 1040 1041 # Because the time to calculate the derivatives scales 1042 # as the number of parameters per processor increases, we 1043 # can do increasingly finer searches in lambda, 1044 # without dominating the total computational time. 1045 # bestl1, bestp1 = self.search_lambda(a, b, bestp, self.lamb) 1046 # if bestp is None or ( bestp1 is not None 1047 # and bestp1.sumsq()<bestp.sumsq() ): 1048 # bestp = bestp1 1049 # self.lamb = bestl1 1050 bestp1, bestrsq1 = self.search_with_update(self.lamb, bestp) 1051 if bestp is None or ( bestp1 is not None 1052 and bestp1.sumsq()<bestp.sumsq() ): 1053 bestp = bestp1 1054 bestrsq = bestrsq1 1055 1056 1057 sys.stderr.flush() 1058 sys.stdout.flush() 1059 if bestp.sumsq() < self.sumsq() + self.dE(): 1060 self.p = bestp 1061 self.lambdir = bestrsq 1062 self.lamb = max(self.lamb * lamb_correct(self.newscale, self.scale), 1e-10) 1063 self.scale = self.newscale 1064 self.steps += 1 1065 self.last = 'step' 1066 else: 1067 die.warn("No Downhill, search. Will anneal.") 1068 self.one_anneal(self.Td())

1069 1070

1071 - def predicted_sumsq_delta(self, last_p):

1072 """The energy change of a step, assuming you're at the minimum.""" 1073 p_step = self.p.p - last_p.p 1074 ps = p_step / self.curv_scale 1075 return Num.dot(ps, Num.matrixmultiply(self.curv, ps))

1076 1077 1078 1079

1080 - def one_anneal(self, T):

1081 """Take one simulated annealing step, as fast as possible. 1082 We send all the processors in random directions, and the one 1083 that reports an acceptable step first is returned. 1084 Returns 1 on success, 0 if no good step could be found 1085 in a reasonable number of tries. 1086 The temperature can also be a function of one argument 1087 (as an alternative to a float), in which case, the temperature 1088 is calculated as T(self).""" 1089 1090 # Allow T to be calculated dynamically. 1091 Ta = _opt_eval(T, self) 1092 1093 attempt = 0 1094 NTRIES = 3 + int(round(math.sqrt(self.np))) * self.sem.nthreads 1095 EPS = 1e-13 1096 lamb = self.lamb / 30.0 + EPS*Num.trace(self.curv) 1097 rv = LockedList() 1098 self.last = 'anneal' 1099 while attempt < NTRIES: 1100 processor = self.sem.acquire() 1101 tmp = threading.Thread(target=anneal_guts, 1102 args=(self.p, rv, Ta, lamb, processor, self)) 1103 tmp.start() 1104 1105 processor = self.sem.acquire() # A thread may have finished. 1106 if len(rv) > 0: 1107 processor.release() 1108 self.p, rlamb = rv.pop(0) 1109 self.lamb = math.sqrt(rlamb * self.lamb) 1110 # THIS FUNCTION CAN EXIT WITH CHILD THREADS STILL RUNNING. 1111 return 1 1112 processor.release() 1113 attempt += 1 1114 lamb *= math.pow(3.0, 1.0/self.sem.nthreads) 1115 die.warn("simulated annealing fails after %d attempts" % attempt) 1116 return 0

1117 1118

1119 - def terminate(self):

1120 """Returns number between 0 and 1 if the optimization seems finished. 1121 Returns -1 if it's clearly not done yet. 1122 This is designed to be called after step(), 1123 so that both self.p and self.last_p refer to the results 1124 of downhill steps. 1125 """ 1126 # The step is small compared to the region 1127 # that was explored for differentiation: 1128 # Unfortunately, this test (sort-of-correctly) indicates 1129 # that optimizations are never done when the chi-squared surface 1130 # is nearly degenerate. Thus, we weaken this first test, 1131 # so termination will eventually happen, but it will be slow 1132 # if the first condition is not really met. 1133 p_step = self.p.p - self.last_p.p 1134 # print "# pstep=", p_step 1135 # print "# explored=", self.explored 1136 step_to_explored = Num.sum(Num.absolute(p_step/self.explored))/p_step.shape[0] 1137 small1 = 1.0/(1 + math.sqrt(step_to_explored)) 1138 if self.verbose: 1139 print "TERM: step_to_explored=", step_to_explored 1140 1141 # The real energy change (downward step) is small compared to the temperature: 1142 delta_E = self.last_p.sumsq() - self.sumsq() 1143 small2 = delta_E < self.Td()/math.sqrt((self.steps+1)*self.np) 1144 if self.verbose: 1145 print "TERM: last_sumsq=", self.last_p.sumsq(), \ 1146 "sumsq=", self.sumsq(), \ 1147 "diff=", self.last_p.sumsq() - self.sumsq() 1148 1149 # The predicted energy change is comparable to the temperature: 1150 predicted = self.predicted_sumsq_delta(self.last_p) 1151 small3 = predicted < 10 * self.Td() 1152 if self.verbose: 1153 print "TERM: predicted improvement=", predicted, "T=", self.Td() 1154 print "TERM: tests:", small1, small2, small3 1155 1156 self.last_p = self.p 1157 if small2 and small3: 1158 return small1 1159 return -1.0

1160 1161

1162 - def covariance_under(self):

1163 """This is an underestimate of covariance. 1164 Adding in the self.lamb term means that it attempts to 1165 correct for the nonlinearity of the problem, 1166 at least along the direction of search. 1167 """ 1168 # c = LA.inverse(self.curv + self.lamb*self.lambdir) 1169 c = LA.inv(self.curv + self.lamb*self.lambdir) 1170 return self.curv_scale[:,Num.NewAxis] * c * self.curv_scale[Num.NewAxis,:]

1171 1172

1173 - def covariance(self):

1174 """Covariance estimate.""" 1175 try: 1176 # c = LA.inverse(self.curv) 1177 c = LA.inv(self.curv) 1178 except LA.LinAlgError: 1179 die.warn("Singular covariance matrix") 1180 return None 1181 return self.curv_scale[:,Num.NewAxis] * c * self.curv_scale[Num.NewAxis,:]

1182 1183

1184 - def run(self, maxsteps=None, misc={}, Ta=None):

1185 """This is the normal top-level routine. 1186 You should feel free to make your own, though.""" 1187 if Ta is None: 1188 Ta = self.Td() 1189 term = 0.0 1190 while term < 3 and (maxsteps is None or self.steps<maxsteps): 1191 self.step() 1192 self.print_ev_diag(vectors=1) 1193 ttmp = self.terminate() 1194 if ttmp < 0: 1195 term = term/2.0 1196 else: 1197 term += self.terminate() 1198 self.print_at(misc=misc) 1199 if Ta > 0: 1200 self.one_anneal(Ta) 1201 self.print_at(misc=misc)

1202 1203

1204 - def print_errbar(self):

1205 log_fd = self._log() 1206 c = self.covariance() 1207 assert c.shape[0] == c.shape[1] 1208 log_fd.write('#var:') 1209 for i in range(c.shape[0]): 1210 log_fd.write(" %6g" % c[i,i]) 1211 log_fd.write('\n') 1212 log_fd.flush()

1213 1214

1215 - def sim_anneal(self, nsteps, Ta=None, misc={}):

1216 """A top-level simulated annealing routine.""" 1217 if Ta is None: 1218 Ta = self.Td() 1219 self.measure() 1220 self.steps = 0 1221 no = 1 + self.np/2 1222 while nsteps is None or self.steps<nsteps: 1223 if self.steps%self.np == self.np-1 : 1224 self.measure() 1225 self.print_errbar() 1226 self.print_ev_diag(vectors=1) 1227 if self.steps%no == no-1 : 1228 self.print_at(misc=misc) 1229 self.one_anneal(Ta) 1230 self.steps += 1 1231 self.print_at(misc=misc)

1232 1233

1234 - def print_at(self, misc={}):

1235 """Prints a file to match the old C++ optimizer output.""" 1236 log_fd = self._log() 1237 log_fd.write("#move:") 1238 p_step = self.p.p - self.last_p.p 1239 for i in range(p_step.shape[0]): 1240 log_fd.write(" %6g" % p_step[i]) 1241 log_fd.write('\n') 1242 o = misc.copy() 1243 o['E'] = self.sumsq() 1244 o['iter'] = self.steps 1245 o['lambda'] = self.lamb 1246 o['last'] = self.last 1247 log_fd.write('# ' + avio.concoct(o) + '\n') 1248 for i in range(self.p.p.shape[0]): 1249 log_fd.write(" %g" % self.p.p[i]) 1250 log_fd.write('\n') 1251 log_fd.flush()

1252 1253

1254 -def test1_fcn(p, args, *d):

1255 return p * (1+Num.arrayrange(p.shape[0]))

1256 1257

1258 -def test1():

1259 """Linear.""" 1260 # p = Num.array([1], Num.Float) 1261 p = Num.array([1, 99, -1, 2, 3, 4, 5, 0, 0, -45, 3, 4, 5, 2, 2, 1, 2], Num.Float) 1262 o = opt(p, 0.0001, test1_fcn) 1263 o.set_nthreads( 1 ) 1264 o.scale = Num.ones(p.shape) * 0.1 1265 o.run() 1266 print "p=", o.p 1267 print "cov=", o.covariance()

1268

1269 -def test1a_fcn(p, args, *d):

1270 pp = Num.array([p]) 1271 x = Num.matrixmultiply(pp, args) 1272 assert x.shape[0] == 1 1273 return x[0]

1274

1275 -def _errsize(a, b):

1276 q = a - b 1277 n = Num.sum(Num.ravel(a)**2) + Num.sum(Num.ravel(b)**2) 1278 s = Num.sum(Num.ravel(q)**2) / (0.5*n) 1279 return math.sqrt(s)

1280

1281 -def test1a():

1282 """Linear.""" 1283 # p = Num.array([1], Num.Float) 1284 p = Num.array([-0.2222, 0.1, 0.5, -1, 5, -0.3], Num.Float) 1285 q = RA.standard_normal((p.shape[0], p.shape[0])) 1286 qq = Num.matrixmultiply(q, Num.transpose(q)) 1287 o = opt(p, 0.0001, test1a_fcn, args=q) 1288 o.scale = Num.exp(RA.standard_normal((p.shape[0],))) 1289 o.set_nthreads( 1 ) 1290 o.run() 1291 print "p=", o.p 1292 print "cov=", o.covariance() 1293 print "cov_under=", o.covariance_under() 1294 # print "Inv cov=", LA.inverse(o.covariance()) 1295 print "Inv cov=", LA.inv(o.covariance()) 1296 print "qq=", qq 1297 # print "Inv qq=", LA.inverse(qq) 1298 print "Inv qq=", LA.inv(qq) 1299 # print "Inv cov_under=", LA.inverse(o.covariance_under()) 1300 print "Inv cov_under=", LA.inv(o.covariance_under()) 1301 # print "cov error=", _errsize(o.covariance(), LA.inverse(qq)) 1302 print "cov error=", _errsize(o.covariance(), LA.inv(qq)) 1303 # print "curv error=", _errsize( LA.inverse(o.covariance()), qq) 1304 print "curv error=", _errsize( LA.inv(o.covariance()), qq) 1305 # ce = o.covariance_under()-LA.inverse(qq) 1306 ce = o.covariance_under()-LA.inv(qq) 1307 print "ce=", ce 1308 # print "cov_under error=", _errsize( o.covariance_under(), LA.inverse(qq)) 1309 print "cov_under error=", _errsize( o.covariance_under(), LA.inv(qq)) 1310 # print "curv_under error=", _errsize( LA.inverse(o.covariance_under()), qq) 1311 print "curv_under error=", _errsize( LA.inv(o.covariance_under()), qq)

1312

1313 -def test2_fcn(p, args, *d):

1314 return Num.sin(p) - Num.cos(p)

1315 1316

1317 -def test2():

1318 """Nonlinear""" 1319 p = Num.array([0, 0.2, 0.4, 0.1, 0.3, 0.5, 0.6, 0.7, 0.8, 0.77, 0.76, 0.33, 0.01], Num.Float) 1320 o = opt(p, 0.0001, test2_fcn) 1321 # o.set_nthreads( 1 ) 1322 o.scale = Num.ones(p.shape) * 0.1 1323 o.run() 1324 print "p=", o.p

1325 1326

1327 -def test3_fcn(p, args, *d):

1328 o = Num.zeros(p.shape, Num.Float) 1329 for i in range(p.shape[0]): 1330 s = 0.0 1331 for j in range(i): 1332 s += p[i] - p[j]*(j+1) 1333 o[i] = s + p[i] 1334 return o

1335

1336 -def test3():

1337 """Linear, but correlated""" 1338 p = Num.array([0, 0.2, 0.4, 0.1, 0.3, 0.5, 0.6, 0.7, 0.8, 0.77, 0.76, 0.33, 0.01], Num.Float) 1339 o = opt(p, 0.0001, test3_fcn) 1340 o.scale = Num.ones(p.shape) * 0.1 1341 o.run() 1342 print "steps=", o.steps 1343 print "p=", o.p

1344 1345 1346

1347 -def test4_fcn(p, args, *d):

1348 if p[0]<0 or p[1]<0: 1349 return None 1350 return [(p[0]-p[1])*100 + p[0], p[1] ]

1351 1352

1353 -def test4():

1354 """Linear, but with a constraint""" 1355 p = Num.array([1, 0.1], Num.Float) 1356 o = opt(p, 0.0001, test4_fcn) 1357 o.set_nthreads( 1 ) 1358 o.run() 1359 print "steps=", o.steps 1360 print "p=", o.p

1361 1362

1363 -def test5_fcn(p, args, *d):

1364 o = Num.zeros(p.shape, Num.Float) 1365 for i in range(p.shape[0]): 1366 s = 0.0 1367 for j in range(i): 1368 s += p[i] - p[j]*(j+1) 1369 o[i] = s + p[i] 1370 if gpkmisc.N_maximum(o) > 30 or gpkmisc.N_minimum(o)<-30: 1371 return None 1372 return Num.exp(o) - 1

1373 1374

1375 -def test5():

1376 """Very Nonlinear""" 1377 p = Num.array([0, 4.2, 1.4, 0.1, 0.8, -1.77, 0.01], 1378 Num.Float) 1379 o = opt(p, 0.0001, test5_fcn) 1380 o.set_nthreads(1) 1381 o.verbose = 1 1382 o.scale = Num.ones(p.shape) * 0.1 1383 try: 1384 o.run() 1385 except NoDownhill, q: 1386 die.warn("No downhill step:" + str(q)) 1387 except: 1388 raise 1389 print "p=", o.p

1390 1391

1392 -def test6_fcn(p, args, *d):

1393 o = Num.zeros((3,), Num.Float) 1394 v = p[0] - p[1]*p[1] 1395 o[0] = v + 0.001*p[0] 1396 o[1] = v + 0.001*p[1] 1397 o[2] = v - 0.001*p[1] 1398 print "test_fcn: p=", p, "o=", o 1399 return o

1400 1401 1402

1403 -def test6():

1404 """Curved valley.""" 1405 p = Num.array([1, 8], Num.Float) 1406 o = opt(p, 0.001, test6_fcn) 1407 o.set_nthreads( 1 ) 1408 o.scale = Num.ones(p.shape) * 0.1 1409 try: 1410 o.run() 1411 except NoDownhill, q: 1412 die.warn("No downhill step:" + str(q)) 1413 except: 1414 raise 1415 print "p=", o.p

1416 1417

1418 -def _test7_fcn(p, args, *d):

1419 if p[0]<0 or p[1]<0: 1420 return None 1421 return [(p[0]-p[1])*100, p[1]+1 ]

1422 1423

1424 -def linconst_min(nparam, param, min):

1425 """Generate a linear constraint to be added onto a list and passed 1426 to linear_constraint(). This constraint expresses that p[param]>=min.""" 1427 1428 vec = Num.zeros((nparam)) 1429 vec[param] = 1. 1430 shift = -1. * min 1431 return vec, shift

1432 1433

1434 -def linconst_max(nparam, param, max):

1435 """Generate a linear constraint to be added onto a list and passed 1436 to linear_constraint(). This constraint expresses that p[param]<=max.""" 1437 1438 vec = Num.zeros((nparam)) 1439 vec[param] = -1. 1440 shift = 1. * max 1441 return vec, shift

1442 1443 1444

1445 -def linear_constraint(op, np, list_of_constraints):

1446 """This constrains a step to lie inside a region bounded by a list_of_constraints. 1447 @param op: parameters at the beginning of the step. 1448 @param np: target parameters. 1449 @param list_of_constraints: C{[ (vector, shift), ... ]} where the allowable region of each 1450 constraint is defined by C{Num.dot(x, vector)+shift >= 0} . 1451 C{Vector} is a C{numpy} array. 1452 """ 1453 # die.dbg("constraint( %s -> %s )" % (str(op), str(np))) 1454 dir = np-op 1455 possibles = list_of_constraints[:] 1456 while len(possibles)>0: 1457 # First, we find the nearest constraint that we will encounter. 1458 eta = None 1459 whichc = None 1460 for i in range(len(possibles)): 1461 cdir, cdist = possibles[i] 1462 distance = Num.dot(op, cdir) + cdist 1463 if distance < 0: 1464 return None # Whoops. Starting point (old_parameters) fails constraint. 1465 progress = Num.dot(cdir, dir) 1466 if progress >= 0: # Will never hit constraint, assuming we're outside already. 1467 # die.dbg("Wrong dir to hit constraint %d: %s %g" % (i, str(cdir), cdist)) 1468 continue 1469 new_eta = -distance / progress 1470 # print "i=", i, "distance=", distance, "progress=", progress, "new eta=", new_eta 1471 if new_eta > 1: # Won't get there on this step... 1472 # die.dbg("Constraint %d: %s %g is too far" % (i, str(cdir), cdist)) 1473 continue 1474 if eta is None or eta > new_eta: 1475 eta = new_eta 1476 whichc = i 1477 if whichc is None: # We won't hit any more constraints. 1478 # die.dbg("No more constraints: %d left" % len(possibles)) 1479 return np 1480 cdir, cdist = possibles.pop(whichc) 1481 # print "Dir=", dir, "eta=", eta, "np=", np 1482 np = op + dir*eta # Move to the intersection. 1483 # die.dbg("Projecting via %d to %s" % (whichc, str(np))) 1484 assert abs(Num.dot(np, cdir)+cdist) < 1e-6, "Should go to edge" 1485 # Future moves will be in the plane of the constraint: 1486 dir -= cdir*Num.dot(dir, cdir)/Num.dot(cdir, cdir) 1487 # die.dbg("Used all constraints") 1488 return np

1489 1490

1491 -def test7():

1492 """Constraint.""" 1493 c = [(Num.array((1, 0), Num.Float), 0), (Num.array((0, 1), Num.Float), 0) ] 1494 p = Num.array([1, 8], Num.Float) 1495 o = opt(p, 0.001, _test7_fcn) 1496 o.set_nthreads( 1 ) 1497 o.constrain = lambda op, np, args, loc=c: linear_constraint(op, np, loc) 1498 o.run() 1499 print "p=", o.p

1500 1501 1502 1503 1504 if __name__ == '__main__': 1505 test1a() 1506

Source Code for Module gmisclib.opt