source: mystic/examples/example11.py @ 855

Revision 855, 4.2 KB checked in by mmckerns, 5 months ago (diff)

updated copyright to 2016

  • Property svn:executable set to *
Line 
1#!/usr/bin/env python
2#
3# Author: Mike McKerns (mmckerns @caltech and @uqfoundation)
4# Copyright (c) 1997-2016 California Institute of Technology.
5# License: 3-clause BSD.  The full license text is available at:
6#  - http://trac.mystic.cacr.caltech.edu/project/mystic/browser/mystic/LICENSE
7"""
8Example:
9    - Solve 8th-order Chebyshev polynomial coefficients with Nelder-Mead.
10    - Callable plot of fitting to Chebyshev polynomial.
11    - Plot (x2) of convergence to Chebyshev polynomial.
12    - Monitor (x2) Chi-Squared for Chebyshev polynomial.
13
14Demonstrates:
15    - standard models
16    - expanded solver interface
17    - parameter bounds constraints
18    - solver interactivity
19    - customized monitors and termination conditions
20"""
21
22# Nelder-Mead Simplex solver
23from mystic.solvers import NelderMeadSimplexSolver
24
25# Chebyshev polynomial and cost function
26from mystic.models.poly import chebyshev8, chebyshev8cost
27from mystic.models.poly import chebyshev8coeffs
28
29# tools
30from mystic.termination import CandidateRelativeTolerance as CRT
31from mystic.monitors import VerboseMonitor, Monitor
32from mystic.tools import getch
33from mystic.math import poly1d
34import pylab
35pylab.ion()
36
37# draw the plot
38def plot_frame(label=None):
39    pylab.close()
40    pylab.title("8th-order Chebyshev coefficient convergence")
41    pylab.xlabel("Nelder-Mead Simplex Solver %s" % label)
42    pylab.ylabel("Chi-Squared")
43    pylab.draw()
44    return
45 
46# plot the polynomial trajectories
47def plot_params(monitor):
48    x = range(len(monitor))
49    y = monitor.y
50    pylab.plot(x,y,'b-')
51    pylab.axis([1,0.5*x[-1],0,y[1]],'k-')
52    pylab.draw()
53    return
54
55# draw the plot
56def plot_exact():
57    pylab.title("fitting 8th-order Chebyshev polynomial coefficients")
58    pylab.xlabel("x")
59    pylab.ylabel("f(x)")
60    import numpy
61    x = numpy.arange(-1.2, 1.2001, 0.01)
62    exact = chebyshev8(x)
63    pylab.plot(x,exact,'b-')
64    pylab.legend(["Exact"])
65    pylab.axis([-1.4,1.4,-2,8],'k-')
66    pylab.draw()
67    return
68 
69# plot the polynomial
70def plot_solution(params,style='y-'):
71    import numpy
72    x = numpy.arange(-1.2, 1.2001, 0.01)
73    f = poly1d(params)
74    y = f(x)
75    pylab.plot(x,y,style)
76    pylab.legend(["Exact","Fitted"])
77    pylab.axis([-1.4,1.4,-2,8],'k-')
78    pylab.draw()
79    return
80
81if __name__ == '__main__':
82
83    print "Nelder-Mead Simplex"
84    print "==================="
85
86    # initial guess
87    import random
88    from mystic.tools import random_seed
89    random_seed(123)
90    ndim = 9
91    x0 = [random.uniform(-5,5) + chebyshev8coeffs[i] for i in range(ndim)]
92
93    # suggest that the user interacts with the solver
94    print "NOTE: while solver is running, press 'Ctrl-C' in console window"
95    getch()
96
97    # draw frame and exact coefficients
98    plot_exact()
99
100    # select parameter bounds constraints
101    from numpy import inf
102    min_bounds = [  0,-1,-300,-1,  0,-1,-100,-inf,-inf]
103    max_bounds = [200, 1,   0, 1,200, 1,   0, inf, inf]
104
105    # configure monitors
106    stepmon = VerboseMonitor(100)
107    evalmon = Monitor()
108
109    # use Nelder-Mead to solve 8th-order Chebyshev coefficients
110    solver = NelderMeadSimplexSolver(ndim)
111    solver.SetInitialPoints(x0)
112    solver.SetEvaluationLimits(generations=999)
113    solver.SetEvaluationMonitor(evalmon)
114    solver.SetGenerationMonitor(stepmon)
115    solver.SetStrictRanges(min_bounds,max_bounds)
116    solver.enable_signal_handler()
117    solver.Solve(chebyshev8cost, termination=CRT(1e-4,1e-4), \
118                 sigint_callback=plot_solution)
119    solution = solver.bestSolution
120
121    # get solved coefficients and Chi-Squared (from solver members)
122    iterations = solver.generations
123    cost = solver.bestEnergy
124    print "Generation %d has best Chi-Squared: %f" % (iterations, cost)
125    print "Solved Coefficients:\n %s\n" % poly1d(solver.bestSolution)
126
127    # compare solution with actual 8th-order Chebyshev coefficients
128    print "Actual Coefficients:\n %s\n" % poly1d(chebyshev8coeffs)
129
130    # plot solution versus exact coefficients
131    plot_solution(solution)
132    getch()
133
134    # plot convergence of coefficients per iteration
135    plot_frame('iterations')
136    plot_params(stepmon)
137    getch()
138
139    # plot convergence of coefficients per function call
140    plot_frame('function calls')
141    plot_params(evalmon)
142    getch()
143
144# end of file
Note: See TracBrowser for help on using the repository browser.