[WIP] Gradient-based protocol optimization

[maxentile]

Proof-of-concept for gradient-based optimization of T.I. protocols that involve coordinated control of more than one parameter.

Shouldn't be merged into main until after #404 , but I'd like to open the PR at this stage in case feedback about how to pose the problem of protocol optimization "differentiably" should flow back up to the design in #404 , #403 , or #420 , and in case this provides any heuristic insights about how to improve handwritten protocols for distance-based decoupling in the meantime.

---

For context, I'd like to point to a paper that influenced my approach, but that I haven't directly implemented here: [Rotskoff, Crooks, Vanden-Eijnden, 2016] "A geometric approach to optimal nonequilibrium control: Minimizing dissipation in nanomagnetic spin systems" https://arxiv.org/abs/1607.07425 .

There, the authors represent protocols as discretized paths in a high-dimensional space of control parameters, estimate the thermodynamic metric tensor at each discretized point along the current path, then perform an efficient string-method-style update of that discretized path to reduce the corresponding protocol "length".

In contrast, here I've done things a bit differently, trying to lean as heavily as possible on autodiff and generic optimization:

• I represent the protocol as an explicit parameterized function that maps a scalar "lambda" to a high-dimensional vector of control parameters (rather than as a collection of high-dimensional "beads on a string"),
• rather than explicitly forming an estimate of a metric tensor (a square matrix) and multiplying it on both sides by the protocol "speed" vector (to get the scalar integrand that appears in the length functional of eq 5), I directly form a differentiable estimate of the stddev of the scalar du_dlam (distributed according to du_dlam(x), x ~ p(x; control_dials(lam; protocol_params))),
• rather than performing string-method-style updates of the locations of protocol "beads," I apply a generic scipy.optimize routine to a differentiable estimate of "protocol quality" as a function of protocol_params.

(Also be aware of a minor notational difference: in this PR "lambda" is the scalar that governs a vector of "control parameters", but in the referenced paper, "lambda" is a vector of control parameters, in turn governed by a different scalar. Also note I sometimes refer to "control parameters" as "control dials" when it helps distinguish from other uses of "parameters".)

---

I'd also like to point to Chapter 2 of Levi Naden's thesis, which described a framework for T.I. protocol optimization, and showed ways to efficiently estimate the variance of du_dlam for new trial protocols based on previous simulation results. Again, there are quite a few implementation differences, but this did strongly influence my approach.

---

For the toy problem constructed here, there are only two control parameters: independently controllable distance offsets for LJ and Coulomb terms. A protocol for this problem is a function specifying how these two control parameters depend on a single governing parameter "lambda."

For applied problems, timemachine exposes many possible control parameters, such as independently controllable distance offsets per particle, and independently interpolated parameters per force field term. (And derivatives of U w.r.t. each!) To take full advantage of this, we would like to be able to efficiently optimize protocols that adjust all of these controllable dials independently.

[maxentile]

Closing for now, may revisit this approach to protocol optimization later.