text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/InterfaceMechanics.jl/blob/master/LICENSE
module InterfaceMechanics

using LinearAlgebra, Parameters, StaticArrays, ForwardDiff

abstract type AbstractInterface end
abstract type AbstractInterfaceState end

@generated function Base.:+(state::T, dstate::T) where {T <: AbstractInterfaceState}
   expr = [:(state.$p+ dstate.$p) for p in fieldnames(T)]
   return :(T($(expr...)))
end

function update!(interface::I) where {I <: AbstractInterface}
    interface.drivers += interface.ddrivers
    interface.parameters += interface.dparameters
    interface.variables = interface.variables_new
    reset!(interface)
end

function reset!(interface::I) where {I <: AbstractInterface}
    interface.ddrivers = typeof(interface.ddrivers)()
    interface.dparameters = typeof(interface.dparameters)()
    interface.variables_new = typeof(interface.variables_new)()
end

function integrate_interface!(interface::I) where {I<:AbstractInterface}
    error("One needs to define how to integrate interface $I!")
end

export integrate_interface!

export update!, reset!

include("couloumb.jl")
export Couloumb, CouloumbDriverState, CouloumbParameterState, CouloumbVariableState
# export hello, domath
#
# hello(who::String) = "Hello, $who"
# domath(x::Number) = x + 5
#
# println("ladataan paketti")
# greet() = print("Hello World!")

end # module


text:
@doc raw"""
    GradientProblem <: Problem
specify a problem for gradient based algorithms.

# Fields
* `M`            – a manifold $\mathcal M$
* `cost` – a function $F\colon\mathcal M\to\mathbb R$ to minimize
* `gradient`     – the gradient $\nabla F\colon\mathcal M
  \to \mathcal T\mathcal M$ of the cost function $F$

# See also
[`gradient_descent`](@ref)
[`GradientDescentOptions`](@ref)

# """
struct GradientProblem{mT<:Manifold,TCost,TGradient} <: Problem
    M::mT
    cost::TCost
    gradient::TGradient
end
"""
    get_gradient(p,x)

evaluate the gradient of a [`GradientProblem`](@ref)`p` at the point `x`.
"""
function get_gradient(p::GradientProblem, x)
    return p.gradient(x)
end
#
# Options
#
"""
    GradientDescentOptions{P,T} <: Options

Describes a Gradient based descent algorithm, with

# Fields
a default value is given in brackets if a parameter can be left out in initialization.

* `x0` – an a point (of type `P`) on a manifold as starting point
* `stopping_criterion` – ([`StopAfterIteration`](@ref)`(100)`) a [`StoppingCriterion`](@ref)
* `stepsize` – ([`ConstantStepsize`](@ref)`(1.)`)a [`Stepsize`](@ref)
* `retraction_method` – (`ExponentialRetraction()`) the rectraction to use, defaults to
  the exponential map

# Constructor

    GradientDescentOptions(x, stop, s [, retr=ExponentialRetraction()])

construct a Gradient Descent Option with the fields and defaults as above

# See also
[`gradient_descent`](@ref), [`GradientProblem`](@ref)
"""
mutable struct GradientDescentOptions{
    P,TStop<:StoppingCriterion,TStepsize<:Stepsize,TRTM<:AbstractRetractionMethod
} <: Options
    x::P
    stop::TStop
    stepsize::TStepsize
    ∇::P
    retraction_method::TRTM
    function GradientDescentOptions{P}(
        initialX::P,
        s::StoppingCriterion=StopAfterIteration(100),
        stepsize::Stepsize=ConstantStepsize(1.0),
        retraction::AbstractRetractionMethod=ExponentialRetraction(),
    ) where {P}
        o = new{P,typeof(s),typeof(stepsize),typeof(retraction)}()
        o.x = initialX
        o.stop = s
        o.retraction_method = retraction
        o.stepsize = stepsize
        return o
    end
end
function GradientDescentOptions(
    x::P,
    stop::StoppingCriterion=StopAfterIteration(100),
    s::Stepsize=ConstantStepsize(1.0),
    retraction::AbstractRetractionMethod=ExponentialRetraction(),
) where {P}
    return GradientDescentOptions{P}(x, stop, s, retraction)
end
#
# Conjugate Gradient Descent
#
"""
    DirectionUpdateRule

A general functor, that handles direction update rules. It's field(s) is usually
only a [`StoreOptionsAction`](@ref) by default initialized to the fields required
for the specific coefficient, but can also be replaced by a (common, global)
individual one that provides these values.
"""
abstract type DirectionUpdateRule end

@doc raw"""
    ConjugateGradientOptions <: Options

specify options for a conjugate gradient descent algoritm, that solves a
[`GradientProblem`].

# Fields
* `x` – the current iterate, a point on a manifold
* `∇` – the current gradient
* `δ` – the current descent direction, i.e. also tangent vector
* `β` – the current update coefficient rule, see .
* `coefficient` – a [`DirectionUpdateRule`](@ref) function to determine the new `β`
* `stepsize` – a [`Stepsize`](@ref) function
* `stop` – a [`StoppingCriterion`](@ref)
* `retraction` – a function to perform a step on the manifold

# See also
[`conjugate_gradient_descent`](@ref), [`GradientProblem`](@ref), [`ArmijoLinesearch`](@ref)
"""
mutable struct ConjugateGradientDescentOptions{
    T,
    TCoeff<:DirectionUpdateRule,
    TStepsize<:Stepsize,
    TStop<:StoppingCriterion,
    TRetr<:AbstractRetractionMethod,
    TVTM<:AbstractVectorTransportMethod,
} <: Options
    x::T
    ∇::T
    δ::T
    β::Float64
    coefficient::TCoeff
    stepsize::TStepsize
    stop::TStop
    retraction_method::TRetr
    vector_transport_method::TVTM
    function ConjugateGradientDescentOptions{T}(
        x0::T,
        sC::StoppingCriterion,
        s::Stepsize,
        dC::DirectionUpdateRule,
        retr::AbstractRetractionMethod=ExponentialRetraction(),
        vtr::AbstractVectorTransportMethod=ParallelTransport(),
    ) where {T}
        o = new{T,typeof(dC),typeof(s),typeof(sC),typeof(retr),typeof(vtr)}()
        o.x = x0
        o.stop = sC
        o.retraction_method = retr
        o.stepsize = s
        o.coefficient = dC
        o.vector_transport_method = vtr
        return o
    end
end
function ConjugateGradientDescentOptions(
    x::T,
    sC::StoppingCriterion,
    s::Stepsize,
    dU::DirectionUpdateRule,
    retr::AbstractRetractionMethod=ExponentialRetraction(),
    vtr::AbstractVectorTransportMethod=ParallelTransport(),
) where {T}
    return ConjugateGradientDescentOptions{T}(x, sC, s, dU, retr, vtr)
end

@doc raw"""
    ConjugateDescentCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the
variables with prequel `Old` based on [^Flethcer1987] adapted to manifolds:

$\beta_k =
\frac{ \lVert \xi_{k+1} \rVert_{x_{k+1}}^2 }
{\langle -\delta_k,\xi_k \rangle_{x_k}}.$

See also [`conjugate_gradient_descent`](@ref)

# Constructor
    ConjugateDescentCoefficient(a::StoreOptionsAction=())

Construct the conjugate descnt coefficient update rule, a new storage is created by default.

[^Flethcer1987]:
    > R. Fletcher, __Practical Methods of Optimization vol. 1: Unconstrained Optimization__
    > John Wiley & Sons, New York, 1987. doi [10.1137/1024028](https://doi.org/10.1137/1024028)
"""
mutable struct ConjugateDescentCoefficient <: DirectionUpdateRule
    storage::StoreOptionsAction
    function ConjugateDescentCoefficient(a::StoreOptionsAction=StoreOptionsAction((:x, :∇)))
        return new(a)
    end
end
function (u::ConjugateDescentCoefficient)(
    p::GradientProblem, o::ConjugateGradientDescentOptions, i
)
    if !all(has_storage.(Ref(u.storage), [:x, :∇]))
        update_storage!(u.storage, o) # if not given store current as old
        return 0.0
    end
    xOld, ∇Old = get_storage.(Ref(u.storage), [:x, :∇])
    update_storage!(u.storage, o)
    return inner(p.M, o.x, o.∇, o.∇) / inner(p.M, xOld, -o.δ, ∇Old)
end

@doc raw"""
    DaiYuanCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the
variables with prequel `Old` based on [^DaiYuan1999]

adapted to manifolds: let $\nu_k = \xi_{k+1} - P_{x_{k+1}\gets x_k}\xi_k$,
where $P_{a\gets b}(\cdot)$ denotes a vector transport from the tangent space at $a$ to $b$.

Then the coefficient reads

````math
\beta_k =
\frac{ \lVert \xi_{k+1} \rVert_{x_{k+1}}^2 }
{\langle P_{x_{k+1}\gets x_k}\delta_k, \nu_k \rangle_{x_{k+1}}}.
````

See also [`conjugate_gradient_descent`](@ref)

# Constructor
    DaiYuanCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=(),
    )

Construct the Dai Yuan coefficient update rule, where the parallel transport is the
default vector transport and a new storage is created by default.

[^DaiYuan1999]:
    > [Y. H. Dai and Y. Yuan, A nonlinear conjugate gradient method with a strong global convergence property,
    > SIAM J. Optim., 10 (1999), pp. 177–182.
    > doi: [10.1137/S1052623497318992](https://doi.org/10.1137/S1052623497318992)
"""
mutable struct DaiYuanCoefficient{TVTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    transport_method::TVTM
    storage::StoreOptionsAction
    function DaiYuanCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=StoreOptionsAction((:x, :∇, :δ)),
    )
        return new{typeof(t)}(t, a)
    end
end
function (u::DaiYuanCoefficient)(p::GradientProblem, o::ConjugateGradientDescentOptions, i)
    if !all(has_storage.(Ref(u.storage), [:x, :∇, :δ]))
        update_storage!(u.storage, o) # if not given store current as old
        return 0.0
    end
    xOld, ∇Old, δOld = get_storage.(Ref(u.storage), [:x, :∇, :δ])
    update_storage!(u.storage, o)

    ∇tr = vector_transport_to(p.M, xOld, ∇Old, o.x, u.transport_method)
    ν = o.∇ - ∇tr #notation y from [HZ06]
    δtr = vector_transport_to(p.M, xOld, δOld, o.x, u.transport_method)
    return inner(p.M, o.x, o.∇, o.∇) / inner(p.M, xOld, δtr, ν)
end

@doc raw"""
    FletcherReevesCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the
variables with prequel `Old` based on [^FletcherReeves1964] adapted to manifolds:

````math
\beta_k =
\frac{\lVert \xi_{k+1}\rVert_{x_{k+1}}^2}{\lVert \xi_{k}\rVert_{x_{k}}^2}.
````

See also [`conjugate_gradient_descent`](@ref)

# Constructor
    FletcherReevesCoefficient(a::StoreOptionsAction=())

Construct the Fletcher Reeves coefficient update rule, a new storage is created by default.

[^FletcherReeves1964]:
    > R. Fletcher and C. Reeves, Function minimization by conjugate gradients,
    > Comput. J., 7 (1964), pp. 149–154.
    > doi: [10.1093/comjnl/7.2.149](http://dx.doi.org/10.1093/comjnl/7.2.149)
"""
mutable struct FletcherReevesCoefficient <: DirectionUpdateRule
    storage::StoreOptionsAction
    FletcherReevesCoefficient(a::StoreOptionsAction=StoreOptionsAction((:x, :∇))) = new(a)
end
function (u::FletcherReevesCoefficient)(
    p::GradientProblem, o::ConjugateGradientDescentOptions, i
)
    if !all(has_storage.(Ref(u.storage), [:x, :∇]))
        update_storage!(u.storage, o) # if not given store current as old
    end
    xOld, ∇Old = get_storage.(Ref(u.storage), [:x, :∇])
    update_storage!(u.storage, o)
    return inner(p.M, o.x, o.∇, o.∇) / inner(p.M, xOld, ∇Old, ∇Old)
end

@doc raw"""
    HagerZhangCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the variables with
prequel `Old` based on [^HagerZhang2005]
adapted to manifolds: let $\nu_k = \xi_{k+1} - P_{x_{k+1}\gets x_k}\xi_k$,
where $P_{a\gets b}(\cdot)$ denotes a vector transport from the tangent space at $a$ to $b$.

````math
\beta_k = \Bigl\langle\nu_k -
\frac{ 2\lVert \nu_k\rVert_{x_{k+1}}^2 }{ \langle P_{x_{k+1}\gets x_k}\delta_k, \nu_k \rangle_{x_{k+1}} }
P_{x_{k+1}\gets x_k}\delta_k,
\frac{\xi_{k+1}}{ \langle P_{x_{k+1}\gets x_k}\delta_k, \nu_k \rangle_{x_{k+1}} }
\Bigr\rangle_{x_{k+1}}.
````

This method includes a numerical stability proposed by those authors.

See also [`conjugate_gradient_descent`](@ref)

# Constructor
    HagerZhangCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=(),
    )

Construct the Hager Zhang coefficient update rule, where the parallel transport is the
default vector transport and a new storage is created by default.

[^HagerZhang2005]:
    > [W. W. Hager and H. Zhang, __A new conjugate gradient method with guaranteed descent and an efficient line search__,
    > SIAM J. Optim, (16), pp. 170-192, 2005.
    > doi: [10.1137/030601880](https://doi.org/10.1137/030601880)
"""
mutable struct HagerZhangCoefficient{TVTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    transport_method::TVTM
    storage::StoreOptionsAction
    function HagerZhangCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=StoreOptionsAction((:x, :∇, :δ)),
    )
        return new{typeof(t)}(t, a)
    end
