text:
SensitivityAlg(args...;kwargs...) = @error("The SensitivtyAlg choice mechanism was completely overhauled. Please consult the local sensitivity documentation for more information")

abstract type AbstractForwardSensitivityAlgorithm{CS,AD,FDT} <: DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT} end
abstract type AbstractAdjointSensitivityAlgorithm{CS,AD,FDT} <: DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT} end
abstract type AbstractSecondOrderSensitivityAlgorithm{CS,AD,FDT} <: DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT} end
abstract type AbstractShadowingSensitivityAlgorithm{CS,AD,FDT} <: DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT} end

"""
ForwardSensitivity{CS,AD,FDT} <: AbstractForwardSensitivityAlgorithm{CS,AD,FDT}

An implementation of continuous forward sensitivity analysis for propagating
derivatives by solving the extended ODE. When used within adjoint differentiation
(i.e. via Zygote), this will cause forward differentiation of the `solve` call
within the reverse-mode automatic differentiation environment.

## Constructor

```julia
function ForwardSensitivity(;
                            chunk_size=0,autodiff=true,
                            diff_type=Val{:central},
                            autojacvec=autodiff,
                            autojacmat=false)
```

## Keyword Arguments

* `autodiff`: Use automatic differentiation in the internal sensitivity algorithm
  computations. Default is `true`.
* `chunk_size`: Chunk size for forward mode differentiation if full Jacobians are
  built (`autojacvec=false` and `autodiff=true`). Default is `0` for automatic
  choice of chunk size.
* `autojacvec`: Calculate the Jacobian-vector product via automatic
  differentiation with special seeding.
* `diff_type`: The method used by FiniteDiff.jl for constructing the Jacobian
  if the full Jacobian is required with `autodiff=false`.

Further details:

- If `autodiff=true` and `autojacvec=true`, then the one chunk `J*v` forward-mode
  directional derivative calculation trick is used to compute the product without
  constructing the Jacobian (via ForwardDiff.jl).
- If `autodiff=false` and `autojacvec=true`, then the numerical direction derivative
  trick `(f(x+epsilon*v)-f(x))/epsilon` is used to compute `J*v` without constructing
  the Jacobian.
- If `autodiff=true` and `autojacvec=false`, then the Jacobian is constructed via
  chunked forward-mode automatic differentiation (via ForwardDiff.jl).
- If `autodiff=false` and `autojacvec=false`, then the Jacobian is constructed via
  finite differences via FiniteDiff.jl.

## SciMLProblem Support

This `sensealg` only supports `ODEProblem`s without callbacks (events).
"""
struct ForwardSensitivity{CS,AD,FDT} <: AbstractForwardSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::Bool
  autojacmat::Bool
