source: DEV_TODO @ 464

Revision 464, 10.2 KB checked in by mmckerns, 5 years ago (diff)

updated citation, todo's from scipy pub

Line 
1#####################################################################
2#   M.M. McKerns, L. Strand, T. Sullivan, A. Fang, M.A.G. Aivazis,
3#   "Building a framework for predictive science", Proceedings of
4#   the 10th Python in Science Conference, (submitted 2011).
5#####################################################################
6
7# the function to be minimized and initial values
8from mystic.models import rosen as my_model
9x0 = [0.8, 1.2, 0.7]
10
11# configure the solver and obtain the solution
12from mystic.solvers import fmin
13solution = fmin(my_model, x0)
14
15#####################################################################
16
17# the function to be minimized and initial values
18from mystic.models import rosen as my_model
19x0 = [0.8, 1.2, 0.7]
20
21# get monitor and termination condition objects
22from mystic.monitors import Monitor, VerboseMonitor
23stepmon = VerboseMonitor(5)
24evalmon = Monitor()
25from mystic.termination import ChangeOverGeneration
26COG = ChangeOverGeneration()
27
28# instantiate and configure the solver
29from mystic.solvers import NelderMeadSimplexSolver
30solver = NelderMeadSimplexSolver(len(x0))
31solver.SetInitialPoints(x0)
32solver.SetGenerationMonitor(stepmon)
33solver.SetEvaluationMonitor(evalmon)
34solver.Solve(my_model, COG)
35
36# obtain the solution
37solution = solver.bestSolution
38
39# obtain diagnostic information
40function_evals = solver.evaluations
41iterations = solver.generations
42cost = solver.bestEnergy
43
44# modify the solver configuration, and continue
45COG = ChangeOverGeneration(tolerance=1e-8)
46solver.Solve(my_model, COG)
47
48# obtain the new solution
49solution = solver.bestSolution
50
51#####################################################################
52
53# a user-provided constraints function
54def constrain(x):
55  x[1] = x[0]
56  return x
57
58# the function to be minimized and the bounds
59from mystic.models import rosen as my_model
60lb = [0.0, 0.0, 0.0]
61ub = [2.0, 2.0, 2.0]
62
63# get termination condition object
64from mystic.termination import ChangeOverGeneration
65COG = ChangeOverGeneration()
66
67# instantiate and configure the solver
68from mystic.solvers import NelderMeadSimplexSolver
69solver = NelderMeadSimplexSolver(len(x0))
70solver.SetRandomInitialPoints(lb, ub)
71solver.SetStrictRanges(lb, ub)
72solver.SetConstraints(constrain)
73solver.Solve(my_model, COG)
74
75# obtain the solution
76solution = solver.bestSolution
77
78#####################################################################
79
80# a user-provided constraints function
81constraints = """
82x2 = x1
83"""
84from mystic.constraints import parse
85constrain = parse(constraints)
86
87#####################################################################
88
89# generate a model from a stock 'model factory'
90from mystic.models.lorentzian import Lorentzian
91lorentz = Lorentzian(coeffs)
92
93# evaluate the model
94y = lorentz(x)
95
96#####################################################################
97
98# a user-provided model function
99def identify(x)
100  return x
101
102# add pathos infrastructure (included in mystic)
103from mystic.tools import modelFactory, Monitor
104evalmon = Monitor()
105my_model = modelFactory(identify, monitor=evalmon)
106
107# evaluate the model
108y = my_model(x)
109
110# evaluate the model with a map function
111from mystic.tools import PythonMap
112my_map = PythonMap()
113z = my_map(my_model, range(10))
114
115#####################################################################
116
117# a user-provided model function
118def identify(x)
119  return x
120
121# cast the model as a distributed service
122from pathos.servers import sshServer
123id = 'foo.caltech.edu:50000:spike42'
124my_proxy = sshServer(identify, server=id)
125
126# evaluate the model via proxy
127y = my_proxy(x)
128
129#####################################################################
130
131# a user-provided model function
132def identify(x)
133  return x
134
135# select and configure a parallel map
136from pathos.maps import ipcPool
137my_map = ipcPool(2, servers=['foo.caltech.edu'])
138
139# evaluate the model in parallel
140z = my_map(identify, range(10))
141
142#####################################################################
143
144# configure and build map
145from pathos.launchers import ipc
146from pathos.strategies import pool
147from pathos.tools import mapFactory
148my_map = mapFactory(launcher=ipc, strategy=pool)
149
150#####################################################################
151
152# establish a tunnel
153from pathos.tunnel import sshTunnel
154uid = 'foo.caltech.edu:12345:tunnel69'
155tunnel_proxy = sshTunnel(uid)
156
157# inspect the ports
158localport = tunnel_proxy.lport
159remoteport = tunnel_proxy.rport
160
161# a user-provided model function
162def identify(x)
163  return x
164
165# cast the model as a distributed service
166from pathos.servers import ipcServer
167id = 'localhost:%s:bug01' % localport
168my_proxy = ipcServer(identify, server=id)
169
170# evaluate the model via tunneled proxy
171y = my_proxy(x)
172
173# disconnect the tunnel
174tunnel_proxy.disconnect()
175
176#####################################################################
177
178# configure and build map
179from pyina.launchers import mpirun
180from pyina.strategies import carddealer as card
181from pyina.tools import mapFactory
182my_map = mapFactory(4, launcher=mpirun, strategy=card)
183
184#####################################################################
185
186# the function to be minimized and the bounds
187from mystic.models import rosen as my_model
188lb = [0.0, 0.0, 0.0]
189ub = [2.0, 2.0, 2.0]
190
191# get termination condition object
192from mystic.termination import ChangeOverGeneration
193COG = ChangeOverGeneration()
194
195# select the parallel launch configuration
196from pyina.maps import MpirunCarddealer
197my_map = MpirunCarddealer(4)
198
199# instantiate and configure the solver
200from mystic.solvers import DifferentialEvolutionSolver
201solver = DifferentialEvolutionSolver(len(lb), 20)
202solver.SetRandomInitialPoints(lb, ub)
203solver.SetStrictRanges(lb, ub)
204solver.SetEvaluationMap(my_map)
205solver.Solve(my_model, COG)
206
207# obtain the solution
208solution = solver.bestSolution
209
210#####################################################################
211
212# the function to be minimized and the bounds
213from mystic.models import rosen as my_model
214lb = [0.0, 0.0, 0.0]
215ub = [2.0, 2.0, 2.0]
216
217# get monitor and termination condition objects
218from mystic.monitors import LoggingMonitor
219stepmon = LoggingMonitor(1, ’log.txt’)
220from mystic.termination import ChangeOverGeneration
221COG = ChangeOverGeneration()
222
223# select the parallel launch configuration
224from pyina.maps import TorqueMpirunCarddealer
225my_map = TorqueMpirunCarddealer(’5:ppn=4’)
226
227# instantiate and configure the nested solver
228from mystic.solvers import PowellDirectionalSolver
229my_solver = PowellDirectionalSolver(len(lb))
230my_solver.SetStrictRanges(lb, ub)
231my_solver.SetEvaluationLimits(50)
232
233# instantiate and configure the outer solver
234from mystic.solvers import BuckshotSolver
235solver = BuckshotSolver(len(lb), 20)
236solver.SetRandomInitialPoints(lb, ub)
237solver.SetGenerationMonitor(stepmon)
238solver.SetNestedSolver(my_solver)
239solver.SetSolverMap(my_map)
240solver.Solve(my_model, COG)
241
242# obtain the solution
243solution = solver.bestSolution
244
245#####################################################################
246
247# prepare a (F(X) - G)**2 a metric
248def costFactory(my_model, my_data):
249  def cost(param):
250
251    # compute the cost
252    return ( my_model(param) - my_data )**2
253
254  return cost
255
256#####################################################################
257'''
258The calculation of the diameter is performed as a nested
259optimization, as shown above for the BuckshotSolver. Each
260inner optimization is a calculation of a component
261suboscillation, using the a global optimizer
262(such as DifferentialEvolutionSolver) and the cost
263metric shown above.
264'''
265
266# prepare a (F(X) - F(X'))**2 cost metric
267def suboscillationFactory(my_model, i):
268  def cost(param):
269
270    # get X and X' (Xi' is appended to X at param[-1])
271    x = param[:-1]
272    x_prime = param[:i] + param[-1:] + param[i+1:-1]
273
274    # compute the suboscillation
275    return -( my_model(x) - my_model(x_prime) )**2
276
277  return cost
278
279#####################################################################
280'''
281Global optimizations used in solving OUQ problems are
282composed in the same manner as shown above for the
283DifferentialEvolutionSolver.
284'''
285
286# OUQ requires bounds in a very specific form...
287# param = [wxi]*nx + [xi]*nx + [wyi]*ny + [yi]*ny + [wzi]*nz + [zi]*nz
288npts = (nx,ny,nz)
289lb = (nx * w_lower) + (nx * x_lower) \
290   + (ny * w_lower) + (ny * y_lower) \
291   + (nz * w_lower) + (nz * z_lower)
292ub = (nx * w_upper) + (nx * x_upper) \
293   + (ny * w_upper) + (ny * y_upper) \
294   + (nz * w_upper) + (nz * z_upper)
295
296from mystic.math.measures import split_param
297from mystic.math.dirac_measure import product_measure
298from mystic.math import almostEqual
299
300# split bounds into weight-only & sample-only
301w_lb, m_lb = split_param(lb, npts)
302w_ub, m_ub = split_param(ub, npts)
303
304# generate constraints function
305def constraints(param):
306  prodmeasure = product_measure()
307  prodmeasure.load(param, npts)
308
309  # impose norm on measures
310  for measure in prodmeasure:
311    if not almostEqual(float(measure.mass), 1.0):
312      measure.normalize()
313
314  # impose expectation on product measure
315  E = float(prodmeasure.get_expect(my_model))
316    if not (E <= float(target_mean + error)) \
317    or not (float(target_mean - error) <= E):
318      prodmeasure.set_expect((target_mean, error), my_model, (m_lb, m_ub))
319
320  # extract weights and positions
321  return prodmeasure.flatten()
322
323# generate maximizing function
324def cost(param):
325  prodmeasure = product_measure()
326  prodmeasure.load(param, npts)
327  return MINMAX * prodmeasure.pof(my_model)
328
329#####################################################################
330"""
331DIRECT:
332* Add more python optimizers: scipy, OpenOpt, PARK (snobfit)
333* Allow for derivative and gradient capture -> use Sow?
334* get "handler" to work in parallel
335* Better 'programmatic' interface for handler
336* Add more options to handler (i.e. toggle_verbosity?, get_cost?, ...?)
337* Allow sigint_callback to take a list (i.e. provide call[i])
338* Add "constraints" to models (design similar to pyre.inventory and validators)
339
340INDIRECT:
341* Build a failure test suite, and a proper test suite
342* Try one of the VTF apps... or Sean's "cain" library
343
344REFERENCE:
345* Look at PARK's rangemap.py for bounds and range mapping
346* Look at PARK's parameter.py, deps.py, expression.py, & assembly.py
347* <-- Find OpenOpt's model & optimizer API -->
348* <-- Find DAKOTA's model & optimizer API -->
349"""
Note: See TracBrowser for help on using the repository browser.