source: mystic/examples/example12.py @ 855

Revision 855, 3.0 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 5th-order polynomial coefficients with Powell's method.
10    - Plot of fitting to 5th-order polynomial.
11
12Demonstrates:
13    - model and cost function construction
14    - minimal solver interface
15"""
16
17# Powell's Directonal solver
18from mystic.solvers import fmin_powell
19
20# cost function factory
21from mystic.forward_model import CostFactory
22
23# tools
24from mystic.tools import getch, random_seed
25random_seed(123)
26import pylab
27pylab.ion()
28
29# build the forward model
30def ForwardPolyFactory(coeff):
31    """generate a 1-D polynomial instance from a list of coefficients"""
32    from numpy import poly1d
33    return poly1d(coeff)
34
35# build the cost function
36def PolyCostFactory(evalpts,datapts,ndim):
37    """generate a cost function instance from evaluation points and data"""
38    from numpy import sum
39    F = CostFactory()
40    F.addModel(ForwardPolyFactory,ndim,"noisy_polynomial")
41    return F.getCostFunction(evalpts=evalpts,observations=datapts, \
42                             metric=lambda x: sum(x*x),sigma=1000.)
43
44# 'data' generators
45def data(params):
46    """generate 'data' from polynomial coefficients"""
47    from numpy import array
48    x = 0.1*(array([range(101)])-50.)[0]
49    fwd = ForwardPolyFactory(params)
50    return x,fwd(x)
51
52def noisy_data(params):
53    """generate noisy data from polynomial coefficients"""
54    from numpy import random
55    x,y = data(params)
56    y = [random.normal(0,1) + i for i in y]
57    return x,y
58
59# draw the plot
60def plot_frame(label=None):
61    pylab.close()
62    pylab.title("fitting noisy 5th-order polynomial coefficients")
63    pylab.xlabel("x")
64    pylab.ylabel("f(x)")
65    pylab.draw()
66    return
67 
68# plot the polynomial
69def plot_data(evalpts,datapts,style='k.'):
70    pylab.plot(evalpts,datapts,'%s' % style)
71    pylab.axis([0,5,0,50],'k-')
72    pylab.draw()
73    return
74   
75def plot_solution(params,style='b-'):
76    x,y = data(params)
77    plot_data(x,y,style)
78    pylab.legend(["Data","Fitted"])
79    pylab.draw()
80    return
81
82
83if __name__ == '__main__':
84
85    print "Powell's Method"
86    print "==============="
87
88    # target and initial guess
89    target = [-1.,4.,-5.,20.,5.]
90    x0     = [-1.,2.,-3.,10.,5.]
91
92    # generate 'observed' data
93    x,datapts = noisy_data(target)
94
95    # plot observed data
96    plot_frame()
97    plot_data(x,datapts)
98
99    # generate cost function
100    costfunction = PolyCostFactory(x,datapts,len(target))
101
102    # use Powell's method to solve 5th-order polynomial coefficients
103    solution = fmin_powell(costfunction,x0)
104
105    # compare solution with actual target 5th-order polynomial coefficients
106    print "\nSolved Coefficients:\n %s\n" % ForwardPolyFactory(solution)
107    print "Target Coefficients:\n %s\n" % ForwardPolyFactory(target)
108 
109    # plot solution versus target coefficients
110    plot_solution(solution)
111    getch()
112
113# end of file
Note: See TracBrowser for help on using the repository browser.