end
Base.@pure function ForwardSensitivity(;
                                       chunk_size=0,autodiff=true,
                                       diff_type=Val{:central},
                                       autojacvec=autodiff,
                                       autojacmat=false)
  autojacvec && autojacmat && error("Choose either Jacobian matrix products or Jacobian vector products,
                                      autojacmat and autojacvec cannot both be true")
  ForwardSensitivity{chunk_size,autodiff,diff_type}(autojacvec,autojacmat)
end

"""
ForwardDiffSensitivity{CS,CTS} <: AbstractForwardSensitivityAlgorithm{CS,Nothing,Nothing}

An implementation of discrete forward sensitivity analysis through ForwardDiff.jl.
When used within adjoint differentiation (i.e. via Zygote), this will cause forward
differentiation of the `solve` call within the reverse-mode automatic differentiation
environment.

## Constructor

```julia
ForwardDiffSensitivity(;chunk_size=0,convert_tspan=nothing)
```

## Keyword Arguments

* `chunk_size`: the chunk size used by ForwardDiff for computing the Jacobian, i.e. the
  number of simultaneous columns computed.
* `convert_tspan`: whether to convert time to also be `Dual` valued. By default this is
  `nothing` which will only convert if callbacks are found. Conversion is required in order
  to accurately differentiate callbacks (hybrid equations).

## SciMLProblem Support

This `sensealg` supports any `SciMLProblem`s, provided that the solver algorithms is
`SciMLBase.isautodifferentiable`. Note that `ForwardDiffSensitivity` can
accurately differentiate code with callbacks only when `convert_tspan=true`.
"""
struct ForwardDiffSensitivity{CS,CTS} <: AbstractForwardSensitivityAlgorithm{CS,Nothing,Nothing}
end
Base.@pure function ForwardDiffSensitivity(;chunk_size=0,convert_tspan=nothing)
  ForwardDiffSensitivity{chunk_size,convert_tspan}()
end

"""
BacksolveAdjoint{CS,AD,FDT,VJP,NOISE} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}

An implementation of adjoint sensitivity analysis using a backwards solution of the ODE.
By default this algorithm will use the values from the forward pass to perturb the
backwards solution to the correct spot, allowing reduced memory (O(1) memory). Checkpointing
stabilization is included for additional numerical stability over the naive implementation.

## Constructor

```julia
BacksolveAdjoint(;chunk_size=0,autodiff=true,
                  diff_type=Val{:central},
                  autojacvec=autodiff,
                  checkpointing=true, noise=true, noisemixing=false)
```

## Keyword Arguments

* `autodiff`: Use automatic differentiation for constructing the Jacobian
  if the Jacobian needs to be constructed.  Defaults to `true`.
* `chunk_size`: Chunk size for forward-mode differentiation if full Jacobians are
  built (`autojacvec=false` and `autodiff=true`). Default is `0` for automatic
  choice of chunk size.
* `diff_type`: The method used by FiniteDiff.jl for constructing the Jacobian
  if the full Jacobian is required with `autodiff=false`.
* `autojacvec`: Calculate the vector-Jacobian product (`J'*v`) via automatic
  differentiation with special seeding. The default is `true`. The total set
  of choices are:
    - `false`: the Jacobian is constructed via FiniteDiff.jl
    - `true`: the Jacobian is constructed via ForwardDiff.jl
    - `TrackerVJP`: Uses Tracker.jl for the vjp.
    - `ZygoteVJP`: Uses Zygote.jl for the vjp.
    - `EnzymeVJP`: Uses Enzyme.jl for the vjp.
    - `ReverseDiffVJP(compile=false)`: Uses ReverseDiff.jl for the vjp. `compile`
      is a boolean for whether to precompile the tape, which should only be done
      if there are no branches (`if` or `while` statements) in the `f` function.
* `checkpointing`: whether checkpointing is enabled for the reverse pass. Defaults
  to `true`.
* `noise`: Calculate the vector-Jacobian product (`J'*v`) of the diffusion term
  of an SDE via automatic differentiation with special seeding. The default is `true`.
  The total set of choices are:
    - `false`: the Jacobian is constructed via FiniteDiff.jl
    - `true`: the Jacobian is constructed via ForwardDiff.jl
    - `DiffEqSensitivity.ZygoteNoise()`: Uses Zygote.jl for the vjp.
    - `DiffEqSensitivity.ReverseDiffNoise(compile=false)`: Uses ReverseDiff.jl for
      the vjp. `compile` is a boolean for whether to precompile the tape, which
      should only be done if there are no branches (`if` or `while` statements) in
      the `f` function.
* `noisemixing`: Handle noise processes that are not of the form `du[i] = f(u[i])`.
  For example, to compute the sensitivities of an SDE with diagonal diffusion
  ```julia
  function g_mixing!(du,u,p,t)
    du[1] = p[3]*u[1] + p[4]*u[2]
    du[2] = p[3]*u[1] + p[4]*u[2]
    nothing
  end
  ```
  correctly, `noisemixing=true` must be enabled. The default is `false`.

For more details on the vjp choices, please consult the sensitivity algorithms
documentation page or the docstrings of the vjp types.

## Applicability of Backsolve and Caution

When `BacksolveAdjoint` is applicable, it is a fast method and requires the least memory.
However, one must be cautious because not all ODEs are stable under backwards integration
by the majority of ODE solvers. An example of such an equation is the Lorenz equation.
Notice that if one solves the Lorenz equation forward and then in reverse with any
adaptive time step and non-reversible integrator, then the backwards solution diverges
from the forward solution. As a quick demonstration:

```julia
using Sundials
function lorenz(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob,Tsit5(),reltol=1e-12,abstol=1e-12)
prob2 = ODEProblem(lorenz,sol[end],(100.0,0.0))
sol = solve(prob,Tsit5(),reltol=1e-12,abstol=1e-12)
@show sol[end]-u0 #[-3.22091, -1.49394, 21.3435]
```

Thus one should check the stability of the backsolve on their type of problem before
enabling this method. Additionally, using checkpointing with backsolve can be a
low memory way to stabilize it.

For more details on this topic, see
[Stiff Neural Ordinary Differential Equations](https://aip.scitation.org/doi/10.1063/5.0060697).

## Checkpointing

To improve the numerical stability of the reverse pass, `BacksolveAdjoint` includes a checkpointing
feature. If `sol.u` is a time series, then whenever a time `sol.t` is hit while reversing, a callback
will replace the reversing ODE portion with `sol.u[i]`. This nudges the solution back onto the appropriate
trajectory and reduces the numerical caused by drift.

## SciMLProblem Support

This `sensealg` only supports `ODEProblem`s, `SDEProblem`s, and `RODEProblem`s. This `sensealg` supports
callback functions (events).

## References

ODE:
 Rackauckas, C. and Ma, Y. and Martensen, J. and Warner, C. and Zubov, K. and Supekar,
 R. and Skinner, D. and Ramadhana, A. and Edelman, A., Universal Differential Equations
 for Scientific Machine Learning,	arXiv:2001.04385

 Hindmarsh, A. C. and Brown, P. N. and Grant, K. E. and Lee, S. L. and Serban, R.
 and Shumaker, D. E. and Woodward, C. S., SUNDIALS: Suite of nonlinear and
 differential/algebraic equation solvers, ACM Transactions on Mathematical
 Software (TOMS), 31, pp:363–396 (2005)

 Chen, R.T.Q. and Rubanova, Y. and Bettencourt, J. and Duvenaud, D. K.,
 Neural ordinary differential equations. In Advances in neural information processing
 systems, pp. 6571–6583 (2018)

 Pontryagin, L. S. and Mishchenko, E.F. and Boltyanskii, V.G. and Gamkrelidze, R.V.
 The mathematical theory of optimal processes. Routledge, (1962)

 Rackauckas, C. and Ma, Y. and Dixit, V. and Guo, X. and Innes, M. and Revels, J.
 and Nyberg, J. and Ivaturi, V., A comparison of automatic differentiation and
 continuous sensitivity analysis for derivatives of differential equation solutions,
 arXiv:1812.01892

DAE:
 Cao, Y. and Li, S. and Petzold, L. and Serban, R., Adjoint sensitivity analysis
 for differential-algebraic equations: The adjoint DAE system and its numerical
 solution, SIAM journal on scientific computing 24 pp: 1076-1089 (2003)

SDE:
 Gobet, E. and Munos, R., Sensitivity Analysis Using Ito-Malliavin Calculus and
 Martingales, and Application to Stochastic Optimal Control,
 SIAM Journal on control and optimization, 43, pp. 1676-1713 (2005)

 Li, X. and Wong, T.-K. L.and Chen, R. T. Q. and Duvenaud, D.,
 Scalable Gradients for Stochastic Differential Equations,
 PMLR 108, pp. 3870-3882 (2020), http://proceedings.mlr.press/v108/li20i.html
"""
struct BacksolveAdjoint{CS,AD,FDT,VJP,NOISE} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::VJP
  checkpointing::Bool
  noise::NOISE
  noisemixing::Bool
end
Base.@pure function BacksolveAdjoint(;chunk_size=0,autodiff=true,
                                      diff_type=Val{:central},
                                      autojacvec=autodiff,
                                      checkpointing=true, noise=true,noisemixing=false)
  BacksolveAdjoint{chunk_size,autodiff,diff_type,typeof(autojacvec),typeof(noise)}(autojacvec,checkpointing,noise,noisemixing)
end

"""
InterpolatingAdjoint{CS,AD,FDT,VJP,NOISE} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}

An implementation of adjoint sensitivity analysis which uses the interpolation of
the forward solution for the reverse solve vector-Jacobian products. By
default it requires a dense solution of the forward pass and will internally
ignore saving arguments during the gradient calculation. When checkpointing is
enabled it will only require the memory to interpolate between checkpoints.

## Constructor

```julia
function InterpolatingAdjoint(;chunk_size=0,autodiff=true,
                               diff_type=Val{:central},
                               autojacvec=autodiff,
                               checkpointing=false, noise=true, noisemixing=false)
```

## Keyword Arguments

* `autodiff`: Use automatic differentiation for constructing the Jacobian
  if the Jacobian needs to be constructed.  Defaults to `true`.
* `chunk_size`: Chunk size for forward-mode differentiation if full Jacobians are
  built (`autojacvec=false` and `autodiff=true`). Default is `0` for automatic
  choice of chunk size.
* `diff_type`: The method used by FiniteDiff.jl for constructing the Jacobian
  if the full Jacobian is required with `autodiff=false`.
* `autojacvec`: Calculate the vector-Jacobian product (`J'*v`) via automatic
  differentiation with special seeding. The default is `true`. The total set
  of choices are:
    - `false`: the Jacobian is constructed via FiniteDiff.jl
    - `true`: the Jacobian is constructed via ForwardDiff.jl
    - `TrackerVJP`: Uses Tracker.jl for the vjp.
    - `ZygoteVJP`: Uses Zygote.jl for the vjp.
    - `EnzymeVJP`: Uses Enzyme.jl for the vjp.
    - `ReverseDiffVJP(compile=false)`: Uses ReverseDiff.jl for the vjp. `compile`
      is a boolean for whether to precompile the tape, which should only be done
      if there are no branches (`if` or `while` statements) in the `f` function.
* `checkpointing`: whether checkpointing is enabled for the reverse pass. Defaults
  to `true`.
* `noise`: Calculate the vector-Jacobian product (`J'*v`) of the diffusion term
  of an SDE via automatic differentiation with special seeding. The default is `true`.
  The total set of choices are:
    - `false`: the Jacobian is constructed via FiniteDiff.jl
    - `true`: the Jacobian is constructed via ForwardDiff.jl
    - `DiffEqSensitivity.ZygoteNoise()`: Uses Zygote.jl for the vjp.
    - `DiffEqSensitivity.ReverseDiffNoise(compile=false)`: Uses ReverseDiff.jl for
      the vjp. `compile` is a boolean for whether to precompile the tape, which
      should only be done if there are no branches (`if` or `while` statements) in
      the `f` function.
* `noisemixing`: Handle noise processes that are not of the form `du[i] = f(u[i])`.
  For example, to compute the sensitivities of an SDE with diagonal diffusion
  ```julia
  function g_mixing!(du,u,p,t)
    du[1] = p[3]*u[1] + p[4]*u[2]
    du[2] = p[3]*u[1] + p[4]*u[2]
    nothing
  end
  ```
  correctly, `noisemixing=true` must be enabled. The default is `false`.

For more details on the vjp choices, please consult the sensitivity algorithms
documentation page or the docstrings of the vjp types.

## Checkpointing

To reduce the memory usage of the reverse pass, `InterpolatingAdjoint` includes a checkpointing
feature. If `sol` is `dense`, checkpointing is ignored and the continuous solution is used for
calculating `u(t)` at arbitrary time points. If `checkpointing=true` and `sol` is not `dense`,
then dense intervals between `sol.t[i]` and `sol.t[i+1]` are reconstructed on-demand for calculating
`u(t)` at arbitrary time points. This reduces the total memory requirement to only the cost of
holding the dense solution over the largest time interval (in terms of number of required steps).
The total compute cost is no more than double the original forward compute cost.

## SciMLProblem Support

This `sensealg` only supports `ODEProblem`s, `SDEProblem`s, and `RODEProblem`s. This `sensealg`
supports callbacks (events).

## References

 Rackauckas, C. and Ma, Y. and Martensen, J. and Warner, C. and Zubov, K. and Supekar,
 R. and Skinner, D. and Ramadhana, A. and Edelman, A., Universal Differential Equations
 for Scientific Machine Learning,	arXiv:2001.04385

 Hindmarsh, A. C. and Brown, P. N. and Grant, K. E. and Lee, S. L. and Serban, R.
 and Shumaker, D. E. and Woodward, C. S., SUNDIALS: Suite of nonlinear and
 differential/algebraic equation solvers, ACM Transactions on Mathematical
 Software (TOMS), 31, pp:363–396 (2005)

 Rackauckas, C. and Ma, Y. and Dixit, V. and Guo, X. and Innes, M. and Revels, J.
 and Nyberg, J. and Ivaturi, V., A comparison of automatic differentiation and
 continuous sensitivity analysis for derivatives of differential equation solutions,
 arXiv:1812.01892
"""
struct InterpolatingAdjoint{CS,AD,FDT,VJP,NOISE} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::VJP
  checkpointing::Bool
  noise::NOISE
  noisemixing::Bool
end
Base.@pure function InterpolatingAdjoint(;chunk_size=0,autodiff=true,
                                         diff_type=Val{:central},
                                         autojacvec=autodiff,
                                         checkpointing=false, noise=true,noisemixing=false)
  InterpolatingAdjoint{chunk_size,autodiff,diff_type,typeof(autojacvec),typeof(noise)}(autojacvec,checkpointing,noise,noisemixing)
end

"""
QuadratureAdjoint{CS,AD,FDT,VJP} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}

An implementation of adjoint sensitivity analysis which develops a full
continuous solution of the reverse solve in order to perform a post-ODE
quadrature. This method requires the the dense solution and will ignore
saving arguments during the gradient calculation. The tolerances in the
constructor control the inner quadrature. The inner quadrature uses a
ReverseDiff vjp if autojacvec, and `compile=false` by default but can
compile the tape under the same circumstances as `ReverseDiffVJP`.

This method is O(n^3 + p) for stiff / implicit equations (as opposed to the
O((n+p)^3) scaling of BacksolveAdjoint and InterpolatingAdjoint), and thus
is much more compute efficient. However, it requires holding a dense reverse
pass and is thus memory intensive.

## Constructor

```julia
function QuadratureAdjoint(;chunk_size=0,autodiff=true,
                            diff_type=Val{:central},
                            autojacvec=autodiff,abstol=1e-6,
                            reltol=1e-3,compile=false)
```

## Keyword Arguments

* `autodiff`: Use automatic differentiation for constructing the Jacobian
  if the Jacobian needs to be constructed.  Defaults to `true`.
* `chunk_size`: Chunk size for forward-mode differentiation if full Jacobians are
  built (`autojacvec=false` and `autodiff=true`). Default is `0` for automatic
  choice of chunk size.
* `diff_type`: The method used by FiniteDiff.jl for constructing the Jacobian
  if the full Jacobian is required with `autodiff=false`.
* `autojacvec`: Calculate the vector-Jacobian product (`J'*v`) via automatic
  differentiation with special seeding. The default is `true`. The total set
  of choices are:
    - `false`: the Jacobian is constructed via FiniteDiff.jl
    - `true`: the Jacobian is constructed via ForwardDiff.jl
    - `TrackerVJP`: Uses Tracker.jl for the vjp.
    - `ZygoteVJP`: Uses Zygote.jl for the vjp.
    - `EnzymeVJP`: Uses Enzyme.jl for the vjp.
    - `ReverseDiffVJP(compile=false)`: Uses ReverseDiff.jl for the vjp. `compile`
      is a boolean for whether to precompile the tape, which should only be done
      if there are no branches (`if` or `while` statements) in the `f` function.
* `abstol`: absolute tolerance for the quadrature calculation
* `reltol`: relative tolerance for the quadrature calculation
* `compile`: whether to compile the vjp calculation for the integrand calculation.
  See `ReverseDiffVJP` for more details.

For more details on the vjp choices, please consult the sensitivity algorithms
documentation page or the docstrings of the vjp types.

## SciMLProblem Support

This `sensealg` only supports `ODEProblem`s. This `sensealg` supports events (callbacks).

## References

 Rackauckas, C. and Ma, Y. and Martensen, J. and Warner, C. and Zubov, K. and Supekar,
 R. and Skinner, D. and Ramadhana, A. and Edelman, A., Universal Differential Equations
 for Scientific Machine Learning,	arXiv:2001.04385

 Hindmarsh, A. C. and Brown, P. N. and Grant, K. E. and Lee, S. L. and Serban, R.
 and Shumaker, D. E. and Woodward, C. S., SUNDIALS: Suite of nonlinear and
 differential/algebraic equation solvers, ACM Transactions on Mathematical
 Software (TOMS), 31, pp:363–396 (2005)

 Rackauckas, C. and Ma, Y. and Dixit, V. and Guo, X. and Innes, M. and Revels, J.
 and Nyberg, J. and Ivaturi, V., A comparison of automatic differentiation and
 continuous sensitivity analysis for derivatives of differential equation solutions,
 arXiv:1812.01892

 Kim, S., Ji, W., Deng, S., Ma, Y., & Rackauckas, C. (2021). Stiff neural ordinary
 differential equations. Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(9), 093122.
"""
struct QuadratureAdjoint{CS,AD,FDT,VJP} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::VJP
  abstol::Float64
  reltol::Float64
  compile::Bool
end
Base.@pure function QuadratureAdjoint(;chunk_size=0,autodiff=true,
                                         diff_type=Val{:central},
                                         autojacvec=autodiff,abstol=1e-6,
                                         reltol=1e-3,compile=false)
  QuadratureAdjoint{chunk_size,autodiff,diff_type,typeof(autojacvec)}(autojacvec,abstol,reltol,compile)
end

"""
TrackerAdjoint <: AbstractAdjointSensitivityAlgorithm{nothing,true,nothing}

An implementation of discrete adjoint sensitivity analysis
using the Tracker.jl tracing-based AD. Supports in-place functions through
an Array of Structs formulation, and supports out of place through struct of
arrays.

## Constructor

```julia
TrackerAdjoint()
```

## SciMLProblem Support

This `sensealg` supports any `DEProblem` if the algorithm is `SciMLBase.isautodifferentiable`
Compatible with a limited subset of `AbstractArray` types for `u0`, including `CuArrays`.

!!! warn

    TrackerAdjoint is incompatible with Stiff ODE solvers using forward-mode automatic
    differentiation for the Jacobians. Thus for example, `TRBDF2()` will error. Instead,
    use `autodiff=false`, i.e. `TRBDF2(autodiff=false)`. This will only remove the
    forward-mode automatic differentiation of the Jacobian construction, not the reverse-mode
    AD usage, and thus performance will still be nearly the same, though Jacobian accuracy
    may suffer which could cause more steps to be required.
"""
struct TrackerAdjoint <: AbstractAdjointSensitivityAlgorithm{nothing,true,nothing} end

"""
ReverseDiffAdjoint <: AbstractAdjointSensitivityAlgorithm{nothing,true,nothing}

An implementation of discrete adjoint sensitivity analysis using the ReverseDiff.jl
tracing-based AD. Supports in-place functions through an Array of Structs formulation,
and supports out of place through struct of arrays.

## Constructor

```julia
ReverseDiffAdjoint()
```

## SciMLProblem Support

This `sensealg` supports any `DEProblem` if the algorithm is `SciMLBase.isautodifferentiable`.
Requires that the state variables are CPU-based `Array` types.
"""
struct ReverseDiffAdjoint <: AbstractAdjointSensitivityAlgorithm{nothing,true,nothing} end

"""
ZygoteAdjoint <: AbstractAdjointSensitivityAlgorithm{nothing,true,nothing}

An implementation of discrete adjoint sensitivity analysis
using the Zygote.jl source-to-source AD directly on the differential equation
solver.

## Constructor

```julia
ZygoteAdjoint()
```

## SciMLProblem Support

Currently fails on almost every solver.
"""
struct ZygoteAdjoint <: AbstractAdjointSensitivityAlgorithm{nothing,true,nothing} end

"""
ForwardLSS{CS,AD,FDT,RType,gType} <: AbstractShadowingSensitivityAlgorithm{CS,AD,FDT}

An implementation of the discrete, forward-mode
[least squares shadowing](https://arxiv.org/abs/1204.0159) (LSS) method. LSS replaces
the ill-conditioned initial value probem (`ODEProblem`) for chaotic systems by a
well-conditioned least-squares problem. This allows for computing sensitivities of
long-time averaged quantities with respect to the parameters of the `ODEProblem`. The
computational cost of LSS scales as (number of states x number of time steps). Converges
to the correct sensitivity at a rate of `T^(-1/2)`, where `T` is the time of the trajectory.
See `NILSS()` and `NILSAS()` for a more efficient non-intrusive formulation.

## Constructor

```julia
ForwardLSS(;
          chunk_size=0,autodiff=true,
          diff_type=Val{:central},
          LSSregularizer=TimeDilation(10.0,0.0,0.0),
          g=nothing)
```

## Keyword Arguments

* `autodiff`: Use automatic differentiation for constructing the Jacobian
  if the Jacobian needs to be constructed.  Defaults to `true`.
* `chunk_size`: Chunk size for forward-mode differentiation if full Jacobians are
  built (`autojacvec=false` and `autodiff=true`). Default is `0` for automatic
  choice of chunk size.
* `diff_type`: The method used by FiniteDiff.jl for constructing the Jacobian
  if the full Jacobian is required with `autodiff=false`.
* `LSSregularizer`: Using `LSSregularizer`, one can choose between three different
  regularization routines. The default choice is `TimeDilation(10.0,0.0,0.0)`.
    - `CosWindowing()`: cos windowing of the time grid, i.e. the time grid (saved
      time steps) is transformed using a cosine.
    - `Cos2Windowing()`: cos^2 windowing of the time grid.
    - `TimeDilation(alpha::Number,t0skip::Number,t1skip::Number)`: Corresponds to
      a time dilation. `alpha` controls the weight. `t0skip` and `t1skip` indicate
      the times truncated at the beginnning and end of the trajectory, respectively.
* `g`: instantaneous objective function of the long-time averaged objective.

## SciMLProblem Support

This `sensealg` only supports `ODEProblem`s. This `sensealg` does not support
events (callbacks). This `sensealg` assumes that the objective is a long-time averaged
quantity and ergodic, i.e. the time evolution of the system behaves qualitatively the
same over infinite time independent of the specified initial conditions, such that only
the sensitivity with respect to the parameters is of interest.

## References

Wang, Q., Hu, R., and Blonigan, P. Least squares shadowing sensitivity analysis of
chaotic limit cycle oscillations. Journal of Computational Physics, 267, 210-224 (2014).

Wang, Q., Convergence of the Least Squares Shadowing Method for Computing Derivative of Ergodic
Averages, SIAM Journal on Numerical Analysis, 52, 156–170 (2014).

Blonigan, P., Gomez, S., Wang, Q., Least Squares Shadowing for sensitivity analysis of turbulent
fluid flows, in: 52nd Aerospace Sciences Meeting, 1–24 (2014).
"""
struct ForwardLSS{CS,AD,FDT,RType,gType} <: AbstractShadowingSensitivityAlgorithm{CS,AD,FDT}
  LSSregularizer::RType
  g::gType
end
Base.@pure function ForwardLSS(;
  chunk_size=0, autodiff=true,
  diff_type=Val{:central},
  LSSregularizer=TimeDilation(10.0,0.0,0.0),
  g=nothing)

  ForwardLSS{chunk_size,autodiff,diff_type,typeof(LSSregularizer),typeof(g)}(LSSregularizer, g)
end

"""
AdjointLSS{CS,AD,FDT,RType,gType} <: AbstractShadowingSensitivityAlgorithm{CS,AD,FDT}

An implementation of the discrete, adjoint-mode
[least square shadowing](https://arxiv.org/abs/1204.0159) method. LSS replaces
the ill-conditioned initial value probem (`ODEProblem`) for chaotic systems by a
well-conditioned least-squares problem. This allows for computing sensitivities of
long-time averaged quantities with respect to the parameters of the `ODEProblem`. The
computational cost of LSS scales as (number of states x number of time steps). Converges
to the correct sensitivity at a rate of `T^(-1/2)`, where `T` is the time of the trajectory.
See `NILSS()` and `NILSAS()` for a more efficient non-intrusive formulation.

## Constructor

```julia
AdjointLSS(;
          chunk_size=0,autodiff=true,
          diff_type=Val{:central},
          LSSRegularizer=CosWindowing(),
          g=nothing)
```

## Keyword Arguments

* `autodiff`: Use automatic differentiation for constructing the Jacobian
  if the Jacobian needs to be constructed.  Defaults to `true`.
* `chunk_size`: Chunk size for forward-mode differentiation if full Jacobians are
  built (`autojacvec=false` and `autodiff=true`). Default is `0` for automatic
  choice of chunk size.
* `diff_type`: The method used by FiniteDiff.jl for constructing the Jacobian
  if the full Jacobian is required with `autodiff=false`.
* `LSSregularizer`: Using `LSSregularizer`, one can choose between different
  regularization routines. The default choice is `TimeDilation(10.0,0.0,0.0)`.
    - `TimeDilation(alpha::Number,t0skip::Number,t1skip::Number)`: Corresponds to
      a time dilation. `alpha` controls the weight. `t0skip` and `t1skip` indicate
      the times truncated at the beginnning and end of the trajectory, respectively.
      The default value for `t0skip` and `t1skip` is `zero(alpha)`.
* `g`: instantaneous objective function of the long-time averaged objective.

## SciMLProblem Support

This `sensealg` only supports `ODEProblem`s. This `sensealg` does not support
events (callbacks). This `sensealg` assumes that the objective is a long-time averaged
quantity and ergodic, i.e. the time evolution of the system behaves qualitatively the
same over infinite time independent of the specified initial conditions, such that only
the sensitivity with respect to the parameters is of interest.

## References

Wang, Q., Hu, R., and Blonigan, P. Least squares shadowing sensitivity analysis of
chaotic limit cycle oscillations. Journal of Computational Physics, 267, 210-224 (2014).
"""
struct AdjointLSS{CS,AD,FDT,RType,gType} <: AbstractShadowingSensitivityAlgorithm{CS,AD,FDT}
  LSSregularizer::RType
  g::gType
end
Base.@pure function AdjointLSS(;
  chunk_size=0, autodiff=true,
  diff_type=Val{:central},
  LSSregularizer=TimeDilation(10.0, 0.0, 0.0),
  g=nothing)
  AdjointLSS{chunk_size,autodiff,diff_type,typeof(LSSregularizer),typeof(g)}(LSSregularizer, g)
end

abstract type AbstractLSSregularizer end
abstract type AbstractCosWindowing <: AbstractLSSregularizer end
struct CosWindowing <: AbstractCosWindowing end
struct Cos2Windowing <: AbstractCosWindowing end

"""
TimeDilation{T1<:Number} <: AbstractLSSregularizer

A regularization method for `LSS`. See `?LSS` for
additional information and other methods.

## Constructor

```julia
TimeDilation(alpha;
          t0skip=zero(alpha),
          t1skip=zero(alpha))
```
"""
struct TimeDilation{T1<:Number} <: AbstractLSSregularizer
  alpha::T1 # alpha: weight of the time dilation term in LSS.
  t0skip::T1
  t1skip::T1
end
function TimeDilation(alpha,t0skip=zero(alpha),t1skip=zero(alpha))
  TimeDilation{typeof(alpha)}(alpha,t0skip,t1skip)
end
"""
struct NILSS{CS,AD,FDT,RNG,nType,gType} <: AbstractShadowingSensitivityAlgorithm{CS,AD,FDT}

An implementation of the forward-mode, continuous
[non-intrusive least squares shadowing](https://arxiv.org/abs/1611.00880) method. `NILSS`
allows for computing sensitivities of long-time averaged quantities with respect to the
parameters of an `ODEProblem` by constraining the computation to the unstable subspace.
`NILSS` employs the continuous-time `ForwardSensitivity` method as tangent solver. To
avoid an exponential blow-up of the (homogenous and inhomogenous) tangent solutions,
the trajectory should be divided into sufficiently small segments, where the tangent solutions
are rescaled on the interfaces. The computational and memory cost of NILSS scale with
the number of unstable (positive) Lyapunov exponents (instead of the number of states as
in the LSS method). `NILSS` avoids the explicit construction of the Jacobian at each time
step and thus should generally be preferred (for large system sizes) over `ForwardLSS`.

## Constructor

```julia
NILSS(nseg, nstep; nus = nothing,
                   rng = Xorshifts.Xoroshiro128Plus(rand(UInt64)),
                   chunk_size=0,autodiff=true,
                   diff_type=Val{:central},
                   autojacvec=autodiff,
                   g=nothing)
```

## Arguments

* `nseg`: Number of segments on full time interval on the attractor.
* `nstep`: number of steps on each segment.

## Keyword Arguments

* `nus`: Dimension of the unstable subspace. Default is `nothing`. `nus` must be
  smaller or equal to the state dimension (`length(u0)`). With the default choice,
  `nus = length(u0) - 1` will be set at compile time.
* `rng`: (Pseudo) random number generator. Used for initializing the homogenous
  tangent states (`w`). Default is `Xorshifts.Xoroshiro128Plus(rand(UInt64))`.
* `autodiff`: Use automatic differentiation in the internal sensitivity algorithm
  computations. Default is `true`.
* `chunk_size`: Chunk size for forward mode differentiation if full Jacobians are
  built (`autojacvec=false` and `autodiff=true`). Default is `0` for automatic
  choice of chunk size.
* `autojacvec`: Calculate the Jacobian-vector product via automatic
  differentiation with special seeding.
* `diff_type`: The method used by FiniteDiff.jl for constructing the Jacobian
  if the full Jacobian is required with `autodiff=false`.
* `g`: instantaneous objective function of the long-time averaged objective.

## SciMLProblem Support

This `sensealg` only supports `ODEProblem`s. This `sensealg` does not support
events (callbacks). This `sensealg` assumes that the objective is a long-time averaged
quantity and ergodic, i.e. the time evolution of the system behaves qualitatively the
same over infinite time independent of the specified initial conditions, such that only
the sensitivity with respect to the parameters is of interest.

## References
Ni, A., Blonigan, P. J., Chater, M., Wang, Q., Zhang, Z., Sensitivity analy-
sis on chaotic dynamical system by Non-Intrusive Least Square Shadowing
(NI-LSS), in: 46th AIAA Fluid Dynamics Conference, AIAA AVIATION Forum (AIAA 2016-4399),
American Institute of Aeronautics and Astronautics, 1–16 (2016).

Ni, A., and Wang, Q. Sensitivity analysis on chaotic dynamical systems by Non-Intrusive
Least Squares Shadowing (NILSS). Journal of Computational Physics 347, 56-77 (2017).
"""
struct NILSS{CS,AD,FDT,RNG,nType,gType} <: AbstractShadowingSensitivityAlgorithm{CS,AD,FDT}
  rng::RNG
  nseg::Int
  nstep::Int
  nus::nType
  autojacvec::Bool
  g::gType
end
Base.@pure function NILSS(nseg, nstep; nus=nothing, rng=Xorshifts.Xoroshiro128Plus(rand(UInt64)),
  chunk_size=0, autodiff=true,
  diff_type=Val{:central},
  autojacvec=autodiff,
  g=nothing
)
  NILSS{chunk_size,autodiff,diff_type,typeof(rng),typeof(nus),typeof(g)}(rng, nseg, nstep, nus, autojacvec, g)
end

"""
NILSAS{CS,AD,FDT,RNG,SENSE,gType} <: AbstractShadowingSensitivityAlgorithm{CS,AD,FDT}

An implementation of the adjoint-mode, continuous
[non-intrusive adjoint least squares shadowing](https://arxiv.org/abs/1801.08674) method.
`NILSAS` allows for computing sensitivities of long-time averaged quantities with respect
to the parameters of an `ODEProblem` by constraining the computation to the unstable subspace.
`NILSAS` employs SciMLSensitivity.jl's continuous adjoint sensitivity methods on each segment
to compute (homogenous and inhomogenous) adjoint solutions. To avoid an exponential blow-up
of the adjoint solutions, the trajectory should be divided into sufficiently small segments,
where the adjoint solutions are rescaled on the interfaces. The computational and memory cost
of NILSAS scale with the number of unstable, adjoint Lyapunov exponents (instead of the number
of states as in the LSS method). `NILSAS` avoids the explicit construction of the Jacobian at
each time step and thus should generally be preferred (for large system sizes) over `AdjointLSS`.
`NILSAS` is favourable over `NILSS` for many parameters because NILSAS computes the gradient
with respect to multiple parameters with negligible additional cost.

## Constructor

```julia
NILSAS(nseg, nstep, M=nothing; rng = Xorshifts.Xoroshiro128Plus(rand(UInt64)),
                                adjoint_sensealg = BacksolveAdjoint(),
                                chunk_size=0,autodiff=true,
                                diff_type=Val{:central},
                                autojacvec=autodiff,
                                g=nothing
                                )
```

## Arguments

* `nseg`: Number of segments on full time interval on the attractor.
* `nstep`: number of steps on each segment.
* `M`: number of homogenous adjoint solutions. This number must be bigger or equal
  than the number of (positive, adjoint) Lyapunov exponents. Default is `nothing`.

## Keyword Arguments

* `rng`: (Pseudo) random number generator. Used for initializing the terminate
  conditions of the homogenous adjoint states (`w`). Default is `Xorshifts.Xoroshiro128Plus(rand(UInt64))`.
* `adjoint_sensealg`: Continuous adjoint sensitivity method to compute homogenous
  and inhomogenous adjoint solutions on each segment. Default is `BacksolveAdjoint()`.
* `autodiff`: Use automatic differentiation for constructing the Jacobian
  if the Jacobian needs to be constructed.  Defaults to `true`.
* `chunk_size`: Chunk size for forward-mode differentiation if full Jacobians are
  built (`autojacvec=false` and `autodiff=true`). Default is `0` for automatic
  choice of chunk size.
* `diff_type`: The method used by FiniteDiff.jl for constructing the Jacobian
  if the full Jacobian is required with `autodiff=false`.
* `autojacvec`: Calculate the vector-Jacobian product (`J'*v`) via automatic
  differentiation with special seeding. The default is `true`. The total set
  of choices are:
    - `false`: the Jacobian is constructed via FiniteDiff.jl
    - `true`: the Jacobian is constructed via ForwardDiff.jl
    - `TrackerVJP`: Uses Tracker.jl for the vjp.
    - `ZygoteVJP`: Uses Zygote.jl for the vjp.
    - `EnzymeVJP`: Uses Enzyme.jl for the vjp.
    - `ReverseDiffVJP(compile=false)`: Uses ReverseDiff.jl for the vjp. `compile`
      is a boolean for whether to precompile the tape, which should only be done
      if there are no branches (`if` or `while` statements) in the `f` function.
* `g`: instantaneous objective function of the long-time averaged objective.

## SciMLProblem Support

This `sensealg` only supports `ODEProblem`s. This `sensealg` does not support
events (callbacks). This `sensealg` assumes that the objective is a long-time averaged
quantity and ergodic, i.e. the time evolution of the system behaves qualitatively the
same over infinite time independent of the specified initial conditions, such that only
the sensitivity with respect to the parameters is of interest.

## References

Ni, A., and Talnikar, C., Adjoint sensitivity analysis on chaotic dynamical systems
by Non-Intrusive Least Squares Adjoint Shadowing (NILSAS). Journal of Computational
Physics 395, 690-709 (2019).
"""
struct NILSAS{CS,AD,FDT,RNG,SENSE,gType} <: AbstractShadowingSensitivityAlgorithm{CS,AD,FDT}
  rng::RNG
  adjoint_sensealg::SENSE
  M::Int
  nseg::Int
  nstep::Int
  autojacvec::Bool
  g::gType
end
Base.@pure function NILSAS(nseg, nstep, M=nothing; rng=Xorshifts.Xoroshiro128Plus(rand(UInt64)),
  adjoint_sensealg=BacksolveAdjoint(),
  chunk_size=0, autodiff=true,
  diff_type=Val{:central},
  autojacvec=autodiff,
  g=nothing
)
  # integer dimension of the unstable subspace
  M === nothing && error("Please provide an `M` with `M >= nus + 1`, where nus is the number of unstable covariant Lyapunov vectors.")

  NILSAS{chunk_size,autodiff,diff_type,typeof(rng),typeof(adjoint_sensealg),typeof(g)}(rng, adjoint_sensealg, M,
    nseg, nstep, autojacvec, g)
end

"""
SteadyStateAdjoint{CS,AD,FDT,VJP,LS} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}

An implementation of the adjoint differentiation of a nonlinear solve. Uses the
implicit function theorem to directly compute the derivative of the solution to
``f(u,p) = 0`` with respect to `p`.

## Constructor

```julia
SteadyStateAdjoint(;chunk_size = 0, autodiff = true,
                    diff_type = Val{:central},
                    autojacvec = autodiff, linsolve = nothing)
```

## Keyword Arguments

## Keyword Arguments

* `autodiff`: Use automatic differentiation for constructing the Jacobian
  if the Jacobian needs to be constructed.  Defaults to `true`.
* `chunk_size`: Chunk size for forward-mode differentiation if full Jacobians are
  built (`autojacvec=false` and `autodiff=true`). Default is `0` for automatic
  choice of chunk size.
* `diff_type`: The method used by FiniteDiff.jl for constructing the Jacobian
  if the full Jacobian is required with `autodiff=false`.
* `autojacvec`: Calculate the vector-Jacobian product (`J'*v`) via automatic
  differentiation with special seeding. The default is `true`. The total set
  of choices are:
    - `false`: the Jacobian is constructed via FiniteDiff.jl
    - `true`: the Jacobian is constructed via ForwardDiff.jl
    - `TrackerVJP`: Uses Tracker.jl for the vjp.
    - `ZygoteVJP`: Uses Zygote.jl for the vjp.
    - `EnzymeVJP`: Uses Enzyme.jl for the vjp.
    - `ReverseDiffVJP(compile=false)`: Uses ReverseDiff.jl for the vjp. `compile`
      is a boolean for whether to precompile the tape, which should only be done
      if there are no branches (`if` or `while` statements) in the `f` function.
* `linsolve`: the linear solver used in the adjoint solve. Defaults to `nothing`,
  which uses a polyalgorithm to attempt to automatically choose an efficient
  algorithm.

For more details on the vjp choices, please consult the sensitivity algorithms
documentation page or the docstrings of the vjp types.

## References

Johnson, S. G., Notes on Adjoint Methods for 18.336, Online at
http://math.mit.edu/stevenj/18.336/adjoint.pdf (2007)
"""
struct SteadyStateAdjoint{CS,AD,FDT,VJP,LS} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::VJP
  linsolve::LS
end

Base.@pure function SteadyStateAdjoint(;chunk_size = 0, autodiff = true, diff_type = Val{:central},
                                        autojacvec = autodiff, linsolve = nothing)
  SteadyStateAdjoint{chunk_size,autodiff,diff_type,typeof(autojacvec),typeof(linsolve)}(autojacvec,linsolve)
end

abstract type VJPChoice end

"""
ZygoteVJP <: VJPChoice

Uses Zygote.jl to compute vector-Jacobian products. Tends to be the fastest VJP method if the
ODE/DAE/SDE/DDE is written with mostly vectorized  functions (like neural networks and other
layers from Flux.jl) and the `f` functions is given out-of-place. If the `f` function is
in-place, then `Zygote.Buffer` arrays are used internally which can greatly reduce the
performance of the VJP method.

## Constructor

```julia
ZygoteVJP(compile=false)
```
"""
struct ZygoteVJP <: VJPChoice end

"""
EnzymeVJP <: VJPChoice

Uses Enzyme.jl to compute vector-Jacobian products. Is the fastest VJP whenever applicable,
though Enzyme.jl currently has low coverage over the Julia programming language, for example
restricting the user's defined `f` function to not do things like require garbage collection
or calls to BLAS/LAPACK. However, mutation is supported, meaning that in-place `f` with
fully mutating non-allocating code will work with Enzyme (provided no high level calls to C
like BLAS/LAPACK are used) and this will be the most efficient adjoint implementation.

## Constructor

```julia
EnzymeVJP(compile=false)
```
"""
struct EnzymeVJP <: VJPChoice end

"""
TrackerVJP <: VJPChoice

Uses Tracker.jl to compute the vector-Jacobian products. If `f` is in-place,
then it uses a array of structs formulation to do scalarized reverse mode,
while if `f` is out-of-place then it uses an array-based reverse mode.

Not as efficient as `ReverseDiffVJP`, but supports GPUs when doing array-based
reverse mode.

## Constructor

```julia
TrackerVJP(compile=false)
```
"""
struct TrackerVJP <: VJPChoice end

"""
ReverseDiffVJP{compile} <: VJPChoice

Uses ReverseDiff.jl to compute the vector-Jacobian products. If `f` is in-place,
then it uses a array of structs formulation to do scalarized reverse mode,
while if `f` is out-of-place then it uses an array-based reverse mode.

Usually the fastest when scalarized operations exist in the f function
(like in scientific machine learning applications like Universal Differential Equations)
and the boolean compilation is enabled (i.e. ReverseDiffVJP(true)), if EnzymeVJP fails on
a given choice of `f`.

Does not support GPUs (CuArrays).

## Constructor

```julia
ReverseDiffVJP(compile=false)
```

## Keyword Arguments

* `compile`: Whether to cache the compilation of the reverse tape. This heavily increases
  the performance of the method but requires that the `f` function of the ODE/DAE/SDE/DDE
  has no branching.
"""
struct ReverseDiffVJP{compile} <: VJPChoice
  ReverseDiffVJP(compile=false) = new{compile}()
end

abstract type NoiseChoice end

"""
ZygoteNoise <: NoiseChoice

Uses Zygote.jl to compute vector-Jacobian products for the noise term (for SDE adjoints only).
Tends to be the fastest VJP method if the ODE/DAE/SDE/DDE is written with mostly vectorized
functions (like neural networks and other layers from Flux.jl) and the `f` functions is given
out-of-place. If the `f` function is in-place, then `Zygote.Buffer` arrays are used
internally which can greatly reduce the performance of the VJP method.

## Constructor

```julia
ZygoteNoise()
```
"""
struct ZygoteNoise <: NoiseChoice end

"""
ReverseDiffNoise{compile} <: NoiseChoice

Uses ReverseDiff.jl to compute the vector-Jacobian products for the noise
term differentiation (for SDE adjoints only). If `f` is in-place,
then it uses a array of structs formulation to do scalarized reverse mode,
while if `f` is out-of-place then it uses an array-based reverse mode.

Usually the fastest when scalarized operations exist in the f function
(like in scientific machine learning applications like Universal Differential Equations)
and the boolean compilation is enabled (i.e. ReverseDiffVJP(true)), if EnzymeVJP fails on
a given choice of `f`.

Does not support GPUs (CuArrays).

## Constructor

```julia
ReverseDiffNoise(compile=false)
```

## Keyword Arguments

* `compile`: Whether to cache the compilation of the reverse tape. This heavily increases
  the performance of the method but requires that the `f` function of the ODE/DAE/SDE/DDE
  has no branching.
"""
struct ReverseDiffNoise{compile} <: NoiseChoice
  ReverseDiffNoise(compile=false) = new{compile}()
end

@inline convert_tspan(::ForwardDiffSensitivity{CS,CTS}) where {CS,CTS} = CTS
@inline convert_tspan(::Any) = nothing
@inline alg_autodiff(alg::DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT}) where {CS,AD,FDT} = AD
@inline get_chunksize(alg::DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT}) where {CS,AD,FDT} = CS
@inline diff_type(alg::DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT}) where {CS,AD,FDT} = FDT
@inline function get_jacvec(alg::DiffEqBase.AbstractSensitivityAlgorithm)
  alg.autojacvec isa Bool ? alg.autojacvec : true
