|
| 1 | +#!/usr/bin/env python3 |
| 2 | +# |
| 3 | +# Copyright 2025 ImpactX contributors |
| 4 | +# Authors: Axel Huebl, Chad Mitchell |
| 5 | +# License: BSD-3-Clause-LBNL |
| 6 | +# |
| 7 | +# -*- coding: utf-8 -*- |
| 8 | + |
| 9 | +import time |
| 10 | + |
| 11 | +import matplotlib.pyplot as plt |
| 12 | +import numpy as np |
| 13 | +from scipy.optimize import minimize |
| 14 | + |
| 15 | +from impactx import my_run |
| 16 | + |
| 17 | +verbose = False |
| 18 | +# mode = "forward" |
| 19 | +mode = "backward" |
| 20 | +inputs_file_beam = "examples/fodo_space_charge/input_fodo_envelope_sc.in" |
| 21 | + |
| 22 | +history = {} |
| 23 | +runtime = {} |
| 24 | +optimizer = None |
| 25 | +grads = False |
| 26 | + |
| 27 | + |
| 28 | +class ProcessTimer: |
| 29 | + def __init__(self): |
| 30 | + self.elapsed_time = 0 |
| 31 | + |
| 32 | + def __enter__(self): |
| 33 | + self.start_time = time.process_time_ns() |
| 34 | + return self |
| 35 | + |
| 36 | + def __exit__(self, exc_type, exc_val, exc_tb): |
| 37 | + end_time = time.process_time_ns() |
| 38 | + self.elapsed_time = end_time - self.start_time |
| 39 | + |
| 40 | + |
| 41 | +def objective(parameters: tuple) -> float: |
| 42 | + """ |
| 43 | + A function that is evaluated by the optimizer. |
| 44 | +
|
| 45 | + Parameters |
| 46 | + ---------- |
| 47 | + parameters: tuple |
| 48 | + quadrupole strengths k of quad 1 and quad 2. |
| 49 | +
|
| 50 | + Returns |
| 51 | + ------- |
| 52 | + The L2 norm of alpha and beta of the beam at the end of the |
| 53 | + simulation. |
| 54 | + """ |
| 55 | + global grads |
| 56 | + |
| 57 | + if verbose: |
| 58 | + print(f"Run objective with parameters={parameters}...") |
| 59 | + |
| 60 | + if not grads: |
| 61 | + use_mode = "gradient-free" |
| 62 | + else: |
| 63 | + use_mode = mode |
| 64 | + |
| 65 | + q1_k, q2_k = parameters |
| 66 | + |
| 67 | + with ProcessTimer() as time_ns: |
| 68 | + values = my_run(q1_k, q2_k, use_mode, inputs_file_beam, verbose) |
| 69 | + # print(f"optimizer (grads={grads}): {time_ns.elapsed_time}ns") |
| 70 | + error = values["error"] |
| 71 | + # print(q1_k, q2_k, error, values) |
| 72 | + |
| 73 | + history[optimizer].append([q1_k, q2_k]) |
| 74 | + runtime[optimizer] = min(time_ns.elapsed_time, runtime[optimizer]) |
| 75 | + |
| 76 | + if verbose: |
| 77 | + print(f"error={error}, q1_k={q1_k}, q2_k={q2_k}") |
| 78 | + |
| 79 | + if np.isnan(error): |
| 80 | + error = 1.0e99 |
| 81 | + |
| 82 | + if grads: |
| 83 | + return (error, [values["derror_dq1_k"], values["derror_dq2_k"]]) |
| 84 | + else: |
| 85 | + return error |
| 86 | + |
| 87 | + |
| 88 | +def optimize_and_plot(opti, jac=None, plot=True): |
| 89 | + global grads, history, optimizer, runtime |
| 90 | + grads = jac == True # noqa |
| 91 | + optimizer = opti |
| 92 | + history[optimizer] = [] |
| 93 | + runtime[optimizer] = 9e99 |
| 94 | + res = minimize( |
| 95 | + objective, |
| 96 | + initial_quad_strengths, |
| 97 | + method=optimizer, |
| 98 | + jac=jac, |
| 99 | + tol=1.0e-8, |
| 100 | + options=options, |
| 101 | + bounds=bounds, |
| 102 | + ) |
| 103 | + ln = len(history[optimizer]) |
| 104 | + print(optimizer, ln) |
| 105 | + q1_k, q2_k = zip(*history[optimizer]) |
| 106 | + if plot: |
| 107 | + opt_str = optimizer |
| 108 | + if optimizer == "CG": |
| 109 | + opt_str = "Conjugate Gradient" |
| 110 | + plt.plot(q1_k, q2_k, label=f"{opt_str}: {ln}x") |
| 111 | + print(" Optimal parameters for k:", res.x) |
| 112 | + print(" L2 norm of alpha & beta at the optimum:", res.fun) |
| 113 | + print(f" Min. runtime: {runtime[optimizer] / ln}ns") |
| 114 | + |
| 115 | + |
| 116 | +if __name__ == "__main__": |
| 117 | + # Initial guess for the quadrople strengths |
| 118 | + initial_quad_strengths = np.array([-85, 80]) |
| 119 | + |
| 120 | + # Bounds for values to test: (min, max) |
| 121 | + positive = (0, None) |
| 122 | + negative = (None, 0) |
| 123 | + bounds = [negative, positive] |
| 124 | + |
| 125 | + # optimizer specific values |
| 126 | + # https://docs.scipy.org/doc/scipy/reference/optimize.minimize-neldermead.html |
| 127 | + # https://docs.scipy.org/doc/scipy/reference/optimize.minimize-lbfgsb.html |
| 128 | + options = { |
| 129 | + "maxiter": 1000, |
| 130 | + } |
| 131 | + |
| 132 | + fig = plt.figure(figsize=(4.5, 2.2)) |
| 133 | + |
| 134 | + # Call the optimizer |
| 135 | + optimize_and_plot("Nelder-Mead", False) |
| 136 | + optimize_and_plot("CG", "2-point", False) |
| 137 | + |
| 138 | + # gradient-based: CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr |
| 139 | + optimize_and_plot("CG", True) |
| 140 | + |
| 141 | + # analytical result: |
| 142 | + # q1_k = -103.12574100336 |
| 143 | + # q2_k = -q1_k |
| 144 | + |
| 145 | + fig.gca().annotate( |
| 146 | + "initial quadrupole strengths", # Annotation text |
| 147 | + xy=(-85, 80), # Point to annotate (arrow points here) |
| 148 | + xytext=(-95.8, 62), # Text position (arrow starts here) |
| 149 | + arrowprops=dict( |
| 150 | + facecolor="black", shrink=0.1, width=1.5, headwidth=4, headlength=6 |
| 151 | + ), |
| 152 | + ) |
| 153 | + fig.gca().annotate( |
| 154 | + "solution", # Annotation text |
| 155 | + xy=(-103.12574100336, 103.12574100336), # Point to annotate (arrow points here) |
| 156 | + xytext=(-106, 113), # Text position (arrow starts here) |
| 157 | + arrowprops=dict( |
| 158 | + facecolor="black", shrink=0.1, width=1.5, headwidth=4, headlength=6 |
| 159 | + ), |
| 160 | + ) |
| 161 | + |
| 162 | + plt.gca().set_xlim(-108, -80) |
| 163 | + plt.gca().set_ylim(60, 120) |
| 164 | + plt.xlabel(r"$k_1$ [1/m$^2$]") |
| 165 | + plt.ylabel(r"$k_2$ [1/m$^2$]") |
| 166 | + plt.legend(loc="upper right") |
| 167 | + plt.tight_layout() |
| 168 | + plt.show() |
0 commit comments