end
function (u::HagerZhangCoefficient)(
    p::P, o::O, i
) where {P<:GradientProblem,O<:ConjugateGradientDescentOptions}
    if !all(has_storage.(Ref(u.storage), [:x, :∇, :δ]))
        update_storage!(u.storage, o) # if not given store current as old
        return 0.0
    end
    xOld, ∇Old, δOld = get_storage.(Ref(u.storage), [:x, :∇, :δ])
    update_storage!(u.storage, o)

    ∇tr = vector_transport_to(p.M, xOld, ∇Old, o.x, u.transport_method)
    ν = o.∇ - ∇tr #notation y from [HZ06]
    δtr = vector_transport_to(p.M, xOld, δOld, o.x, u.transport_method)
    denom = inner(p.M, o.x, δtr, ν)
    νknormsq = inner(p.M, o.x, ν, ν)
    β = inner(p.M, o.x, ν, o.∇) / denom - 2 * νknormsq * inner(p.M, o.x, δtr, o.∇) / denom^2
    # Numerical stability from Manopt / Hager-Zhang paper
    ξn = norm(p.M, o.x, o.∇)
    η = -1 / (ξn * min(0.01, norm(p.M, xOld, ∇Old)))
    return max(β, η)
end

@doc raw"""
    HeestenesStiefelCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the
variables with prequel `Old` based on [^HeestensStiefel1952]

adapted to manifolds as follows: let $\nu_k = \xi_{k+1} - P_{x_{k+1}\gets x_k}\xi_k$.
Then the update reads

````math
\beta_k = \frac{\langle \xi_{k+1}, \nu_k \rangle_{x_{k+1}} }
    { \langle P_{x_{k+1}\gets x_k} \delta_k, \nu_k\rangle_{x_{k+1}} },
````
where $P_{a\gets b}(\cdot)$ denotes a vector transport from the tangent space at $a$ to $b$.

# Constructor
    HeestenesStiefelCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=()
    )

Construct the Heestens Stiefel coefficient update rule, where the parallel transport is the
default vector transport and a new storage is created by default.

See also [`conjugate_gradient_descent`](@ref)

[^HeestensStiefel1952]:
    > M.R. Hestenes, E.L. Stiefel, Methods of conjugate gradients for solving linear systems,
    > J. Research Nat. Bur. Standards, 49 (1952), pp. 409–436.
    > doi: [10.6028/jres.049.044](http://dx.doi.org/10.6028/jres.049.044)
"""
mutable struct HeestenesStiefelCoefficient{TVTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    transport_method::TVTM
    storage::StoreOptionsAction
    function HeestenesStiefelCoefficient(
        transport_method::AbstractVectorTransportMethod=ParallelTransport(),
        storage_action::StoreOptionsAction=StoreOptionsAction((:x, :∇, :δ)),
    )
        return new{typeof(transport_method)}(transport_method, storage_action)
    end
end
function (u::HeestenesStiefelCoefficient)(
    p::GradientProblem, o::ConjugateGradientDescentOptions, i
)
    if !all(has_storage.(Ref(u.storage), [:x, :∇, :δ]))
        update_storage!(u.storage, o) # if not given store current as old
        return 0.0
    end
    xOld, ∇Old, δOld = get_storage.(Ref(u.storage), [:x, :∇, :δ])
    update_storage!(u.storage, o)
    ∇tr = vector_transport_to(p.M, xOld, ∇Old, o.x, u.transport_method)
    δtr = vector_transport_to(p.M, xOld, δOld, o.x, u.transport_method)
    ν = o.∇ - ∇tr #notation from [HZ06]
    β = inner(p.M, o.x, o.∇, ν) / inner(p.M, o.x, δtr, ν)
    return max(0, β)
end

@doc raw"""
    LiuStoreyCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the
variables with prequel `Old` based on [^LuiStorey1991]
adapted to manifolds: let $\nu_k = \xi_{k+1} - P_{x_{k+1}\gets x_k}\xi_k$,
where $P_{a\gets b}(\cdot)$ denotes a vector transport from the tangent space at $a$ to $b$.

Then the coefficient reads

````math
\beta_k = -
\frac{ \langle \xi_{k+1},\nu_k \rangle_{x_{k+1}} }
{\langle \delta_k,\xi_k \rangle_{x_k}}.
````

See also [`conjugate_gradient_descent`](@ref)

# Constructor
    LiuStoreyCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=()
    )

Construct the Lui Storey coefficient update rule, where the parallel transport is the
default vector transport and a new storage is created by default.

[^LuiStorey1991]:
    > [Y. Liu and C. Storey, Efficient generalized conjugate gradient algorithms, Part 1: Theory
    > J. Optim. Theory Appl., 69 (1991), pp. 129–137.
    > doi: [10.1007/BF00940464](https://doi.org/10.1007/BF00940464)
"""
mutable struct LiuStoreyCoefficient{TVTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    transport_method::TVTM
    storage::StoreOptionsAction
    function LiuStoreyCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=StoreOptionsAction((:x, :∇, :δ)),
    )
        return new{typeof(t)}(t, a)
    end
end
function (u::LiuStoreyCoefficient)(
    p::GradientProblem, o::ConjugateGradientDescentOptions, i
)
    if !all(has_storage.(Ref(u.storage), [:x, :∇, :δ]))
        update_storage!(u.storage, o) # if not given store current as old
    end
    xOld, ∇Old, δOld = get_storage.(Ref(u.storage), [:x, :∇, :δ])
    update_storage!(u.storage, o)
    ∇tr = vector_transport_to(p.M, xOld, ∇Old, o.x, u.transport_method)
    ν = o.∇ - ∇tr # notation y from [HZ06]
    return inner(p.M, o.x, o.∇, ν) / inner(p.M, xOld, -δOld, ∇Old)
end

@doc raw"""
    PolakRibiereCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the
variables with prequel `Old` based on [^PolakRibiere1969][^Polyak1969]

adapted to manifolds: let $\nu_k = \xi_{k+1} - P_{x_{k+1}\gets x_k}\xi_k$,
where $P_{a\gets b}(\cdot)$ denotes a vector transport from the tangent space at $a$ to $b$.

Then the update reads
````math
\beta_k =
\frac{ \langle \xi_{k+1}, \nu_k \rangle_{x_{k+1}} }
{\lVert \xi_k \rVert_{x_k}^2 }.
````

# Constructor
    PolakRibiereCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=()
    )

Construct the PolakRibiere coefficient update rule, where the parallel transport is the
default vector transport and a new storage is created by default.

See also [`conjugate_gradient_descent`](@ref)

[^PolakRibiere1969]:
    > E. Polak, G. Ribiere, Note sur la convergence de méthodes de directions conjuguées
    > ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique, Tome 3 (1969) no. R1, p. 35-43,
    > url: [http://www.numdam.org/item/?id=M2AN_1969__3_1_35_0](http://www.numdam.org/item/?id=M2AN_1969__3_1_35_0)

[^Polyak1969]:
    > B. T. Polyak, The conjugate gradient method in extreme problems,
    > USSR Comp. Math. Math. Phys., 9 (1969), pp. 94–112.
    > doi: [10.1016/0041-5553(69)90035-4](https://doi.org/10.1016/0041-5553(69)90035-4)
"""
mutable struct PolakRibiereCoefficient{TVTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    transport_method::TVTM
    storage::StoreOptionsAction
    function PolakRibiereCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=StoreOptionsAction((:x, :∇)),
    )
        return new{typeof(t)}(t, a)
    end
end
function (u::PolakRibiereCoefficient)(
    p::GradientProblem, o::ConjugateGradientDescentOptions, i
)
    if !all(has_storage.(Ref(u.storage), [:x, :∇]))
        update_storage!(u.storage, o) # if not given store current as old
    end
    xOld, ∇Old = get_storage.(Ref(u.storage), [:x, :∇])
    update_storage!(u.storage, o)

    ∇tr = vector_transport_to(p.M, xOld, ∇Old, o.x, u.transport_method)
    ν = o.∇ - ∇tr
    β = inner(p.M, o.x, o.∇, ν) / inner(p.M, xOld, ∇Old, ∇Old)
    return max(0, β)
end

@doc raw"""
    SteepestDirectionUpdateRule <: DirectionUpdateRule

The simplest rule to update is to have no influence of the last direction and
hence return an update $\beta = 0$ for all [`ConjugateGradientDescentOptions`](@ref)` o`

See also [`conjugate_gradient_descent`](@ref)
"""
struct SteepestDirectionUpdateRule <: DirectionUpdateRule end
function (u::SteepestDirectionUpdateRule)(
    p::GradientProblem, o::ConjugateGradientDescentOptions, i
)
    return 0.0
end

@doc raw"""
    DebugGradient <: DebugAction

debug for the gradient evaluated at the current iterate

# Constructors
    DebugGradient([long=false,p=print])

display the short (`false`) or long (`true`) default text for the gradient.

    DebugGradient(prefix[, p=print])

display the a `prefix` in front of the gradient.
"""
mutable struct DebugGradient <: DebugAction
    io::IO
    prefix::String
    function DebugGradient(long::Bool=false, io::IO=stdout)
        return new(io, long ? "Gradient: " : "∇F(x):")
    end
    DebugGradient(prefix::String, io::IO=stdout) = new(io, prefix)
end
function (d::DebugGradient)(::GradientProblem, o::GradientDescentOptions, i::Int)
    print(d.io, (i >= 0) ? d.prefix * "" * string(o.∇) : "")
    return nothing
end

@doc raw"""
    DebugGradientNorm <: DebugAction

debug for gradient evaluated at the current iterate.

# Constructors
    DebugGradientNorm([long=false,p=print])

display the short (`false`) or long (`true`) default text for the gradient norm.

    DebugGradientNorm(prefix[, p=print])

display the a `prefix` in front of the gradientnorm.
"""
mutable struct DebugGradientNorm <: DebugAction
    io::IO
    prefix::String
    function DebugGradientNorm(long::Bool=false, io::IO=stdout)
        return new(io, long ? "Norm of the Gradient: " : "|∇F(x)|:")
    end
    DebugGradientNorm(prefix::String, io::IO=stdout) = new(io, prefix)
end
function (d::DebugGradientNorm)(p::GradientProblem, o::GradientDescentOptions, i::Int)
    print(d.io, (i >= 0) ? d.prefix * "$(norm(p.M,o.x,o.∇))" : "")
    return nothing
end

@doc raw"""
    DebugStepsize <: DebugAction

debug for the current step size.

# Constructors
    DebugStepsize([long=false,p=print])

display the short (`false`) or long (`true`) default text for the step size.

    DebugStepsize(prefix[, p=print])

display the a `prefix` in front of the step size.
"""
mutable struct DebugStepsize <: DebugAction
    io::IO
    prefix::String
    function DebugStepsize(long::Bool=false, io::IO=stdout)
        return new(io, long ? "step size:" : "s:")
    end
    DebugStepsize(prefix::String, io::IO=stdout) = new(io, prefix)
end
function (d::DebugStepsize)(p::GradientProblem, o::GradientDescentOptions, i::Int)
    print(d.io, (i > 0) ? d.prefix * "$(get_last_stepsize(p,o,i))" : "")
    return nothing
end

#
# Records
#
@doc raw"""
    RecordGradient <: RecordAction

record the gradient evaluated at the current iterate

# Constructors
    RecordGradient(ξ)

initialize the [`RecordAction`](@ref) to the corresponding type of the tangent vector.
"""
mutable struct RecordGradient{T} <: RecordAction
    recorded_values::Array{T,1}
    RecordGradient{T}() where {T} = new(Array{T,1}())
end
RecordGradient(ξ::T) where {T} = RecordGradient{T}()
function (r::RecordGradient{T})(
    ::GradientProblem, o::GradientDescentOptions, i::Int
) where {T}
    return record_or_reset!(r, o.∇, i)
end

@doc raw"""
    RecordGradientNorm <: RecordAction

record the norm of the current gradient
"""
mutable struct RecordGradientNorm <: RecordAction
    recorded_values::Array{Float64,1}
    RecordGradientNorm() = new(Array{Float64,1}())
end
function (r::RecordGradientNorm)(
    p::P, o::O, i::Int
) where {P<:GradientProblem,O<:GradientDescentOptions}
    return record_or_reset!(r, norm(p.M, o.x, o.∇), i)
end

@doc raw"""
    RecordStepsize <: RecordAction

record the step size
"""
mutable struct RecordStepsize <: RecordAction
    recorded_values::Array{Float64,1}
    RecordStepsize() = new(Array{Float64,1}())
end
function (r::RecordStepsize)(
    p::P, o::O, i::Int
) where {P<:GradientProblem,O<:GradientDescentOptions}
    return record_or_reset!(r, get_last_stepsize(p, o, i), i)
end

ext: jl


text:
function test_apply_rhs()
    grid = generate_grid(Quadrilateral, (20, 20))
    dim = 2
    ip = Lagrange{dim,RefCube,1}()
    qr = QuadratureRule{dim,RefCube}(2)
    cellvalues = CellScalarValues(qr, ip)

    dh = DofHandler(grid)
    push!(dh, :u, 1)
    close!(dh)

    K = create_sparsity_pattern(dh)

    ch = ConstraintHandler(dh)

    ∂Ω = union(getfaceset.((grid,), ["left", "right"])...)
    dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0)
    add!(ch, dbc);

    ∂Ω = union(getfaceset.((grid,), ["top", "bottom"])...)
    dbc = Dirichlet(:u, ∂Ω, (x, t) -> 2)
    add!(ch, dbc);

    close!(ch)
    update!(ch, 0.0);

    function doassemble!(
        cellvalues::CellScalarValues{dim},
        K::SparseMatrixCSC,
        dh::DofHandler,
    ) where {dim}

        n_basefuncs = getnbasefunctions(cellvalues)
        Ke = zeros(n_basefuncs, n_basefuncs)
        fe = zeros(n_basefuncs)

        f = zeros(ndofs(dh))
        assembler = start_assemble(K, f)

        @inbounds for cell in CellIterator(dh)
            fill!(Ke, 0)
            fill!(fe, 0)

            reinit!(cellvalues, cell)

            for q_point = 1:getnquadpoints(cellvalues)
                dΩ = getdetJdV(cellvalues, q_point)

                for i = 1:n_basefuncs
                    v = shape_value(cellvalues, q_point, i)
                    ∇v = shape_gradient(cellvalues, q_point, i)
                    fe[i] += v * dΩ
                    for j = 1:n_basefuncs
                        ∇u = shape_gradient(cellvalues, q_point, j)
                        Ke[i, j] += (∇v ⋅ ∇u) * dΩ
                    end
                end
            end

            assemble!(assembler, celldofs(cell), fe, Ke)
        end
        return K, f
    end

    K, f = doassemble!(cellvalues, K, dh)
    A = create_sparsity_pattern(dh)
    A, g = doassemble!(cellvalues, A, dh)
    rhsdata = get_rhs_data(ch, A)

    apply!(K, f, ch)
    apply!(A, ch) # need to apply bcs to A once
    apply_rhs!(rhsdata, g, ch)
    u₁ = K \ f
    u₂ = A \ g
    return u₁, u₂
