Reference

implicit(solve, residual, x, p=(); drdy=drdy_forward, lsolve=linear_solve)

Make implicit function AD compatible (specifically with ForwardDiff and ReverseDiff).

Arguments

• solve::function: y = solve(x, p). Solve implicit function returning state variable y, for input variables x, and fixed parameters p.
• residual::function: Either r = residual(y, x, p) or in-place residual!(r, y, x, p). Set residual r (scalar or vector), given state y (scalar or vector), variables x (vector) and fixed parameters p (tuple).
• x::vector{float}: evaluation point.
• p::tuple: fixed parameters. default is empty tuple.
• drdy::function: drdy(residual, y, x, p). Provide (or compute yourself): ∂ri/∂yj. Default is forward mode AD.
• lsolve::function: lsolve(A, b). Linear solve A x = b (where A is computed in drdy and b is computed in jvp, or it solves A^T x = c where c is computed in vjp). Default is backslash operator.

implicit_linear(A, b; lsolve=linear_solve, Af=factorize)

Make implicit function AD compatible (specifically with ForwardDiff and ReverseDiff). This version is for linear equations Ay = b

Arguments

• A::matrix, b::vector: components of linear system $A y = b$
• lsolve::function: lsolve(A, b). Function to solve the linear system, default is backslash operator.
• Af::factorization: An optional factorization of A, useful to override default factorize, or if multiple linear solves will be performed with same A matrix.
• mmul: Function to compute $A*y$ for a vector $y$. Defaults to the julia multipy operator.

apply_factorization(A, factfun)

Apply a matrix factorization to the primal portion of a dual number. Avoids user from needing to add ForwardDiff as a dependency. Afactorization = factfun(A)

implicit_eigval(A, B, eigsolve)

Make eigenvalue problems AD compatible with ForwardDiff and ReverseDiff

Arguments

