#
source:
DEV_TODO
@
464

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

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 |

8 | from mystic.models import rosen as my_model |

9 | x0 = [0.8, 1.2, 0.7] |

10 | |

11 | # configure the solver and obtain the solution |

12 | from mystic.solvers import fmin |

13 | solution = fmin(my_model, x0) |

14 | |

15 | ##################################################################### |

16 | |

17 | # the function to be minimized and initial values |

18 | from mystic.models import rosen as my_model |

19 | x0 = [0.8, 1.2, 0.7] |

20 | |

21 | # get monitor and termination condition objects |

22 | from mystic.monitors import Monitor, VerboseMonitor |

23 | stepmon = VerboseMonitor(5) |

24 | evalmon = Monitor() |

25 | from mystic.termination import ChangeOverGeneration |

26 | COG = ChangeOverGeneration() |

27 | |

28 | # instantiate and configure the solver |

29 | from mystic.solvers import NelderMeadSimplexSolver |

30 | solver = NelderMeadSimplexSolver(len(x0)) |

31 | solver.SetInitialPoints(x0) |

32 | solver.SetGenerationMonitor(stepmon) |

33 | solver.SetEvaluationMonitor(evalmon) |

34 | solver.Solve(my_model, COG) |

35 | |

36 | # obtain the solution |

37 | solution = solver.bestSolution |

38 | |

39 | # obtain diagnostic information |

40 | function_evals = solver.evaluations |

41 | iterations = solver.generations |

42 | cost = solver.bestEnergy |

43 | |

44 | # modify the solver configuration, and continue |

45 | COG = ChangeOverGeneration(tolerance=1e-8) |

46 | solver.Solve(my_model, COG) |

47 | |

48 | # obtain the new solution |

49 | solution = solver.bestSolution |

50 | |

51 | ##################################################################### |

52 | |

53 | # a user-provided constraints function |

54 | def constrain(x): |

55 | x[1] = x[0] |

56 | return x |

57 | |

58 | # the function to be minimized and the bounds |

59 | from mystic.models import rosen as my_model |

60 | lb = [0.0, 0.0, 0.0] |

61 | ub = [2.0, 2.0, 2.0] |

62 | |

63 | # get termination condition object |

64 | from mystic.termination import ChangeOverGeneration |

65 | COG = ChangeOverGeneration() |

66 | |

67 | # instantiate and configure the solver |

68 | from mystic.solvers import NelderMeadSimplexSolver |

69 | solver = NelderMeadSimplexSolver(len(x0)) |

70 | solver.SetRandomInitialPoints(lb, ub) |

71 | solver.SetStrictRanges(lb, ub) |

72 | solver.SetConstraints(constrain) |

73 | solver.Solve(my_model, COG) |

74 | |

75 | # obtain the solution |

76 | solution = solver.bestSolution |

77 | |

78 | ##################################################################### |

79 | |

80 | # a user-provided constraints function |

81 | constraints = """ |

82 | x2 = x1 |

83 | """ |

84 | from mystic.constraints import parse |

85 | constrain = parse(constraints) |

86 | |

87 | ##################################################################### |

88 | |

89 | # generate a model from a stock 'model factory' |

90 | from mystic.models.lorentzian import Lorentzian |

91 | lorentz = Lorentzian(coeffs) |

92 | |

93 | # evaluate the model |

94 | y = lorentz(x) |

95 | |

96 | ##################################################################### |

97 | |

98 | # a user-provided model function |

99 | def identify(x) |

100 | return x |

101 | |

102 | # add pathos infrastructure (included in mystic) |

103 | from mystic.tools import modelFactory, Monitor |

104 | evalmon = Monitor() |

105 | my_model = modelFactory(identify, monitor=evalmon) |

106 | |

107 | # evaluate the model |

108 | y = my_model(x) |

109 | |

110 | # evaluate the model with a map function |

111 | from mystic.tools import PythonMap |

112 | my_map = PythonMap() |

113 | z = my_map(my_model, range(10)) |

114 | |

115 | ##################################################################### |

116 | |

117 | # a user-provided model function |

118 | def identify(x) |

119 | return x |

120 | |

121 | # cast the model as a distributed service |

122 | from pathos.servers import sshServer |

123 | id = 'foo.caltech.edu:50000:spike42' |

124 | my_proxy = sshServer(identify, server=id) |

125 | |

126 | # evaluate the model via proxy |

127 | y = my_proxy(x) |

128 | |

129 | ##################################################################### |

130 | |

131 | # a user-provided model function |

132 | def identify(x) |

133 | return x |

134 | |

135 | # select and configure a parallel map |

136 | from pathos.maps import ipcPool |

137 | my_map = ipcPool(2, servers=['foo.caltech.edu']) |

138 | |

139 | # evaluate the model in parallel |

140 | z = my_map(identify, range(10)) |

141 | |

142 | ##################################################################### |

143 | |

144 | # configure and build map |

145 | from pathos.launchers import ipc |

146 | from pathos.strategies import pool |

147 | from pathos.tools import mapFactory |

148 | my_map = mapFactory(launcher=ipc, strategy=pool) |

149 | |

150 | ##################################################################### |

151 | |

152 | # establish a tunnel |

153 | from pathos.tunnel import sshTunnel |

154 | uid = 'foo.caltech.edu:12345:tunnel69' |

155 | tunnel_proxy = sshTunnel(uid) |

156 | |

157 | # inspect the ports |

158 | localport = tunnel_proxy.lport |

159 | remoteport = tunnel_proxy.rport |

160 | |