end
@inline function get_jacmat(alg::DiffEqBase.AbstractSensitivityAlgorithm)
  alg.autojacmat isa Bool ? alg.autojacmat : true
end
@inline ischeckpointing(alg::DiffEqBase.AbstractSensitivityAlgorithm, sol=nothing) = false
@inline ischeckpointing(alg::InterpolatingAdjoint) = alg.checkpointing
@inline ischeckpointing(alg::InterpolatingAdjoint, sol) = alg.checkpointing || !sol.dense
@inline ischeckpointing(alg::BacksolveAdjoint, sol=nothing) = alg.checkpointing

@inline isnoise(alg::DiffEqBase.AbstractSensitivityAlgorithm) = false
@inline isnoise(alg::InterpolatingAdjoint) = alg.noise
@inline isnoise(alg::BacksolveAdjoint) = alg.noise

@inline isnoisemixing(alg::DiffEqBase.AbstractSensitivityAlgorithm) = false
@inline isnoisemixing(alg::InterpolatingAdjoint) = alg.noisemixing
@inline isnoisemixing(alg::BacksolveAdjoint) = alg.noisemixing

@inline compile_tape(vjp::ReverseDiffVJP{compile}) where compile = compile
@inline compile_tape(noise::ReverseDiffNoise{compile}) where compile = compile
@inline compile_tape(autojacvec::Bool) = false
@inline compile_tape(sensealg::QuadratureAdjoint) = sensealg.compile

