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 | """ |
---|
8 | Example: |
---|
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 | |
---|
14 | Demonstrates: |
---|
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 |
---|
23 | from mystic.solvers import NelderMeadSimplexSolver |
---|
24 | |
---|
25 | # Chebyshev polynomial and cost function |
---|
26 | from mystic.models.poly import chebyshev8, chebyshev8cost |
---|
27 | from mystic.models.poly import chebyshev8coeffs |
---|
28 | |
---|
29 | # tools |
---|
30 | from mystic.termination import CandidateRelativeTolerance as CRT |
---|
31 | from mystic.monitors import VerboseMonitor, Monitor |
---|
32 | from mystic.tools import getch |
---|
33 | from mystic.math import poly1d |
---|
34 | import pylab |
---|
35 | pylab.ion() |
---|
36 | |
---|
37 | # draw the plot |
---|
38 | def 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 |
---|
47 | def 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 |
---|
56 | def 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 |
---|
70 | def 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 | |
---|
81 | if __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 |
---|