161 | # a user-provided model function |

162 | def identify(x) |

163 | return x |

164 | |

165 | # cast the model as a distributed service |

166 | from pathos.servers import ipcServer |

167 | id = 'localhost:%s:bug01' % localport |

168 | my_proxy = ipcServer(identify, server=id) |

169 | |

170 | # evaluate the model via tunneled proxy |

171 | y = my_proxy(x) |

172 | |

173 | # disconnect the tunnel |

174 | tunnel_proxy.disconnect() |

175 | |

176 | ##################################################################### |

177 | |

178 | # configure and build map |

179 | from pyina.launchers import mpirun |

180 | from pyina.strategies import carddealer as card |

181 | from pyina.tools import mapFactory |

182 | my_map = mapFactory(4, launcher=mpirun, strategy=card) |

183 | |

184 | ##################################################################### |

185 | |

186 | # the function to be minimized and the bounds |

187 | from mystic.models import rosen as my_model |

188 | lb = [0.0, 0.0, 0.0] |

189 | ub = [2.0, 2.0, 2.0] |

190 | |

191 | # get termination condition object |

192 | from mystic.termination import ChangeOverGeneration |

193 | COG = ChangeOverGeneration() |

194 | |

195 | # select the parallel launch configuration |

196 | from pyina.maps import MpirunCarddealer |

197 | my_map = MpirunCarddealer(4) |

198 | |

199 | # instantiate and configure the solver |

200 | from mystic.solvers import DifferentialEvolutionSolver |

201 | solver = DifferentialEvolutionSolver(len(lb), 20) |

202 | solver.SetRandomInitialPoints(lb, ub) |

203 | solver.SetStrictRanges(lb, ub) |

204 | solver.SetEvaluationMap(my_map) |

205 | solver.Solve(my_model, COG) |

206 | |

207 | # obtain the solution |

208 | solution = solver.bestSolution |

209 | |

210 | ##################################################################### |

211 | |

212 | # the function to be minimized and the bounds |

213 | from mystic.models import rosen as my_model |

214 | lb = [0.0, 0.0, 0.0] |

215 | ub = [2.0, 2.0, 2.0] |

216 | |

217 | # get monitor and termination condition objects |

218 | from mystic.monitors import LoggingMonitor |

219 | stepmon = LoggingMonitor(1, âlog.txtâ) |

220 | from mystic.termination import ChangeOverGeneration |

221 | COG = ChangeOverGeneration() |

222 | |

223 | # select the parallel launch configuration |

224 | from pyina.maps import TorqueMpirunCarddealer |

225 | my_map = TorqueMpirunCarddealer(â5:ppn=4â) |

226 | |

227 | # instantiate and configure the nested solver |

228 | from mystic.solvers import PowellDirectionalSolver |

229 | my_solver = PowellDirectionalSolver(len(lb)) |

230 | my_solver.SetStrictRanges(lb, ub) |

231 | my_solver.SetEvaluationLimits(50) |

232 | |

233 | # instantiate and configure the outer solver |

234 | from mystic.solvers import BuckshotSolver |

235 | solver = BuckshotSolver(len(lb), 20) |

236 | solver.SetRandomInitialPoints(lb, ub) |

237 | solver.SetGenerationMonitor(stepmon) |

238 | solver.SetNestedSolver(my_solver) |

239 | solver.SetSolverMap(my_map) |

240 | solver.Solve(my_model, COG) |

241 | |

242 | # obtain the solution |

243 | solution = solver.bestSolution |

244 | |

245 | ##################################################################### |

246 | |

247 | # prepare a (F(X) - G)**2 a metric |

248 | def 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 | ''' |

258 | The calculation of the diameter is performed as a nested |

259 | optimization, as shown above for the BuckshotSolver. Each |

260 | inner optimization is a calculation of a component |

261 | suboscillation, using the a global optimizer |

262 | (such as DifferentialEvolutionSolver) and the cost |

263 | metric shown above. |

264 | ''' |

265 | |

266 | # prepare a (F(X) - F(X'))**2 cost metric |

267 | def 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 | ''' |

281 | Global optimizations used in solving OUQ problems are |

282 | composed in the same manner as shown above for the |

283 | DifferentialEvolutionSolver. |

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 |

288 | npts = (nx,ny,nz) |

289 | lb = (nx * w_lower) + (nx * x_lower) \ |

290 | + (ny * w_lower) + (ny * y_lower) \ |

291 | + (nz * w_lower) + (nz * z_lower) |

292 | ub = (nx * w_upper) + (nx * x_upper) \ |

293 | + (ny * w_upper) + (ny * y_upper) \ |

294 | + (nz * w_upper) + (nz * z_upper) |

295 | |

296 | from mystic.math.measures import split_param |

297 | from mystic.math.dirac_measure import product_measure |

298 | from mystic.math import almostEqual |

299 | |

300 | # split bounds into weight-only & sample-only |

301 | w_lb, m_lb = split_param(lb, npts) |

302 | w_ub, m_ub = split_param(ub, npts) |

303 | |

304 | # generate constraints function |

305 | def 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 |

324 | def cost(param): |

325 | prodmeasure = product_measure() |

326 | prodmeasure.load(param, npts) |

327 | return MINMAX * prodmeasure.pof(my_model) |

328 | |

329 | ##################################################################### |

330 | """ |

331 | DIRECT: |

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 | |

340 | INDIRECT: |

341 | * Build a failure test suite, and a proper test suite |

342 | * Try one of the VTF apps... or Sean's "cain" library |

343 | |

344 | REFERENCE: |

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.