"""
ForwardDiffOverAdjoint{A} <: AbstractSecondOrderSensitivityAlgorithm{nothing,true,nothing}

ForwardDiff.jl over a choice of `sensealg` method for the adjoint.

## Constructor

```julia
ForwardDiffOverAdjoint(sensealg)
```

## SciMLProblem Support

This supports any SciMLProblem that the `sensealg` choice supports, provided the solver algorithm
is `SciMLBase.isautodifferentiable`.

## References

Hindmarsh, A. C. and Brown, P. N. and Grant, K. E. and Lee, S. L. and Serban, R.
and Shumaker, D. E. and Woodward, C. S., SUNDIALS: Suite of nonlinear and
differential/algebraic equation solvers, ACM Transactions on Mathematical
Software (TOMS), 31, pp:363–396 (2005)
"""
struct ForwardDiffOverAdjoint{A} <: AbstractSecondOrderSensitivityAlgorithm{nothing,true,nothing}
  adjalg::A
end


text:
@symbol_func function K_CS(cur_reactor::AbstractReactor)
  cur_K = cur_reactor.k

  cur_K *= cur_reactor.B_CS

  cur_K *= sqrt( K_LP(cur_reactor) )

  cur_K /= sqrt(1e5)

  cur_K
end

ext: jl


text:
export SpUri, SpId, SpCategoryId, SpUserId, SpUrl, SpPlaylistId, SpAlbumId,
       SpArtistId, SpShowId, SpEpisodeId
"""
All web API arguments are strings, but types
`SpUri`, `SpId`, `CategoryId`, `SpUserId`, `SpUrl`
aids in picking default values. [format](
https://developer.spotify.com/documentation/web-api/#spotify-uris-and-ids)

PARAMETER  | DESCRIPTION | VALUE
:----------:| :---------- |:------------:
SpUri       |The resource identifier that you can enter, for example, in the |
            | Spotify Desktop client’s search box to locate an artist,|
            | album, or track. To find a Spotify URI simply right-|
            | click (on Windows) or Ctrl-Click (on a Mac) on |spotify:track:
            | the artist’s or album’s or track’s name.| 6rqhFgbbKwnb9MLmUQDhG6|
SpId        |The base-62 identifier that you can find at the end of the |
            | Spotify URI (see above) for an artist, track, album, |
            | playlist, etc. Unlike a Spotify URI, a Spotify ID does|
            | not clearly identify the type of resource; that information is |
            | provided elsewhere in the call.| 6rqhFgbbKwnb9MLmUQDhG6
SpCategoryId|The unique string identifying the Spotify category.| party|
SpUserId    |The unique string identifying the Spotify user that you can |
            | find at the end of the Spotify URI for the user. The ID |
            | of the current user can be obtained via the Web API endpoint.| wizzler
SpUrl       |An HTML link that opens a track, album, app, playlist or other |
            | Spotify resource in a Spotify client (which  client |
            | is determined by the user’s device and |  http://open.spotify.com/
            | account settings at play.spotify.com. |   track/6rqhFgbbKwnb9MLmUQDhG6

"""
SpUri, SpId, SpCategoryId, SpUserId, SpUrl, SpPlaylistId, SpAlbumId, SpArtistId, SpShowId, SpEpisodeId

mutable struct SpUri
    s::String
    SpUri(s) = isuri(s) ? new(s) : error("must be 'spotify:track:<base 62 string>")
