================== Sequence Optimizer ================== The menu selection :menuselection:`Pulse Sequence --> Sequence Optimizer...` opens the Optimization window. This window provides tools to optimize the currently loaded pulse sequence so that it performs a desired state transformation or implements a target unitary gate. .. drops:: I1z :width: 99% :sequence: Standard/PulDel100msec :class: border :window: OptimWindow :caption: The Sequence Optimization window. :link: The optimizer modifies only the pulse elements of the sequence (RF amplitudes and, optionally, durations), leaving delay periods unchanged. All optimization runs in real time, so the DROPS display updates as the algorithm iterates. Overview ======== A typical optimization workflow: 1. Load or create a pulse sequence. 2. Set the :guilabel:`Initial State` (shown at the bottom of the window -- this is the initial density operator, e.g. ``I1z``). 3. Enter a :guilabel:`Target State` expression (e.g. ``I1x``). 4. Select an :ref:`algorithm ` from the dropdown. 5. Choose a :ref:`cost function ` (PP or UR). 6. Click :guilabel:`Run` to start iterating, or :guilabel:`Step` for a single iteration. 7. Watch the cost decrease in the cost history plot. 8. When satisfied, :guilabel:`Dump` the result for external use. Buttons ======= :guilabel:`Step` Perform a single optimization iteration. If the optimizer is currently running, clicking :guilabel:`Step` stops it first. :guilabel:`Run` Start continuous optimization. The algorithm iterates repeatedly until convergence or until :guilabel:`Run` is clicked again to stop. The DROPS display updates each iteration so you can watch the pulse sequence evolve. :guilabel:`Reset` Reset the sequence to its original un-optimized parameters. This clears the cost history and step counter and recompiles the pulse sequence. :guilabel:`Dump` Export the complete optimization state to a file. See :ref:`Export Format` below. Status Display ============== The optimization window shows three status values at the top: :guilabel:`Cost` The current best cost value (lower is better). A perfectly optimized sequence has a cost near ``-1.0`` for normalized targets. :guilabel:`|f'|` The maximum gradient norm across all frames. When this value approaches zero, the optimizer has found a stationary point. :guilabel:`Step` The iteration counter. Below the parameters, a multi-line status field shows algorithm-specific information (e.g. line-search progress, convergence criteria). At the bottom, a **cost history** plot shows how the cost evolved over iterations. Cost Functions ============== The cost function measures how well the optimized sequence achieves the desired goal. Each ensemble group can use a different cost function, selected from the dropdown next to the target expression. Point-to-Point (PP) -------------------- .. math:: C_\text{PP} = -2\,\text{Re}\langle \rho_T \,|\, \rho(T) \rangle The PP cost measures the overlap between the target density operator :math:`\rho_T` and the final state :math:`\rho(T)` at the end of the pulse sequence. Use PP for **state-to-state transfer** tasks, for example preparing a specific quantum state from thermal equilibrium. The :guilabel:`Target State` field accepts any operator expression (e.g. ``I1x``, ``I1xI2y``, ``0.5*I1z + 0.5*I2z``). Universal Rotation (UR) ----------------------- .. math:: C_\text{UR} = -\frac{|\text{Tr}(C^\dagger U)|^2}{N^2} The UR cost measures how well the cumulative propagator :math:`U` matches a target unitary gate generated by :math:`C`, **up to a global phase**. Because of the modulus, :math:`U` and :math:`e^{i\phi} C` are treated as equivalent for any phase :math:`\phi`. This is the phase-insensitive (SU) form of the GRAPE quality factor; see [Khaneja2005]_ and [deFouquieres2011]_. The :guilabel:`Target State` field specifies the generator of the target unitary (e.g. ``pi/2*I1x`` for a 90-degree rotation about x). .. note:: In practice the phase-insensitive form has been found to converge less reliably than the phase-sensitive form below; prefer **UR-PSU** for new optimizations. Universal Rotation, phase-sensitive (UR-PSU) -------------------------------------------- .. math:: C_\text{UR-PSU} = -\frac{\text{Re}\,\text{Tr}(C^\dagger U)}{N} The phase-sensitive form drops the modulus from the original UR cost. The minimum :math:`C_\text{UR-PSU} = -1` is reached when :math:`U = C` exactly (global phase included); :math:`U = -C` gives :math:`+1`. This is the phase-sensitive (PSU) form of the GRAPE quality factor recommended in [Khaneja2005]_. The corresponding gradient is .. math:: \frac{\partial C_\text{UR-PSU}}{\partial u_k(m)} = -\frac{1}{N}\,\text{Re}\, \text{Tr}\!\left(\Lambda_{m+1}^{\dagger}\, \frac{\partial U_m}{\partial u_k(m)}\,P_m\right) where :math:`P_m = U_{m-1}\cdots U_0` is the forward propagator and :math:`\Lambda_{m+1} = U_{m+1}^{\dagger}\cdots U_{M-1}^{\dagger}\,C` is the backward co-state. Use UR (or UR-PSU) for **unitary gate design** — for example, designing a broadband excitation pulse or a spin-lock sequence. .. note:: The PP cost function uses the target as a density operator, while the UR / UR-PSU cost functions interpret it as the generator of a unitary rotation. Switching cost type without updating the target expression may give unexpected results. Optimization Algorithms ======================= Select the algorithm from the :guilabel:`Optimizer` dropdown at the top of the window. Each algorithm exposes its own parameters below the dropdown. Gradient Descent ----------------- Simple first-order gradient descent with optional Armijo backtracking line search. This is the default algorithm and a good starting point for most problems. Gradients of the cost with respect to the controls are computed by reverse-mode automatic differentiation, equivalent to the GRAPE recursion of [Khaneja2005]_. Parameters: :guilabel:`Linesearch` (checkbox, default: on) When enabled, the step size is determined automatically by an Armijo backtracking line search. When disabled, the step size is fixed at the :guilabel:`alpha` value. :guilabel:`alpha` (slider, default: 0.01) The step size parameter. When line search is enabled, this value is updated each iteration to show the step size that was selected. When line search is disabled, this is the fixed step size -- larger values converge faster but may overshoot. Convergence criteria: - Gradient norm < :math:`10^{-6}`, or - Cost improvement < :math:`10^{-8}` with decreasing cost. Conjugate Gradient ------------------ Nonlinear conjugate gradient method using the Polak-Ribiere-Polyak (PR+) formula with Powell restart. This typically converges faster than gradient descent, especially for problems with many control parameters. The algorithm maintains a conjugate search direction that combines the current gradient with information from previous iterations. A Powell restart resets to steepest descent when the search directions lose conjugacy. Parameters are the same as gradient descent (:guilabel:`Linesearch` and :guilabel:`alpha`). Interior Point -------------- An interior-point method using the `Ipopt `_ library [WaechterBiegler2006]_ with L-BFGS Hessian approximation. This is a second-order method that handles constraints natively and often converges in fewer iterations than first-order methods. Application of quasi-Newton methods to NMR pulse design is discussed in [deFouquieres2011]_. .. note:: Interior Point is only available when SpinDrops is built with Ipopt support (``SPIND_DISABLE_OPTIM=OFF``). The algorithm enforces RF amplitude magnitude constraints directly (rather than by projection as in GD/CG), so it naturally respects hardware limits. It supports warm-starting across multiple :guilabel:`Step` clicks. Optimize Time ============= The :guilabel:`Optimize time` checkbox enables optimization of the duration of each pulse segment in addition to the RF amplitudes. When enabled, the optimizer can adjust how long each time step lasts, subject to duration bounds defined in the pulse sequence. When disabled (the default), only RF amplitudes are optimized and all timings remain fixed. Ensemble Optimization ===================== When multiple :doc:`ensemble ` groups are defined, the optimizer window shows a separate :guilabel:`Target State` field and :guilabel:`Cost Type` dropdown for each group. This allows different groups to have different target states and even different cost functions. The optimizer finds a **single set of controls** (one pulse sequence) that works well across all ensemble points simultaneously. The overall cost is the worst (maximum) cost across all frames and groups. This enables designing pulse sequences that are **robust** against parameter variations (e.g. B1 inhomogeneity or chemical shift offsets), or that simultaneously achieve multiple transfer goals for different ensemble members. Classic examples of such broadband / B1-robust pulses include the BURBOP family [Skinner2003]_ [Kobzar2008]_ and broadband universal rotations more generally [Kobzar2012]_. See :ref:`Ensemble and Optimization` in the manual for more details on configuring ensemble groups. Control Constraints =================== During optimization, the controls are subject to the following constraints: **RF amplitude pairs** For each pair of {x, y} RF channels, the magnitude is constrained to :math:`\sqrt{x^2 + y^2} \leq c_\text{max}`, where :math:`c_\text{max}` is the maximum RF amplitude defined in the pulse sequence. Gradient descent and conjugate gradient enforce this by projection after each step; interior point handles it as a native nonlinear constraint. **Offset / drift channels** Clamped to ``[controlmin, controlmax]``. **Time / duration** When :ref:`Optimize Time` is enabled, each segment duration is clamped to ``[controlmin, controlmax]`` as defined in the sequence. Export Format ============= The :guilabel:`Dump` button exports the optimization state to a text file (default name ``optim.txt``). The file uses a simple ``key=value`` format compatible with `Octave `_ and Matlab. For each frame (ensemble point) in the simulation, the file contains: **cost** The current cost value for this frame. **step** The iteration counter. **params** The optimizer parameters as a JSON object (e.g. ``{"Linesearch": true, "alpha": 0.01}``). **costhist** The cost at each iteration, formatted as a row vector. **x** The control matrix. Each column is one time step; each row is a control parameter. The first rows are pairs of {x, y} RF channels corresponding to the RF operators in the experiment. After the RF channels come any user-defined :ref:`Hamiltonian channels `, and the last row is always the duration (in seconds) of that time step. **dx** The gradient of the cost function with respect to the controls ``x``. At a true optimum, these values should be close to zero. Non-zero values at convergence may indicate a local minimum. **xmin** Lower bounds on the control values. **xmax** Upper bounds on the control values. **xon** A boolean matrix indicating which controls are active (i.e. included in the optimization). Pulse elements are active; delays are not. **acqN** If the pulse sequence contains acquisition elements (either a table-sequence :ref:`Table Acq ` or :ref:`go `), each acquisition is included separately. **acq** The total accumulated acquisition across all frames. Optimizer References ==================== .. [Khaneja2005] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, S. J. Glaser, *Optimal Control of Coupled Spin Dynamics: Design of NMR Pulse Sequences by Gradient Ascent Algorithms*, J. Magn. Reson. **172**, 296-305 (2005). `doi:10.1016/j.jmr.2004.11.004 `_ .. [deFouquieres2011] P. de Fouquieres, S. G. Schirmer, S. J. Glaser, I. Kuprov, *Second order gradient ascent pulse engineering*, J. Magn. Reson. **212**, 412-417 (2011). `doi:10.1016/j.jmr.2011.07.023 `_ .. [Skinner2003] T. E. Skinner, T. O. Reiss, B. Luy, N. Khaneja, S. J. Glaser, *Application of optimal control theory to the design of broadband excitation pulses for high-resolution NMR*, J. Magn. Reson. **163**, 8-15 (2003). `doi:10.1016/S1090-7807(03)00153-8 `_ .. [Kobzar2008] K. Kobzar, B. Luy, N. Khaneja, S. J. Glaser, *Pattern pulses: design of arbitrary excitation profiles as a function of pulse amplitude and offset*, J. Magn. Reson. **194**, 58-66 (2008). `doi:10.1016/j.jmr.2008.05.023 `_ .. [Kobzar2012] K. Kobzar, S. Ehni, T. E. Skinner, S. J. Glaser, B. Luy, *Exploring the limits of broadband 90° and 180° universal rotation pulses*, J. Magn. Reson. **225**, 142-160 (2012). `doi:10.1016/j.jmr.2012.09.013 `_ .. [WaechterBiegler2006] A. Wächter, L. T. Biegler, *On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming*, Mathematical Programming **106**, 25-57 (2006). `doi:10.1007/s10107-004-0559-y `_