source: branches/decorate/test_timed_monitor.py @ 855

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

updated copyright to 2016

Line 
1#!/usr/bin/env python
2#
3# Author: Mike McKerns (mmckerns @caltech and @uqfoundation)
4# Copyright (c) 2012-2016 California Institute of Technology.
5# License: 3-clause BSD.  The full license text is available at:
6#  - http://mmckerns.github.io/project/mystic/browser/mystic/LICENSE
7
8#######################################################################
9# scaling and mpi info; also optimizer configuration parameters
10# hard-wired: use DE solver, don't use mpi, F-F' calculation
11#######################################################################
12scale = 1.0
13#XXX: <mpi config goes here>
14
15npop = 20
16maxiter = 1000
17maxfun = 1e+6
18convergence_tol = 1e-4
19crossover = 0.9
20percent_change = 0.9
21
22
23#######################################################################
24# the model function
25#######################################################################
26from surrogate import marc_surr as model
27from surrogate import ballistic_limit as limit
28
29
30#######################################################################
31# the subdiameter calculation
32#######################################################################
33def costFactory(i):
34  """a cost factory for the cost function"""
35
36  def cost(rv):
37    """compute the diameter as a calculation of cost
38
39  Input:
40    - rv -- 1-d array of model parameters
41
42  Output:
43    - diameter -- scale * | F(x) - F(x')|**2
44    """
45
46    # prepare x and xprime
47    params = rv[:-1]                         #XXX: assumes Xi' is at rv[-1]
48    params_prime = rv[:i]+rv[-1:]+rv[i+1:-1] #XXX: assumes Xi' is at rv[-1]
49
50    # get the F(x) response
51    Fx = model(params)
52
53    # get the F(x') response
54    Fxp = model(params_prime)
55
56    # compute diameter
57    return -scale * (Fx - Fxp)**2
58
59  return cost
60
61
62#######################################################################
63# the differential evolution optimizer
64#######################################################################
65def optimize(cost,lb,ub):
66  from mystic.solvers import DifferentialEvolutionSolver2
67  from mystic.termination import CandidateRelativeTolerance as CRT
68  from mystic.strategy import Best1Exp
69  from mystic.monitors import VerboseMonitor, Monitor
70  from mystic.tools import getch, random_seed
71
72  random_seed(123)
73
74 #stepmon = VerboseMonitor(100)
75 #stepmon = Monitor()
76  from timer import TimedMonitor
77  stepmon = TimedMonitor()
78  stepmon.timer.stamp = False
79  evalmon = Monitor()
80
81  ndim = len(lb) # [(1 + RVend) - RVstart] + 1
82
83  solver = DifferentialEvolutionSolver2(ndim,npop)
84  solver.SetRandomInitialPoints(min=lb,max=ub)
85  solver.SetStrictRanges(min=lb,max=ub)
86  solver.SetEvaluationLimits(maxiter,maxfun)
87  solver.SetEvaluationMonitor(evalmon)
88  solver.SetGenerationMonitor(stepmon)
89
90  tol = convergence_tol
91  solver.Solve(cost,termination=CRT(tol,tol),strategy=Best1Exp, \
92               CrossProbability=crossover,ScalingFactor=percent_change)
93
94  print "solved: %s" % solver.bestSolution
95  print "generations: %s" % len(stepmon.x)
96  print "time: %s seconds" % stepmon.timer._t
97  print "-"*70
98  diameter_squared = -solver.bestEnergy / scale  #XXX: scale != 0
99  func_evals = solver.evaluations
100  return diameter_squared, func_evals
101
102
103#######################################################################
104# loop over model parameters to calculate concentration of measure
105#######################################################################
106def UQ(start,end,lower,upper):
107  diameters = []
108  function_evaluations = []
109  total_func_evals = 0
110  total_diameter = 0.0
111
112  for i in range(start,end+1):
113    lb = lower + [lower[i]]
114    ub = upper + [upper[i]]
115 
116    #construct cost function and run optimizer
117    cost = costFactory(i)
118    subdiameter, func_evals = optimize(cost,lb,ub) #XXX: no initial conditions
119
120    function_evaluations.append(func_evals)
121    diameters.append(subdiameter)
122
123    total_func_evals += function_evaluations[-1]
124    total_diameter += diameters[-1]
125
126  print "subdiameters (squared): %s" % diameters
127  print "diameter (squared): %s" % total_diameter
128  print "func_evals: %s => %s" % (function_evaluations, total_func_evals)
129
130  return total_diameter
131
132
133#######################################################################
134# rank, bounds, and restart information
135#######################################################################
136if __name__ == '__main__':
137  from math import sqrt
138
139  function_name = "marc_surr"
140  lower_bounds = [60.0, 0.0, 2.1]
141  upper_bounds = [105.0, 30.0, 2.8]
142# h = thickness = [60,105]
143# a = obliquity = [0,30]
144# v = speed = [2.1,2.8]
145
146  RVstart = 0; RVend = 2
147  RVmax = len(lower_bounds) - 1
148
149  # when not a random variable, set the value to the lower bound
150  for i in range(0,RVstart):
151    upper_bounds[i] = lower_bounds[i]
152  for i in range(RVend+1,RVmax+1):
153    upper_bounds[i] = lower_bounds[i]
154
155  lbounds = lower_bounds[RVstart:1+RVend]
156  ubounds = upper_bounds[RVstart:1+RVend]
157
158  print "...SETTINGS..."
159  print "npop = %s" % npop
160  print "maxiter = %s" % maxiter
161  print "maxfun = %s" % maxfun
162  print "convergence_tol = %s" % convergence_tol
163  print "crossover = %s" % crossover
164  print "percent_change = %s" % percent_change
165  print "..............\n\n"
166
167  print " model: f(x) = %s(x)" % function_name
168  param_string = "["
169  for i in range(RVmax+1):
170    param_string += "'x%s'" % str(i+1)
171    if i == (RVmax):
172      param_string += "]"
173    else:
174      param_string += ", "
175
176  print " parameters: %s" % param_string
177  print "  varying 'xi', with i = %s" % range(RVstart+1,RVend+2)
178  print " lower bounds: %s" % lower_bounds
179  print " upper bounds: %s" % upper_bounds
180# print " ..."
181  try: model.load()
182  except: pass
183  diameter = UQ(RVstart,RVend,lower_bounds,upper_bounds)
184  try: model.dump()
185  except: pass
186  try: print model.info()
187  except: pass
188
189# EOF
Note: See TracBrowser for help on using the repository browser.