end
SpUri() = SpUri("spotify:track:6rqhFgbbKwnb9MLmUQDhG6")

mutable struct SpId
    s::String
    SpId(s) =  isid(s) ? new(s) : error("must be <base 62 string>")
end
SpId() = SpId("6rqhFgbbKwnb9MLmUQDhG6")

mutable struct SpCategoryId
    s::String
end
SpCategoryId() = SpCategoryId("party")

mutable struct SpUserId
    s::String
end
SpUserId() = get_user_name() == "" ? throw(".ini file has no user name") : SpUserId(get_user_name())

mutable struct SpUrl
    s::String
    SpUrl(s) = isurl(s) ? new(s) : error("must be an url")
end
SpUrl() =  SpUrl("http://open.spotify.com/track/6rqhFgbbKwnb9MLmUQDhG6")

mutable struct SpPlaylistId
    s::String
    SpPlaylistId(s) = isid(s) ? new(s) : error("must be <base 62 string> playlist ID")
end
SpPlaylistId() = SpPlaylistId("37i9dQZF1E4vUblDJbCkV3")

mutable struct SpAlbumId
    s::String
    SpAlbumId(s) = isid(s) ? new(s) : error("must be <base 62 string> album ID")
end
SpAlbumId() = SpAlbumId("5XgEM5g3xWEwL4Zr6UjoLo")

mutable struct SpArtistId
    s::String
    SpArtistId(s) = isid(s) ? new(s) : error("must be <base 62 string> artist ID")
end
SpArtistId() = SpArtistId("0YC192cP3KPCRWx8zr8MfZ")

mutable struct SpShowId
    s::String
    SpShowId(s) = isid(s) ? new(s) : error("must be <base 62 string> show ID")
end
SpShowId() = SpShowId("2MAi0BvDc6GTFvKFPXnkCL")

mutable struct SpEpisodeId
    s::String
    SpEpisodeId(s) = isid(s) ? new(s) : error("must be <base 62 string> episode ID")
end
SpEpisodeId() = SpEpisodeId("512ojhOuo1ktJprKbVcKyQ")

# Verify if the url, id and uri have the correct structure
isurl(s) = !isnothing(match(r"(http(s)?:\/\/.)?(www\.)?[-a-zA-Z0-9@:%._\+~#=]{2,256}\.[a-z]{2,6}\b([-a-zA-Z0-9@:%_\+.~#?&//=]*)", s))
isid(s) = !isnothing(match(r"\b[a-zA-Z0-9]{22}", s))
isuri(s) = count(==( ':'), s) == 2 && length(s) - findlast(':', s) == 22

# These types should work with string interpolation using $.
# So we need to alter their undecorated representations.
# For argument inspection, we don't want to hide them not
# being ordinary strings, so we add color. The color will
# not be used in contexts which can't display color (i.e.
# when used as an url)
show(io::IO, s::SpUri) =        printstyled(io, s.s; color = :blue)
show(io::IO, s::SpId) =         printstyled(io, s.s; color = 176)
show(io::IO, s::SpCategoryId) = printstyled(io, s.s; color = 172)
show(io::IO, s::SpUrl) =        printstyled(io, s.s; color = :blue, bold = true)
show(io::IO, s::SpUserId) =     printstyled(io, s.s; color = :green)
show(io::IO, s::SpPlaylistId) =     printstyled(io, s.s; color = :light_red)
show(io::IO, s::SpAlbumId) =     printstyled(io, s.s; color = :light_yellow)
show(io::IO, s::SpArtistId) =     printstyled(io, s.s; color = :light_magenta)
show(io::IO, s::SpShowId) =     printstyled(io, s.s; color = :light_cyan)
show(io::IO, s::SpEpisodeId) =     printstyled(io, s.s; color = :yellow)

# Add quotes "" around the String representation
show(io::IO,  ::MIME"text/plain", s::SpUri) =       printstyled(io, '"' * s.s * '"'; color = :blue)
show(io::IO, ::MIME"text/plain", s::SpId) =         printstyled(io, '"' * s.s * '"'; color = 176)
show(io::IO, ::MIME"text/plain", s::SpCategoryId) = printstyled(io, '"' * s.s * '"'; color = 172)
show(io::IO, ::MIME"text/plain", s::SpUrl) =        printstyled(io, '"' * s.s * '"'; color = :blue, bold = true)
show(io::IO, ::MIME"text/plain", s::SpUserId) =     printstyled(io, '"' * s.s * '"'; color = :green)
show(io::IO, ::MIME"text/plain", s::SpPlaylistId) =     printstyled(io, '"' * s.s * '"'; color = :light_red)
show(io::IO, ::MIME"text/plain", s::SpAlbumId) =     printstyled(io, '"' * s.s * '"'; color = :light_yellow)
show(io::IO, ::MIME"text/plain", s::SpArtistId) =     printstyled(io, '"' * s.s * '"'; color = :light_magenta)
show(io::IO, ::MIME"text/plain", s::SpShowId) =     printstyled(io, '"' * s.s * '"'; color = :light_cyan)
show(io::IO, ::MIME"text/plain", s::SpEpisodeId) =     printstyled(io, '"' * s.s * '"'; color = :yellow)

show(io::IO, v::Vector{SpUri}) =        show_vector(io, v, "", "")
show(io::IO, v::Vector{SpId}) =         show_vector(io, v, "", "")
show(io::IO, v::Vector{SpCategoryId}) = show_vector(io, v, "", "")
show(io::IO, v::Vector{SpUrl}) =        show_vector(io, v, "", "")
show(io::IO, v::Vector{SpUserId}) =     show_vector(io, v, "", "")
show(io::IO, v::Vector{SpPlaylistId}) = show_vector(io, v, "", "")
show(io::IO, v::Vector{SpAlbumId}) =    show_vector(io, v, "", "")
show(io::IO, v::Vector{SpArtistId}) =   show_vector(io, v, "", "")
show(io::IO, v::Vector{SpShowId}) =     show_vector(io, v, "", "")
show(io::IO, v::Vector{SpEpisodeId}) =  show_vector(io, v, "", "")

typeinfo_implicit(::Type{SpUri}) = true
typeinfo_implicit(::Type{SpId}) = true
typeinfo_implicit(::Type{SpCategoryId}) = true
typeinfo_implicit(::Type{SpUrl}) = true
typeinfo_implicit(::Type{SpUserId}) = true
typeinfo_implicit(::Type{SpPlaylistId}) = true
typeinfo_implicit(::Type{SpAlbumId}) = true
typeinfo_implicit(::Type{SpArtistId}) = true
typeinfo_implicit(::Type{SpShowId}) = true
typeinfo_implicit(::Type{SpEpisodeId}) = true


text:
using GLVisualize, GeometryTypes, Reactive, GLAbstraction, Colors

if !isdefined(:runtests)
    window = glscreen()
    timesignal = loop(linspace(0f0, 1f0, 360))
end
const n1 = 30
positions = rand(Point2f0, n1).*1000f0

xy_data(x,y,i) = Float32(sin(y/2f0/i)+cos(i*x))
const n2 = 128
# a distance field is a Matrix{Float32} array, which encodes the distance to
# the border of a filled shape. Positive numbers are inside the shape, 0 is the
# border and negative number are outside.
# this is basically how we render text, since you can do anti aliasing very
# nicely when you know the distance to the border.
# For text we use one big texture and specify uv coordinates into this big texture
# for every particle. How this is done can be seen in example
# partices/sprites/image_texture_atlas.jl

dfield = map(timesignal) do t
    tpi = (2pi*t)+0.2
    Float32[xy_data(x,y,tpi)+0.5f0 for x=1:n2, y=1:n2]
end
Base.rand(m::MersenneTwister, ::Type{U8}) = U8(rand(m, UInt8))
Base.rand{T <: Colorant}(m::MersenneTwister, ::Type{T}) = T(ntuple(x->rand(m, eltype(T)), Val{length(T)})...)

distfield = visualize((DISTANCEFIELD, positions),
    stroke_width=4f0,
    scale=Vec2f0(120),
    stroke_color=rand(RGBA{Float32}, n1),
    color=rand(RGBA{Float32}, n1),
    distancefield=dfield
)
_view(distfield, window)

if !isdefined(:runtests)
    renderloop(window)
end


text:
include("auto_encoder_task.jl")

# ::: NETWORK
networks = Dict(
    "500_300_150" => ([
        DenseLayer(INPUT_SIZE, 500),
        ActivationLayer(500, relu, Δrelu),
        DenseLayer(500, 300),
        ActivationLayer(300, relu, Δrelu),
        DenseLayer(300, 150),
        ActivationLayer(150, relu, Δrelu),
        DenseLayer(150, LATENT_SIZE),
        ActivationLayer(LATENT_SIZE, diminish, Δdiminish)
        ],[
        DenseLayer(LATENT_SIZE, 150),
        ActivationLayer(150, relu, Δrelu),
        DenseLayer(150, 300),
        ActivationLayer(300, relu, Δrelu),
        DenseLayer(300, 500),
        ActivationLayer(500, relu, Δrelu),
        DenseLayer(500, INPUT_SIZE),
        ActivationLayer(INPUT_SIZE, diminish, Δdiminish)
        ]),
   "600_400_200_100" => ([
        DenseLayer(INPUT_SIZE, 600),
        ActivationLayer(600, relu, Δrelu),
        DenseLayer(600, 400),
        ActivationLayer(400, relu, Δrelu),
        DenseLayer(400, 200),
        ActivationLayer(200, relu, Δrelu),
        DenseLayer(200, 100),
        ActivationLayer(100, relu, Δrelu),
        DenseLayer(100, LATENT_SIZE),
        ActivationLayer(LATENT_SIZE, diminish, Δdiminish)
        ],[
        DenseLayer(LATENT_SIZE, 100),
        ActivationLayer(100, relu, Δrelu),
        DenseLayer(100, 200),
        ActivationLayer(200, relu, Δrelu),
        DenseLayer(200, 400),
        ActivationLayer(400, relu, Δrelu),
        DenseLayer(400, 600),
        ActivationLayer(600, relu, Δrelu),
        DenseLayer(600, INPUT_SIZE),
        ActivationLayer(INPUT_SIZE, diminish, Δdiminish)
        ]),
    "500_500_250_250_125_125" => ([
        DenseLayer(INPUT_SIZE, 500),
        ActivationLayer(500, relu, Δrelu),
        DenseLayer(500, 500),
        ActivationLayer(500, relu, Δrelu),
        DenseLayer(500, 250),
        ActivationLayer(250, relu, Δrelu),
        DenseLayer(250, 250),
        ActivationLayer(250, relu, Δrelu),
        DenseLayer(250, 125),
        ActivationLayer(125, relu, Δrelu),
        DenseLayer(125, 125),
        ActivationLayer(125, relu, Δrelu),
        DenseLayer(125, LATENT_SIZE),
        ActivationLayer(LATENT_SIZE, diminish, Δdiminish)
       ],[
        DenseLayer(LATENT_SIZE, 125),
        ActivationLayer(125, relu, Δrelu),
        DenseLayer(125, 125),
        ActivationLayer(125, relu, Δrelu),
        DenseLayer(125, 250),
        ActivationLayer(250, relu, Δrelu),
        DenseLayer(250, 250),
        ActivationLayer(250, relu, Δrelu),
        DenseLayer(250, 500),
        ActivationLayer(500, relu, Δrelu),
        DenseLayer(500, 500),
        ActivationLayer(500, relu, Δrelu),
        DenseLayer(500, INPUT_SIZE),
        ActivationLayer(INPUT_SIZE, diminish, Δdiminish)
        ]),
    "500" => ([
        DenseLayer(INPUT_SIZE, 500),
        ActivationLayer(500, relu, Δrelu),
        DenseLayer(500, LATENT_SIZE),
        ActivationLayer(LATENT_SIZE, diminish, Δdiminish)
       ],[
        DenseLayer(LATENT_SIZE, 500),
        ActivationLayer(500, relu, Δrelu),
        DenseLayer(500, INPUT_SIZE),
        ActivationLayer(INPUT_SIZE, diminish, Δdiminish)
        ]))

