Changeset 104
- Timestamp:
- 02/20/09 18:36:22 (7 years ago)
- Files:
-
- 1 added
- 1 deleted
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
DEV_NOTES
r103 r104 5 5 */*/*py: 6 6 - move "circle" to mystic.models, and move more of "mogi" to mystic.models 7 mystic/mystic/ fmin_powell.py:7 mystic/mystic/scipy_optimize.py: 8 8 - derive an AbstractSolver class for all solvers to inherit from 9 - add scipy solver "Brent", to _remove_ the scipy dependency (*) 10 - migrate all contents into scipy_optimize.py 11 - clean up & edit inline documentation 12 - better way to deal with callback, direc, maxfun, maxiter (?) 13 mystic/mystic/scipy_optimize.py: 14 - add the "__all__" scipy trick 9 - better way to deal with direc, maxfun, maxiter (?) 15 10 mystic/mystic/differential_evolution.py: 16 11 - add ability to choose initial point -
mystic/examples/test_rosenbrock3.py
r103 r104 29 29 30 30 #from scipy.optimize import fmin_powell 31 from mystic. fmin_powellimport fmin_powell, PowellDirectionalSolver31 from mystic.scipy_optimize import fmin_powell, PowellDirectionalSolver 32 32 #print fmin_powell(rosen,x0,retall=0,full_output=0) 33 33 solver = PowellDirectionalSolver(len(x0)) -
mystic/mystic/Make.mm
r103 r104 29 29 differential_evolution.py \ 30 30 scipy_optimize.py \ 31 fmin_powell.py \31 _scipy060optimize.py \ 32 32 termination.py \ 33 33 strategy.py \ -
mystic/mystic/__init__.py
r103 r104 16 16 # solvers 17 17 import differential_evolution, scipy_optimize 18 import fmin_powell #FIXME: move to scipy_optimize after remove scipy dependency19 18 20 19 # strategies, termination conditions -
mystic/mystic/differential_evolution.py
r102 r104 42 42 43 43 """ 44 __all__ = ['DifferentialEvolutionSolver','DifferentialEvolutionSolver2'] 44 45 45 46 from mystic.tools import Null, wrap_function -
mystic/mystic/scipy_optimize.py
r103 r104 4 4 # (derives from optimize.py module by Travis E. Oliphant) 5 5 # 6 # adapted from scipy.optimize.fmin, scipy version 0.4.86 # adapted scipy.optimize.fmin (from scipy version 0.4.8) 7 7 # by Patrick Hung, Caltech. 8 8 # 9 # adapted to class (& added bounds) 10 # updated to scipy version 0.6.0 9 # adapted from function to class (& added bounds) 10 # adapted scipy.optimize.fmin_powell 11 # updated solvers to scipy version 0.6.0 11 12 # by mmckerns@caltech.edu 12 13 13 14 """ 14 Algorithm adapted from scipy.optimize,15 Algorithms adapted from scipy.optimize, 15 16 16 17 fmin --- Nelder-Mead Simplex algorithm. … … 19 20 -- StepMonitor 20 21 22 fmin_powell --- Powell's Directional Search method. 23 Takes two additional input args 24 -- EvaluationMonitor 25 -- StepMonitor 26 21 27 """ 28 __all__ = ['NelderMeadSimplexSolver','PowellDirectionalSolver', 29 'fmin','fmin_powell'] 30 22 31 23 32 from mystic.tools import Null, wrap_function … … 27 36 from numpy import atleast_1d, eye, zeros, shape, \ 28 37 asarray, absolute, sqrt, Inf, asfarray 29 from numpy import clip 38 from numpy import clip, squeeze 30 39 31 40 abs = absolute 41 42 from _scipy060optimize import brent #XXX: local copy to avoid dependency! 43 32 44 33 45 class NelderMeadSimplexSolver(object): … … 176 188 ExtraArgs -- extra arguments for func. 177 189 178 """ 190 Further Inputs: 191 192 callback -- an optional user-supplied function to call after each 193 iteration. It is called as callback(xk), where xk is the 194 current parameter vector 195 disp -- non-zero to print convergence messages. 196 radius -- percentage change for initial simplex values 197 198 """ 179 199 # set arg names to scipy.optimize.fmin names; set fixed inputs 180 200 x0 = self.population[0] … … 185 205 callback=None #user-supplied function, called after each step 186 206 radius=0.05 #percentage change for initial simplex values 207 if kwds.has_key('callback'): callback = kwds['callback'] 187 208 if kwds.has_key('disp'): disp = kwds['disp'] 188 209 if kwds.has_key('radius'): radius = kwds['radius'] 189 if kwds.has_key('callback'): callback = kwds['callback']190 210 #------------------------------------------------------------- 191 211 … … 431 451 return retlist 432 452 453 ############################################################################ 454 455 def _linesearch_powell(func, p, xi, tol=1e-3): 456 # line-search algorithm using fminbound 457 # find the minimium of the function 458 # func(x0+ alpha*direc) 459 def myfunc(alpha): 460 return func(p + alpha * xi) 461 alpha_min, fret, iter, num = brent(myfunc, full_output=1, tol=tol) 462 xi = alpha_min*xi 463 return squeeze(fret), p+xi, xi 464 465 466 class PowellDirectionalSolver(NelderMeadSimplexSolver): #FIXME: not a simplex solver 467 """ 468 Powell Direction Search optimization adapted from scipy.optimize.fmin_powell. 469 """ 470 471 def __init__(self, dim): 472 """ 473 Takes one initial input: 474 dim -- dimensionality of the problem 475 """ 476 NelderMeadSimplexSolver.__init__(self,dim) 477 self._direc = None #FIXME: this is the easy way to return 'direc'... 478 #FIXME: NO SIMPLEX, so maybe best not to inherit __init__? 479 480 481 def Solve(self, func, termination, 482 maxiter=None, maxfun=None, sigint_callback=None, 483 EvaluationMonitor=Null, StepMonitor=Null, ExtraArgs=(), **kwds): 484 """Minimize a function using modified Powell's method. 485 486 Description: 487 488 Uses a modified Powell Directional Search algorithm to find 489 the minimum of function of one or more variables. 490 491 Inputs: 492 493 func -- the Python function or method to be minimized. 494 termination -- callable object providing termination conditions. 495 496 Additional Inputs: 497 498 maxiter -- the maximum number of iterations to perform. 499 maxfun -- the maximum number of function evaluations. 500 sigint_callback -- callback function for signal handler. 501 EvaluationMonitor -- a callable object that will be passed x, fval 502 whenever the cost function is evaluated. 503 StepMonitor -- a callable object that will be passed x, fval 504 after the end of a simplex iteration. 505 ExtraArgs -- extra arguments for func. 506 507 Further Inputs: 508 509 callback -- an optional user-supplied function to call after each 510 iteration. It is called as callback(xk), where xk is the 511 current parameter vector 512 direc -- initial direction set 513 xtol -- line-search error tolerance. 514 disp -- non-zero to print convergence messages. 515 516 """ 517 518 # set arg names to scipy.optimize.fmin_powell names; set fixed inputs 519 x0 = self.population[0] 520 args = ExtraArgs 521 full_output=1 #non-zero if fval and warnflag outputs are desired. 522 disp=0 #non-zero to print convergence messages. 523 retall=0 #non-zero to return all steps 524 direc=None 525 callback=None #user-supplied function, called after each step 526 xtol=1e-4 #line-search error tolerance 527 if kwds.has_key('callback'): callback = kwds['callback'] 528 if kwds.has_key('direc'): direc = kwds['direc'] #XXX: best interface? 529 if kwds.has_key('xtol'): xtol = kwds['xtol'] 530 if kwds.has_key('disp'): disp = kwds['disp'] 531 #------------------------------------------------------------- 532 533 import signal 534 import mystic.termination as detools 535 detools.EARLYEXIT = False 536 537 fcalls, func = wrap_function(func, args, EvaluationMonitor) 538 if self._useStrictRange: 539 x0 = self._setGuessWithinRangeBoundary(x0) 540 func = wrap_bounds(func, self._strictMin, self._strictMax) 541 542 def handler(signum, frame): 543 import inspect 544 print inspect.getframeinfo(frame) 545 print inspect.trace() 546 while 1: 547 s = raw_input(\ 548 """ 549 550 Enter sense switch. 551 552 sol: Write current best solution. 553 cont: Continue calculation. 554 call: Executes sigint_callback [%s]. 555 exit: Exits with current best solution. 556 557 >>> """ % sigint_callback) 558 if s.lower() == 'sol': 559 print "sw1." 560 print self.bestSolution 561 elif s.lower() == 'cont': 562 return 563 elif s.lower() == 'call': 564 # sigint call_back 565 if sigint_callback is not None: 566 sigint_callback(self.bestSolution) 567 elif s.lower() == 'exit': 568 detools.EARLYEXIT = True 569 return 570 else: 571 print "unknown option : %s ", s 572 573 self.signal_handler = handler 574 575 if self._handle_sigint: signal.signal(signal.SIGINT, self.signal_handler) 576 #------------------------------------------------------------- 577 578 x = asfarray(x0).flatten() 579 if retall: 580 allvecs = [x] 581 N = len(x) #XXX: this should be equal to self.nDim 582 rank = len(x.shape) 583 if not -1 < rank < 2: 584 raise ValueError, "Initial guess must be a scalar or rank-1 sequence." 585 if maxiter is None: 586 maxiter = N * 1000 587 if maxfun is None: 588 maxfun = N * 1000 589 self._maxiter = maxiter #XXX: better to just copy the code? 590 self._maxfun = maxfun #XXX: better to just copy the code? 591 592 if direc is None: 593 direc = eye(N, dtype=float) 594 else: 595 direc = asarray(direc, dtype=float) 596 fval = squeeze(func(x)) 597 x1 = x.copy() 598 599 self._direc = direc #XXX: instead, use a monitor? 600 self.bestSolution = x 601 self.bestEnergy = fval 602 self.population = [x] #XXX: pointless, if simplex not used 603 self.popEnergy = [fval] #XXX: pointless, if simplex not used 604 self.energy_history.append(self.bestEnergy) 605 StepMonitor(x,fval) # get initial values 606 607 iter = 0; 608 ilist = range(N) 609 610 CONTINUE = True 611 while CONTINUE: 612 fx = fval 613 bigind = 0 614 delta = 0.0 615 for i in ilist: 616 direc1 = direc[i] 617 fx2 = fval 618 fval, x, direc1 = _linesearch_powell(func, x, direc1, tol=xtol*100) 619 if (fx2 - fval) > delta: 620 delta = fx2 - fval 621 bigind = i 622 623 iter += 1 624 if callback is not None: 625 callback(x) 626 if retall: 627 allvecs.append(x) 628 629 self.energy_history.append(fval) #XXX: the 'best' for now... 630 if detools.EARLYEXIT or termination(self): CONTINUE = False #break 631 elif fcalls[0] >= maxfun: CONTINUE = False #break 632 elif iter >= maxiter: CONTINUE = False #break 633 634 else: # Construct the extrapolated point 635 direc1 = x - x1 636 x2 = 2*x - x1 637 x1 = x.copy() 638 fx2 = squeeze(func(x2)) 639 640 if (fx > fx2): 641 t = 2.0*(fx+fx2-2.0*fval) 642 temp = (fx-fval-delta) 643 t *= temp*temp 644 temp = fx-fx2 645 t -= delta*temp*temp 646 if t < 0.0: 647 fval, x, direc1 = _linesearch_powell(func, x, direc1, tol=xtol*100) 648 direc[bigind] = direc[-1] 649 direc[-1] = direc1 650 651 self.energy_history[-1] = fval #...update to 'best' energy 652 653 self._direc = direc #XXX: instead, use a monitor? 654 self.bestSolution = x 655 self.bestEnergy = fval 656 self.population = [x] #XXX: pointless, if simplex not used 657 self.popEnergy = [fval] #XXX: pointless, if simplex not used 658 StepMonitor(x,fval) # get ith values; #XXX: should be [x],[fval] ? 659 660 self.generations = iter 661 signal.signal(signal.SIGINT,signal.default_int_handler) 662 663 # code below here is dead, unless disp!=0 664 warnflag = 0 665 666 if fcalls[0] >= maxfun: 667 warnflag = 1 668 if disp: 669 print "Warning: Maximum number of function evaluations has "\ 670 "been exceeded." 671 elif iter >= maxiter: 672 warnflag = 2 673 if disp: 674 print "Warning: Maximum number of iterations has been exceeded" 675 else: 676 if disp: 677 print "Optimization terminated successfully." 678 print " Current function value: %f" % fval 679 print " Iterations: %d" % iter 680 print " Function evaluations: %d" % fcalls[0] 681 682 x = squeeze(x) 683 684 if full_output: 685 retlist = x, fval, direc, iter, fcalls[0], warnflag 686 if retall: 687 retlist += (allvecs,) 688 else: 689 retlist = x 690 if retall: 691 retlist = (x, allvecs) 692 693 return #retlist 694 695 696 def fmin_powell(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None, 697 maxfun=None, full_output=0, disp=1, retall=0, callback=None, 698 direc=None): 699 """fmin_powell using the 'original' scipy.optimize.fmin_powell interface""" 700 701 #FIXME: need to resolve "direc" 702 # - should just pass 'direc', and then hands-off ? How return it ? 703 704 from mystic.tools import Sow 705 stepmon = Sow() 706 evalmon = Sow() 707 from mystic.termination import NormalizedChangeOverGeneration as NCOG 708 709 solver = PowellDirectionalSolver(len(x0)) 710 solver.SetInitialPoints(x0) 711 #solver.enable_signal_handler() 712 solver.Solve(func,termination=NCOG(ftol),\ 713 maxiter=maxiter,maxfun=maxfun,\ 714 EvaluationMonitor=evalmon,StepMonitor=stepmon,\ 715 xtol=xtol, callback=callback, \ 716 disp=disp, direc=direc) #XXX: last two lines use **kwds 717 solution = solver.Solution() 718 719 # code below here pushes output to scipy.optimize.fmin_powell interface 720 #x = list(solver.bestSolution) 721 x = solver.bestSolution 722 fval = solver.bestEnergy 723 warnflag = 0 724 fcalls = len(evalmon.x) 725 iterations = len(stepmon.x) - 1 726 allvecs = [] 727 for i in range(len(stepmon.x)): 728 #allvecs.append(list(stepmon.x[i])) 729 allvecs.append(stepmon.x[i]) 730 direc = solver._direc #FIXME: better way to get direc from Solve() ? 731 732 if fcalls >= solver._maxfun: 733 warnflag = 1 734 elif iterations >= solver._maxiter: 735 warnflag = 2 736 737 x = squeeze(x) #FIXME: write squeezed x to stepmon instead? 738 739 if full_output: 740 retlist = x, fval, direc, iterations, fcalls, warnflag 741 if retall: 742 retlist += (allvecs,) 743 else: 744 retlist = x 745 if retall: 746 retlist = (x, allvecs) 747 748 return retlist 749 433 750 434 751 if __name__=='__main__':
Note: See TracChangeset
for help on using the changeset viewer.