1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
#!/usr/bin/env python
#
# Author: Mike McKerns (mmckerns @caltech and @uqfoundation)
# Copyright (c) 1997-2016 California Institute of Technology.
# Copyright (c) 2016-2024 The Uncertainty Quantification Foundation.
# License: 3-clause BSD. The full license text is available at:
# - https://github.com/uqfoundation/mystic/blob/master/LICENSE
"""
Example:
- Solve 5th-order polynomial coefficients with Powell's method.
- Plot of fitting to 5th-order polynomial.
Demonstrates:
- model and cost function construction
- minimal solver interface
"""
# Powell's Directonal solver
from mystic.solvers import fmin_powell
# cost function factory
from mystic.forward_model import CostFactory
# tools
from mystic.tools import getch, random_seed
random_seed(123)
import matplotlib.pyplot as plt
plt.ion()
# build the forward model
def ForwardPolyFactory(coeff):
"""generate a 1-D polynomial instance from a list of coefficients"""
from numpy import poly1d
return poly1d(coeff)
# build the cost function
def PolyCostFactory(evalpts,datapts,ndim):
"""generate a cost function instance from evaluation points and data"""
from numpy import sum
F = CostFactory()
F.addModel(ForwardPolyFactory,ndim,"noisy_polynomial")
return F.getCostFunction(evalpts=evalpts,observations=datapts, \
metric=lambda x: sum(x*x),sigma=1000.)
# 'data' generators
def data(params):
"""generate 'data' from polynomial coefficients"""
from numpy import array
x = 0.1*(array([list(range(101))])-50.)[0]
fwd = ForwardPolyFactory(params)
return x,fwd(x)
def noisy_data(params):
"""generate noisy data from polynomial coefficients"""
from numpy import random
x,y = data(params)
y = [random.normal(0,1) + i for i in y]
return x,y
# draw the plot
def plot_frame(label=None):
plt.close()
plt.title("fitting noisy 5th-order polynomial coefficients")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.draw()
plt.pause(0.001)
return
# plot the polynomial
def plot_data(evalpts,datapts,style='k.'):
plt.plot(evalpts,datapts,'%s' % style)
plt.axis([0,5,0,50])#,'k-')
plt.draw()
plt.pause(0.001)
return
def plot_solution(params,style='b-'):
x,y = data(params)
plot_data(x,y,style)
plt.legend(["Data","Fitted"])
plt.draw()
plt.pause(0.001)
return
if __name__ == '__main__':
print("Powell's Method")
print("===============")
# target and initial guess
target = [-1.,4.,-5.,20.,5.]
x0 = [-1.,2.,-3.,10.,5.]
# generate 'observed' data
x,datapts = noisy_data(target)
# plot observed data
plot_frame()
plot_data(x,datapts)
# generate cost function
costfunction = PolyCostFactory(x,datapts,len(target))
# use Powell's method to solve 5th-order polynomial coefficients
solution = fmin_powell(costfunction,x0)
# compare solution with actual target 5th-order polynomial coefficients
print("\nSolved Coefficients:\n %s\n" % ForwardPolyFactory(solution))
print("Target Coefficients:\n %s\n" % ForwardPolyFactory(target))
# plot solution versus target coefficients
plot_solution(solution)
getch()
# end of file
|