# ::: TRAINING
runtime_plotting = false
plot_intervall = 20
begin
    tm1 = [NaN for _ in 1:TRAIN_EPISODES];
    tm2 = [NaN for _ in 1:TRAIN_EPISODES];
    tm3 = [NaN for _ in 1:TRAIN_EPISODES];
    tm4 = [NaN for _ in 1:TRAIN_EPISODES];

    println("---"); @time for e in 1:TRAIN_EPISODES
        t = time()

        z1 = forward!(networks["500_300_150"][1], train_data_x[e])
        z2 = forward!(networks["600_400_200_100"][1], train_data_x[e])
        z3 = forward!(networks["500_500_250_250_125_125"][1], train_data_x[e])
        z4 = forward!(networks["500"][1], train_data_x[e])

        y1 = forward!(networks["500_300_150"][2], z1)
        y2 = forward!(networks["600_400_200_100"][2], z2)
        y3 = forward!(networks["500_500_250_250_125_125"][2], z3)
        y4 = forward!(networks["500"][2], z4)

        dloss1 = -Δmse(train_data_x[e], y1)
        dloss2 = -Δmse(train_data_x[e], y2)
        dloss3 = -Δmse(train_data_x[e], y3)
        dloss4 = -Δmse(train_data_x[e], y4)

        dz1 = backward!(reverse(networks["500_300_150"][2]), dloss1, opt)
        _ = backward!(reverse(networks["500_300_150"][1]), dz1, opt)
        dz2 = backward!(reverse(networks["600_400_200_100"][2]), dloss2, opt)
        _ = backward!(reverse(networks["600_400_200_100"][1]), dz2, opt)
        dz3 = backward!(reverse(networks["500_500_250_250_125_125"][2]), dloss3, opt)
        _ = backward!(reverse(networks["500_500_250_250_125_125"][1]), dz3, opt)
        dz4 = backward!(reverse(networks["500"][2]), dloss4, opt)
        _ = backward!(reverse(networks["500"][1]), dz4, opt)

        tm1[e] = mse(train_data_x[e], y1)
        tm2[e] = mse(train_data_x[e], y2)
        tm3[e] = mse(train_data_x[e], y3)
        tm4[e] = mse(train_data_x[e], y4)

        # if runtime_plotting && e % plot_intervall == 0
        #     p = plot(FC_train_metric)
        #     display(p)
        # end

        # println("episode = $e : loss = $loss : dt = $(time() - t)")
    end
end

begin
    gr = 30
    plot([sum(tm1[i:i+gr])/gr for i in 1:TRAIN_EPISODES-gr], ylabel="MSE (lr=$LEARNING_RATE, batch size=$BATCH_SIZE)", xlabel="training iterations", label="FC [500_300_150]")
    plot!([sum(tm2[i:i+gr])/gr for i in 1:TRAIN_EPISODES-gr], xlabel="training iterations", label="FC [600_400_200_100]")
    plot!([sum(tm3[i:i+gr])/gr for i in 1:TRAIN_EPISODES-gr], xlabel="training iterations", label="FC [500_500_250_250_125_125]")
    plot!([sum(tm4[i:i+gr])/gr for i in 1:TRAIN_EPISODES-gr], xlabel="training iterations", label="FC [500]")
end

generate_images(10, "/Users/BBM2H16AHO/Documents/programming/DNN-Framework/auto_encoder/images/relu_sigmoid_test")


text:
module TestRoxygenize
	using Roxygen

	roxygenize("test", "output")
end


text:
abstract type CompiledDomain <: Domain end

get_source(::CompiledDomain) = error("Not implemented.")

function generate_domain_type(domain::Domain, state::State)
    name = pddl_to_type_name(get_name(domain))
    if domain isa AbstractedDomain
        domain_type = gensym("CompiledAbstracted" * name * "Domain")
    else
        domain_type = gensym("Compiled" * name * "Domain")
    end
    domain_typedef = :(struct $domain_type <: CompiledDomain end)
    domain_method_defs = generate_domain_methods(domain, domain_type)
    return (domain_type, domain_typedef, domain_method_defs)
end

function generate_domain_methods(domain::Domain, domain_type::Symbol)
    get_name_def = :(get_name(::$domain_type) = $(QuoteNode(get_name(domain))))
    get_source_def = :(get_source(::$domain_type) = $(QuoteNode(domain)))
    get_requirements_def = :(get_requirements(::$domain_type) =
            $(QuoteNode((; get_requirements(domain)...))))
    get_typetree_def = :(get_typetree(::$domain_type) =
        $(QuoteNode((; get_typetree(domain)...))))
    get_constants_def = :(get_constants(::$domain_type) =
        $(QuoteNode(Tuple(get_constants(domain)))))
    get_constypes_def = :(get_constypes(domain::$domain_type) =
        get_constypes(get_source(domain)))
    get_predicates_def = :(get_predicates(::$domain_type) =
        $(QuoteNode((; get_predicates(domain)...))))
    get_functions_def = :(get_functions(::$domain_type) =
        $(QuoteNode((; get_functions(domain)...))))
    get_funcdefs_def = :(get_funcdefs(::$domain_type) =
        $(QuoteNode((; get_funcdefs(domain)...))))
    get_fluents_def = :(get_fluents(::$domain_type) =
        $(QuoteNode((; get_fluents(domain)...))))
    get_axioms_def = :(get_axioms(::$domain_type) =
        $(QuoteNode((; get_axioms(domain)...))))
    domain_defs = Expr(:block, get_name_def, get_source_def,
        get_requirements_def, get_typetree_def, get_constants_def,
        get_constypes_def, get_predicates_def, get_functions_def,
        get_funcdefs_def, get_fluents_def, get_axioms_def)
end


text:
using Logging

using StructuralIdentifiability

#QWWC
logger = Logging.SimpleLogger(stdout, Logging.Debug)
global_logger(logger)

ode = @ODEmodel(
    x'(t) = a * (y(t) - x(t)) + y(t) * z(t),
    y'(t) = b * (x(t) + y(t)) - x(t) * z(t),
    z'(t) = - c * z(t) - d * w(t) + x(t) * y(t),
    w'(t) = e * z(t) - f * w(t) + x(t) * y(t),
    g(t) = x(t)
)

println(assess_identifiability(ode))


text:
using Complementarity
using Test

@info("-------[Testing Complementarity/NLsolve]------------------------------------------")

@testset "mcp test 0 with NLsolve example" begin

    lb = [0., 0., 0., 0.]
    ub = [Inf, Inf, Inf, Inf]
    x0 = [1.25, 0., 0., 0.5]

    m = nothing
    m = MCPModel()

    @variable(m, x[1:4] >= 0)
    @mapping(m, F1, 3*x[1]^2+2*x[1]*x[2]+2*x[2]^2+x[3]+3*x[4]-6)
    @mapping(m, F2, 2*x[1]^2+x[1]+x[2]^2+3*x[3]+2*x[4]-2)
    @mapping(m, F3, 3*x[1]^2+x[1]*x[2]+2*x[2]^2+2*x[3]+3*x[4]-1)
    @mapping(m, F4, x[1]^2+3*x[2]^2+2*x[3]+3*x[4]-3)
    @complementarity(m, [F1 F2 F3 F4], x)

    set_start_value.(x, x0)

    status = solveMCP(m, solver=:NLsolve)
    @show status

    z = result_value.(x)
    # Fz = [result_value(F1) result_value(F2) result_value(F3) result_value(F4)]

    @show z
    # @show Fz

    @test isapprox(z[1], 1.22474, atol=1e-4)
    @test isapprox(z[2], 0.0, atol=1e-4)
    @test isapprox(z[3], -1.378e-19, atol=1e-4)
    @test isapprox(z[4], 0.5, atol=1e-4)

    # @test isapprox(Fz[1], -1.26298e-9, atol=1e-4)
    # @test isapprox(Fz[2], 3.22474, atol=1e-4)
    # @test isapprox(Fz[3], 5.0, atol=1e-4)
    # @test isapprox(Fz[4], 3.62723e-11, atol=1e-4)

    set_start_value.(x, x0)
    status = solveMCP(m, solver=:PATH)
    @show status

    z = result_value.(x)
    # Fz = [result_value(F1) result_value(F2) result_value(F3) result_value(F4)]

    @show z
    # @show Fz

    @test isapprox(z[1], 1.22474, atol=1e-4)
    @test isapprox(z[2], 0.0, atol=1e-4)
    @test isapprox(z[3], -1.378e-19, atol=1e-4)
    @test isapprox(z[4], 0.5, atol=1e-4)

    # @test isapprox(Fz[1], -1.26298e-9, atol=1e-4)
    # @test isapprox(Fz[2], 3.22474, atol=1e-4)
    # @test isapprox(Fz[3], 5.0, atol=1e-4)
    # @test isapprox(Fz[4], 3.62723e-11, atol=1e-4)

    set_start_value.(x, x0)
    status = solveMCP(m, solver=:NLsolve)
    @show status

    z = result_value.(x)
    # Fz = [result_value(F1) result_value(F2) result_value(F3) result_value(F4)]

    @show z
    # @show Fz

    @test isapprox(z[1], 1.22474, atol=1e-4)
    @test isapprox(z[2], 0.0, atol=1e-4)
    @test isapprox(z[3], -1.378e-19, atol=1e-4)
    @test isapprox(z[4], 0.5, atol=1e-4)

    # @test isapprox(Fz[1], -1.26298e-9, atol=1e-4)
    # @test isapprox(Fz[2], 3.22474, atol=1e-4)
    # @test isapprox(Fz[3], 5.0, atol=1e-4)
    # @test isapprox(Fz[4], 3.62723e-11, atol=1e-4)

end

println("------------------------------------------------------------------")

@testset "mcp test 1 with NLsolve" begin
    m = nothing
    m = MCPModel()

    M = [0  0 -1 -1 ;
         0  0  1 -2 ;
         1 -1  2 -2 ;
         1  2 -2  4 ]

    q = [2; 2; -2; -6]

    lb = zeros(4)
    ub = Inf*ones(4)

    items = 1:4

    # @variable(m, lb[i] <= x[i in items] <= ub[i])
    @variable(m, x[i in items] >= 0)
    @mapping(m, F[i in items], sum(M[i,j]*x[j] for j in items) + q[i])
    @complementarity(m, F, x)

    status = solveMCP(m, solver=:NLsolve)
    @show status

    z = result_value.(x)
    # Fz = result_value(F)

    @show z
    # @show Fz

    @test isapprox(z[1], 2.8, atol=1e-4)
    @test isapprox(z[2], 0.0, atol=1e-4)
    @test isapprox(z[3], 0.8, atol=1e-4)
    @test isapprox(z[4], 1.2, atol=1e-4)
end

println("------------------------------------------------------------------")

@testset "mcp test 2 with NLsolve" begin

    m = nothing
    m = MCPModel()

    M = [0  0 -1 -1 ;
         0  0  1 -2 ;
         1 -1  2 -2 ;
         1  2 -2  4 ]

    q = [2; 2; -2; -6]

    lb = zeros(4)
    ub = Inf*ones(4)

    items = 1:4

    @variable(m, lb[i] <= x[i in items] <= ub[i])
    # @variable(m, x[i in items] >= 0)
    @mapping(m, F[i in items], sum(M[i,j]*x[j] for j in items) + q[i])
    @complementarity(m, F, x)

    status = solveMCP(m, solver=:NLsolve)
    @show status

    z = result_value.(x)
    # Fz = result_value(F) # currently produces an error

    @show z
    # @show Fz

    @test isapprox(z[1], 2.8, atol=1e-4)
    @test isapprox(z[2], 0.0, atol=1e-4)
    @test isapprox(z[3], 0.8, atol=1e-4)
    @test isapprox(z[4], 1.2, atol=1e-4)
end

println("------------------------------------------------------------------")

@testset "mcp test 3 with NLsolve" begin

    m = nothing
    m = MCPModel()

    M = [0  0 -1 -1 ;
         0  0  1 -2 ;
         1 -1  2 -2 ;
         1  2 -2  4 ]

    q = [2; 2; -2; -6]

    lb = zeros(4)
    ub = Inf*ones(4)

    @variable(m, lb[i] <= myvariablename[i in 1:4] <= ub[i])
    @mapping(m, myconst[i=1:4], sum(M[i,j]*myvariablename[j] for j in 1:4) + q[i])
    @complementarity(m, myconst, myvariablename)

    status = solveMCP(m, solver=:NLsolve)
    @show status

    z = result_value.(myvariablename)
    # Fz = result_value(myconst)

    @show z
    # @show Fz

    @test isapprox(z[1], 2.8, atol=1e-4)
    @test isapprox(z[2], 0.0, atol=1e-4)
    @test isapprox(z[3], 0.8, atol=1e-4)
    @test isapprox(z[4], 1.2, atol=1e-4)

end

println("------------------------------------------------------------------")

@testset "mcp test 4 with NLsolve" begin

    m = nothing
    m = MCPModel()

    lb = zeros(4)
    ub = Inf*ones(4)
    items = 1:4
    @variable(m, lb[i] <= x[i in items] <= ub[i])

    @mapping(m, F1, 3*x[1]^2 + 2*x[1]*x[2] + 2*x[2]^2 + x[3] + 3*x[4] -6)
    @mapping(m, F2, 2*x[1]^2 + x[1] + x[2]^2 + 3*x[3] + 2*x[4] -2)
    @mapping(m, F3, 3*x[1]^2 + x[1]*x[2] + 2*x[2]^2 + 2*x[3] + 3*x[4] -1)
    @mapping(m, F4, x[1]^2 + 3*x[2]^2 + 2*x[3] + 3*x[4] - 3)

    @complementarity(m, F1, x[1])
    @complementarity(m, F2, x[2])
    @complementarity(m, F3, x[3])
    @complementarity(m, F4, x[4])

    set_start_value(x[1], 1.25)
    set_start_value(x[2], 0.)
    set_start_value(x[3], 0.)
    set_start_value(x[4], 0.5)

    status = solveMCP(m, solver=:NLsolve)
    @show status

    z = result_value.(x)
    # Fz = [result_value(F1), result_value(F2), result_value(F3), result_value(F4)]

    @show z
    # @show Fz

    @test isapprox(z[1], 1.22474, atol=1e-4)
    @test isapprox(z[2], 0.0, atol=1e-4)
    @test isapprox(z[3], 0.0, atol=1e-4)
    @test isapprox(z[4], 0.5, atol=1e-4)

    set_start_value(x[1], 1.19)
    set_start_value(x[2], 0.2)
    set_start_value(x[3], 0.1)
    set_start_value(x[4], 0.49)

    status = solveMCP(m, solver=:PATH)
    @show status

    z = result_value.(x)
    # Fz = [result_value(F1), result_value(F2), result_value(F3), result_value(F4)]

    @show z
    # @show Fz

    @test isapprox(z[1], 1.22474, atol=1e-4)
    @test isapprox(z[2], 0.0, atol=1e-4)
    @test isapprox(z[3], 0.0, atol=1e-4)
    @test isapprox(z[4], 0.5, atol=1e-4)

    set_start_value(x[1], 1.40)
    set_start_value(x[2], 0.1)
    set_start_value(x[3], 0.2)
    set_start_value(x[4], 0.43)

    status = solveMCP(m, solver=:NLsolve)
    @show status

    z = result_value.(x)
    # Fz = [result_value(F1), result_value(F2), result_value(F3), result_value(F4)]

    @show z
    # @show Fz

    @test isapprox(z[1], 1.22474, atol=1e-4)
    @test isapprox(z[2], 0.0, atol=1e-4)
    @test isapprox(z[3], 0.0, atol=1e-4)
    @test isapprox(z[4], 0.5, atol=1e-4)