end

u1, u2 = test_apply_rhs()
@test u1 == u2


text:
module StatsdClient
# A simple statsd client written in Julia
# Usage:
# > using StatsdClient
# > server = StatsdClient.Statsd()
# Statsd Server: 127.0.0.1 Port: 8125
# You can also specify the server address and port number
# (port defaults to 8125)
# > sdc = StatsdClient.Statsd("0.0.0.0",8125)
# Statsd Server: 0.0.0.0 Port: 8125
# > StatsdClient.increment(sdc,"test.me")
# Please refer to the readme for more examples.
# Note: Requires Julia 0.3 with commit sha 6585e3de1b or later.

export increment, decrement, count, gauge, timing

type Statsd
    server_address::IPv4
    server_port::Int
    send_msg::Function
    Statsd(ipv4::IPv4,port) = new(ipv4,port,_make_send(ipv4,port))
end
Statsd() = Statsd(IPv4(127,0,0,1),8125)
Statsd(ip::String,port=8125) = Statsd(getaddrinfo(ip),port)

import Base.show
show(io::IO,client::Statsd) = print(io,string("Statsd Server: ",
                                              client.server_address,
                                              " Port: ",
                                              client.server_port))

function _make_send(ip,port)
    sock = UdpSocket()
    Base.bind(sock,ip,0)
    Base.setopt(sock,enable_broadcast=1)
    function _send(data,sample_rate::Union(Number,Nothing)=nothing)
        if is(sample_rate,nothing)
            send(sock,ip,port,data)
        elseif (0 <= sample_rate <= 1) &&
               rand() < sample_rate
            send(sock,ip,port,string(data,"|@",sample_rate))
        end
    end
end

increment(cl::Statsd,metric,sample_rate=nothing) = count(cl,metric,1,sample_rate)

decrement(cl::Statsd,metric,sample_rate=nothing) = count(cl,metric,-1,sample_rate)

count(cl::Statsd,metric,value,
      sample_rate=nothing) = cl.send_msg(string(metric,":",value,"|c"),sample_rate)

timing(cl::Statsd,metric,value,
       sample_rate=nothing) = cl.send_msg(string(metric,":",value,"|ms"),sample_rate)

gauge(cl::Statsd,metric,value,
      sample_rate=nothing) = cl.send_msg(string(metric,":",value,"|g"),sample_rate)

set(cl::Statsd,metric,value,
    sample_rate=nothing) = cl.send_msg(string(metric,":",value,"|s"),sample_rate)

# versions of the above that accept multiple metrics at once
for f in [:count, :timing, :gauge, :set]
    @eval begin
        function $(f)(cl::Statsd, metrics::Union(AbstractArray, Tuple), value,
                      sample_rate=nothing)
            for m in metrics
                $(f)(cl, m, value, sample_rate)
            end
        end
    end
end

end


text:
using Interact,
      Plots,
      PotentialCalculation

"""
    min_distance(mpp::MoleculePairPotential, points)

Gives minimum distance of molecules in `mpp` on given points.
Used to help plotting.
"""
function min_distance(mpp::MoleculePairPotential, points)
    l1 = length(mpp.mol1)
    l2 = length(mpp.mol2)
    return map(x-> minimum(distances(x[1:l1],x[l1+1:l1+l2])), points)
end

"""
    min_distance(points, c1::AbstractCluster, c2::AbstractCluster)

Minimum without needing fitted potential.
"""
function min_distance(points, c1::AbstractCluster, c2::AbstractCluster)
    l1 = length(c1)
    l2 = length(c2)
    return map(x-> minimum(distances(x[1:l1],x[l1+1:l1+l2])), points)
end

"""
    plot_potential(points, mpp::MoleculePairPotential; emax=100, unit="cm^-1",
                        leg=false, figsize=(800,400), font=font(20))

Plots potential

# Arguments
- `points`  : array of [`Cluster`](@ref) where potential is plotted
- `mpp::MoleculePairPotential`  : potential

# Keywords
- `emax=100`  : maximum energy in plot - cut all values with energy greater
- `unit="cm^-1"`        : energy unit
- `leg=false`           : draw legend
- `figsize=(800,400)`   : size of picture
- `font=font(20)`       : font size
"""
function plot_potential(points, mpp::MoleculePairPotential; emax=100, unit="cm^-1",
                        leg=false, figsize=(800,400), font=font(20))
    E = energy_to.(mpp.(points), unit)
    r = min_distance(mpp, points)
    i = E .< emax
    plot(r[i], E[i],
         xlabel = "Distance   [Å]",
         ylabel = "Energy   ["*unit*"]",
         label = "Fitted", size=figsize, tickfont=font)
end

"""
    plot_potential(pes::Dict; emax=100, unit="cm^-1", figsize=(800,400), font=font(20))

Plots calculated potential

# Arguments
- `pes::Dict`     : PES in Dict from that can be get with `load_data_file` command

# Keywords
- `emax=100`      : maximum energy in plot - cut all values with energy greater
- `unit="cm^-1"`        : energy unit
- `figsize=(800,400)`   : size of picture
- `font=font(20)`       : font size
"""
function plot_potential(pes::Dict; emax=100, unit="cm^-1", figsize=(800,400), font=font(20))
    @assert haskey(pes, "Energy")
    @assert haskey(pes, "Points")
    @assert haskey(pes, "cluster1")
    @assert haskey(pes, "cluster2")
    E = energy_to.(pes["Energy"], unit)
    r = min_distance(pes["Points"], pes["cluster1"], pes["cluster2"])
    s = size(E, 2)
    plt = @manipulate for col in slider(1:s, label="Collumn")
        i = E[:,col] .< emax
        plot(r[:,col][i], E[:,col][i];
              size=figsize,
              font=font,
              leg=false,
              seriestype=:scatter,
              xlabel = "Distance   [Å]",
              ylabel = "Energy   [$unit]")
    end
    plt
end

function plot_potential(pes::Dict, i::Int; emax=100, unit="cm^-1", figsize=(800,400), font=font(20))
    @assert haskey(pes, "Energy")
    @assert haskey(pes, "Points")
    @assert haskey(pes, "cluster1")
    @assert haskey(pes, "cluster2")
    E = energy_to.(pes["Energy"][:,i], unit)
    r = min_distance(pes["Points"][:,i], pes["cluster1"], pes["cluster2"])
    j = E .< emax
    plt=plot(r[j], E[j];
          size=figsize,
          font=font,
          leg=false,
          seriestype=:scatter,
          xlabel = "Distance   [Å]",
          ylabel = "Energy   [$unit]")

    plt
end

"""
    plot_compare(points, energy, mpp::MoleculePairPotential...; emax=100, unit="cm^-1",
                      leg=false, size=(800,400), font=font(20))

Compares fitted energy to calculated one.

# Arguments
- `points`  : array of [`Cluster`](@ref) where potential is plotted
- `energy`  : calculated energy
- `mpp::MoleculePairPotential...`  : potentials to be fitted

# Keywords
- `emax=100`  : maximum energy in plot - cut all values with energy greater
- `unit="cm^-1"`  : energy unit
- `leg=false`     : draw legend
- `size=(800,400)`   : size of picture
- `font=font(20)`    : font size
"""
function plot_compare(points, energy, mpp::MoleculePairPotential...; emax=100, unit="cm^-1",
                      leg=false, figsize=(800,400), font=font(20))
    Ecal = energy_to.(energy, unit)
    r = min_distance(mpp[1], points)
    i = Ecal .< emax
    out=plot(r[i], Ecal[i], xlabel="Min Distance  [Å]", ylabel="Energy   [$unit]",  label="Calculated",
         leg=leg, size=figsize, tickfont=font, seriestype=:scatter)

    for x in mpp
        Efit = energy_to.(x.(points), unit)
        plot!(r[i], Efit[i], label="Fit", tickfont=font)
    end
    out
end

"""
visualize_points(points; stdout=devnull, stderr=devnull, command="vmd")

Visualize using external program.
"""
function visualize_points(points; stdout=devnull, stderr=devnull, command="vmd")
    fname = tempname()*".xyz"
    open(fname,"w+") do f
        for x in points
            print_xyz(f, x)
        end
    end
    run(pipeline(`$(command) $(fname)`, stdout=stdout, stderr=stderr))
    rm(fname)
end

"""
scan_compare(points,energy, mppe...; emax=100, unit="cm^-1",
                      leg=false, figsize=(800,400))

Use [`Interact`](https://juliagizmos.github.io/Interact.jl/stable/) to view
potentials ineractively on different points.

Visualization is done on collumn vise.

# Arguments
- `points`  : array of points, first dimension is displayd while second can be chosen
- `energy`  : array of reference energy, first dimension is displayd while second can be chosen
- `mppe...` : [`MoleculePairPotential`](@ref) which are plotted

# Keywords
- `emax=100`  : maximum energy in plot - cut all values with energy greater
- `unit="cm^-1"`  : energy unit
- `leg=false`     : draw legend
- `figsize=(800,400)`   : size of picture
"""
function scan_compare(points,energy, mppe...; emax=100, unit="cm^-1",
                      leg=false, figsize=(800,400))
    @assert size(points) == size(energy) "points and energy need to have same size"
    e = energy_from(emax,unit)
    s=size(points)
    units = Observable(["cm^-1", "kcal/mol", "kJ/mol", "eV", "K", "hartree"])
    wdg = dropdown(units, label="Energy unit")
    display(wdg)

    plt = @manipulate for x in wdg, col in slider(1:s[2], label="Collumn")
        plot_compare(points[:,col],energy[:,col], mppe..., emax=energy_to(e,x), unit=x, leg=leg, figsize=figsize)
    end
    return plt
end


text:
###########################################################################
#
#  Random-Walk Metropolis (RWM)
#
#  Parameters:
#    -scale: for scaling the jumps
#    -tuner: for tuning the scale parameter
#
###########################################################################

export RWM

println("Loading RWM(scale, tuner) sampler")

###########################################################################
#                  RWM specific 'tuners'
###########################################################################
abstract RWMTuner <: MCMCTuner

###########################################################################
#                  RWM type
###########################################################################

immutable RWM <: MCMCSampler
  scale::Float64
  tuner::Union(Nothing, RWMTuner)

  function RWM(x::Float64, t::Union(Nothing, RWMTuner))
    @assert x>0 "scale should be > 0"
    new(x, t)
  end
end
RWM() = RWM(1., nothing)
RWM(x::Float64) = RWM(x, nothing)
RWM(t::RWMTuner) = RWM(1., t)
RWM(;scale::Float64=1.0, tuner::Union(Nothing, RWMTuner)=nothing) = RWM(scale, tuner)

###########################################################################
#                  RWM task
###########################################################################

# RWM sampling
function SamplerTask(model::MCMCModel, sampler::RWM, runner::MCMCRunner)
	local pars, proposedPars
	local logTarget, proposedLogTarget
  local scale

	# hook inside Task to allow remote resetting
	task_local_storage(:reset, (resetPars::Vector{Float64}) -> (pars = copy(resetPars); logTarget = model.eval(pars)))

	# initialization
	scale = model.scale .* sampler.scale  # rescale model scale by sampler scale
	pars = copy(model.init)
	logTarget = model.eval(pars)
  @assert isfinite(logTarget) "Initial values out of model support, try other values"

	# main loop
	while true
		proposedPars = pars + randn(model.size) .* scale
		proposedLogTarget = model.eval(proposedPars)

	  ratio = proposedLogTarget-logTarget
		if ratio > 0 || (ratio > log(rand()))
			ms = MCMCSample(proposedPars, proposedLogTarget, pars, logTarget, {"accept" => true})
	    produce(ms)
	    pars, logTarget = copy(proposedPars), copy(proposedLogTarget)
	  else
			ms = MCMCSample(pars, logTarget, pars, logTarget, {"accept" => false})
      produce(ms)
    end
	end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "TempestModel"