• A::matrix, B::matrix: generlized eigenvalue problem. A v = λ B v (B is identity for standard eigenvalue problem)
• eigsolve::function: λ, V, U = eigsolve(A, B). Function to solve the eigenvalue problem. λ is a vector containing eigenvalues. V is a matrix whose columns are the corresponding eigenvectors (i.e., λ[i] corresponds to V[:, i]). U is a matrix whose columns contain the left eigenvectors (u^H A = λ u^H B) The left eigenvectors must be in the same order as the right ones (i.e., U' * B * V must be diagonal). U can be provided with any normalization as we normalize internally s.t. U' * B * V = I If A and B are symmetric/Hermitian then U = V.

Returns

• λ::vector: eigenvalues and their derivatives. (Currently only eigenvalue derivatives are provided. not eigenvectors)

explicit_unsteady(initialize, onestep, solve, t, xd, xc, p=(); cache=nothing)

Make reverse-mode efficient for explicit ODEs Builds tape over each time step separately and analytically propagates between, rather than recording a long tape over the entire time sequence.

Arguments:

• initialize::function: Return initial state variables. y0 = initialize(t0, xd, xc0, p). May or may not depend on t0 (initial time), xd (design variables), xc0 (initial control variables), p (fixed parameters)
• onestep::function: define how states are updated (assuming one-step methods). y = onestep(yprev, t, tprev, xd, xci, p) or in-place onestep!(y, yprev, t, tprev, xd, xci, p). Set the next set of state variables y, given the previous state yprev, current time t, previous time tprev, design variables xd, current control variables xc, and fixed parameters p.
• t::vector{float}: time steps that simulation runs across
• xd::vector{float}: design variables, don't change in time, but do change from one run to the next (otherwise they would be parameters)
• xc::matrix{float}, size nxc x nt: control variables. xc[:, i] are the control variables at the ith time step.
• p::tuple: fixed parameters, i.e., they are always constant and so do not affect derivatives. Default is empty tuple.

Keyword Arguments:

• cache=nothing: see explicit_unsteady_cache. If computing derivatives more than once, you should compute the cache beforehand the save for later iterations. Otherwise, it will be created internally.

explicit_unsteady_cache(initialize, onestep!, ny, nxd, nxc, p=(); compile=false)

Initialize arrays and functions needed for explicit_unsteady reverse mode

Arguments

• initialize::function: Return initial state variables. y0 = initialize(t0, xd, xc0, p). May or may not depend on t0 (initial time), xd (design variables), xc0 (initial control variables), p (fixed parameters)
• onestep::function: define how states are updated (assuming one-step methods). y = onestep(yprev, t, tprev, xd, xci, p) or in-place onestep!(y, yprev, t, tprev, xd, xci, p). Set the next set of state variables y, given the previous state yprev, current time t, previous time tprev, design variables xd, current control variables for that time step xci, and fixed parameters p.
• ny::int: number of states
• nxd::int: number of design variables
• nxc::int: number of control variables (at a given time step)
• p::tuple: fixed parameters, i.e., they are always constant and so do not affect derivatives. Default is empty tuple.
• compile::bool: indicates whether the tape for the function onestep can be prerecorded. Will be much faster but should only be true if onestep does not contain any branches. Otherwise, ReverseDiff may return incorrect gradients.

implicit_unsteady(initialize, onestep, residual, t, xd, xc, p=(); cache=nothing, drdy=drdy_forward, lsolve=linear_solve)

Make reverse-mode efficient for implicit ODEs (solvers compatible with AD).

Arguments:

• initialize::function: Return initial state variables. y0 = initialize(t0, xd, xc0, p). May or may not depend on t0 (initial time), xd (design variables), xc0 (initial control variables), p (fixed parameters)
• onestep::function: define how states are updated (assuming one-step methods) including any solver. y = onestep(yprev, t, tprev, xd, xci, p) or in-place onestep!(y, yprev, t, tprev, xd, xci, p). Set the next set of state variables y, given the previous state yprev, current time t, previous time tprev, design variables xd, current control variables for that time step xci, and fixed parameters p.
• residual::function: define residual that is solved in onestep. Either r = residual(y, yprev, t, tprev, xd, xci, p) where variables are same as above or residual!(r, y, yprev, t, tprev, xd, xci, p).
• t::vector{float}: time steps that simulation runs across
• xd::vector{float}: design variables, don't change in time, but do change from one run to the next (otherwise they would be parameters)
• xc::matrix{float}, size nxc x nt: control variables. xc[:, i] are the control variables at the ith time step.
• p::tuple: fixed parameters, i.e., they are always constant and so do not affect derivatives. Default is empty tuple.

Keyword Arguments:

• cache=nothing: see implicit_unsteady_cache. If computing derivatives more than once, you should compute the cache beforehand the save for later iterations. Otherwise, it will be created internally.
• drdy: drdy(residual, r, y, yprev, t, tprev, xd, xci, p). Provide (or compute yourself): ∂ri/∂yj. Default is forward mode AD.
• lsolve::function: lsolve(A, b). Linear solve A x = b (where A is computed in drdy and b is computed in jvp, or it solves A^T x = c where c is computed in vjp). Default is backslash operator.

implicit_unsteady_cache(initialize, residual, ny, nxd, nxc, p=(); compile=false)

Initialize arrays and functions needed for implicit_unsteady reverse mode

Arguments

• initialize::function: Return initial state variables. y0 = initialize(t0, xd, xc0, p). May or may not depend on t0 (initial time), xd (design variables), xc0 (initial control variables), p (fixed parameters)
• residual::function: define residual that is solved in onestep. Either r = residual(y, yprev, t, tprev, xd, xci, p) where variables are same as above residual!(r, y, yprev, t, tprev, xd, xci, p).
• ny::int: number of states
• nxd::int: number of design variables
• nxc::int: number of control variables (at a given time step)
• p::tuple: fixed parameters, i.e., they are always constant and so do not affect derivatives. Default is empty tuple.
• compile::bool: indicates whether the tape for the function onestep can be prerecorded. Will be much faster but should only be true if onestep does not contain any branches. Otherwise, ReverseDiff may return incorrect gradients.

provide_rule(func, x, p=(); mode="ffd", jacobian=nothing, jvp=nothing, vjp=nothing)

Provide partials rather than rely on AD. For cases where AD is not available or to provide your own rule, or to use mixed mode, etc.

Arguments

• func::function, x::vector{float}, p::tuple: function of form: y = func(x, p), where x are variables and p are fixed parameters.
• mode::string:
  • "ffd": forward finite difference
  • "cfd": central finite difference
  • "cs": complex step
  • "jacobian": provide your own jacobian (see jacobian function)
  • "vp": provide jacobian vector product and vector jacobian product (see jvp and vjp)
• jacobian::function: only used if mode="jacobian". J = jacobian(x, p) provide Jij = dyi / dx_j
• jvp::function: only used if mode="vp" and in forward mode. ydot = jvp(x, p, v) provide Jacobian vector product J*v
• vjp::function: only used if mode="vp" and in revese mode. xbar = vjp(x, p, v) provide vector Jacobian product v'*J

derivativesetup(func, x, p, ad, compiletape=false)

Set up a derivative function to make it easier to repeatedly differentiate (e.g., in an optimization). Primary use case is to pass derivatives to Python or another language.

Arguments

• func::Function: The function to differentiate. Must have signature f = func(x, p) where x are the variables to differentiate with respect to and p are additional parameters.
• x::AbstractVector: a typical input vector, used to size the derivative arrays.
• p::Any: Additional parameters to pass to func. If this input changes, a new function must be setup.
• ad::String: The type of automatic differentiation to use. Options are:
  • "fjacobian": Forward mode AD to compute the full Jacobian df/dx.
  • "rjacobian": Reverse mode AD to compute the full Jacobian df/dx.
  • "jvp": Forward mode AD to compute the Jacobian-vector product df/dx * v.
  • "vjp": Reverse mode AD to compute the vector-Jacobian product v' * df/dx.
• compiletape::Bool: (optional, default=false) If using "rjacobian", whether to compile the tape for faster execution (assumes no branching behavior).

Returns

A function that computes the requested derivative. The signature depends on the type of derivative:

• For "fjacobian" and "rjacobian": df(J, x) where J is a preallocated array to hold the Jacobian.
• For "jvp": djvp(ydot, x, xdot) where ydot is a preallocated array to hold the output and xdot is the input direction vector.
• For "vjp": dvjp(xbar, x, ybar) where xbar is a preallocated array to hold the output and ybar is the input direction vector.