end


text:
using LinearAlgebra
using Printf

function getInverse(U::ITensor,s::ITensor,V::ITensor)::ITensor
    as=LinearAlgebra.diag(array(s))
    asinv=Vector(as)
    for ia in 1:length(asinv) asinv[ia]=1.0/asinv[ia] end
    sinv=diagITensor(eltype(asinv),asinv,inds(s))
    Minv=dag(V)*dag(sinv)*dag(U)
    return Minv
end

#
# Solve RL=M*RL_prime we only need to solve the last few columns to the right of
# the unit matrix inside RL_prime
#
function SolveRLprime(RL::ITensor,RL_prime::ITensor,U::ITensor,s::ITensor,V::ITensor,iq::Index,im::Index,ms::matrix_state)::ITensor
    Minv=getInverse(U,s,V)
    ic1=noncommonind(RL_prime,im) #RL_prime col index
    ic2=noncommonind(RL      ,iq) #RL       col index
    @assert ic1==ic2
    @assert dim(iq)==dim(im)
    #eps=1e-14
    #@show "RL="
    #pprint(iq,RL,ic1,eps)
    #@show "RL_prime="
    #pprint(im,RL_prime,ic1,eps)
    Dm,Dc=dim(im),dim(ic1)
    @assert Dm>2
    @assert Dc>Dm
    imm=filterinds(Minv,tags=tags(im))[1]
    iqm=filterinds(Minv,tags=tags(iq))[1]
    @assert dim(imm)==dim(im)-2
    @assert dim(iqm)==dim(iq)-2
    icm=Index(Dc-Dm,tags(ic1))
    R2=ITensor(0.0,iqm,icm)
    j1_offset=1
    if ms.ul==lower
        j2_offset=Dm-1
    else #upper
        j2_offset=1
    end
    #@show ms Dm j2_offset

    for j2 in eachindval(icm)
        for j1 in eachindval(iqm)
            R2[j1,j2]=RL[iq=>j1.second+j1_offset,ic1=>j2.second+j2_offset]
        end
    end
    # @show "R2="
    # pprint(iqm,R2,icm,1e-14)
    R2_prime=Minv*R2

    #@show inds(R2_prime)
    #pprint(imm,R2_prime,icm,eps)
    for j2 in eachindval(icm)
        for j1 in eachindval(imm)
            RL_prime[im=>j1.second+j1_offset,ic1=>j2.second+j2_offset]=R2_prime[j1,j2]
        end
    end
    #pprint(im,RL_prime,ic1,eps)

    return RL_prime
end
#
#  Compress one site
#
function truncate(W::ITensor,ul::reg_form;kwargs...)::Tuple{ITensor,ITensor,bond_spectrum}
    d,n,r,c=parse_links(W) # W[l=$(n-1)l=$n]=W[r,c]
    lr::orth_type=get(kwargs, :orth, right)
    ms=matrix_state(ul,lr)
    eps=1e-14 #relax used for testing upper/lower/regular-form etc.
    # establish some tag strings then depend on lr.
    (tsvd,tuv,tln) = lr==left ? ("qx","u","l=$n") : ("m","v","l=$(n-1)")

#
# Block repecting QR/QL/LQ/RQ factorization.  RL=L or R for upper and lower.
# here we purposely turn off rank reavealing feature (epsrr=0.0) to (mostly) avoid
# horizontal rectangular RL matricies which are hard to handle accurately.
#
    Q,RL,lq=block_qx(W,ul;epsrr=0.0,kwargs...) #left Q[r,qx], RL[qx,c] - right RL[r,qx] Q[qx,c]
    ITensors.@debug_check begin
        if order(Q)==4
            @assert is_canonical(Q,ms,eps)
            @assert is_regular_form(Q,ul,eps)
        end
    end
#
#  Factor RL=M*L' (left/lower) = L'*M (right/lower) = M*R' (left/upper) = R'*M (right/upper)
#
    c=noncommonind(RL,lq) #if size changed the old c is not lnger valid
    M,RL_prime,im,RLnz=getM(RL,ms,eps) #left M[lq,im] RL_prime[im,c] - right RL_prime[r,im] M[im,lq]
#
#  At last we can svd and compress M using epsSVD as the cutoff.
#
    min_cutoff=1e-12
    cutoff=get(kwargs,:cutoff,min_cutoff)
    cutoff=cutoff<min_cutoff ? min_cutoff : cutoff
    kwargs=add_or_replace(kwargs,:cutoff,cutoff) #make sure cutoff is not zero.
    isvd=findinds(M,tsvd)[1] #decide the left index
    U,s,V=svd(M,isvd;kwargs...) # ns sing. values survive compression
    ns=dim(inds(s)[1])

    spectrum=bond_spectrum(s,n)
    #@show diag(array(s))
    #
    #  If RL is rectangular we need to solve RL=M*RL_prime for RL_prime
    #  since we know UsV anyway we can calculate M^-1=dag(V)*1.0/s*dag(U)
    #
    if ns>0 && RLnz
        @show "fixing RL_prime"
        RL_prime=SolveRLprime(RL,RL_prime,U,s,V,lq,im,ms)
    end
    Mplus=grow(M,lq,im)
    D=RL-Mplus*RL_prime
    if norm(D)>get(kwargs, :cutoff, 1e-14)
        @printf "High normD(D)=%.1e min(s)=%.1e \n" norm(D) Base.min(diag(array(s))...)
    end

    luv=Index(ns+2,"Link,$tuv") #link for expanded U,US,V,sV matricies.
    if lr==left
        ITensors.@debug_check begin
            @assert is_upper_lower(lq,RL      ,c,ul,eps)
        end
        RL=grow(s*V,luv,im)*RL_prime #RL[l=n,u] dim ns+2 x Dw2
        W=Q*grow(U,lq,luv) #W[l=n-1,u]
    else # right
        ITensors.@debug_check begin
            @assert is_upper_lower(r,RL      ,lq,ms.ul,eps)
        end
        RL=RL_prime*grow(U*s,im,luv) #RL[l=n-1,v] dim Dw1 x ns+2
        W=grow(V,lq,luv)*Q #W[l=n-1,v]
    end
    replacetags!(RL,tuv,tln) #RL[l=n,l=n] sames tags, different id's and possibly diff dimensions.
    replacetags!(W ,tuv,tln) #W[l=n-1,l=n]
    ITensors.@debug_check begin
        @assert is_regular_form(W,ul,eps)
        @assert is_canonical(W,ms,eps)
    end
    return W,RL,spectrum
end

@doc """
    truncate!(H::MPO)

Compress an MPO using block respecting SVD techniques as described in
> *Daniel E. Parker, Xiangyu Cao, and Michael P. Zaletel Phys. Rev. B 102, 035147*

# Arguments
- `H` MPO for decomposition. If H is not already in the correct canonical form for compression, it will automatically be put into the correct form prior to compression.

# Keywords
- `orth::orth_type = left` : choose `left` or `right` canonical form for the final output.
- `cutoff::Float64 = 1e-14` : Using a `cutoff` allows the SVD algorithm to truncate as many states as possible while still ensuring a certain accuracy.
- `maxdim::Int64` : If the number of singular values exceeds `maxdim`, only the largest `maxdim` will be retained.
- `mindim::Int64` : At least `mindim` singular values will be retained, even if some fall below the cutoff

# Example
```julia
julia> using ITensors
julia> using ITensorMPOCompression
julia> N=10; #10 sites
julia> NNN=7; #Include up to 7th nearest neighbour interactions
julia> sites = siteinds("S=1/2",N);
#
# This makes H directly, bypassing autoMPO.  (AutoMPO is too smart for this
# demo, it makes maximally reduced MPOs right out of the box!)
#
julia> H=make_transIsing_MPO(sites,NNN);
#
#  Make sure we have a regular form or truncate! won't work.
#
julia> is_lower_regular_form(H)==true
true

#
#  Now we can truncate with defaults of left orthogonal cutoff=1e-14.
#  truncate! returns the spectrum of singular values at each bond.  The largest
#  singular values are remaining well under control.  i.e. no sign of divergences.
#
julia> truncate!(H)
9-element Vector{bond_spectrum}:
 bond_spectrum([0.307], 1)
 bond_spectrum([0.354, 0.035], 2)
 bond_spectrum([0.375, 0.045, 0.021], 3)
 bond_spectrum([0.385, 0.044, 0.026, 0.018], 4)
 bond_spectrum([0.388, 0.043, 0.031, 0.019, 0.001], 5)
 bond_spectrum([0.385, 0.044, 0.026, 0.018], 6)
 bond_spectrum([0.375, 0.045, 0.021], 7)
 bond_spectrum([0.354, 0.035], 8)
 bond_spectrum([0.307], 9)

julia> pprint(H[2])
I 0 0 0
S S S 0
0 S S I

#
#  We can see that bond dimensions have been drastically reduced.
#
julia> get_Dw(H)
9-element Vector{Int64}: 3 4 5 6 7 6 5 4 3

julia> is_lower_regular_form(H)==true
true

julia> is_orthogonal(H,left)==true
true

```
"""
function truncate!(H::MPO;kwargs...)::bond_spectrums
    #@printf "---- start compress ----\n"
    #
    # decide left/right and upper/lower
    #
    lr::orth_type=get(kwargs, :orth, left)
    kwargs=add_or_replace(kwargs,:orth,lr) #if lr is not yet in kwargs, we need to stuff in there
    (bl,bu)=detect_regular_form(H)
    if !(bl || bu)
        throw(ErrorException("truncate!(H::MPO), H must be in either lower or upper regular form"))
    end
    @assert !(bl && bu)
    ul::reg_form = bl ? lower : upper #if both bl and bu are true then something is seriously wrong
    #
    # Now check if H required orthogonalization
    #
    ms=matrix_state(ul,lr)
    if !is_canonical(H,mirror(ms))
        epsrr=get(kwargs,:epsrr,1e-12)
        epsrr= epsrr==0 ? 1e-12 : epsrr #0.0 not allwed here.
        orthogonalize!(H,epsrr=epsrr,orth=mirror(lr)) #TODO why fail if spec ul here??
    end
    N=length(H)
    ss=bond_spectrums(undef,N-1)
    if lr==left
        rng=1:1:N-1 #sweep left to right
        link_offest=0
    else #right
        rng=N:-1:2 #sweep right to left
        link_offest=-1
    end
    for n in rng
        nn=n+rng.step #index to neighbour
        W,RL,s=truncate(H[n],ul;kwargs...)
        H[n]=W
        H[nn]=RL*H[nn]
        is_regular_form(H[nn],ms.ul)
        ss[n+link_offest]=s
    end
    return ss
end


text:
using PowerSystems
using JuMP

base_dir = string(dirname(dirname(pathof(PowerSystems))))
println(joinpath(base_dir,"data/data_5bus.jl"))
include(joinpath(base_dir,"data/data_5bus.jl"))

generators_hg = [
    HydroFix("HydroFix",true,nodes5[2],
        TechHydro(60.0, 15.0, (min = 0.0, max = 60.0), nothing, nothing, nothing, nothing),
        TimeSeries.TimeArray(DayAhead,solar_ts_DA)
    ),
    HydroCurtailment("HydroCurtailment",true,nodes5[3],
        TechHydro(60.0, 10.0, (min = 0.0, max = 60.0), nothing, nothing, (up = 10.0, down = 10.0), nothing),
        100.0,TimeSeries.TimeArray(DayAhead,wind_ts_DA) )
]

m = Model()
sys5b = PowerSystem(nodes5, append!(generators5, generators_hg), loads5_DA, branches5, nothing,  100.0)

test_hy = [d for d in sys5b.generators.hydro if !isa(d, PowerSystems.HydroFix)] # Filter StaticLoads Out
phg, inyection_array = PowerSimulations.activepowervariables(m, devices_netinjection,  test_hy, sys5.time_periods)

true


text:
# interface