version = v"0.1"
sources = [
    GitSource("https://github.com/paullric/tempestmodel",
	      "4ecf5146f3100fb24d36156a2f99433f049f0f66"),
    DirectorySource("./bundled")
]

script = raw"""
mkdir -vp $bindir
mkdir -vp $libdir

cd ${WORKSPACE}/srcdir/tempestmodel*

atomic_patch -p1 ../patches/mk_defs.patch
atomic_patch -p1 ../patches/mk_system_macosx.patch
# linux system fallback
atomic_patch -p1 ../patches/mk_system_agri.patch
atomic_patch -p1 ../patches/netcdf_cpp.patch

export CXXFLAGS="-I${includedir}"
export LDFLAGS="-L${libdir}"

make NETCDF_ROOT=$prefix -j${nproc} all

cp -v test/nonhydro_sphere/ScharMountainSphereTest ${bindir}
cp -v test/nonhydro_sphere/BaroclinicWaveJWTest ${bindir}
cp -v test/nonhydro_sphere/StationaryMountainFlowTest ${bindir}
cp -v test/nonhydro_sphere/MountainWaveSphereTest ${bindir}
cp -v test/nonhydro_sphere/HeldSuarezTest ${bindir}
cp -v test/nonhydro_sphere/InertiaGravityWaveTest ${bindir}
cp -v test/nonhydro_sphere/BaroclinicWaveUMJSTest ${bindir}
cp -v test/nonhydro_sphere/MountainRossby3DTest ${bindir}
cp -v test/nonhydro_sphere/BaldaufGravityWaveTest ${bindir}

cp -v test/nonhydro_xz/RobertBubbleCartesianTest ${bindir}
cp -v test/nonhydro_xz/ShearJetMtnWave2DCartesianTest ${bindir}
cp -v test/nonhydro_xz/ThermalBubbleCartesian3DTest ${bindir}
cp -v test/nonhydro_xz/Baroclinic3DCartesianTest ${bindir}
cp -v test/nonhydro_xz/DensityCurrentCartesianTest ${bindir}
cp -v test/nonhydro_xz/NonHydroMountainCartesianTest ${bindir}
cp -v test/nonhydro_xz/HydrostaticMountainCartesianTest ${bindir}
cp -v test/nonhydro_xz/ThermalBubbleCartesianTest ${bindir}
cp -v test/nonhydro_xz/Baroclinic3DCartesianRidgeTest ${bindir}
cp -v test/nonhydro_xz/ScharMountainCartesianTest ${bindir}
cp -v test/nonhydro_xz/InertialGravityCartesianXZTest ${bindir}

cp -v test/shallowwater_sphere/BarotropicInstabilityTest ${bindir}
cp -v test/shallowwater_sphere/MountainRossbyTest ${bindir}
cp -v test/shallowwater_sphere/RossbyHaurwitzWaveTest ${bindir}
"""

# Note: We are restricted to the platforms that NetCDF supports, the library is Unix only
platforms = [
    Platform("x86_64", "linux"),
    Platform("x86_64", "macos"),
]
platforms = expand_cxxstring_abis(platforms)

products = [
    # non-hydro sphere test cases
    ExecutableProduct("ScharMountainSphereTest", :ScharMountainSphereTest_exe),
    ExecutableProduct("BaroclinicWaveJWTest", :BaroclinicWaveJWTest_exe),
    ExecutableProduct("StationaryMountainFlowTest", :StationaryMountainFlowTest_exe),
    ExecutableProduct("MountainWaveSphereTest", :MountainWaveSphereTest_exe),
    ExecutableProduct("HeldSuarezTest", :HeldSuarezTest_exe),
    ExecutableProduct("InertiaGravityWaveTest", :InertiaGravityWaveTest_exe),
    ExecutableProduct("BaroclinicWaveUMJSTest", :BaroclinicWaveUMJSTest_exe),
    ExecutableProduct("MountainRossby3DTest", :MountainRossby3DTest_exe),
    ExecutableProduct("BaldaufGravityWaveTest", :BaldaufGravityWaveTest_exe),

    # non-hydro xz test cases
    ExecutableProduct("RobertBubbleCartesianTest", :RobertBubbleCartesianTest_exe),
    ExecutableProduct("ShearJetMtnWave2DCartesianTest", :ShearJetMtnWave2DCartesianTest_exe),
    ExecutableProduct("ThermalBubbleCartesian3DTest", :ThermalBubbleCartesian3DTest_exe),
    ExecutableProduct("Baroclinic3DCartesianTest", :Baroclinic3DCartesianTest_exe),
    ExecutableProduct("DensityCurrentCartesianTest", :DensityCurrentCartesianTest_exe),
    ExecutableProduct("NonHydroMountainCartesianTest", :NonHydroMountainCartesianTest_exe),
    ExecutableProduct("HydrostaticMountainCartesianTest", :HydrostaticMountainCartesianTest_exe),
    ExecutableProduct("ThermalBubbleCartesianTest", :ThermalBubbleCartesianTest_exe),
    ExecutableProduct("Baroclinic3DCartesianRidgeTest", :Baroclinic3DCartesianRidgeTest_exe),
    ExecutableProduct("ScharMountainCartesianTest", :ScharMountainCartesianTest_exe),
    ExecutableProduct("InertialGravityCartesianXZTest", :InertialGravityCartesianXZTest_exe),

    # shallow water test cases
    ExecutableProduct("BarotropicInstabilityTest", :BarotropicInstabilityTest_exe),
    ExecutableProduct("MountainRossbyTest", :MountainRossbyTest_exe),
    ExecutableProduct("RossbyHaurwitzWaveTest", :RossbyHaurwitzWaveTest_exe),
]