@inline getindex(b::AbstractVector, K::BlockIndex{1}) = b[Block(K.I[1])][K.α[1]]
@inline getindex(b::AbstractArray{T,N}, K::BlockIndex{N}) where {T,N} =
    b[block(K)][K.α...]
@inline getindex(b::AbstractArray, K::BlockIndex{1}, J::BlockIndex{1}...) =
    b[BlockIndex(tuple(K, J...))]

@inline getindex(b::AbstractArray{T,N}, K::BlockIndexRange{N}) where {T,N} = b[block(K)][K.indices...]
@inline getindex(b::LayoutArray{T,N}, K::BlockIndexRange{N}) where {T,N} = b[block(K)][K.indices...]
@inline getindex(b::LayoutArray{T,1}, K::BlockIndexRange{1}) where {T} = b[block(K)][K.indices...]

function findblockindex(b::AbstractVector, k::Integer)
    K = findblock(b, k)
    K[searchsortedfirst(b[K], k)] # guaranteed to be in range
end

function _BlockedUnitRange end

"""
    BlockedUnitRange

is an `AbstractUnitRange{Int}` that has been divided
into blocks, and is used to represent axes of block arrays.
Construction is typically via `blockrange` which converts
a vector of block lengths to a `BlockedUnitRange`.
```jldoctest; setup = quote using BlockArrays end
julia> blockedrange([2,2,3])
3-blocked 7-element BlockedUnitRange{Vector{Int64}}:
 1
 2
 ─
 3
 4
 ─
 5
 6
 7
```
"""
struct BlockedUnitRange{CS} <: AbstractUnitRange{Int}
    first::Int
    lasts::CS
    global _BlockedUnitRange(f, cs::CS) where CS = new{CS}(f, cs)
end

const DefaultBlockAxis = BlockedUnitRange{Vector{Int}}

@inline _BlockedUnitRange(cs) = _BlockedUnitRange(1,cs)

BlockedUnitRange(::BlockedUnitRange) = throw(ArgumentError("Forbidden due to ambiguity"))
_blocklengths2blocklasts(blocks) = cumsum(blocks) # extra level to allow changing default cumsum behaviour
@inline blockedrange(blocks::Union{Tuple,AbstractVector}) = _BlockedUnitRange(_blocklengths2blocklasts(blocks))

@inline blockfirsts(a::BlockedUnitRange) = [a.first; @views(a.lasts[1:end-1]) .+ 1]
@inline blocklasts(a::BlockedUnitRange) = a.lasts

_diff(a::AbstractVector) = diff(a)
_diff(a::Tuple) = diff(collect(a))
@inline blocklengths(a::BlockedUnitRange) = isempty(a.lasts) ? _diff(a.lasts) : [first(a.lasts)-a.first+1; _diff(a.lasts)]

length(a::BlockedUnitRange) = isempty(a.lasts) ? 0 : Integer(last(a.lasts)-a.first+1)

"""
   blockisequal(a::AbstractUnitRange{Int}, b::AbstractUnitRange{Int})

returns true if a and b have the same block structure.
"""
blockisequal(a::AbstractUnitRange{Int}, b::AbstractUnitRange{Int}) = first(a) == first(b) && blocklasts(a) == blocklasts(b)
blockisequal(a, b, c, d...) = blockisequal(a,b) && blockisequal(b,c,d...)
"""
   blockisequal(a::Tuple, b::Tuple)

returns true if `all(blockisequal.(a,b))`` is true.
"""
blockisequal(a::Tuple, b::Tuple) = all(blockisequal.(a, b))

Base.convert(::Type{BlockedUnitRange}, axis::BlockedUnitRange) = axis
Base.convert(::Type{BlockedUnitRange}, axis::AbstractUnitRange{Int}) = _BlockedUnitRange(first(axis),[last(axis)])
Base.convert(::Type{BlockedUnitRange}, axis::Base.Slice) = _BlockedUnitRange(first(axis),[last(axis)])
Base.convert(::Type{BlockedUnitRange}, axis::Base.IdentityUnitRange) = _BlockedUnitRange(first(axis),[last(axis)])
Base.convert(::Type{BlockedUnitRange{CS}}, axis::BlockedUnitRange{CS}) where CS = axis
Base.convert(::Type{BlockedUnitRange{CS}}, axis::BlockedUnitRange) where CS = _BlockedUnitRange(first(axis), convert(CS, blocklasts(axis)))
Base.convert(::Type{BlockedUnitRange{CS}}, axis::AbstractUnitRange{Int}) where CS = convert(BlockedUnitRange{CS}, convert(BlockedUnitRange, axis))

Base.unitrange(b::BlockedUnitRange) = first(b):last(b)

"""
    blockaxes(A)

Return the tuple of valid block indices for array `A`.
```jldoctest; setup = quote using BlockArrays end
julia> A = BlockArray([1,2,3],[2,1])
2-blocked 3-element BlockVector{Int64}:
 1
 2
 ─
 3

julia> blockaxes(A)[1]
2-element BlockRange{1, Tuple{Base.OneTo{Int64}}}:
 Block(1)
 Block(2)
```
"""
blockaxes(b::BlockedUnitRange) = _blockaxes(b.lasts)
_blockaxes(b::AbstractVector) = (Block.(axes(b,1)),)
_blockaxes(b::Tuple) = (Block.(Base.OneTo(length(b))),)
blockaxes(b) = blockaxes.(axes(b), 1)

"""
    blockaxes(A, d)

Return the valid range of block indices for array `A` along dimension `d`.
```jldoctest; setup = quote using BlockArrays end
julia> A = BlockArray([1,2,3],[2,1])
2-blocked 3-element BlockVector{Int64}:
 1
 2
 ─
 3

julia> blockaxes(A,1)
2-element BlockRange{1, Tuple{Base.OneTo{Int64}}}:
 Block(1)
 Block(2)
```
"""
function blockaxes(A::AbstractArray{T,N}, d) where {T,N}
    @_inline_meta
    d::Integer <= N ? blockaxes(A)[d] : Base.OneTo(1)
end

"""
    blocksize(A)
    blocksize(A,i)

Return the tuple of the number of blocks along each
dimension. See also size and blocksizes.
```jldoctest; setup = quote using BlockArrays end
julia> A = BlockArray(ones(3,3),[2,1],[1,1,1])
2×3-blocked 3×3 BlockMatrix{Float64}:
 1.0  │  1.0  │  1.0
 1.0  │  1.0  │  1.0
 ─────┼───────┼─────
 1.0  │  1.0  │  1.0

julia> blocksize(A)
(2, 3)

julia> blocksize(A,2)
3
```
"""
blocksize(A) = map(length, blockaxes(A))
blocksize(A,i) = length(blockaxes(A,i))
blocklength(t) = (@_inline_meta; prod(blocksize(t)))

"""
    blocksizes(A)
    blocksizes(A,i)

Return the tuple of the sizes of blocks along each
dimension. See also size and blocksize.
```jldoctest; setup = quote using BlockArrays end
julia> A = BlockArray(ones(3,3),[2,1],[1,1,1])
2×3-blocked 3×3 BlockMatrix{Float64}:
 1.0  │  1.0  │  1.0
 1.0  │  1.0  │  1.0
 ─────┼───────┼─────
 1.0  │  1.0  │  1.0

julia> blocksizes(A)
([2, 1], [1, 1, 1])

julia> blocksizes(A,2)
3-element Vector{Int64}:
 1
 1
 1
```
"""
blocksizes(A) = map(blocklengths, axes(A))
blocksizes(A,i) = blocklengths(axes(A,i))

axes(b::BlockedUnitRange) = (_BlockedUnitRange(blocklasts(b) .- (first(b)-1)),)
unsafe_indices(b::BlockedUnitRange) = axes(b)
first(b::BlockedUnitRange) = b.first
last(b::BlockedUnitRange) = isempty(blocklasts(b)) ? first(b)-1 : last(blocklasts(b))

function getindex(b::BlockedUnitRange, K::Block{1})
    k = Integer(K)
    bax = blockaxes(b,1)
    cs = blocklasts(b)
    @boundscheck K in bax || throw(BlockBoundsError(b, k))
    S = first(bax)
    K == S && return first(b):first(cs)
    return cs[k-1]+1:cs[k]
end

function getindex(b::BlockedUnitRange, KR::BlockRange{1})
    cs = blocklasts(b)
    isempty(KR) && return _BlockedUnitRange(1,cs[1:0])
    K,J = first(KR),last(KR)
    k,j = Integer(K),Integer(J)
    bax = blockaxes(b,1)
    @boundscheck K in bax || throw(BlockBoundsError(b,K))
    @boundscheck J in bax || throw(BlockBoundsError(b,J))
    K == first(bax) && return _BlockedUnitRange(first(b),cs[k:j])
    _BlockedUnitRange(cs[k-1]+1,cs[k:j])
end

function getindex(b::BlockedUnitRange, KR::BlockRange{1,Tuple{Base.OneTo{Int}}})
    cs = blocklasts(b)
    isempty(KR) && return _BlockedUnitRange(1,cs[Base.OneTo(0)])
    J = last(KR)
    j = Integer(J)
    bax = blockaxes(b,1)
    @boundscheck J in bax || throw(BlockBoundsError(b,J))
    _BlockedUnitRange(first(b),cs[Base.OneTo(j)])
end

_searchsortedfirst(a::AbstractVector, k) = searchsortedfirst(a, k)
function _searchsortedfirst(a::Tuple, k)
    k ≤ first(a) && return 1
    1+_searchsortedfirst(tail(a), k)
end
_searchsortedfirst(a::Tuple{}, k) = 1

function findblock(b::BlockedUnitRange, k::Integer)
    @boundscheck k in b || throw(BoundsError(b,k))
    Block(_searchsortedfirst(blocklasts(b), k))
end

Base.dataids(b::BlockedUnitRange) = Base.dataids(blocklasts(b))

###
# BlockedUnitRange interface
###
Base.checkindex(::Type{Bool}, b::BlockRange, K::Int) = checkindex(Bool, Int.(b), K)
Base.checkindex(::Type{Bool}, b::AbstractUnitRange{Int}, K::Block{1}) = checkindex(Bool, blockaxes(b,1), Int(K))

function getindex(b::AbstractUnitRange{Int}, K::Block{1})
    @boundscheck K == Block(1) || throw(BlockBoundsError(b, K))
    b
end

function getindex(b::AbstractUnitRange{Int}, K::BlockRange)
    @boundscheck K == Block.(1:1) || throw(BlockBoundsError(b, K))
    b
end

blockaxes(b::AbstractUnitRange{Int}) = (Block.(Base.OneTo(1)),)

function findblock(b::AbstractUnitRange{Int}, k::Integer)
    @boundscheck k in axes(b,1) || throw(BoundsError(b,k))
    Block(1)
end

"""
   blockfirsts(a::AbstractUnitRange{Int})

returns the first index of each block of `a`.
"""
blockfirsts(a::AbstractUnitRange{Int}) = Ones{Int}(1)
"""
   blocklasts(a::AbstractUnitRange{Int})

returns the last index of each block of `a`.
"""
blocklasts(a::AbstractUnitRange{Int}) = Fill(length(a),1)
"""
   blocklengths(a::AbstractUnitRange{Int})

returns the length of each block of `a`.
"""
blocklengths(a::AbstractUnitRange) = blocklasts(a) .- blockfirsts(a) .+ 1

Base.summary(a::BlockedUnitRange) = _block_summary(a)
Base.summary(io::IO, a::BlockedUnitRange) =  _block_summary(io, a)

###
# Slice{<:BlockedUnitRange}
###

Base.axes(S::Base.Slice{<:BlockedUnitRange}) = (S.indices,)
Base.unsafe_indices(S::Base.Slice{<:BlockedUnitRange}) = (S.indices,)
Base.axes1(S::Base.Slice{<:BlockedUnitRange}) = S.indices
blockaxes(S::Base.Slice) = blockaxes(S.indices)
getindex(S::Base.Slice, b::Block{1}) = S.indices[b]
getindex(S::Base.Slice, b::BlockRange{1}) = S.indices[b]

# This supports broadcasting with infinite block arrays
Base.BroadcastStyle(::Type{BlockedUnitRange{R}}) where R = Base.BroadcastStyle(R)

###
# Special Fill/Range cases
#
# We want to use lazy types when possible
###

_blocklengths2blocklasts(blocks::AbstractRange) = RangeCumsum(blocks)
function blockfirsts(a::BlockedUnitRange{Base.OneTo{Int}})
    a.first == 1 || error("Offset axes not supported")
    Base.OneTo{Int}(length(a.lasts))
end
function blocklengths(a::BlockedUnitRange{Base.OneTo{Int}})
    a.first == 1 || error("Offset axes not supported")
    Ones{Int}(length(a.lasts))
end
function blockfirsts(a::BlockedUnitRange{<:AbstractRange})
    st = step(a.lasts)
    a.first == 1 || error("Offset axes not supported")
    @assert first(a.lasts)-a.first+1 == st
    range(1; step=st, length=length(a.lasts))
end
function blocklengths(a::BlockedUnitRange{<:AbstractRange})
    st = step(a.lasts)
    a.first == 1 || error("Offset axes not supported")
    @assert first(a.lasts)-a.first+1 == st
    Fill(st,length(a.lasts))
end

# TODO: Remove

function _last end
function _length end