dependencies = [
    Dependency("MKL_jll"),
    Dependency("MPICH_jll"),
    Dependency("NetCDF_jll"; compat="400.701.400 - 400.799"),
    Dependency("HDF5_jll"),
    # `MbedTLS_jll` is an indirect dependency through NetCDF, we need to specify
    # a compatible build version for this to work.
    BuildDependency(PackageSpec(; name="MbedTLS_jll", version=v"2.24.0")),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    julia_compat="1.6",
)

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
# Test complete Lanczos factorization
@testset "Complete Lanczos factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr) # tests fail miserably for cgs and mgs
            A = rand(T,(n,n))
            v = rand(T,(n,))
            A = (A+A')
            iter = LanczosIterator(wrapop(A), wrapvec(v), orth)
            fact = initialize(iter)
            while length(fact) < n
                expand!(iter, fact; verbosity = 1)
            end

            V = hcat(unwrapvec.(basis(fact))...)
            H = rayleighquotient(fact)
            @test normres(fact) < 10*n*eps(real(T))
            @test V'*V ≈ I
            @test A*V ≈ V*H

            @inferred initialize!(iter, deepcopy(fact); verbosity = 1)
            states = collect(Iterators.take(iter, n)) # collect tests size and eltype?
            @test rayleighquotient(last(states)) ≈ H
        end
    end
end

# Test complete Arnoldi factorization
@testset "Complete Arnoldi factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs, mgs, cgs2, mgs2, cgsr, mgsr)
            A = rand(T,(n,n))
            v = rand(T,(n,))
            iter = ArnoldiIterator(wrapop(A), wrapvec(v), orth)
            fact = initialize(iter)
            while length(fact) < n
                expand!(iter, fact)
            end

            V = hcat(unwrapvec.(basis(fact))...)
            H = rayleighquotient(fact)
            factor = (orth == cgs || orth == mgs ? 100 : 10)
            @test normres(fact) < factor*n*eps(real(T))
            @test V'*V ≈ I
            @test A*V ≈ V*H

            @inferred initialize!(iter, deepcopy(fact); verbosity = 1)
            states = collect(Iterators.take(iter, n)) # collect tests size and eltype?
            @test rayleighquotient(last(states)) ≈ H
        end
    end
end

# Test complete Arnoldi factorization
@testset "Complete Golub-Kahan-Lanczos factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T,(n,n))
            v = A*rand(T,(n,)) # ensure v is in column space of A
            iter = GKLIterator(wrapop(A), wrapvec(v), orth)
            fact = initialize(iter)
            while length(fact) < n
                expand!(iter, fact)
            end

            U = hcat(unwrapvec.(basis(fact, :U))...)
            V = hcat(unwrapvec.(basis(fact, :V))...)
            B = rayleighquotient(fact)
            @test normres(fact) < 10*n*eps(real(T))
            @test U'*U ≈ I
            @test V'*V ≈ I
            @test A*V ≈ U*B
            @test A'*U ≈ V*B'
        end
    end
end

# Test incomplete Lanczos factorization
@testset "Incomplete Lanczos factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64, Complex{Int})
        @testset for orth in (cgs2, mgs2, cgsr, mgsr) # tests fail miserably for cgs and mgs
            if T == Complex{Int}
                A = rand(-100:100, (N, N)) + im * rand(-100:100, (N, N))
                v = rand(-100:100, (N,))
            else
                A = rand(T,(N,N))
                v = rand(T,(N,))
            end
            A = (A+A')
            iter = @inferred LanczosIterator(wrapop(A), wrapvec(v), orth)
            krylovdim = n
            fact = @inferred initialize(iter)
            while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim
                @inferred expand!(iter, fact)

                V = hcat(unwrapvec.(basis(fact))...)
                H = rayleighquotient(fact)
                r = unwrapvec(residual(fact))
                β = normres(fact)
                e = rayleighextension(fact)
                @test V'*V ≈ I
                @test norm(r) ≈ β
                @test A*V ≈ V*H + r*e'
            end

            fact = @inferred shrink!(fact, div(n,2))
            V = hcat(unwrapvec.(@inferred basis(fact))...)
            H = @inferred rayleighquotient(fact)
            r = unwrapvec(@inferred residual(fact))
            β = @inferred normres(fact)
            e = @inferred rayleighextension(fact)
            @test V'*V ≈ I
            @test norm(r) ≈ β
            @test A*V ≈ V*H + r*e'
        end
    end
end

# Test incomplete Arnoldi factorization
@testset "Incomplete Arnoldi factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64, Complex{Int})
        @testset for orth in (cgs, mgs, cgs2, mgs2, cgsr, mgsr)
            if T == Complex{Int}
                A = rand(-100:100, (N, N)) + im * rand(-100:100, (N, N))
                v = rand(-100:100, (N,))
            else
                A = rand(T,(N,N))
                v = rand(T,(N,))
            end
            iter = @inferred ArnoldiIterator(wrapop(A), wrapvec(v), orth)
            krylovdim = 3*n
            fact = @inferred initialize(iter)
            while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim
                @inferred expand!(iter, fact)

                V = hcat(unwrapvec.(basis(fact))...)
                H = rayleighquotient(fact)
                r = unwrapvec(residual(fact))
                β = normres(fact)
                e = rayleighextension(fact)
                @test V'*V ≈ I
                @test norm(r) ≈ β
                @test A*V ≈ V*H + r*e'
            end

            fact = @inferred shrink!(fact, div(n,2))
            V = hcat(unwrapvec.(@inferred basis(fact))...)
            H = @inferred rayleighquotient(fact)
            r = unwrapvec(@inferred residual(fact))
            β = @inferred normres(fact)
            e = @inferred rayleighextension(fact)
            @test V'*V ≈ I
            @test norm(r) ≈ β
            @test A*V ≈ V*H + r*e'
        end
    end
end

# Test incomplete Arnoldi factorization
@testset "Incomplete GKL factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64, Complex{Int})
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            if T == Complex{Int}
                A = rand(-100:100, (N, N)) + im * rand(-100:100, (N, N))
                v = rand(-100:100, (N,))
            else
                A = rand(T,(N,N))
                v = rand(T,(N,))
            end
            iter = @inferred GKLIterator(wrapop(A), wrapvec(v), orth)
            krylovdim = 3*n
            fact = @inferred initialize(iter)
            while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim
                @inferred expand!(iter, fact)

                U = hcat(unwrapvec.(basis(fact, :U))...)
                V = hcat(unwrapvec.(basis(fact, :V))...)
                B = rayleighquotient(fact)
                r = unwrapvec(residual(fact))
                β = normres(fact)
                e = rayleighextension(fact)
                @test U'*U ≈ I
                @test V'*V ≈ I
                @test norm(r) ≈ β
                @test A*V ≈ U*B + r*e'
                @test A'*U ≈ V*B'
            end

            fact = @inferred shrink!(fact, div(n,2))
            U = hcat(unwrapvec.(@inferred basis(fact, :U))...)
            V = hcat(unwrapvec.(@inferred basis(fact, :V))...)
            B = @inferred rayleighquotient(fact)
            r = unwrapvec(@inferred residual(fact))
            β = @inferred normres(fact)
            e = @inferred rayleighextension(fact)
            @test U'*U ≈ I
            @test V'*V ≈ I
            @test norm(r) ≈ β
            @test A*V ≈ U*B + r*e'
            @test A'*U ≈ V*B'
        end
    end
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Grammatical Evolution for the ACASX problem.
Example usage: config=configure(ACASX_GE,"normal","nvn_dasc"); acasx_ge(;config...)
"""
module ACASX_GE

export configure, acasx_ge

import Compat.ASCIIString
using ExprSearch.GE
using Datasets
using RLESUtils, FileUtils, ArrayUtils, Configure, LogSystems, Loggers
import RLESTypes.SymbolTable

using GrammarExpts
using ACASXProblem, DerivTreeVis
import Configure.configure

const CONFIGDIR = joinpath(dirname(@__FILE__), "..", "config")
const RESULTDIR = joinpath(dirname(@__FILE__), "..", "..", "..", "results")

configure(::Type{Val{:ACASX_GE}}, configs::AbstractString...) = configure_path(CONFIGDIR, configs...)

#nmacs vs nonnmacs
"""
Example call:
config=configure(ACASX_GE, "nvn_dasc", "normal")
acasx_ge(; config...)
"""
function acasx_ge(;outdir::AbstractString=joinpath(RESULTDIR, "./ACASX_GE"),
                  seed=1,
                  logfileroot::AbstractString="acasx_ge_log",

                  runtype::Symbol=:nmacs_vs_nonnmacs,
                  data::AbstractString="dasc",
                  manuals::AbstractString="",
                  clusterdataname::AbstractString="",

                  genome_size::Int64=20,
                  pop_size::Int64=50,
                  maxwraps::Int64=0,
                  top_keep::Float64=0.25,
                  top_seed::Float64=0.5,
                  rand_frac::Float64=0.25,
                  prob_mutation::Float64=0.2,
                  mutation_rate::Float64=0.2,
                  defaultcode::Union{Symbol,Expr}=:(eval(false)),
                  maxiterations::Int64=3,

                  limit_members::Int64=30,
                  vis::Bool=true)
    srand(seed)
    mkpath(outdir)

    problem = ACASXClustering(runtype, data, manuals, clusterdataname)

    logsys = GE.logsystem()
    empty_listeners!(logsys)
    send_to!(STDOUT, logsys, ["verbose1", "current_best_print", "result"])
    logs = TaggedDFLogger()
    send_to!(logs, logsys, ["code", "computeinfo", "current_best", "elapsed_cpu_s", "fitness",
        "fitness5", "parameters", "result"])

    ge_params = GEESParams(genome_size, pop_size, maxwraps, top_keep, top_seed, rand_frac,
        prob_mutation, mutation_rate, defaultcode, maxiterations, logsys;
        userargs=SymbolTable(:ids=>collect(1:length(problem.Dl))))

    result = exprsearch(ge_params, problem)

    #manually push! extra info to log
    push!(logs, "parameters", ["seed", seed])
    push!(logs, "parameters", ["runtype", runtype])
    push!(logs, "parameters", ["data", data])
    add_folder!(logs, "expression", [ASCIIString, ASCIIString, ASCIIString],
        ["raw", "pretty", "natural"])
    push!(logs, "expression", [string(result.expr), pretty_string(result.tree, FMT_PRETTY),
        pretty_string(result.tree, FMT_NATURAL, true)])

    #save log
    outfile = joinpath(outdir, "$(logfileroot).txt")
    save_log(outfile, logs)

    if vis
        derivtreevis(get_derivtree(result), joinpath(outdir, "$(logfileroot)_derivtreevis"))
    end

    textfile(joinpath(outdir, "summary.txt"), "mcts", seed=seed, maxiterations=maxiterations,
           fitness=result.fitness, expr=string(result.expr))

    result
end

end #module


text:
#=

Authors: Chase Coleman, Thomas Sargent

=#

using QuantEcon

#---------------------------------------------------------------------#
# Set up parameter values and LQ matrices
# Remember state is x_t = [1, y_{1, t}, y_{2, t}] and
# control is u_{i, t} = [y_{i, t+1} - y_{i, t}]
#---------------------------------------------------------------------#
a0 = 10.
a1 = 1.
beta = 1.
d = .5

a = eye(3)
b1 = [0.; 1.; 0.]
b2 = [0.; 0.; 1.]

r1 = [a0 0. 0.
      0. -a1 -a1 / 2.
      0  -a1 / 2. 0.]

r2 = [a0 0. 0.
      0. 0. -a1/2.
      0 -a1/2. -a1]

q1 = [-.5 * d]
q2 = [-.5 * d]

#---------------------------------------------------------------------#
# Solve using QE's nnash function
#---------------------------------------------------------------------#

f1, f2, p1, p2 = nnash(a, b1, b2, r1, r2, q1, q2, 0., 0., 0.,
                       0., 0., 0., tol=1e-8, max_iter=1000)


text:
## ROW INDEXING

function Base.getindex(task::SupervisedTask, r)
    X = selectrows(task.X, r)
    y = selectrows(task.y, r)
    is_probabilistic = task.is_probabilistic
    input_scitypes =  NamedTuple{schema(X).names}(schema(X).scitypes)
    target = task.target
    input_scitype = scitype(X)
    target_scitype = scitype(y)
    return SupervisedTask(X,
                          y,
                          is_probabilistic,
                          input_scitypes,
                          target,
                          input_scitype,
                          target_scitype)
end

## ROWS SHUFFLING

function Random.shuffle!(rng::AbstractRNG, task::SupervisedTask)
    rows = shuffle!(rng, Vector(1:nrows(task)))
    task.X = selectrows(task.X, rows)
    task.y = selectrows(task.y, rows)
    return task
end

function Random.shuffle!(task::SupervisedTask)
    rows = shuffle!(Vector(1:nrows(task)))
    task.X = selectrows(task.X, rows)
    task.y = selectrows(task.y, rows)
    return task
end

Random.shuffle(rng::AbstractRNG, task::SupervisedTask) =
    task[shuffle!(rng, Vector(1:nrows(task)))]
Random.shuffle(task::SupervisedTask) = task[shuffle!(Vector(1:nrows(task)))]

## TASK CONSTRUCORS WITH OPTIONAL TYPE COERCION

"""
    task = supervised(data=nothing,
                      types=Dict(),
                      target=nothing,
                      ignore=Symbol[],
                      is_probabilistic=false,
                      verbosity=1)

Construct a supervised learning task with input features `X` and
target `y`, where: `y` is the column vector from `data` named
`target`, if this is a single symbol, or a vector of tuples, if
`target` is a vector; `X` consists of all remaining columns of `data`
not named in `ignore`, and is a table unless it has only one column, in
which case it is a vector.

The data types of elements in a column of `data` named as a key of the
dictionary `types` are coerced to have a scientific type given by the
corresponding value. Possible values are `Continuous`, `Multiclass`,
`OrderedFactor` and `Count`. So, for example,
`types=Dict(:x1=>Count)` means elements of the column of `data` named
`:x1` will be coerced into integers (whose scitypes are always `Count`).

    task = supervised(X, y;
                      is_probabilistic=false,
                      verbosity=1)

A more customizable constructor, this returns a supervised learning
task with input features `X` and target `y`, where: `X` is a table or
vector (univariate inputs), while `y` must be a vector whose elements
are scalars, or tuples scalars (of constant length for ordinary
multivariate predictions, and of variable length for sequence
prediction). Table rows must correspond to patterns and columns to
features. Type coercion is not available for this constructor (but see
also [`coerce`](@ref)).

    X, y = task()

Returns the input `X` and target `y` of the task, also available as
`task.X` and `task.y`.

"""
supervised(; data=nothing, types=Dict(), kwargs...) =
	    SupervisedTask(; data = coerce(data, types), kwargs...)
supervised(X, y; kwargs...) = SupervisedTask(X, y; kwargs...)
"""
    task = unsupervised(data=nothing, types=Dict(), ignore=Symbol[], verbosity=1)

Construct an unsupervised learning task with given input `data`, which
should be a table or, in the case of univariate inputs, a single
vector.

The data types of elements in a column of `data` named as a key of the
dictionary `types` are coerced to have a scientific type given by the
corresponding value. Possible values are `Continuous`, `Multiclass`,
`OrderedFactor` and `Count`. So, for example,
`types=Dict(:x1=>Count)` means elements of the column of `data` named
`:x1` will be coerced into integers (whose scitypes are always `Count`).

Rows of `data` must correspond to patterns and columns to
features. Columns in `data` whose names appear in `ignore` are
ignored.

    X = task()

Return the input data in form to be used in models.

See also [`scitype`](@ref), [`scitype_union`](@ref)

"""
unsupervised(; data=nothing, types=Dict(), kwargs...) =
	    UnsupervisedTask(; data = coerce(types, data), kwargs...)

## MODEL MATCHING

models(task::SupervisedTask) = models() do model
    model.is_supervised &&
        task.target_scitype <: model.target_scitype &&
        task.input_scitype <: model.input_scitype &&
        ((task.is_probabilistic &&
          model.prediction_type == :probabilistic) ||
         (!(task.is_probabilistic) &&
          model.prediction_type == :deterministic))
end

models(task::UnsupervisedTask) = models() do model
    !(model.is_supervised) &&
        task.input_scitype <: model.input_scitype
end


text:
# This file is part of Kpax3. License is MIT.

abstract type KData end

"""
# Genetic data

## Description

DNA data and its metadata.

## Fields

* `data` Multiple sequence alignment (MSA) encoded as a binary (UInt8) matrix
* `id` units' ids
* `ref` reference sequence, i.e. a vector of the same length of the original
sequences storing the values of homogeneous sites. SNPs are instead represented
by a value of 29
* `val` vector with unique values per MSA site
* `key` vector with indices of each value

## Details

Let `n` be the total number of units and `ml` be the total number of unique
values observed at SNP `l`. Define m = m1 + ... + mL, where L is the total
number of SNPs.

`data` is a `m`-by-`n` indicator matrix, i.e. `data[j, i]` is `1` if unit `i`
possesses value `j`, `0` otherwise.

The value associated with column `j` can be obtained by `val[j]` while the SNP
position by `findall(ref == 29)[key[j]]`.

## References

Pessia A., Grad Y., Cobey S., Puranen J. S. and Corander J. (2015). K-Pax2:
Bayesian identification of cluster-defining amino acid positions in large
sequence datasets. _Microbial Genomics_ **1**(1).
<http://dx.doi.org/10.1099/mgen.0.000025>.
"""
struct NucleotideData <: KData
  data::Matrix{UInt8}
  id::Vector{String}
  ref::Vector{UInt8}
  val::Vector{UInt8}
  key::Vector{Int}
end

function NucleotideData(settings::KSettings)
  missval = if (length(settings.miss) == 1) && (settings.miss[1] == 0x00)
              0x00
            else
              UInt8('?')
            end

  (data, id, ref) = readfasta(settings.ifile, false, settings.miss, settings.l,
                              settings.verbose, settings.verbosestep)
  (bindata, val, key) = categorical2binary(data, UInt8(127), missval)

  NucleotideData(bindata, id, ref, val, key)
end

"""
# Genetic data

## Description

Amino acid data and its metadata.

## Fields

* `data` Multiple sequence alignment (MSA) encoded as a binary (UInt8) matrix
* `id` units' ids
* `ref` reference sequence, i.e. a vector of the same length of the original
sequences storing the values of homogeneous sites. SNPs are instead represented
by a value of 29
* `val` vector with unique values per MSA site
* `key` vector with indices of each value

## Details

Let `n` be the total number of units and `ml` be the total number of unique
values observed at SNP `l`. Define m = m1 + ... + mL, where L is the total
number of SNPs.

`data` is a `m`-by-`n` indicator matrix, i.e. `data[j, i]` is `1` if unit `i`
possesses value `j`, `0` otherwise.

The value associated with column `j` can be obtained by `val[j]` while the SNP
position by `findall(ref == 29)[key[j]]`.

## References

Pessia A., Grad Y., Cobey S., Puranen J. S. and Corander J. (2015). K-Pax2:
Bayesian identification of cluster-defining amino acid positions in large
sequence datasets. _Microbial Genomics_ **1**(1).
<http://dx.doi.org/10.1099/mgen.0.000025>.
"""
struct AminoAcidData{T} <: KData
  # generalize AminoAcidData to represent generic categorical data. Code must
  # be refactored
  data::Matrix{UInt8}
  id::Vector{String}
  ref::Vector{T}
  val::Vector{T}
  key::Vector{Int}
end

function AminoAcidData(settings::KSettings)
  missval = if (length(settings.miss) == 1) && (settings.miss[1] == 0x00)
              0x00
            else
              UInt8('?')
            end

  (data, id, ref) = readfasta(settings.ifile, true, settings.miss, settings.l,
                              settings.verbose, settings.verbosestep)
  (bindata, val, key) = categorical2binary(data, UInt8(127), missval)

  AminoAcidData(bindata, id, ref, val, key)
end

"""
# CSV data

## Description

Generic data and its metadata.

## Fields

* `data` original data matrix encoded as a binary (UInt8) matrix
* `id` units' ids
* `ref` reference observation, i.e. a vector of the same length of the original
observations storing the values of homogeneous sites. Polymorphisms are instead
represented by the string "."
* `val` vector with unique values per dataset column
* `key` vector with indices of each value

## Details

Let `n` be the total number of units and `ml` be the total number of unique
values observed at polymorphic column `l`. Define m = m1 + ... + mL, where L is
the total number of polymorphic columns.

`data` is a `m`-by-`n` indicator matrix, i.e. `data[j, i]` is `1` if unit `i`
possesses value `j`, `0` otherwise.

The value associated with column `j` can be obtained by `val[j]` while the
polymorphic position by `findall(ref == ".")[key[j]]`.

## References

Pessia A., Grad Y., Cobey S., Puranen J. S. and Corander J. (2015). K-Pax2:
Bayesian identification of cluster-defining amino acid positions in large
sequence datasets. _Microbial Genomics_ **1**(1).
<http://dx.doi.org/10.1099/mgen.0.000025>.
"""
function CategoricalData(settings::KSettings)
  (data, id, ref) = readdata(settings.ifile, ',', settings.misscsv,
                             settings.verbose, settings.verbosestep)
  (bindata, val, key) = categorical2binary(data, "")

  #CategoricalData(bindata, id, ref, val, key)
  # return an AminoAcidData object until the code has been refactored
  AminoAcidData(bindata, id, ref, val, key)
end

#=
struct CategoricalData <: KData
  data::Matrix{UInt8}
  id::Vector{String}
  ref::Vector{String}
  val::Vector{String}
  key::Vector{Int}
end
=#


text:
abstract type AbstractGP{L<:Likelihood,I<:Inference,T<:Real,V<:AbstractArray{T}} end

function Random.rand!(model::AbstractGP,A::DenseArray{T},X::AbstractArray{T}) where {T<:Real}
    rand!(MvNormal(predict_f(model,X,covf=true,fullcov=true)...),A)
end

function Random.rand(model::AbstractGP,X::AbstractArray{T},n::Int=1) where {T<:Real}
    if model.nLatent == 1
        rand!(model,Array{T}(undef,size(X,1),n),X)
    else
        @error "Sampling not implemented for multiple output GPs"
    end
end


text:
# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp

"""
  gen_grad_func(vi, spl, model)

Generate a function that takes `θ` and returns log-joint probabilty and gradient at `θ`, using Turing-related variables `vi`, `spl` and `model`.
"""
function gen_grad_func(vi, spl, model)

  grad_func(θ::T) where {T<:Union{Vector,SubArray}} = begin

    if ADBACKEND == :forward_diff
      vi[spl] = θ
      grad = gradient(vi, model, spl)
    elseif ADBACKEND == :reverse_diff
      grad = gradient_r(θ, vi, model, spl)
    end

    return getlogp(vi), grad

  end

  return grad_func

end

"""
  gen_lj_func(vi, spl, model)

Generate a function that takes `θ` and returns log-joint probabilty at `θ`, using Turing-related variables `vi`, `spl` and `model`.
"""
function gen_lj_func(vi, spl, model)

  lj_func(θ) = begin

    vi[spl][:] = θ[:]

    return runmodel(model, vi, spl).logp

  end

  return lj_func

end

"""
  gen_rev_func(vi, spl)

Generate a function that takes `x_old` and `lp_old` and resets the values in `vi` for `spl` as `x_old` and the log-joint probabilty as `lp_old`, using local variables `vi` and `spl` in the scope where `gen_rev_func` is called.
"""
function gen_rev_func(vi, spl)

  rev_func(θ_old::T, old_logp::R) where {T<:Union{Vector,SubArray},R<:Real} = begin

    if ADBACKEND == :forward_diff
      vi[spl] = θ_old
    elseif ADBACKEND == :reverse_diff
      vi_spl = vi[spl]
      for i = 1:length(θ_old)
        vi_spl[i] = θ_old[i]
      end
    end
    setlogp!(vi, old_logp)

  end

  return rev_func

end

"""
  gen_log_func(spl)

Generate a function that takes no argument and performs logging for number of leapfrog steps used, using the local variable `spl` in the scope where `gen_log_func` is called.
"""
function gen_log_func(spl)

  log_func() = begin

    spl.info[:lf_num] += 1
    spl.info[:total_lf_num] += 1  # record leapfrog num

  end

  return log_func

end

runmodel(model::Function, vi::VarInfo, spl::Union{Nothing,Sampler}) = begin
  @debug "run model..."
  setlogp!(vi, zero(Real))
  if spl != nothing spl.info[:total_eval_num] += 1 end
  # model(vi=vi, sampler=spl) # run model
  Base.invokelatest(model, vi, spl)
end

sample_momentum(vi::VarInfo, spl::Sampler) = begin
  @debug "sampling momentum..."
  # randn(length(getranges(vi, spl))) ./ spl.info[:wum][:stds]

  d = length(getranges(vi, spl))
  stds = spl.info[:wum][:stds]

  return _sample_momentum(d, stds)

end

function _sample_momentum(d::Int, stds::Vector)

  return randn(d) ./ stds

end

# Leapfrog step
# NOTE: leapfrog() doesn't change θ in place!
leapfrog(_θ::T, p::Vector{Float64}, τ::Int, ϵ::Float64,
          model::Function, vi::VarInfo, spl::Sampler) where {T<:Union{Vector,SubArray}} = begin

  θ = realpart(_θ)

  lp_grad_func = gen_grad_func(vi, spl, model)
  rev_func = gen_rev_func(vi, spl)
  log_func = gen_log_func(spl)

  return _leapfrog(θ, p, τ, ϵ, lp_grad_func; rev_func=rev_func, log_func=log_func)

end

function _leapfrog(θ::T, p::Vector{Float64}, τ::Int, ϵ::Float64, lp_grad_func::Function;
                   rev_func=nothing, log_func=nothing) where {T<:Union{Vector,SubArray}}

  old_logp, grad = lp_grad_func(θ)
  verifygrad(grad) || (return θ, p, 0)

  τ_valid = 0
  for t in 1:τ
    # NOTE: we dont need copy here becase arr += another_arr
    #       doesn't change arr in-place
    p_old = p; θ_old = copy(θ)

    p -= ϵ .* grad / 2
    θ += ϵ .* p  # full step for state

    if log_func != nothing log_func() end

    old_logp, grad = lp_grad_func(θ)
    if ~verifygrad(grad)
      if rev_func != nothing rev_func(θ_old, old_logp) end
      θ = θ_old; p = p_old; break
    end

    p -= ϵ * grad / 2

    τ_valid += 1
  end

  θ, p, τ_valid
end

# Compute Hamiltonian
find_H(p::Vector, model::Function, vi::VarInfo, spl::Sampler) = begin

  # Old code
  # # NOTE: getlogp(vi) = 0 means the current vals[end] hasn't been used at all.
  # #       This can be a result of link/invlink (where expand! is used)
  # if getlogp(vi) == 0 vi = runmodel(model, vi, spl) end
  #
  # p_orig = p .* spl.info[:wum][:stds]
  #
  # H = dot(p_orig, p_orig) / 2 + realpart(-getlogp(vi))
  # if isnan(H) H = Inf else H end
  #
  # H

  logpdf_func_float = gen_lj_func(vi, spl, model)

  return _find_H(vi[spl], p, logpdf_func_float, spl.info[:wum][:stds])

end

function _find_H(theta::T, p::Vector, logpdf_func_float::Function, stds::Vector) where {T<:Union{Vector,SubArray}}

  lp = logpdf_func_float(theta)

  return _find_H(theta, p, lp, stds)

end

function _find_H(theta::T, p::Vector, lp::Real, stds::Vector) where {T<:Union{Vector,SubArray}}

  p_orig = p .* stds

  H = 0.5 * dot(p_orig, p_orig) + (-lp)

  H = isnan(H) ? Inf : H

  return H

end

# TODO: remove used Turing-wrapper functions

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/base_hmc.hpp
find_good_eps(model::Function, vi::VarInfo, spl::Sampler{T}) where T = begin
  println("[Turing] looking for good initial eps...")
  ϵ = 0.1

  p = sample_momentum(vi, spl)
  H0 = find_H(p, model, vi, spl)

  θ = realpart(vi[spl])
  θ_prime, p_prime, τ = leapfrog(θ, p, 1, ϵ, model, vi, spl)
  h = τ == 0 ? Inf : find_H(p_prime, model, vi, spl)

  delta_H = H0 - h
  direction = delta_H > log(0.8) ? 1 : -1

  iter_num = 1

  # Heuristically find optimal ϵ
  while (iter_num <= 12)

    p = sample_momentum(vi, spl)
    H0 = find_H(p, model, vi, spl)

    θ_prime, p_prime, τ = leapfrog(θ, p, 1, ϵ, model, vi, spl)
    h = τ == 0 ? Inf : find_H(p_prime, model, vi, spl)
    @debug "direction = $direction, h = $h"

    delta_H = H0 - h

    if ((direction == 1) && !(delta_H > log(0.8)))
      break;
    elseif ((direction == -1) && !(delta_H < log(0.8)))
      break;
    else
      ϵ = direction == 1 ? 2.0 * ϵ : 0.5 * ϵ
    end

    iter_num += 1
  end

  while h == Inf  # revert if the last change is too big
    ϵ = ϵ / 2               # safe is more important than large
    θ_prime, p_prime, τ = leapfrog(θ, p, 1, ϵ, model, vi, spl)
    h = τ == 0 ? Inf : find_H(p_prime, model, vi, spl)
  end
  println("\r[$T] found initial ϵ: ", ϵ)
  ϵ
end

"""
  mh_accept(H, H_new)

Peform MH accept criteria. Returns a boolean for whether or not accept and the acceptance ratio in log space.

"""
function mh_accept(H, H_new; log_proposal_ratio=nothing)

  logα = min(0, -(H_new - H))

  u = rand()
  logu = log(u)

  if log_proposal_ratio == nothing
    is_accept = (logu + H_new < min(H_new, H))
  else
    is_accept = (logu + H_new < H + log_proposal_ratio)
  end
  return is_accept, logα

end


text:
struct FileLocInfo
    filepath::String
    line::Int
    # 0 if unknown
    column::Int
    # The line at which the current context starts, 0 if unknown
    defline::Int
end

struct BufferLocInfo
    data::String
    line::Int
    # 0 if unknown
    column::Int
    defline::Int
end

function loc_for_fname(file::String, line::Integer, defline::Integer)
    if startswith(file, "REPL[")
        hist_idx = parse(Int,string(file)[6:end-1])
        isdefined(Base, :active_repl) || return nothing, ""
        hp = Base.active_repl.interface.modes[1].hist
        return BufferLocInfo(hp.history[hp.start_idx+hist_idx], line, 0, defline)
    else
        for path in SEARCH_PATH
            fullpath = joinpath(path,string(file))
            if isfile(fullpath)
                return FileLocInfo(fullpath, line, 0, defline)
            end
        end
    end
    return nothing
end

function locinfo(frame::Frame)
    if frame.framecode.scope isa Method
        meth = frame.framecode.scope
        file, def_line = JuliaInterpreter.whereis(meth)
        _, current_line = JuliaInterpreter.whereis(frame)
        return loc_for_fname(file, current_line, def_line)
    else
        println("not yet implemented")
    end
end

# Used for the tests
const _print_full_path = Ref(true)

function locdesc(frame::Frame)
    sprint() do io
        if frame.framecode.scope isa Method
            meth = frame.framecode.scope
            argnames = frame.framecode.src.slotnames[2:meth.nargs]
            spectypes = Any[Any for i=1:length(argnames)]
            print(io, meth.name,'(')
            first = true
            for (argname, argT) in zip(argnames, spectypes)
                first || print(io, ", ")
                first = false
                print(io, argname)
                !(argT === Any) && print(io, "::", argT)
            end
            print(io, ") at ", _print_full_path[] ? meth.file : basename(String(meth.file)), ":",meth.line)
        else
            println("not yet implemented")
        end
    end
end

"""
Determine the offsets in the source code to print, based on the offset of the
currently highlighted part of the code, and the start and stop line of the
entire function.
"""
function compute_source_offsets(code::String, offset::Integer, startline::Integer, stopline::Integer; file::SourceFile = SourceFile(code))
    offsetline = compute_line(file, offset)
    if offsetline - 3 > length(file.offsets) || startline > length(file.offsets)
        return -1, -1
    end
    startoffset = max(file.offsets[max(offsetline-3,1)], file.offsets[startline])
    stopoffset = lastindex(code)-1
    if offsetline + 3 < lastindex(file.offsets)
        stopoffset = min(stopoffset, file.offsets[offsetline + 3]-1)
    end
    if stopline + 1 < lastindex(file.offsets)
        stopoffset = min(stopoffset, file.offsets[stopline + 1]-1)
    end
    startoffset, stopoffset
end

function print_sourcecode(io::IO, code::String, line::Integer, defline::Integer; file::SourceFile = SourceFile(code))
    startoffset, stopoffset = compute_source_offsets(code, file.offsets[line], defline, line+3; file=file)

    if startoffset == -1
        printstyled(io, "Line out of file range (bad debug info?)", color=:bold)
        return
    end

    # Compute necessary data for line numbering
    startline = compute_line(file, startoffset)
    stopline = compute_line(file, stopoffset)
    current_line = line
    stoplinelength = length(string(stopline))

    code = split(code[(startoffset+1):(stopoffset+1)],'\n')
    lineno = startline

    if !isempty(code) && isempty(code[end])
        pop!(code)
    end

    for textline in code
        printstyled(io,
            string(lineno, " "^(stoplinelength-length(lineno)+1));
            color = lineno == current_line ? :yellow : :bold)
        println(io, textline)
        lineno += 1
    end
    println(io)
end


text:
module Pkg

export Git, Dir, GitHub, Types, Reqs, Cache, Read, Query, Resolve, Write, Generate, Entry
export dir, init, rm, add, available, installed, status, clone, checkout,
       release, fix, update, resolve, register, tag, publish, generate, test

const DEFAULT_META = "git://github.com/JuliaLang/METADATA.jl"
const META_BRANCH = "metadata-v2"

for file in split("git dir github types reqs cache read query resolve write generate entry")
    include("pkg/$file.jl")
end
const cd = Dir.cd

dir(path...) = Dir.path(path...)
init(meta::String=DEFAULT_META, branch::String=META_BRANCH) = Dir.init(meta,branch)

edit() = cd(Entry.edit)
rm(pkg::String) = cd(Entry.rm,pkg)
add(pkg::String, vers::VersionNumber...) = cd(Entry.add,pkg,vers...)

available() = cd(Entry.available)
available(pkg::String) = cd(Entry.available,pkg)

installed() = cd(Entry.installed)
installed(pkg::String) = cd(Entry.installed,pkg)

status(io::IO=STDOUT) = cd(Entry.status,io)

clone(url_or_pkg::String) = cd(Entry.clone,url_or_pkg)
clone(url::String, pkg::String) = cd(Entry.clone,url,pkg)

checkout(pkg::String, branch::String="master"; merge::Bool=true, pull::Bool=true) =
    cd(Entry.checkout,pkg,branch,merge,pull)

free(pkg::String) = cd(Entry.free,pkg)

pin(pkg::String) = cd(Entry.pin,pkg)
pin(pkg::String, ver::VersionNumber) = cd(Entry.pin,pkg,ver)

update() = cd(Entry.update,META_BRANCH)
resolve() = cd(Entry.resolve)

register(pkg::String) = cd(Entry.register,pkg)
register(pkg::String, url::String) = cd(Entry.register,pkg,url)

tag(pkg::String, sym::Symbol=:patch) = cd(Entry.tag,pkg,sym)
tag(pkg::String, sym::Symbol, commit::String) = cd(Entry.tag,pkg,sym,false,commit)

tag(pkg::String, ver::VersionNumber; force::Bool=false) = cd(Entry.tag,pkg,ver,force)
tag(pkg::String, ver::VersionNumber, commit::String; force::Bool=false) =
	cd(Entry.tag,pkg,ver,force,commit)

submit(pkg::String) = cd(Entry.submit,pkg)
submit(pkg::String, commit::String) = cd(Entry.submit,pkg,commit)

publish() = cd(Entry.publish,META_BRANCH)

build() = cd(Entry.build)
build(pkgs::String...) = cd(Entry.build,[pkgs...])

generate(pkg::String, license::String; force::Bool=false) =
	cd(Generate.package,pkg,license,force=force)

test() = cd(Entry.test)
test(pkgs::String...) = cd(Entry.test,String[pkgs...])

@deprecate release free
@deprecate fixup build
@deprecate fix pin

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using eFEMpart, JLD

using LinearAlgebra # for condition number calculation

function main()
  function computeNorms(N,order)
    start = time()
    # load mesh
    mesh = squareMesh([-2,2,0,1],N,order)

    OperatorType = :AdvDiffAS

    Pe = 2.3
    windX(x,y) = -x^4*y^2 + 5.9; windY(x,y) = x^3*y^3
    wx = [windX.(mesh.xy[i].x,mesh.xy[i].y) for i=1:length(mesh.xy)]
    wy = [windY.(mesh.xy[i].x,mesh.xy[i].y) for i=1:length(mesh.xy)]
    param  = AdvDiffParam(wx,wy,Pe*ones(length(mesh.xy)))

    dNodes = Dirichlet(:left,:right,:top,:bottom)
    dBCf = Dirichlet((x,y) -> x^4*y^4)
    ff   = Forcing((x,y) -> 4*5.9*x^3*y^4 - 4/Pe*x^2*y^2*(4*x^2 + 3*y^2))
    Nodes = [dNodes]
    bcfun = [dBCf,ff]

    prob = Problem(mesh,Nodes,bcfun,OperatorType)
    sol = solve(prob,mesh,param)

    elapsed = time() - start

    # compute condition number of operator matrix
    #LinOp = GenerateSystem(mesh,prob,param)
    #ApplyBC!(LinOp,mesh,prob,param,OperatorType)
    #κ = cond(Array(LinOp.Op),2)
    κ = 0.5

    # compute mesh h
    h = hCalc(mesh)

    # generate exact solution
    Uexact(x,y) = dBCf.f(x,y)
    UexactArr = [Uexact.(mesh.xy[i].x,mesh.xy[i].y) for i=1:length(mesh.xy)]

    # compute difference
    err = sol.u - UexactArr

    # compute Li norm
    L1   = DomainNorm(mesh.xy,mesh.cm,err;normID="1")
    L2   = DomainNorm(mesh.xy,mesh.cm,err;normID="2")
    Linf = DomainNorm(mesh.xy,mesh.cm,err;normID="Inf")

    return κ,h,L1,L2,Linf,elapsed
  end

  Narr    = [8,8,16,32]#,64]#,128]#,256]
  order   = [1,2]
  N       = length(Narr)
  κarr    = zeros(N,2)
  harr    = zeros(N,2)
  L1arr   = zeros(N,2)
  L2arr   = zeros(N,2)
  LInfarr = zeros(N,2)
  timearr = zeros(N,2)

  for i=1:N
    for j=1:2
      n = Narr[i]
      o = order[j]
      κarr[i,j],harr[i,j],L1arr[i,j],L2arr[i,j],LInfarr[i,j],timearr[i,j] = computeNorms(n,o)
      println("completed N=$(n), order=$(o)")
    end
  end

  κarr = κarr[2:end,:]
  harr = harr[2:end,:]
  L1arr = L1arr[2:end,:]
  L2arr = L2arr[2:end,:]
  LInfarr = LInfarr[2:end,:]
  timearr = timearr[2:end,:]

  # export output to *.jld file
  run(`mkdir -p data`)

  save("data/TEMP_advdiffAS_validation.jld", "harr",harr[:,1],
                                          "o1L1arr", L1arr[:,1],
                                          "o1L2arr", L2arr[:,1],
                                          "o1LInfarr",LInfarr[:,1],
                                          "o1timearr",timearr[:,1],
                                          "o2L1arr", L1arr[:,2],
                                          "o2L2arr", L2arr[:,2],
                                          "o2LInfarr",LInfarr[:,2],
                                          "o2timearr",timearr[:,2],
                                          "κarr",κarr[:,2])

  nothing
end

main()

max_stars_repo_name: pseastham/eFEMpart


text:
using SpecialFunctions: erfcinv

@doc raw"""
The Morlet wavelet.
```math
 g(t) = \sqrt{\frac{\delta}{2\pi}}\exp(2\pi\,i\,t-t^2\,\delta)\,,
```
where δ specifies the time-frequency resolution of the wavelet.
"""
struct MorletWavelet <: GenericContinuousWavelet
    δ::Float64
    norm::Float64
    σ::Float64

    @doc raw"""
        MorletWavelet(δ::Real; ϵ::Real=1e-3)

    Constructs a new Morlet wavelet, whith the given `dff` parameter specifying the time-frequency
    resolution. The optional keyword arguent ϵ specifies the cutoff at which the kernel evaluation
    gets truncated. It is defined as the fraction of total power loss of the mother wavelet. Smaller
    values of ϵ will increase the precision of the wavelet transform on the cost of longer kernels
    leading to slower convolutions.
    """
    function MorletWavelet(δ::Real; ϵ::Real=1e-3)
        new(δ, δ/2/π, sqrt(2)*erfcinv(2*ϵ));
    end
end

function eval_analysis(wav::MorletWavelet, t::Float64)
    exp(-t^2 * wav.δ/2 + 2π*t*1im) * wav.norm;
end

eval_synthesis(wav::MorletWavelet, t::Float64) = eval_analysis(wav, t);

function eval_repkern(wav::MorletWavelet, a::Float64, b::Float64)
    a2p1 = (a^2 + 1); am1 = (a - 1);
    d = wav.δ; d2 = d*d;
    c = d2/sqrt(2*π*d*a2p1);
    re = -(d2*b^2 + 4*π^2*am1^2)/(2*d*a2p1);
    im = 2*π*b*(1 + am1/a2p1)/a;
    c * exp(re + im*1im);
end

function cutoff_time(wav::MorletWavelet)
    # 99.9% power at scale 1
    wav.σ/sqrt(wav.δ);
end;

function cutoff_freq(wav::MorletWavelet)
    # 99.9% power at scale 1
    1 + wav.σ*sqrt(wav.δ);
end;


text:
using OrdinaryDiffEq, DifferentialEquations

tspan = (0.0f0, 200.0f0)

"""
FitzHugh Nagumo oscillator model from:
    https://en.wikipedia.org/wiki/FitzHugh%E2%80%93Nagumo_model
"""
function fhn_system(du, u, p, t;
                     U_func_time=U_func_time_trivial,
                     U_func_space=U_func_space_trivial)
    Ft = U_func_time(t)
    Fs = U_func_space(u)
    a, b, I_external, τ = p
    v, ω = u

    du[1] = v - (v^3)/3 - ω + I_external
    du[2] = (v + a - b*ω) / τ

    # Add controllers in
    du .+= Ft .+ Fs
end

# Generate data, with a time component
p = [0.7, 0.8, 0.5, 12.5]
u0 = [1.0, 0.0]
ts = range(tspan[1], tspan[2], length=5001)

function solve_fhn_system(;U_func_time=U_func_time_trivial,
                              U_func_space=U_func_space_trivial)
    prob = ODEProblem((du, u, p, t)->fhn_system(du, u, p, t,
                                U_func_time=U_func_time,
                                U_func_space=U_func_space),
                    u0, tspan, p)
    sol = solve(prob, Tsit5(), saveat=ts);
    return sol
end

#####
##### Turing function for Bayesian parameter estimation
#####
# Unforced version for use with Turing.jl
# function vdp_system(u, p, t)
#     u = convert.(eltype(p),u)
#     beta, gamma = p
#     S, I, R = u
#
#     N = sum([S,I,R])
#     tmp = beta*I*S/N
#     du = [-tmp; tmp - gamma*I; gamma*I]
#     return du
# end

#####
##### True model in SINDy syntax
#####
a, b, I, τ = p
#      x    y   c   xx xy yy xxx  xxy xyy yyy
 A = [[1   -1   I   0  0  0  -1/3 0   0   0];
      [1/τ -b/τ a/τ 0  0  0  0    0   0   0]]
n = size(A, 1)
sindy_library = Dict(
    "cross_terms"=>[2,3],
    "constant"=>nothing
);
core_dyn_true = sindyModel(ts, A,
                            convert_string2function(sindy_library),
                            ["v", "ω"])

# Actually export
export ts, solve_fhn_system, fhn_system, core_dyn_true


text:
@testset "6.4.7 (d hyper)^m (a+b (c coth)^n)^p" begin
    (a, b, c, d, x, ) = @variables a b c d x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*sinh(e+f*x)^m*(a+b*coth(e+f*x)^n)^p=#

    #= ::Title:: =#
    #=Integrands*of*the*form*cosh(e+f*x)^m*(a+b*coth(e+f*x)^n)^p=#

    #= ::Title:: =#
    #=Integrands*of*the*form*coth(e+f*x)^m*(a+b*coth(e+f*x)^n)^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*coth(e+f*x)^m*(a+b*coth(e+f*x)^2)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*coth(e+f*x)^m*(a+b*coth(e+f*x)^2)^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [(a + b*coth(c + d*x)^2)^5, x, 4, (a + b)^5*x - (b*(5*a^4 + 10*a^3*b + 10*a^2*b^2 + 5*a*b^3 + b^4)*coth(c + d*x))/d - (b^2*(10*a^3 + 10*a^2*b + 5*a*b^2 + b^3)*coth(c + d*x)^3)/(3*d) - (b^3*(10*a^2 + 5*a*b + b^2)*coth(c + d*x)^5)/(5*d) - (b^4*(5*a + b)*coth(c + d*x)^7)/(7*d) - (b^5*coth(c + d*x)^9)/(9*d)]
    @test_int [(a + b*coth(c + d*x)^2)^4, x, 4, (a + b)^4*x - (b*(2*a + b)*(2*a^2 + 2*a*b + b^2)*coth(c + d*x))/d - (b^2*(6*a^2 + 4*a*b + b^2)*coth(c + d*x)^3)/(3*d) - (b^3*(4*a + b)*coth(c + d*x)^5)/(5*d) - (b^4*coth(c + d*x)^7)/(7*d)]
    @test_int [(a + b*coth(c + d*x)^2)^3, x, 4, (a + b)^3*x - (b*(3*a^2 + 3*a*b + b^2)*coth(c + d*x))/d - (b^2*(3*a + b)*coth(c + d*x)^3)/(3*d) - (b^3*coth(c + d*x)^5)/(5*d)]
    @test_int [(a + b*coth(c + d*x)^2)^2, x, 4, (a + b)^2*x - (b*(2*a + b)*coth(c + d*x))/d - (b^2*coth(c + d*x)^3)/(3*d)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [1/(a + b*coth(c + d*x)^2),x, 3, x/(a + b) - (sqrt(b)*atan((sqrt(a)*tanh(c + d*x))/sqrt(b)))/(sqrt(a)*(a + b)*d)]
    @test_int [1/(a + b*coth(c + d*x)^2)^2,x, 5, x/(a + b)^2 - (sqrt(b)*(3*a + b)*atan((sqrt(a)*tanh(c + d*x))/sqrt(b)))/(2*a^(3/2)*(a + b)^2*d) + (b*coth(c + d*x))/(2*a*(a + b)*d*(a + b*coth(c + d*x)^2))]
    @test_int [1/(a + b*coth(c + d*x)^2)^3,x, 6, x/(a + b)^3 - (sqrt(b)*(15*a^2 + 10*a*b + 3*b^2)*atan((sqrt(a)*tanh(c + d*x))/sqrt(b)))/(8*a^(5/2)*(a + b)^3*d) + (b*coth(c + d*x))/(4*a*(a + b)*d*(a + b*coth(c + d*x)^2)^2) + (b*(7*a + 3*b)*coth(c + d*x))/(8*a^2*(a + b)^2*d*(a + b*coth(c + d*x)^2))]
    @test_int [1/(a + b*coth(c + d*x)^2)^4,x, 7, x/(a + b)^4 - (sqrt(b)*(35*a^3 + 35*a^2*b + 21*a*b^2 + 5*b^3)*atan((sqrt(a)*tanh(c + d*x))/sqrt(b)))/(16*a^(7/2)*(a + b)^4*d) + (b*coth(c + d*x))/(6*a*(a + b)*d*(a + b*coth(c + d*x)^2)^3) + (b*(11*a + 5*b)*coth(c + d*x))/(24*a^2*(a + b)^2*d*(a + b*coth(c + d*x)^2)^2) + (b*(19*a^2 + 16*a*b + 5*b^2)*coth(c + d*x))/(16*a^3*(a + b)^3*d*(a + b*coth(c + d*x)^2))]

    @test_int [1/(1 - 2*coth(x)^2), x, 3, -x + sqrt(2)*atanh(tanh(x)/sqrt(2))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*coth(e+f*x)^m*(a+b*coth(e+f*x)^2)^(p/2)*when*a+b=0=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [sqrt(1 - coth(x)^2), x, 3, asin(coth(x))]
    @test_int [sqrt(-1 + coth(x)^2), x, 4, -atanh(coth(x)/sqrt(csch(x)^2))]

    @test_int [(1 - coth(x)^2)^(3/2), x, 4, (1/2)*asin(coth(x)) + (1/2)*coth(x)*sqrt(-csch(x)^2)]
    @test_int [(-1 + coth(x)^2)^(3/2), x, 5, (1/2)*atanh(coth(x)/sqrt(csch(x)^2)) - (1/2)*coth(x)*sqrt(csch(x)^2)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [1/sqrt(1 - coth(x)^2), x, 3, coth(x)/sqrt(-csch(x)^2)]
    @test_int [1/sqrt(-1 + coth(x)^2), x, 3, coth(x)/sqrt(csch(x)^2)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*coth(e+f*x)^m*(a+b*coth(e+f*x)^2)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [coth(x)^3*sqrt(a + b*coth(x)^2), x, 6, sqrt(a + b)*atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b)) - sqrt(a + b*coth(x)^2) - (a + b*coth(x)^2)^(3/2)/(3*b)]
    @test_int [coth(x)^2*sqrt(a + b*coth(x)^2), x, 7, -(((a + 2*b)*atanh((sqrt(b)*coth(x))/sqrt(a + b*coth(x)^2)))/(2*sqrt(b))) + sqrt(a + b)*atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2)) - (1/2)*coth(x)*sqrt(a + b*coth(x)^2)]
    @test_int [coth(x)^1*sqrt(a + b*coth(x)^2), x, 5, sqrt(a + b)*atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b)) - sqrt(a + b*coth(x)^2)]
    @test_int [coth(x)^0*sqrt(a + b*coth(x)^2), x, 6, (-sqrt(b))*atanh((sqrt(b)*coth(x))/sqrt(a + b*coth(x)^2)) + sqrt(a + b)*atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2))]
    @test_int [tanh(x)^1*sqrt(a + b*coth(x)^2), x, 7, (-sqrt(a))*atanh(sqrt(a + b*coth(x)^2)/sqrt(a)) + sqrt(a + b)*atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b))]
    @test_int [tanh(x)^2*sqrt(a + b*coth(x)^2), x, 5, sqrt(a + b)*atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2)) - sqrt(a + b*coth(x)^2)*tanh(x)]

    @test_int [coth(x)^3*(a + b*coth(x)^2)^(3/2), x, 7, (a + b)^(3/2)*atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b)) - (a + b)*sqrt(a + b*coth(x)^2) - (1/3)*(a + b*coth(x)^2)^(3/2) - (a + b*coth(x)^2)^(5/2)/(5*b)]
    @test_int [coth(x)^2*(a + b*coth(x)^2)^(3/2), x, 8, -(((3*a^2 + 12*a*b + 8*b^2)*atanh((sqrt(b)*coth(x))/sqrt(a + b*coth(x)^2)))/(8*sqrt(b))) + (a + b)^(3/2)*atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2)) - (1/8)*(5*a + 4*b)*coth(x)*sqrt(a + b*coth(x)^2) - (1/4)*b*coth(x)^3*sqrt(a + b*coth(x)^2)]
    @test_int [coth(x)^1*(a + b*coth(x)^2)^(3/2), x, 6, (a + b)^(3/2)*atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b)) - (a + b)*sqrt(a + b*coth(x)^2) - (1/3)*(a + b*coth(x)^2)^(3/2)]
    @test_int [coth(x)^0*(a + b*coth(x)^2)^(3/2), x, 7, (-(1/2))*sqrt(b)*(3*a + 2*b)*atanh((sqrt(b)*coth(x))/sqrt(a + b*coth(x)^2)) + (a + b)^(3/2)*atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2)) - (1/2)*b*coth(x)*sqrt(a + b*coth(x)^2)]
    @test_int [tanh(x)^1*(a + b*coth(x)^2)^(3/2), x, 8, (-a^(3/2))*atanh(sqrt(a + b*coth(x)^2)/sqrt(a)) + (a + b)^(3/2)*atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b)) - b*sqrt(a + b*coth(x)^2)]
    @test_int [tanh(x)^2*(a + b*coth(x)^2)^(3/2), x, 7, (-b^(3/2))*atanh((sqrt(b)*coth(x))/sqrt(a + b*coth(x)^2)) + (a + b)^(3/2)*atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2)) - a*sqrt(a + b*coth(x)^2)*tanh(x)]

    @test_int [sqrt(1 + coth(x)^2), x, 5, -asinh(coth(x)) + sqrt(2)*atanh((sqrt(2)*coth(x))/sqrt(1 + coth(x)^2))]
    @test_int [sqrt(-1 - coth(x)^2), x, 6, atan(coth(x)/sqrt(-1 - coth(x)^2)) - sqrt(2)*atan((sqrt(2)*coth(x))/sqrt(-1 - coth(x)^2))]

    @test_int [(1 + coth(x)^2)^(3/2), x, 6, (-(5/2))*asinh(coth(x)) + 2*sqrt(2)*atanh((sqrt(2)*coth(x))/sqrt(1 + coth(x)^2)) - (1/2)*coth(x)*sqrt(1 + coth(x)^2)]
    @test_int [(-1 - coth(x)^2)^(3/2), x, 7, (-(5/2))*atan(coth(x)/sqrt(-1 - coth(x)^2)) + 2*sqrt(2)*atan((sqrt(2)*coth(x))/sqrt(-1 - coth(x)^2)) + (1/2)*coth(x)*sqrt(-1 - coth(x)^2)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [coth(x)^3/sqrt(a + b*coth(x)^2), x, 5, atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b))/sqrt(a + b) - sqrt(a + b*coth(x)^2)/b]
    @test_int [coth(x)^2/sqrt(a + b*coth(x)^2), x, 6, -(atanh((sqrt(b)*coth(x))/sqrt(a + b*coth(x)^2))/sqrt(b)) + atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2))/sqrt(a + b)]
    @test_int [coth(x)^1/sqrt(a + b*coth(x)^2), x, 4, atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b))/sqrt(a + b)]
    @test_int [coth(x)^0/sqrt(a + b*coth(x)^2), x, 3, atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2))/sqrt(a + b)]
    @test_int [tanh(x)^1/sqrt(a + b*coth(x)^2), x, 7, -(atanh(sqrt(a + b*coth(x)^2)/sqrt(a))/sqrt(a)) + atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b))/sqrt(a + b)]
    @test_int [tanh(x)^2/sqrt(a + b*coth(x)^2), x, 5, atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2))/sqrt(a + b) - (sqrt(a + b*coth(x)^2)*tanh(x))/a]

    @test_int [coth(x)^3/(a + b*coth(x)^2)^(3/2), x, 5, atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b))/(a + b)^(3/2) + a/(b*(a + b)*sqrt(a + b*coth(x)^2))]
    @test_int [coth(x)^2/(a + b*coth(x)^2)^(3/2), x, 4, atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2))/(a + b)^(3/2) - coth(x)/((a + b)*sqrt(a + b*coth(x)^2))]
    @test_int [coth(x)/(a + b*coth(x)^2)^(3/2), x, 5, atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b))/(a + b)^(3/2) - 1/((a + b)*sqrt(a + b*coth(x)^2))]
    @test_int [tanh(x)/(a + b*coth(x)^2)^(3/2), x, 8, -(atanh(sqrt(a + b*coth(x)^2)/sqrt(a))/a^(3/2)) + atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b))/(a + b)^(3/2) + b/(a*(a + b)*sqrt(a + b*coth(x)^2))]
    @test_int [tanh(x)^2/(a + b*coth(x)^2)^(3/2), x, 6, atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2))/(a + b)^(3/2) + (b*tanh(x))/(a*(a + b)*sqrt(a + b*coth(x)^2)) - ((a + 2*b)*sqrt(a + b*coth(x)^2)*tanh(x))/(a^2*(a + b))]

    @test_int [coth(x)^3/(a + b*coth(x)^2)^(5/2), x, 6, atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b))/(a + b)^(5/2) + a/(3*b*(a + b)*(a + b*coth(x)^2)^(3/2)) - 1/((a + b)^2*sqrt(a + b*coth(x)^2))]
    @test_int [coth(x)^2/(a + b*coth(x)^2)^(5/2), x, 6, atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2))/(a + b)^(5/2) - coth(x)/(3*(a + b)*(a + b*coth(x)^2)^(3/2)) - ((2*a - b)*coth(x))/(3*a*(a + b)^2*sqrt(a + b*coth(x)^2))]
    @test_int [coth(x)/(a + b*coth(x)^2)^(5/2), x, 6, atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b))/(a + b)^(5/2) - 1/(3*(a + b)*(a + b*coth(x)^2)^(3/2)) - 1/((a + b)^2*sqrt(a + b*coth(x)^2))]
    @test_int [tanh(x)/(a + b*coth(x)^2)^(5/2), x, 9, -(atanh(sqrt(a + b*coth(x)^2)/sqrt(a))/a^(5/2)) + atanh(sqrt(a + b*coth(x)^2)/sqrt(a + b))/(a + b)^(5/2) + b/(3*a*(a + b)*(a + b*coth(x)^2)^(3/2)) + (b*(2*a + b))/(a^2*(a + b)^2*sqrt(a + b*coth(x)^2))]
    @test_int [tanh(x)^2/(a + b*coth(x)^2)^(5/2), x, 7, atanh((sqrt(a + b)*coth(x))/sqrt(a + b*coth(x)^2))/(a + b)^(5/2) + (b*tanh(x))/(3*a*(a + b)*(a + b*coth(x)^2)^(3/2)) + (b*(7*a + 4*b)*tanh(x))/(3*a^2*(a + b)^2*sqrt(a + b*coth(x)^2)) - ((3*a + 2*b)*(a + 4*b)*sqrt(a + b*coth(x)^2)*tanh(x))/(3*a^3*(a + b)^2)]

    @test_int [1/sqrt(1 + coth(x)^2), x, 3, atanh((sqrt(2)*coth(x))/sqrt(1 + coth(x)^2))/sqrt(2)]
    @test_int [1/sqrt(-1 - coth(x)^2), x, 3, atan((sqrt(2)*coth(x))/sqrt(-1 - coth(x)^2))/sqrt(2)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*coth(e+f*x)^m*(a+b*coth(e+f*x)^3)^p=#

    @test_int [1/(1 + coth(x)^3), x, 6, x/2 - (2*atan((1 - 2*coth(x))/sqrt(3)))/(3*sqrt(3)) - 1/(6*(1 + coth(x)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*coth(e+f*x)^m*(a+b*coth(e+f*x)^4)^p=#

    @test_int [coth(x)*(a + b*coth(x)^4)^(1/2), x, 8, (-(1/2))*sqrt(b)*atanh((sqrt(b)*coth(x)^2)/sqrt(a + b*coth(x)^4)) + (1/2)*sqrt(a + b)*atanh((a + b*coth(x)^2)/(sqrt(a + b)*sqrt(a + b*coth(x)^4))) - (1/2)*sqrt(a + b*coth(x)^4)]
    @test_int [coth(x)/(a + b*coth(x)^4)^(1/2), x, 4, atanh((a + b*coth(x)^2)/(sqrt(a + b)*sqrt(a + b*coth(x)^4)))/(2*sqrt(a + b))]
    @test_int [coth(x)/(a + b*coth(x)^4)^(3/2), x, 6, atanh((a + b*coth(x)^2)/(sqrt(a + b)*sqrt(a + b*coth(x)^4)))/(2*(a + b)^(3/2)) - (a - b*coth(x)^2)/(2*a*(a + b)*sqrt(a + b*coth(x)^4))]

    #= ::Section:: =#
    #=Integrands*of*the*form*coth(e+f*x)^m*(a+b*coth(e+f*x)^n)^p=#
    end

max_stars_repo_path: test/6 Hyperbolic functions/6.4 Hyperbolic cotangent/6.4.7 (d hyper)^m (a+b (c coth)^n)^p.jl